1. Introduction
Since the early work of Sandström (Reference Sandström1908) on marine glacial discharges in Norwegian fjords and his pioneering work on ocean circulation through differential input of fresh/cold and salty/warm water, attempts to link differential heating (Rossby Reference Rossby1965, Reference Rossby1998) and/or salt and fresh water input to a deep meridional circulation capable of driving the world's ocean circulation led to series of results predicting that differential buoyancy forcing on a horizontal surface alone could not explain the deep water cycle that, over a millennial time scale, overturns the global ocean (Defant Reference Defant1961). Paparella & Young (Reference Paparella and Young2002) (hereafter, PY) derived for horizontal convection (HC) a bound on the amount of kinetic energy dissipation 
$\epsilon$ and argued that in the limit of infinitely large Rayleigh numbers with constant Prandtl number 
$\epsilon$, properly normalised, vanishes. This is different from Rayleigh–Bénard convection (RBC), where it is expected to reach a finite value in what is known as the ‘ultimate’ regime. Since the publication of PY's anti-turbulence theorem, as is known, HC has been shown to undergo turbulence (see e.g. Scotti & White Reference Scotti and White2011; Gayen, Griffiths & Hughes Reference Gayen, Griffiths and Hughes2014), although with a different landscape in the Prandtl number–Rayleigh number phase space (Shishkina & Wagner Reference Shishkina and Wagner2016). This is in part due to the rate at which energy is dissipated which, so far, has been shown to consistently follow laminar-type scaling laws with respect to the magnitude of the forcing (Rossby Reference Rossby1998; Hughes et al. Reference Hughes, Griffiths, Mullarney and Peterson2007; Shishkina & Wagner Reference Shishkina and Wagner2016), as long as the buoyancy gradient is unidirectional (Griffiths & Gayen Reference Griffiths and Gayen2015). This is despite the fact that the flow is unstable with respect to two- and three-dimensional perturbations (Gayen et al. Reference Gayen, Griffiths and Hughes2014; Passaggia, Scotti & White Reference Passaggia, Scotti and White2017). An analogue of Rayleigh–Bénard theory was proposed and applied to HC to characterise the regime diagram of laminar and turbulent regimes (Shishkina, Grossmann & Lohse Reference Shishkina, Grossmann and Lohse2016). The 
$Pr=O(0.1\unicode{x2013}10)$ regime has been investigated numerically (Shishkina & Wagner Reference Shishkina and Wagner2016; Reiter & Shishkina Reference Reiter and Shishkina2020; Tsai et al. Reference Tsai, Hussam, King and Sheard2020; Constantinou et al. Reference Constantinou, Rocha, Smith and Young2023). In this work, we investigate numerically the low-Prandtl-number region of this regime diagram. We aim to investigate which limiting regimes are effectively observed and where turbulent HC starts to appear. According to the recent work of Shishkina et al. (Reference Shishkina, Grossmann and Lohse2016), the transition to the limiting turbulent regime, appearing at sufficiently high Rayleigh numbers, should be observed first at low Prandtl numbers, and we show here that this is indeed the case.
Griffiths & Gayen (Reference Griffiths and Gayen2015) considered the problem of HC forced by spatially periodic forcing. Their result showed that HC becomes turbulent in the core. Their forcing, localised on a length scale smaller than the depth of the domain and with fluctuations in both horizontal directions, shows turbulence throughout the domain, a regime transition to a dominant domain-scale circulation and a region of logarithmic velocity in the boundary layer (BL). The same geometry was further analysed by Rosevear, Gayen & Griffiths (Reference Rosevear, Gayen and Griffiths2017) where they observed that the non-dimensional heat flux, the Nusselt number, had a steeper scaling with respect to the Rayleigh number than the (laminar) Rossby (Reference Rossby1965) scaling. Their scaling analysis suggests that for deep enough domains, the flow is fully driven by the interior core of the flow, located under the surface BLs. One interesting fact is that despite the existence of a log layer in their direct numerical simulation (DNS), they did not observe log-type corrections in the scaling for the heat transfer. This is relatively surprising since it is now well established in RBC that heat transfers are buffered through the log layer (Grossmann & Lohse Reference Grossmann and Lohse2011; Ahlers et al. Reference Ahlers, Bodenschatz, Funfschilling, Grossmann, He, Lohse, Stevens and Verzicco2012; Ahlers, Bodenschatz & He Reference Ahlers, Bodenschatz and He2014).
Shishkina & Wagner (Reference Shishkina and Wagner2016) report a similar exponent in the case of large Prandtl numbers and low Rayleigh numbers. Their study shows that when the BL extends to the bottom of the domain, HC is highly effective in transporting heat. In addition, they report the dependence on Prandtl numbers for both their newly identified regime, denoted by 
$I^*_l$, and the Rossby regime, denoted by 
$I_l$. In our analysis, we follow the same nomenclature in an attempt to unify all known results. Note that these results are also observed experimentally in the companion paper (Passaggia, Cohen & Scotti Reference Passaggia, Cohen and Scotti2024).
Recently, Reiter & Shishkina (Reference Reiter and Shishkina2020) analysed classical and symmetrical HC in Rayleigh numbers up to 
$10^{12}$ and three Prandtl numbers. They found that for large Rayleigh numbers at 
$10^{11}$ and large-aspect-ratio domains, the plume detaches and exhibits low-frequency oscillations, whereas the Nusselt exhibits locally a steeper scaling. In this paper, we confirm these results in a different set-up and extend the Rayleigh number range by three orders of magnitude, up to 
$10^{15}$.
Although the ratio of viscosity to heat diffusion, taken here as the Prandtl number (
$Pr$), is 
$O(1)$ or greater in atmospheric and oceanic applications, HC at low Prandtl numbers has interesting geophysical applications, such as in the highly thermally conductive part of Earth's mantle. Although much attention has been devoted to RBC for the outer core dynamics, it is only very recently that HC has attracted the attention of planetary scientists (Alboussiere, Deguen & Melzani Reference Alboussiere, Deguen and Melzani2010). For example, Takehiro (Reference Takehiro2011) suggested that it could be a potential mechanism to drive zonal heat and momentum near the inner core through the Joule effect due to the Earth's magnetic field. At the edge of the inner core of Earth, horizontal regions of thermally stable (crystallising) and unstable (melting) may explain the zonal asymmetry of the inner core (Alboussiere et al. Reference Alboussiere, Deguen and Melzani2010). However, very little is known about the properties of the turbulent horizontal flows generated in these regions, and HC appears to be an interesting candidate to model such flows.
In this study, we use DNS to investigate how the Reynolds number (
$Re$) and the Nusselt number (
$Nu$) depend on the Rayleigh number (
$Ra$) and the Prandtl number (
$Pr$) in turbulent HC at low to intermediate 
$Pr$ for values characteristic of convection in gases (
$0.1< Pr<1$; Taylor, Bauer & McEligot Reference Taylor, Bauer and McEligot1988; Roche et al. Reference Roche, Castaing, Chabaud and Hébral2002), and liquid metals where 
$Pr= {O}(10^{-2})$ (Takehiro Reference Takehiro2011). The results agree with the scaling power laws recently derived by Shishkina et al. (Reference Shishkina, Grossmann and Lohse2016) (hereafter, SGL) based on the original work of (Grossmann & Lohse Reference Grossmann and Lohse2000) (hereafter, GL) on RBC, and the numerical simulations of Takehiro (Reference Takehiro2011). Furthermore, we provide evidence that the regime observed by Rosevear et al. (Reference Rosevear, Gayen and Griffiths2017) generalises to HC over a monotonic temperature profile with a turbulent log layer which indeed acts as a buffer to heat transfer (i.e. when compared with homogeneous turbulence) and slightly decreases the exponent reported previously. We also provide a connection between the GL theory and the plume-driven dynamics derived by Hughes et al. (Reference Hughes, Griffiths, Mullarney and Peterson2007).
Our simulations cover the laminar regime 
$I_l$ (Shishkina & Wagner Reference Shishkina and Wagner2016), the high-
$Pr$ laminar regime 
$I^*_l$ recently reported by Shishkina & Wagner (Reference Shishkina and Wagner2016) and a new low-
$Pr$ turbulent regime, named 
$II_l$, which is a new turbulent limiting regime for HC. We also observe the plume-dominated flow regime of Hughes et al. (Reference Hughes, Griffiths, Mullarney and Peterson2007), which we call 
$I^+_u$ according to the nomenclature introduced in SGL. We also report the existence of the turbulent interior-dominated regime 
$IV_u$ at high Rayleigh numbers amended with the appropriate log-type corrections. An important contribution of our work is that these results agree and extend known HC regimes diagrams at low 
$Ra$ (Hughes & Griffiths Reference Hughes and Griffiths2008) to fit within the theoretical prediction of SGL.
Finally, we further explore the relation between the Reynolds number characterising the magnitude of the overturning flow and turbulent dissipation. This analysis allows for condensing this complex regime transitions diagram into a more traditional laminar to transitional to and turbulent regime. In analogy with RBC (Lohse & Shishkina Reference Lohse and Shishkina2023), the turbulent regime is further divided into a classical regime, characterised by a normalised dissipation which still depends on the Reynolds number, and a ultimate regime where the normalised dissipation leans towards a constant, thereby achieving complete similarity with respect to the Reynolds number. We further confirm that a ultimate turbulent regime cannot be observed for the Prandtl numbers considered in this work for the present values of 
$Ra$ which are too small to attain a such regime. This assumption is theoretically justified by the fact that the Richardson number (the ratio of stratification to shear) in the stably stratified layer remains asymptotically above the threshold for instability, suggesting that for 
$Pr= {O}(1)$ or less, the layer under the warming boundary remains laminar.
Similarly to Shishkina et al. (Reference Shishkina, Grossmann and Lohse2016), we exploit the idea that in turbulent thermal convection, the thermal and viscous dissipation rates spatially and temporally averaged are determined to leading order by their bulk or BL contributions.
The rest of the paper is organised as follows. In § 2 we provide the mathematical background for the problem, followed in § 3 by a review of the previously known laminar and turbulent HC regimes. Section 4 describes the numerical simulations; § 5 presents the results and the scaling analysis; § 6 considers the question of whether ultimate turbulence is possible in HC; and, finally, in § 7 we present conclusions. The high-Prandtl-number regime is considered in the companion paper (Passaggia et al. Reference Passaggia, Cohen and Scotti2024).
2. Problem description
Here, we consider the problem of convection in the Boussinesq limit, where the buoyancy difference 
$(\rho _{max}-\rho _{min})/\rho _{min}$ imposed on the upper boundary, normalised by the bulk density 
$\rho _{0}$, is small. In this limit, the equations describing the motion of the flow are
$$\begin{gather} \frac{{\rm D} \boldsymbol{u}}{{\rm D} t} ={-}\boldsymbol{\nabla} p + b\boldsymbol{e_z} + \left(\frac{Pr}{Ra}\right)^{1/2}\nabla^2 \boldsymbol{u}, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}$$
$$\begin{gather}\frac{{\rm D} b}{{\rm D} t} = (Pr Ra)^{{-}1/2}\nabla^2 b, \end{gather}$$
where 
${\rm D}/{\rm D} t$ denotes the material derivative, 
$\boldsymbol {u}=(u,v,w)^{\rm T}$ is the velocity vector, 
$b=-g(\rho -\rho _{min})/\rho _{min}$ is the buoyancy, 
$g$ is the acceleration of gravity along the vertical unit vector 
$\boldsymbol {e}_z$ and 
$p$ is the hydrodynamic pressure.
The equations are non-dimensionalised using the length of the tank 
$L$ and the free-fall velocity 
$\sqrt {{\rm \Delta} L}$, where 
${\rm \Delta}$ is the buoyancy difference imposed along the upper surface (note that, unlike RBC, it is the length of the tank, rather than its depth, that is dynamically relevant). The Prandtl number is given by 
$Pr=\nu /\kappa$ where 
$\nu$ and 
$\kappa$ are the viscosity and diffusivity of the stratifying agent, respectively. The Rayleigh number is defined as 
$Ra={\rm \Delta} L^3/(\nu \kappa )$ The computational domain, shown in figure 1, is a parallelepiped with aspect ratio 
$\varGamma =L/H=4$ with dimensions 
$[L,W,H]=[1,1/8,1/4]$ where 
$W$ is the width of the computational domain (Scotti & White Reference Scotti and White2011). The parameter landscape of HC is wider than traditional RBC. In addition to selecting a geometry and velocity boundary conditions along the boundary of the domain, an infinite number of possibilities exist in principle to impose the buoyancy forcing at the surface. Regarding the latter, previous numerical studies have used a linear profile, a sinusoidal profile or a step function (see table 1). Here, we impose a buoyancy profile along the upper surface 
$z=H$, where 
$H$ is the height of the domain, that 
$b(x)|_{z=H}=(1+\tanh (10x/L))/2$, which is a smoothed version of the sharp profile used in previous calculations (Shishkina & Wagner Reference Shishkina and Wagner2016; Passaggia et al. Reference Passaggia, Scotti and White2017). This is necessary to keep the numerical simulations stable. It does not appear that there is a consensus within the literature with regard to the velocity boundary conditions, with numerical studies employing either no-slip or free-slip boundary conditions or even a mixture of both (table 1). We use free-slip conditions along top, bottom, left and right boundaries, and periodic conditions along the spanwise direction. Previous studies suggest that the flow is not particularly sensitive to the lateral and bottom boundary conditions (Beardsley & Festa Reference Beardsley and Festa1972; Rossby Reference Rossby1998; Chiu-Webster et al. Reference Chiu-Webster, Hinch and Lister2008). Geophysically inspired studies tend to use free-slip along the surface where the forcing is applied, which is arguably more appropriate as a boundary condition for the ocean than a no-slip condition. Moreover, free-slip conditions sidestep the numerical difficulties involved in resolving no-slip BLs and allow us to efficiently explore Rayleigh numbers as high as 
$1.92\times 10^{15}$ for a wide range of Prandtl numbers. In what follows, we define integral values to be linked with the control parameters 
$Ra$ and 
$Pr$: the magnitude of the large-scale flow is used to define the Reynolds number 
$Re={(\langle {\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {u}}\rangle )}^{1/2}L/\nu$ where 
$\langle {\cdot }\rangle$ denotes the spatiotemporal average over the computational domain; the Péclet number 
$Pe=RePr$ is defined similarly. We define the Nusselt number as the ratio of the buoyancy flux averaged over the stabilising side of the domain over the same quantity obtained for a purely conducting configuration, that is 
$Nu=\overline {\partial b/\partial z}|_{z=H}/\varPhi _c$, where the over bar indicates a horizontal average of the left half of the domain and 
$\varPhi _c=\overline {\partial b_c/\partial z}|_{z=H}$ is the average gradient in the purely conducting case (i.e. when 
$Ra<10^3$; Siggers, Kerswell & Balmforth Reference Siggers, Kerswell and Balmforth2004). Note that other definitions have been considered (Rocha et al. Reference Rocha2020a). For the geometry considered here, its non-dimensional value was found numerically to be 
$\varPhi _c=0.53$.
Schematic of the present showing the different length scales of the thermal BL 
$\lambda _b$ and the kinetic BL 
$\lambda _u$ together with the full depth 
$H$ and the large overturning scale 
$h$ used for the theoretical prediction.
Forcing and velocity boundary conditions of previous HC studies.
3. Previously known regimes of laminar and turbulent HC
In this section, we review the existing scaling laws for heat and momentum exchanges in HC. The areas of the Rayleigh–Prandtl landscape previously explored using DNSs and experiments are reported, together with the sections investigated in this paper and its companion.
3.1. Rossby's (Reference Rossby1965) laminar scaling
Rossby (Reference Rossby1965) explored HC induced by differential heating in a parallelepipedic container with an aspect ratio 
$L/H=2.5$ and derived a scaling law relating the Nusselt number as a function of the Rayleigh number. Rossby collapsed temperature measurements from his experiments with a scaling relationship between the thickness of the BL and the streamfunction (see p. 13 in Rossby Reference Rossby1965). Taking the curl of (2.1a), neglecting the nonlinear terms, and defining the streamfunction 
$\psi$ such that 
$\psi _x=-w$ and 
$\psi _z=u$, the two-dimensional Navier–Stokes equations reduce to
$$\begin{gather} \left(\frac{Pr}{Ra}\right)^{1/2}\nabla^4 \psi = \partial_x b, \end{gather}$$
$$\begin{gather}-\partial_x\psi\partial_zb = (Pr Ra)^{{-}1/2}\nabla^2 b. \end{gather}$$
Near the conducting wall, the flow is governed by the laminar BL whose thickness is defined by 
$\lambda$ and (3.1a), (3.1b) reduce to leading order to
\begin{equation} \left(\frac{Pr}{Ra}\right)^{1/2} \frac{\psi}{\lambda^4} \sim \frac{\varDelta}{L} \quad \mbox{and} \quad \frac{\psi \varDelta}{\lambda L} \sim (Pr Ra)^{{-}1/2} \frac{\varDelta}{\lambda^2}. \end{equation}Combining these equations, Rossby obtained the relation
\begin{equation} \lambda \sim L Ra^{{-}1/5}. \end{equation}
What Rossby failed to recognise in his original work was that the thickness of the BL 
$\lambda$ could be defined using either the thermal BL thickness, denoted by the subscript 
$_b$ or the momentum BL denoted by the subscript 
$_u$, a point that becomes important when the viscosity and diffusivity are not close. While this has no implication for the 
$Ra$-dependence as shown in the next subsection, the Prandtl-number dependence may not be predicted accurately for different values of the Rayleigh number 
$Ra$ and varying 
$Pr$.
3.2. Paparella & Young (Reference Paparella and Young2002) inequality
HC and RBC are both closed systems driven by buoyancy fluxes imposed at the boundaries. However, whereas in the latter the total dissipation in the flow is proportional to the Nusselt number, in the former it is possible to exactly relate it to the Rayleigh number via (Paparella & Young Reference Paparella and Young2002)
\begin{equation} \overline{\epsilon_u} = B(\varGamma/2)\nu^{3}L^{{-}4}RaPr^{{-}2}, \end{equation}
where 
$0< B\leq 1$, which, once combined with the original idea of Rossby, opens possibilities for relating the dissipation in the BL or the core with the heat transfer coefficient near the horizontal boundary.
One interesting consequence is that as 
$Ra$ increases while keeping 
$Pr$ and 
$\varGamma$ constant, the flow becomes progressively confined under the conducting boundary. This effect is also known as the anti-turbulence theorem and implies that beyond a certain point, the overturning depth scale becomes
\begin{equation} h< H, \end{equation}and a zone of stratified fluid nearly at rest will form along the bottom of the tank.
Although such regimes were only observed in DNSs of laminar HC (Ilicak & Vallis Reference Ilicak and Vallis2012) at high 
$Pr$ and were theoretically predicted by Chiu-Webster et al. (Reference Chiu-Webster, Hinch and Lister2008) for the same range of parameters, experiments by Wang & Huang (Reference Wang and Huang2005) show the onset of such behaviour at intermediate 
$Pr$ and relatively low 
$Ra$.
3.3. Shishkina & Wagner (Reference Shishkina and Wagner2016) revisited laminar regime 
$I_l$
Rossby's laminar regime was recast in Shishkina & Wagner (Reference Shishkina and Wagner2016) to obtain a more accurate prediction for the Prandtl number dependence. The idea is to start with the steady thermal BL equation, which is obtained from (2.1c) and write an advection–diffusion balance in the BL
\begin{equation} u b_x + v b_z = \kappa b_{zz}. \end{equation}
The dominant terms in this expression reduce to 
$U\varDelta /L = \kappa \varDelta /\lambda _b^2$ where 
$\lambda _b$ is the thickness of the thermal BL, which scales as 
$\lambda _b/L \sim Nu^{-1}$. Combining the above reduces to the well-known thermal–laminar BL scaling
\begin{equation} Nu=Re^{1/2}Pr^{1/2}, \end{equation}
provides a relation tying 
$Nu$, 
$Re$ and 
$Pr$. In laminar regimes, the buoyancy variance is essentially concentrated in the BL and reads
\begin{equation} \overline{\epsilon_{b, {BL}}} \sim \kappa \frac{\varDelta^{2}}{\lambda_{b}^{2}} \frac{\lambda_{b}}{H}=\kappa \frac{\varDelta^{2}}{H^{2}} \frac{\lambda_{u}}{\lambda_{b}} {Re}^{1/2}, \end{equation}
where the dependence on the aspect ratio 
$\varGamma$ was omitted. Since the thickness of the laminar momentum BL scales as 
$\lambda _u/H \sim Re^{-1/2}$, the scaling of the mean dissipation in the particular case of a laminar BL (Landau & Lifschitz Reference Landau and Lifschitz1987; Grossmann & Lohse Reference Grossmann and Lohse2000) is
\begin{equation} \overline{\epsilon_{u,BL}}\sim\nu\frac{U^2}{\lambda^2_u}\frac{\lambda_u}{H}=\nu^3H^{{-}4}Re^{5/2}. \end{equation}Combining (3.7), (3.4) and (3.9), one recovers the laminar scaling (Rossby Reference Rossby1965, Reference Rossby1998; Gayen et al. Reference Gayen, Griffiths and Hughes2014; Shishkina et al. Reference Shishkina, Grossmann and Lohse2016)
$$\begin{gather} Re \sim Ra^{2/5}Pr^{{-}4/5}, \end{gather}$$
$$\begin{gather}Nu \sim Ra^{1/5}Pr^{1/10}. \end{gather}$$
By analogy to the notation in the GL theory for RBC, this scaling regime is denoted as 
$I_l$, where the subscript 
$l$ stands for low-
$Pr$ fluids.
3.4. Hughes et al.'s (Reference Hughes, Griffiths, Mullarney and Peterson2007) laminar BL/turbulent plume regime $I^+_u$
Increasing 
$Ra$ and for intermediate 
$Pr$, the momentum BL becomes progressively thinner, whereas the thermal boundary remains relatively thick in comparison. In this case, it is the thermal BL that drives the dynamics and leads to a turbulent plume, detached from the bottom (see figure 2b). This particular case was theorised by Hughes et al. (Reference Hughes, Griffiths, Mullarney and Peterson2007) with a plume model inside a filling box. Here, we recast their model according to the SGL theory (i.e. see the plume model definition (2.15)–(2.20) in Hughes et al. Reference Hughes, Griffiths, Mullarney and Peterson2007) and the dissipation in the BL is balanced by the ratio between the thermal and the kinetic BL 
$\lambda _b/\lambda _u$ which reads
\begin{equation} {\overline{\epsilon_{u}}_{HGMP} \sim\nu\frac{U^2}{\lambda^2_u} \frac{\lambda_u}{H}\frac{\lambda_b}{\lambda_u}= \nu^{3}H^{{-}4}Re^{5/2}Pr^{{-}1/2}, } \end{equation}where the dissipation now scales with the thickness of the thermal, rather than the momentum layer. Now, combining (3.7), (3.4) and (3.11), the heat and momentum exchanges are given by
$$\begin{gather} Re\sim Ra^{2/5}Pr^{{-}3/5}, \end{gather}$$
$$\begin{gather}Nu\sim Ra^{1/5}Pr^{1/5}, \end{gather}$$
which we denote as 
$I^+_u$ since it has a similar power laws for 
$Ra$ to the 
$I_l$ regime but different 
$Pr$ dependence. Note that this scaling is partially supported in the experiments of Mullarney, Griffiths & Hughes (Reference Mullarney, Griffiths and Hughes2004) and Wang & Huang (Reference Wang and Huang2005), as well as in the DNSs of Gayen et al. (Reference Gayen, Griffiths and Hughes2014).
Snapshot of iso-contours of 
$\varLambda _2=-Pr^{-1}/8$ criteria (blue) and buoyancy 
$b$ (background) at (a) 
$Ra=6.4\times 10^{13}$, 
$Pr=0.01$ showing the new regime (
$II^*_l$) and (b) 
$Ra=6.4\times 10^{13}$, 
$Pr=1$ corresponding to the Hughes regime 
$I^+_u$.
3.5. Shishkina & Wagner (Reference Shishkina and Wagner2016) laminar regime $I^*_l$
At low 
$Ra$ and for large 
$Pr$ and/or large aspect ratio 
$\varGamma$, the BL thickness 
$\lambda _u$ saturates and reaches the depth of the domain that gives 
$\lambda _u\approx h=H$ and (3.4) becomes equivalent to the dissipation in a pressured-driven laminar flow in a channel which writes
\begin{equation} \overline{\epsilon_{u}}_{SW} \sim\nu\frac{U^2}{H^2}= \nu^{3}H^{{-}4}Re^{2}. \end{equation}Combining (3.7), (3.4) and (3.13), one obtains the laminar scaling derived Shishkina & Wagner (Reference Shishkina and Wagner2016)
$$\begin{gather} Re\sim Ra^{1/2}Pr^{{-}1}, \end{gather}$$
$$\begin{gather}Nu\sim Ra^{1/4}Pr^{0}, \end{gather}$$
denoted as 
$I^*_l$ and first observed by Beardsley & Festa (Reference Beardsley and Festa1972) in their numerical simulations. Note that Rossby also observed a steeper scaling than 
$Nu\sim Ra^{1/5}$ in his numerical simulations for low 
$Ra$ (see p. 245 in Rossby Reference Rossby1998). More recent numerical simulations at infinite Prandtl numbers (Ramme & Hansen Reference Ramme and Hansen2019) suggest a similar scaling. However, it should be noted that the observation of this scaling is confined to less than a decade of Rayleigh numbers. It is important to stress that in this regime, the circulation is assumed to span the entire domain, and this particular aspect will be further explored in the present and companion paper.
3.6. Turbulent regimes and associated bounds
Most of the existing work on HC highlighted laminar-type flows, dominated by the BL behaviour, except for an analogue of HC Griffiths & Gayen (Reference Griffiths and Gayen2015) and Rosevear et al. (Reference Rosevear, Gayen and Griffiths2017). They considered a spatially periodic forcing at the conducting boundary with a short wavelength compared with the depth of the domain. In this particular set-up, Rosevear et al. (Reference Rosevear, Gayen and Griffiths2017) were able to show that 
$Nu\sim Ra^{1/4}$ with a turbulent core driven by inertia. They were also the first to report turbulent BLs with a log-type layer developing along the conducting wall. This feature is important because it is a necessary condition for turbulent convection scaling to arise in RBC (Grossmann & Lohse Reference Grossmann and Lohse2000, Reference Grossmann and Lohse2011).
Siggers et al. (Reference Siggers, Kerswell and Balmforth2004) theorised a fully turbulent regime based on the assumption that the BLs and thus the Nusselt number scales as 
$Nu\sim Ra^{1/3}$. Recent work by Rocha et al. refined their original results and showed that short oscillations with large amplitudes of the forcing boundary were needed to observe such a regime when 
$Ra\rightarrow \infty$. Indeed, if we recast the PY bound on dissipation with (3.4) we obtain a bound on the Richardson number in the thermal layer under the warm side
\begin{align} Ri &\equiv \frac{\partial_z b}{(\partial_z U)^2} \sim \frac{\varDelta}{\lambda_b}\left(\frac{\lambda_u}{U}\right)^2 \sim {\rm \Delta} L^{{-}1}\frac{L}{H}\frac{H}{\lambda_b}L^2\frac{\lambda_u^2}{L^2} \frac{\nu^2}{L^2U^2}{L^2}{\nu^{{-}2}}\nonumber\\&\quad \sim \varGamma Ra Pr^{{-}1} Nu Re^{{-}2}\frac{\lambda_u^2}{L^2}. \end{align}
For a fully turbulent regime, the BLs must scale as 
$Ra^{-1/3}$, thus 
$Nu \sim Ra^{1/3}$ (see Siggers et al. (Reference Siggers, Kerswell and Balmforth2004) and Rocha et al. (Reference Rocha, Bossy, Llewellyn-Smith and Young2020b) for an improved estimate). Thus,
\begin{equation} Ri\sim \varGamma Pr^{{-}1}Ra^{2/3}Re^{{-}2}. \end{equation}
For 
$Ri$ to become asymptotically small (at fixed 
$Pr$), we must have 
$Re\sim Ra^\alpha$, with 
$\alpha >1/3$. However, under these assumptions, the normalised dissipation 
$\overline {\epsilon _u} L/U^3< Ra^{1-3\alpha }$ (see (3.4)) would decay asymptotically, violating the hypothesis of fully turbulent flow. This means that the stability of the shear layer is, to leading order, Rayleigh number independent.
This behaviour was observed by Whitehead & Wang (Reference Whitehead and Wang2008) in their laboratory experiments, where a forcing was added. In summary, the Richardson number cannot decrease with increasing 
$Ra$, and therefore there are no reasons to believe that increasing 
$Ra$ will eventually lead to the development of an unstable layer.
Of course, this conclusion is only valid for non-rotating HC and would not hold in the case where rotation is taken into account (Barkan, Winters & Smith Reference Barkan, Winters and Smith2013; Vreugdenhil, Griffiths & Gayen Reference Vreugdenhil, Griffiths and Gayen2017).
4. Numerical calculations
The Navier–Stokes equations are solved numerically on a Cartesian grid, stretched near the upper boundary, using a standard second order in space and time projection method.
Snapshots of the flow are shown in figure 2(a,b) and by mean of the 
$\varLambda _2$ criteria defined by the second largest eigenvalue of the matrix 
$S^2+\varOmega ^2$ where 
$S$ and 
$\varOmega$ are the symmetric and anti-symmetric parts of the velocity gradient tensor, respectively.
Simulations cover the range 
$Ra=[6.4\times 10^5,1.92\times 10^{15}]$ and 
$0.002\leq Pr\leq 2$. For 
$Ra<10^8$ and 
$0.5\leq Pr\leq 2$, the computed flows are steady (Shishkina & Wagner Reference Shishkina and Wagner2016; Passaggia et al. Reference Passaggia, Scotti and White2017). With increasing 
$Ra$ and/or decreasing values of 
$Pr$, the flows become increasingly unsteady, leading to turbulence (as shown in figure 2c) and the mesh is refined in order to resolve the Kolmogorov length scale. In the case of homogeneous turbulence, the Kolmogorov length scale is given by 
$\eta \approx (\nu ^3/\epsilon _u)^{1/4}$. Using the PY inequality, an approximation yields 
$\eta /L\approx Pr^{1/2}/(\varGamma B Ra)^{1/4} \gtrsim 10^{-4}$ for the highest value of 
$Ra$ and the lowest 
$Pr$ considered in this work. These estimates are valid for homogeneous turbulence. In the present case, a substantial amount of dissipation is located in the BL, where the mesh is refined up to 
${\rm \Delta} z/L = 10^{-4}$ in the vertical direction and 
${\rm \Delta} y/L = 10^{-3}$ in the horizontal direction and should ensure numerical convergence. Mesh sizes are reported in figure 3(b) in the 
$(Ra,Pr)$ plane along with the different regimes reported later in this paper. Note that the turbulence in HC for moderate values of 
$Pr$ is confined to a narrow region located below the cooling / heavy boundary consisting of the plume and the BL where the fluid is statically unstable (cf. figure 2) (Scotti & White Reference Scotti and White2011; Gayen et al. Reference Gayen, Griffiths and Hughes2014). Decreasing the value of 
$Pr$ increases the volume of fluid subject to turbulence (see figure 2) and decreases the depth of the circulation.
Sketch of the phase diagram in the 
$(Ra, Pr)$ plane for the laminar regimes 
$I_l$ and 
$I^*_l$ together with the turbulent scaling 
$II_{l}$ with the conducted DNS. The yellow stripes show the transition from 
$I^*_l$ to 
$I_{l}$, and 
$I_l$ to 
$II_{l}$, with a slope 
$Pr\approx Ra^{1/2}$. The transition from 
$II^*_l$ to 
$I^+_u$ with slope 
$Pr\approx Ra^{-1}$. The symbols reflect the computational meshes in 
$(x,y,z)$, used in the DNS: 
$512\times 256\times 256$ (circle), 
$1024\times 384\times 128$ (squares) and 
$2048\times 256\times 256$ (triangles). The values (
$\alpha,\beta$) in each region provide the exponents 
$Nu\sim Ra^{\alpha }Pr^{\beta }$ measured in the DNS and derived theoretically, with the exception of the exponents in the region 
$IV_u$, which are empirically derived from the DNS. The transition between regimes is given as follows: from 
$I^*_l$ to 
$I_l$ is given by 
$Pr \sim 3.0\times 10^{-5} Ra^{1/2}$, from 
$I_l$ to 
$II^*_l$ by 
$Pr \sim 3.0\times 10^{-6} Ra^{1/2}$, from 
$I^+_u$ to 
$IV_u$ by 
$Pr \sim 3.0\times 10^{-8} Ra^{1/2}$ and from 
$II^*_l$ to 
$I^+_u$ by 
$Pr \sim 3.0\times 10^{11} Ra^{-1}$. The black dashed line shows when the BL becomes turbulent and follows from 
$II^*_l$ to 
$IV_u$ by 
$Ra \approx 3.0\times 10^{12}$. The red dashed line shows the transition from 
$I^+_u$ to 
$IV_u$ but the scaling is 
$Ra$ dependent and was approximated from the data obtained in this study.
5. Results and scaling analysis
The regimes observed in our numerical simulations are summarised in figure 3 together with the exponents computed fitting power laws to the data (figure 4a–d)). The colours shown in figure 3 correspond to the colours shown in figure 4(a–d) and the two new regimes, labelled 
$II^*_l$ in purple and cyan whereas 
$IV_u$ is shown in brown in figure 3 for turbulent scaling laws. Table 2 summarizes the scaling exponents and numerical prefactors for each regime.
(a,c) 
$Ra$ dependencies and (b,d) 
$Pr$ dependencies of (a,b) the Nusselt number and (c,d) the Reynolds number, as obtained in the DNS for (a,c) 
$Pr=1$ (squares), 
$Pr=0.1$ (circles), 
$Pr=0.01$ (triangles) and for (b,d) 
$Ra=6.4\times 10^{10}$ (diamonds), 
$Ra=1.92\times 10^{12}$ (pentagons) and 
$Ra=1.92\times 10^{14}$ (triangles). The DNS results support the scaling in the regime 
$I_l$ (solid lines) (3.10a) and (3.10b), transition to 
$II^*_l$ (dotted lines) (5.2a), (5.2b) and (5.12a), (5.12b) transition to 
$I^+_u$ (dotted lines) (3.12a) and (3.12b). The black squares in (a–d) are from the numerical simulations of Shishkina & Wagner (Reference Shishkina and Wagner2016), with 
$Pr=0.1$ in (a,c) and 
$Ra=10^{10}$ in (b,d).
List of scaling laws with prefactors measured for the Nusselt and Reynolds number dependencies across different regimes. The scaling exponents decorated with an asterisk are estimated from the numerical simulations. All other exponents are derived from theoretical considerations.
The dependence of 
$Nu$ and 
$Re$ with respect to 
$Ra$ and 
$Pr$ is summarised in figure 4(a–d). The Nusselt number obeys a scaling law 
$Nu\sim Ra^{\alpha }$ (see figure 4a) with the exponent 
$\alpha$ depending on 
$Ra$ and 
$Pr$ as follows:
(i)
$\alpha =1/4$, the enhanced laminar scaling, for low $Ra$ and higher $Pr$;(ii)
$\alpha =1/5$, the classical laminar scaling, for small $Ra$;(iii)
$\alpha =1/6$, for small $Pr$;(iv)
$\alpha =1/5$, the entrainment-type regime, at high $Ra$ and not too small $Pr$;(v)
$\alpha \approx {{0.225}}$ for large $Ra$ and small $Pr$.
We observe the laminar scaling 
$Re\sim Ra^{\gamma }$ with 
$\gamma =1/2$ (Shishkina & Wagner Reference Shishkina and Wagner2016) and 
$\gamma =2/5$ (Rossby Reference Rossby1965). At higher 
$Ra$, the new scaling 
$\gamma =1/3$ is also observed, reverting back to 
$\gamma =2/5$ (Hughes et al. Reference Hughes, Griffiths, Mullarney and Peterson2007) (figure 4c). Similarly, when 
$Pr<1$ and 
$Ra$ is fixed, we observe a scaling relationship 
$Nu\sim Pr^{\beta }$ with:
(i)
$\beta =0$ for higher $Pr$ and low $Ra$ (see Shishkina & Wagner Reference Shishkina and Wagner2016);(ii)
$\beta =1/10$ for $Ra<10^{11}$ (see Shishkina & Wagner Reference Shishkina and Wagner2016);(iii)
$\beta =7/24$ at low $Pr$;(iv)
$\beta =1/5$ for $Ra>5\times 10^{11}$ (see Hughes et al. Reference Hughes, Griffiths, Mullarney and Peterson2007);(v)
$\beta \approx 0.417$ at large $Ra$ and low $Pr$.
The Reynolds number dependence 
$Re \sim Pr^{\delta }$ with 
$\delta =-2/3$ for the smaller value of 
$Pr$, then 
$\delta =-1$ for 
$10^{-2}\lesssim Pr\lesssim 0.2$ changes to 
$\delta =-4/5$ with increasing 
$Ra$ regardless of 
$Pr$, and increases at high 
$Ra$ to the HGMP scaling 
$\delta =-3/5$ for the larger values of 
$Ra$ (figure 4d) considered.
Despite the difference in boundary conditions, our simulations recover the 
$I^*_l$ and 
$I_l$ regions in the phase diagram as found by Shishkina et al. (Reference Shishkina, Grossmann and Lohse2016) and the 
$I^+_u$ regime described by Hughes & Griffiths (Reference Hughes and Griffiths2008), albeit with different boundaries (figures 4(a,b) and 5(b)). For example, at 
$Pr=0.1$ and 
$Ra>10^8$, our simulations have already entered the 
$II^*_l$ regime, whereas in Shishkina & Wagner (Reference Shishkina and Wagner2016) they are still in the 
$I_l$ regime (figure 4a). This suggests that the scaling exponents are robust with respect to the boundary conditions. The scaling regimes that occur in the lower regions of the parameter space, with 
$\alpha =1/6, \beta =7/24,$ and 
$\alpha \approx {{0.225}}, \beta \approx 0.417$, are new to this study and are the focus of the following section.
(a,c) 
$Ra$ dependencies and (b,d) 
$Pr$ dependencies of (a,b) 
$NuRe^{-1/2}$ and (c,d) 
$L^4\nu ^{-3}\overline {\epsilon _u}Ra^{-1}$, as obtained in the DNS for (a,c) 
$Pr=1$ (squares), 
$Pr=0.1$ (circles), 
$Pr=0.01$ (triangles) and for (b,d) 
$Ra=10^9$ (diamonds), 
$Ra=2\times 10^{10}$ (pentagons) and 
$Ra=1.92\times 10^{14}$ (triangles). The upper figures support the estimates in (3.7) and (5.12b), whereas the lower figures confirm (3.4). The black squares in (a,b) are from the simulations of Shishkina & Wagner (Reference Shishkina and Wagner2016) at 
$Ra=10^{10}$.
5.1. The low-Prandtl-number core-driven flow $II_l$
5.1.1. The Shishkina et al. (Reference Shishkina, Grossmann and Lohse2016) regime
In the low-Prandtl-number regime, the flow transitions from Rossby's 
$I_l$ regime to the 
$II_l$ regime as 
$Ra$ increases. A snapshot of the flow (figure 2a) shows that, unlike the Rossby regime, the flow here is turbulent in the core. In this low-Prandtl-number regime, the buoyancy flux provided through the boundary is large but the thermal and momentum BLs remain thick, hence laminar, and (3.7) still holds. With decreasing 
$Pr$ and/or increasing 
$Ra$, the bulk dynamics dominates dissipation with a large-scale overturning flow occupying the entire domain and whose horizontal length scale is 
$L$. In this case, it is the large-scale velocity 
$U$ which drives the dissipation of kinetic energy and the latter is given by
\begin{equation} \overline{\epsilon_u} \sim \nu^{3}L^{{-}4}Re^3. \end{equation}
From (3.7), (3.4) and (5.1), it follows that low-
$Pr$ HC exhibits dependencies of the form
$$\begin{gather} Re\sim Ra^{1/3}Pr^{{-}2/3}, \end{gather}$$
$$\begin{gather}Nu\sim Ra^{1/6}Pr^{1/6}, \end{gather}$$
where this scaling regime is denoted as 
$II_l$ (see figure 3(b) and Shishkina et al. Reference Shishkina, Grossmann and Lohse2016). Note that these scalings are only observed for the Rayleigh number dependence, but the Prandtl number dependence is clearly underestimated for both 
$Nu$ and 
$Re$.
The BL scaling observed so far was are consistent with
\begin{equation} Nu\sim Pe^{1/2}. \end{equation}In the following, we show that the balance in the boundary has to be modified to take into account either core-size modifications or turbulent BL effects.
5.1.2. Modification induced by the variable turbulent depth $h$ for the $II_l$ regime
In the low-Prandtl-number regime (
$Pr < 10^{-1}$), turbulence is confined between the plume and the left part of the domain, under the statically unstable BL whose depth is indicated by 
$h< H$. The PY constraint on dissipation can be used to relate the depth 
$h$ occupied by the turbulent core to the Reynolds number 
$Re$. The inbalance between the dissipation in the bulk and the BL is the first occurrence of a turbulent regime subject to two regions with different dissipation rates at low Prandtl numbers. At low values of 
$Pr$, the thermal BL providing the available energy drives the dynamics near the forcing boundary, and its dissipation rate is given in (3.9). Dissipation in the core provides a higher dissipation rate, as given by (5.1) which therefore dissipates energy at a rate that grows faster than the supply by the forcing boundary. As a consequence, the bulk size 
$h$ must decrease with respect to the full depth of the domain 
$H$, to account for this effect.
Chiu-Webster et al. (Reference Chiu-Webster, Hinch and Lister2008) made similar observations in the case of an infinite Prandtl number. In their case, the core dissipates less than the BL, leaving the Nusselt number scaling independent of the core dynamics.
Here, we demonstrate that the turbulent core modifies the Prandtl number dependence by relating the dissipation and buoyancy variance in both the BL and the bulk. The idea is to introduce the ratio 
$h/H$ in the equations and to provide the Rayleigh and Prandtl number dependencies and correct for the above estimate.
We first provide a measure for 
$h/H$. It should first be stated that the variable turbulent depth 
$h$ is not uniform across the entire domain. The depth 
$h$ is somewhat more representative of the region where mixing occurs which is clearly defined by the area 
$\mathcal {A}$ where the turbulent dissipation 
$\epsilon _u > (Re/2)^3$. In other words, the area where mixing takes place is defined as
\begin{equation} \mathcal{A} = \int_{S} {\mathbb 1}|_{\epsilon_u^{1/3} > Re/2}\,\mbox{d}S = h^2, \end{equation}
and is related to the turbulent depth 
$h$, whereas 
${\mathbb 1}$ is the indicator function. It is also important to highlight that depending on the type regime, mixing occurs at different locations. For instance, in the low-Prandtl-number core-driven flow 
$II_l$, the region where mixing occurs is essentially located in the plume and the core. On the other hand, for the turbulent BL regime at large Rayleigh numbers 
$IV_u$, mixing is dominant in both the plume and the BL. Measurements of 
$h$ for the three cases varying the Prandtl number at constant Rayleigh numbers are shown in figure 6 where all measurements scale as 
$h/H \sim Pr^{1/4}$ when appropriately rescaled with the Reynolds number.
From a theoretical point of view, the bulk size can be determined assuming that the mixing efficiency in the bulk and the mixing efficiency in the BL follow similar scaling laws. This is justified since both bulk and BLs are driven by an excess of available potential energy (Scotti & White Reference Scotti and White2014) from which it follows that
\begin{equation} \left(\frac{\overline{\epsilon_{b,bulk}}}{\overline{{\epsilon}_{b,BL}}}\right) \sim \left(\frac{\overline{\epsilon_{u,bulk}}}{\overline{{\epsilon}_{u,BL}}}\right). \end{equation}
We now substitute the expressions for dissipation and the buoyancy variance in the different regions. For 
$\overline {\epsilon _{u,bulk}}$, 
$\overline {{\epsilon }_{u,BL}}$ and 
$\overline {{\epsilon }_{b,BL}}$, the values are expressed in terms of the Reynolds and Prandtl numbers in (5.1), (3.9) and (3.8), respectively. The buoyancy variance in the bulk is
\begin{equation} \overline{\epsilon_{b,bulk}} \sim \frac{U \varDelta^2}{h}\frac{h-\lambda_b}{h} \sim \frac{U \varDelta^2}{h} \sim \frac{\kappa\varDelta^2}{h^2} Pe. \end{equation}
Recalling that 
$\lambda _b/L\sim Nu^{-1} \sim {\rm Pe}^{1/2}$ and that 
$\lambda _u/H\sim Re^{-1/2}$, the balance in (5.5) therefore reduces to
\begin{equation} \frac{L^{{-}4} Re^3}{H^{{-}4}Re^{5/2}} \sim \frac{H^2 Pe}{h^2 Pe^{1/2}}. \end{equation}Rearranging the terms in the above balance leads to provides the following scaling
\begin{equation} \frac h H \sim (\varGamma^{{-}4} {Re^{1/2}}{Pe^{{-}1/2}})^{{-}1/2} \sim \varGamma^2 Pr^{1/4}. \end{equation}
This scaling is shown in figure 6 (squares). This decrease in bulk size is depicted in figure 7(a–c) where the plume depth behaves according to the above scaling. The reduction in 
$h$ also implies that in this transition regime between 
$I_l$ and 
$II_l$ or 
$I^+_u$ and 
$II_l$, (5.1) has to take into account the fact that now 
$Re$ is not based on 
$H$ but on 
$h$, the overturning depth. With respect to this, the overturning may be rescaled so that
\begin{equation} \overline{\epsilon_u} \sim \frac{\nu^{3}}{H^{4}}Re^3\left(\frac h H\right)^3, \end{equation}
to take into account the reduction of the bulk size 
$h$ as 
$Pr$ decreases. The above argument can also be expressed through the heat transport in the laminar BL where the bulk modification, similarly to (5.9) is now tied to the amount of heat transport such that
\begin{equation} U \frac {\rm \Delta} L \sim \frac{\kappa\varDelta}{\lambda_b^2}\left(\frac H h\right). \end{equation}
Substituting (5.8) into (5.9), the modified dissipation leads to a new dependence on the Prandtl number for the regime 
$II_l$ such that
$$\begin{gather} Re \sim Ra^{1/3}Pr^{{-}2/3}\left(\frac{h}{H}\right)^{{-}1} \sim Ra^{1/3}Pr^{{-}11/12}, \end{gather}$$
$$\begin{gather}Nu \sim Re^{1/2} Pr^{1/2}\left(\frac{h}{H}\right) \sim Re^{1/2} Pr^{3/4} , \end{gather}$$
which is confirmed by our DNS (see figures 4(d) and 5(b)). Combining (5.11a) and (5.11b) provides a correction for this 
$Pr$ transition in the 
$II_l$ regime
$$\begin{gather} Re\sim Ra^{1/3}Pr^{{-}11/12}, \end{gather}$$
$$\begin{gather}Nu\sim Ra^{1/6}Pr^{7/24}, \end{gather}$$
found for 
$Pr\lesssim 0.2$ (see figure 4b,d). Therefore, the advection–diffusion balance in the BL (3.7) is modified according to
\begin{equation} Nu \sim Pe^{1/2}Pr^{1/4} \sim Re^{1/2}Pr^{3/4} , \end{equation}
which now takes into account variable depth effects 
$h$, decoupled from the domain's depth 
$H$. This particular scaling is the last controlled by the laminar BL as 
$Ra$ increases. As Shishkina et al. (Reference Shishkina, Emran, Grossmann and Lohse2017) noted, further increasing the Rayleigh number may eventually lead to core-driven dynamics, as was previously observed by Griffiths & Gayen (Reference Griffiths and Gayen2015) for a somewhat different boundary condition. In the following subsection, we report the transition to a similar regime, which marks the shift to the limiting regime of horizontal turbulent convection for large 
$Ra$. In particular, we take advantage of the same scaling analysis for the variable size of the turbulent core with depth 
$h$ to recover the Prandtl number dependencies.
Time-averaged iso-contours of the streamfunction 
$\psi =[-0.1, -0.075, -0.05, -0.025, 0, 2, 4, 6, 8, 10, 12,$ 
$14, 16]/8\times 10^{-4}$ for 
$Ra=1.92\times 10^{15}$: (a) 
$Pr=1$, (b) 
$Pr=0.1$ and (c) 
$Pr=0.01$ showing the narrowing of the circulation as 
$Pr$ decreases and the circulation within the core narrowing beneath the differentially heated surface at 
$z\varGamma =1$.
5.2. The limiting turbulent BL regime at large Rayleigh numbers $IV_u$
The modification of the depth of the circulation at high Rayleigh numbers was first investigated by Griffiths & Gayen (Reference Griffiths and Gayen2015) and Rosevear et al. (Reference Rosevear, Gayen and Griffiths2017) who considered a periodic forcing at the surface and very small aspect ratios 
$\varGamma$. They showed that a laminar-type scaling for the dissipation in the bulk was responsible for the transition to a core-dominated turbulence regime. The scaling obtained from the vorticity equation is related to the previous analysis and considers a balance between the dissipation in the core (the interior) scaling as 
$Re\sim Ra^{1/2}Pr^{-1/2}$ and balances the dissipation in the BL, obtained from (3.7) BL. This is in contrast to the scaling observed in the previous regime (
$II_l$). This is shown in figures 4 and 5 where the dependencies of 
$Ra$ and 
$Pr$ do not support the above scaling. There is a clear departure from the laminar scaling obtained in (3.7) suggesting that turbulent BLs, characterised by a log-type profile and modified heat transfer coefficients (Grossmann & Lohse Reference Grossmann and Lohse2011), may be expected. Note that such BL profiles were already observed in Rosevear et al. (Reference Rosevear, Gayen and Griffiths2017) but the latter did not affect the heat- or the momentum-transfer scaling obtained in their analysis. Here, we report different results and show that log-type profiles do influence heat and momentum transfers, in a similar way to results obtained in RBC (Grossmann & Lohse Reference Grossmann and Lohse2011; van der Poel et al. Reference van der Poel, Ostilla-Mónico, Verzicco, Grossmann and Lohse2015).
5.2.1. Incompatibility with a laminar BL scaling
Low-Prandtl-number flows are particularly interesting with respect to the study of the transition to the limiting regime in HC since the transition to a turbulent-dominated flow is first observed in the 
$II_l$ regime where the dissipation of kinetic energy driving the dynamics scales as 
${\epsilon _u \sim \nu ^{3}L^{-4}Re^3}$.
Therefore, increasing 
$Ra$ may eventually lead to turbulent BLs as well. However, increasing 
$Ra$ leads to thinner BLs and a thinner bulk. Following this logic, once turbulence is triggered in the BL, the thermal BL becomes embedded into the momentum one and the BL profiles exhibit log-type profiles. As recently observed in Reiter & Shishkina (Reference Reiter and Shishkina2020), the velocity of the flow, which carries the temperature in the bulk, reduces from 
$U$ to 
$U(\lambda _{b}/\lambda _u)$ and the buoyancy variance dissipation rate (Shishkina et al. Reference Shishkina, Grossmann and Lohse2016) becomes
\begin{equation} \overline{\epsilon_{b,bulk}}\sim(\varGamma/2)\kappa\varDelta^2 h^{{-}2}{Pr Re^{3/2}Nu^{{-}1}}, \end{equation}
whereas the total volume 
$V$ has a buoyancy variance of
\begin{equation} \overline{\epsilon_{b,V}}=(\varGamma/2) \kappa\varDelta^2 L^{{-}2}Nu. \end{equation}Equating (5.14) and (5.15) gives the following expression for the Nusselt number
\begin{equation} Nu\sim \varGamma Re^{3/4}Pr^{1/2}\left(\frac h H\right). \end{equation}
With increasing 
$Ra$, the bulk dynamics are driven by the large-scale overturning flow whose length scale is 
$h$ as confirmed by our simulations in the previous subsection. In this case, the dissipation of kinetic energy in the bulk is essentially dependent on the large-scale velocity 
$U$ and the Reynolds number is again modified from 
$Re$ based on 
$H$ to 
$Re$ based on 
$h$ such that
\begin{equation} \overline{\epsilon_{u,bulk}} \sim \frac{\nu^{3}}{H^{4}}Re^3\left(\frac h H\right)^3, \end{equation}whereas the Rayleigh and Prandtl number dependencies of the bulk now yield
\begin{equation} \left(\frac{\overline{\epsilon_{b,bulk}}}{\overline{{\epsilon}_{b,BL}}}\right) \sim \left(\frac{\overline{\epsilon_{u,bulk}}}{\overline{{\epsilon}_{u,BL}}}\right) \equiv \frac h H \sim Re^{{-}1/8} Pr^{1/4}. \end{equation}
In the above, the buoyancy variance in the bulk still follows 
${\overline {\epsilon _{b,bulk}}\sim \kappa \varDelta ^2 h^{-2}{\rm Pe}}$, but in the BL, the buoyancy variance follows (5.16) and 
${\overline {{\epsilon }_{b,BL}}\sim \kappa \varDelta ^2 h^{-2}Re^{3/4}Pr^{1/2}}$ (see Reiter & Shishkina Reference Reiter and Shishkina2020). The turbulent kinetic energy dissipation is given by (5.17) and reads 
${\overline {\epsilon _{u,bulk}}\sim \nu ^3 L^{-4}Re^3}$ whereas in the BL, dissipation is bounded by the stably stratified layer 
${\overline {\epsilon _{u,BL}}\sim \nu ^3 L^{-4}Re^{5/2}}$. Combining (5.14) and (5.17), together with the bulk reduction effect 
$(h/H)$ for the Nusselt number as in (5.11b), the relation for the Nusselt number now reads
\begin{equation} Nu\sim Re^{3/4}Pr^{1/2}\left(\frac h H\right) \sim Re^{5/8}Pr^{3/4}, \end{equation}which agrees with the results shown in figure 5(a,b). Combining (5.14), (3.4), (5.17a) and (5.17b), one obtains
$$\begin{gather} Re \sim Ra^{8/21}Pr^{{-}22/21}, \end{gather}$$
$$\begin{gather}Nu \sim Ra^{5/21}Pr^{17/96}, \end{gather}$$
which slightly overestimates the 
$Ra$ number dependence for the Nusselt number with respect to 
$Ra$ (i.e. 
${\approx }0.225$ from the DNS vs 
${\approx }0.238$ from the theory) but clearly underestimates the Prandtl number dependence (
${\approx }0.41$ for the DNS vs 
${\approx }0.17$ for the theory) in that particular region of the 
$(Ra,Pr)$ plane and is denoted as 
$IV_u$ in figures 3(a,b) and 5(a,b).
This scaling analysis shows that the effect of turbulence must play a role. In particular, the dissipation scaling has to take into account the turbulent BL characteristics and log-type corrections have to be eventually reintroduced to accurately predict both the Prandtl and Reynolds number dependencies observed from the simulations.
5.2.2. A fully turbulent BL scaling
The above scaling overestimates the heat flux, since the observed exponent 
$Nu\sim Ra^{0.225}$ is close to the exponent 
$Nu\sim Ra^{5/21}\sim Ra^{0.238}$ but smaller and the turbulent 
$Nu\sim Ra^{1/4}$ scaling (figure 4a) and suggests a new correction for the dissipation in what may appear as turbulent thermal and momentum BLs. In addition, both Prandtl number dependencies found in the numerical simulations are not predicted by (5.20a) and (5.20b) which suggests that turbulence plays a non-trivial role in the above scaling. For small 
$Pr$ and large 
$Ra$, the statically unstable BLs become turbulent and, for decreasing 
$Pr$, the Reynolds number increases (figure 4c) which causes the BL to transition to turbulence. The dissipation of kinetic energy may be split between a viscous sublayer 
$\overline {\epsilon _{vs}}$ and a logarithmic layer 
$\overline {\epsilon _{ll}}$ such that 
$\overline {\epsilon _{u,BL}}=\overline {\epsilon _{vs}}+\overline {\epsilon _{ll}}$ (Grossmann & Lohse Reference Grossmann and Lohse2011). In the log layer, the dissipation reads
\begin{equation} {{\epsilon_{u,ll}}(z) = \frac{u_*^3}{C_{\kappa,u} z}}, \end{equation}
where 
${C_{\kappa,u}\approx 0.4}$ is the Von Kármán constant. The last equation may also be rearranged as
\begin{equation} {{\epsilon_{u,ll}}(z) = \nu^3 L^{{-}3} \frac{Re^3}{C_{\kappa,u} z} \left(\frac{u_*}{U}\right)^3}. \end{equation}The mean kinetic energy dissipation in the log layer can therefore be obtained by integrating the above equation so that
\begin{equation} \overline{\epsilon_{u,ll}} = \frac{1}{h/2}\int_{z^*}^{h/2}\epsilon_{u,ll}(z)\,\mbox{d}z, \end{equation}
where 
$z_* = \nu /u_*$ and 
$h/2$ corresponds to the outer edge of the logarithmic zone. The dissipation in the log layer thus acts as a buffer to heat exchanges and induces a log-type correction denoted as 
$\mathcal {L}({\cdot })$, which depends on the Reynolds number
\begin{equation} \overline{\epsilon_{u,ll}} := \nu^{3}L^{{-}4}Re^3\left(\frac{u_*}{U}\right)^3 \left(\frac{h}{H}\right)^3 \frac{2}{C_{\kappa,u}}\log\left(Re\frac{u_*}{U} \frac{1}{2}\frac{h}{H}\right), \end{equation}
where the length scale is now 
$L$ (i.e. the length of the domain), and 
$u_*^2=\overline {u'w'}$ is the typical velocity fluctuation scale (Grossmann & Lohse Reference Grossmann and Lohse2011). In fact, in the present simulations, the friction directly on the wall is zero because of the free-slip boundary condition. Therefore, the friction velocity 
$u_*$ does not refer to the wall shear stress but to the turbulent shear stresses, induced by the turbulent fluctuations of the buoyancy flux 
$\overline {w'b'}$ originating from the statically unstable buoyancy profile in this region near the wall. The logarithmic profiles for both velocity and buoyancy are shown in figure 8(a,b). For all turbulent profiles, a logarithmic region can be observed for the velocity profile. In addition, the profile seems to be self-similar, at least for the three profiles reported in the figure. The source of turbulent stresses here is suggested by the large logarithmic profile, measured for two decades in figure 8(b) for the buoyancy profile at 
$Ra = 1.92\times 10^{15}$ and 
$Pr=10^{-2}$. The buoyancy variance in the BL can also be expressed using a similar analogy. As shown in figure 8(b), the thermal layer shows a log-type layer, which is a common feature of turbulent statically stable and unstable BLs
\begin{equation} {\epsilon_{b,ll}}(z) = {Pr_t}\frac{\kappa b_*^2}{{C_{\kappa,b}} z^2}, \end{equation}
where 
$b_*= \overline {w'b'}/u_*$. The fluctuating velocity 
$u_*$ can be connected to the outer velocity 
$U$ or 
$Re$, and similarly for the buoyancy variance via
\begin{equation} \frac{u_*}{U} = \frac{{C_{\kappa,u}}}{\log\left(Re\dfrac{u_*}{U} \dfrac{1}{\theta}\right)} \quad \mbox{and} \quad \frac{b_*}{\varDelta} = {\frac{1}{Pr_t}} \frac{{C_{\kappa,b}}}{\log \left(Re\dfrac{u_*}{U}\dfrac{1}{\theta}\right)}. \end{equation}(a) Mean turbulent BL profiles and (b) mean turbulent buoyancy profiles measured at 
$x=-0.325$, rescaled using the slip velocity at the wall 
$u_0=u(x=-0.75,z=H)$ and the maximum velocity at this particular 
$x$ location. Note that we are not in the presence of a free-slip-type turbulent BL but we recover a log-type layer as shown by the dashed line for the highest (lowest) values of 
$Ra$ (
$Pr$) and, thus, the highest 
$Re$. See the main text for a discussion of the origin of the log layer.
In the above, we assumed that the flux Richardson number 
$R_f=\overline {w'b'}/(\overline {u'w'} \partial u/\partial z)$ is constant, which should hold provided that the structure of the turbulent BL remains self-similar with 
$Ra$ and 
$Pr$ (see figure 9a,b). In addition, we assume that the thermal Kolmogorov constant 
$C_{\kappa,b}\approx C_{\kappa,u}\approx 0.4$ which is consistent with the Monin–Obukhov BL theory (Landau & Lifschitz Reference Landau and Lifschitz1987). Further, the turbulent Prandtl number 
$Pr_t\approx {O}(1)$ which is consistent with finite values of the flux Richardson number and the value of the mixing efficiency reported in HC at large Rayleigh numbers (Scotti & White Reference Scotti and White2011; Gayen et al. Reference Gayen, Griffiths and Hughes2014). The empirical constant 
$\theta$ depends on the geometry of the system, along a plate 
$\theta$ is empirically found to be equal to 
$0.13$ (see Grossmann & Lohse Reference Grossmann and Lohse2011). Again, the mean buoyancy variance can be integrated so that
\begin{equation} \overline{\epsilon_{b,ll}} = \frac{1}{h/2}\int_{z^*}^{h/2}\epsilon_{b,ll}(z)\,\mbox{d}z, \end{equation}which reduces to
\begin{equation} \overline{\epsilon_{b,ll}} = {Pr_t}\frac{2 \kappa b^2_*}{{C_{\kappa,b}} L^2} Re\frac{u_*}{U}\frac{h}{H}, \end{equation}and using (5.26a,b), the expression becomes
\begin{equation} \overline{\epsilon_{b,ll}} = \kappa \varDelta^2 L^{{-}2} Re\left(\frac{u_*}{U}\right) \left(\frac{h}{H}\right)\frac{2{Pr_t}}{C_{\kappa,b} }{\log\left(Re\frac{u_*}{U} \frac{1}{\theta}\frac{h}{H}\right)^{{-}2}}. \end{equation}Dependencies of 
$Nu$ (red) and 
$Re$ (blue) with respect to 
$Ra$ (a) and 
$Pr$ (b) showing the modification and weak variations of both 
$Nu$ and 
$Re$ in the regime considered here. The continuous lines show the predictions from (5.34a,b) whereas the dashed line shows the measurements from figure 4(a,d).
In the above expressions, the unknown ratio 
$u_*/U$ can be computed using Lambert's 
$W$-function where 
$u_*/U=C_\kappa /W(Re C_\kappa /\theta )$ (Grossmann & Lohse Reference Grossmann and Lohse2011). The dissipation in the turbulent regime is thus modified from (5.17a) and becomes
\begin{equation} {\overline{\epsilon_{u,ll}} \sim \nu^{3}L^{{-}4}Re^3\left(\frac{h}{H}\right)^3{\mathcal{L}(Re)}}, \end{equation}
where 
${\mathcal {L}(Re)}$ is given by
\begin{equation} \mathcal{L}(Re)\equiv C_{\kappa,u} \left(\frac{u_*}{U}\right)^3\log \left(Re\frac{u_*}{U}\frac{1}{\theta}\frac{h}{H} \right). \end{equation}
However, the bulk still dissipates at a rate expressed in (5.20a) above (see figure 4c,d). We can think of this correction as a decrease in heat transfer through the BL, which is responsible for an even faster modification of the relative depth of the recirculation region 
$(h/H)$. The latter is estimated using (5.18) which writes
\begin{equation} \left(\frac{\overline{\epsilon_{b,bulk}}}{\overline{{\epsilon}_{b,BL}}}\right) \sim \left(\frac{\overline{\epsilon_{u,bulk}}}{\overline{{\epsilon}_{u,BL}}}\right)\equiv \frac h H \sim Re^{{-}1/8} Pr^{1/4}. \end{equation}Thus, the scaling for the Nusselt number becomes
\begin{equation} Nu=Re^{3/4}Pr^{1/2}\left(\frac h H\right)\sim Re^{5/8}Pr^{3/4}\mathcal{L}(Re)^{1/2}, \end{equation}and from (5.17b), (5.14), (3.4) and (5.30) it follows that
$$\begin{gather} Re \sim Ra^{8/21}Pr^{{-}22/21}\mathcal{L}(Re)^{{-}8/21}, \end{gather}$$
$$\begin{gather}Nu \sim Ra^{5/21}Pr^{17/96}\mathcal{L}(Re)^{{-}5/27}. \end{gather}$$ These scalings are verified in figure 9(a,b) for both the Nusselt number and the Reynolds number with respect to both 
$Ra$. The Prandtl number dependence prediction is also improved compared with (5.17a,b). The Prandtl number exponent for the Reynolds number is found at 
$Re \sim Pr^{-1}$ for the DNS whereas the log-corrected exponent is 
$Re \sim Pr^{-1.075}$. The Nusselt number dependence is 
$Nu \sim Pr^{0.41}$ in the DNS while the log-corrected scaling provides 
$Nu \sim Pr^{0.2}$ which hints at a possible Prandtl number dependence due to plume dynamics, as observed in RBC (Grossmann & Lohse Reference Grossmann and Lohse2011; Ni, Huang & Xia Reference Ni, Huang and Xia2011).
To the best of the authors’ knowledge, this is the first time such log-type corrections have been applied and verified for the Prandtl number dependence. The logarithmic region of the time-averaged turbulent BLs is shown in figure 8 for two different Prandtl numbers and two different Rayleigh numbers in the turbulent regime, supporting our assumption and analysis. Thus, these scaling laws provide evidence for a new regime in HC, which can be considered as the limiting regime in HC. It is worth noting that as 
$Ra$ further increases, the 
$\alpha$ exponent progressively reaches the value of 
$5/21$ where log corrections become less significant. Note that this exponent also agrees very well with the recent study (Reiter & Shishkina Reference Reiter and Shishkina2020). The transition from the 
$II_l$ regime to the 
$IV_u$ regime in figure 3 is marked by the vertical black dotted line and is obtained by matching the Reynolds number between each region. After equating (5.34a) with (5.12), the transition is found very close to a constant at 
$Ra \approx 3\times 10^{12}$.
The next subsection investigates whether this last transition can be thought of in terms of turbulent regimes. As one expects to see a transition to the limiting regime of convection, the question therefore arises whether turbulence in the present flow is dependent or not on viscosity and if we are reaching the ultimate regime where dissipation is not dependent on viscosity. Note that a viscous-independent turbulent regime would make the present result relevant for large-scale geophysical applications.
6. ‘Ultimate’ or ‘classical’ turbulence?
The picture drawn in the previous section provides an understanding of the heat and momentum transition leading to the ultimate regime of HC. Although the above picture is convincing, with successive scaling transitions in agreement with what was observed previously in RBC, it remains complex, with two control parameters and at least five different regimes spanning ten orders of magnitude in 
$Ra$.
We propose to re-analyse our data in the framework of Kolmogorov's zeroth law of turbulence (Frisch Reference Frisch1995), and provide a valuable tool to diagnose the state of turbulence observed and how these regimes and their transitions may be further analysed.
We introduce the Kolmogorov number, obtained by rescaling the PY constraint on dissipation using 
$U$ and 
$L$, giving
\begin{equation} \frac{\epsilon L}{U^3}\equiv Ko = Re^{{-}3} Ra Pr^{{-}2}. \end{equation}
The scaling law relating 
$Ko$ with respect to 
$Re$ (or 
$Pe$) provides a new way to analyse whether the flow is laminar, transitional or in a state of ‘classical’ or ‘ultimate’ turbulence. In the laminar case, the Kolmogorov number and, hence, the dissipation is solely controlled by the vertical gradient of the velocity where 
$Ko\sim Re^{-1}$. In contrast, ultimate turbulence achieves complete similarity for parameters that contain viscosity (Vassilicos Reference Vassilicos2015), that is, we should expect 
$Ko={\rm Const}$. An intermediate regime 
$Ko\sim Re^{-1/2}$ can also be expected if the BLs dominate dissipation, as the latter may not achieve complete similarity.
The evolution of 
$Ko$ with respect to both 
$Re$ and 
$Pe$ is shown in figure 10, where for laminar flows we obtain, for all 
$Pr$, 
$Ko\sim Re^{-1}$. Across the 
$I_l$ and 
$I^*_l$ regimes, the dependence of 
$Ko$ exhibits a 
$Re^{1}$ transition, where dissipation is enhanced throughout this core-driven mixing regime. A 
$Ko\sim Re^{-1/2}$ type regime is then recovered for all the values of 
$Pr$ in the 
$IV_u$ regime (see figure 10a). The same observation can be made for the scaling concerning 
$Pe$ where the same conclusions arise (see figure 10b). Variations with respect to the Prandtl number follow the same rationale. Transitions between regimes varying the Prandtl number are found at constant 
$Ko$ (see figure 10c), whereas in the 
$II_l$ and 
$II^*_l$ regimes, variations with respect to 
$Pr$ occur for 
$Ko\sim Re^{-1}$. At higher 
$Ra$ and in the 
$IV_u$ regime, conclusions are difficult to draw, but we may expect 
$Ko\sim Re^{-1/2}$. The transition to the ‘ultimate’ regime of turbulent convection in natural HC does not seem to be attainable even at extremely high 
$Ra$. Ultimately, the bound on the Richardson number under the stable (warming) boundary proves to be what prevents HC in the 
$Pr=O(1)$ or lower regimes from achieving the ultimate regime. Adding active forcing in this same region or, for instance, through radiative heat transfer (Lepot, Aumaître & Gallet Reference Lepot, Aumaître and Gallet2018) may allow an early transition to the ultimate regime to be observed. Such strategies may allow the realisation of the 
$IV_l$ regime predicted by Siggers et al. (Reference Siggers, Kerswell and Balmforth2004) and Rocha et al. (Reference Rocha, Bossy, Llewellyn-Smith and Young2020b).
(a) 
$Re$ dependencies and (b) 
$Pe$ dependencies of the Kolmogorov number 
$Ko$ for variations with respect to 
$Ra$. Same for (c,d) but for variations with respect to 
$Pr$ (refer to figure 3 for colour code).
7. Conclusions
We have reported evidence of two new turbulent regimes in HC based on scaling arguments at low Prandtl numbers. More precisely, we first highlighted the regimes that are known as limiting regimes. For asymptotically small Prandtl numbers, we have a regime where the core is driven by turbulence but where the BL remains laminar and we name this regime 
$II_l$ following the nomenclature of Shishkina et al. (Reference Shishkina, Grossmann and Lohse2016). The second regime is characterised by both a turbulent core and a transition within the BL to a state which has characteristics proper of turbulence even though lacking a clearly discernible log layer. It is possible that this embryonic logarithmic layer will, at higher Rayleigh numbers, expand to form a canonical layer, but at present we cannot state this with confidence. Following SGL's nomenclature, we call this the 
$IV_u$ regime.
Our results support and integrate previous evidence from Shishkina & Wagner (Reference Shishkina and Wagner2016) and the model of Hughes et al. (Reference Hughes, Griffiths, Mullarney and Peterson2007) in the SGL theory of HC.
In the 
$II_l$ regime, we have observed a new scaling, where the modification of the turbulent bulk size modifies the dependence of the Prandtl number for both the Reynolds and Nusselt numbers scaling. This reduction in the size of the bulk has been found to be essentially Prandtl number dependent where the bulk decreases in size when 
$Pr$ decreases.
The transition to the turbulent limiting regime, denoted as 
$IV_u$, has also been observed. In this particular regime, the flow in the core becomes turbulent, whereas the statistics of the geometric invariants in the BL assume turbulent characteristics, although a fully developed logarithmic is not observed. Similar to the work Rosevear et al. (Reference Rosevear, Gayen and Griffiths2017), the turbulent region is essentially located below the horizontal forcing and progressively clusters beneath as both 
$Ra$ increases and 
$Pr$ decreases. According to Shishkina et al. (Reference Shishkina, Grossmann and Lohse2016), this last regime marks the final transition at large 
$Ra$. Log corrections allow recovery of the correct dependencies 
$Nu$ and 
$Re$ with respect to 
$Ra$ and improved estimates for 
$Pr$. However, more work is needed to understand the exact dependence of the Nusselt number on the Prandtl number, which may be attributed to plume dynamics (Grossmann & Lohse Reference Grossmann and Lohse2011; Ni et al. Reference Ni, Huang and Xia2011).
It is important to stress that our calculations were performed with free-slip conditions. We confirm that the previously observed scaling laws in the region 
$I^*_l, I_l$ and 
$I^+_u$, which were reported under no-slip conditions, are observed even with free-slip conditions, albeit with shifted boundaries. Whether the newly reported regimes are robust with respect to the boundary conditions on velocity remains an open question. We also propose a new analysis, based on the Kolmogorov number 
$Ko$, a rescaled dissipation rate. The analysis confirms that the flow transitions from laminar to classical turbulence, but also shows that the flow never transitions to ultimate turbulence, which would be akin to the 
$IV_l$ regime of Siggers et al. (Reference Siggers, Kerswell and Balmforth2004), Shishkina et al. (Reference Shishkina, Grossmann and Lohse2016) and Rocha et al. (Reference Rocha, Bossy, Llewellyn-Smith and Young2020b).
The ultimate regime 
$IV_l$, if it exists, has thus yet to be observed. It is therefore of particular interest to study new types of HC where the turbulence can be strong enough to get rid of the effect of BLs and trigger purely inertial, core-driven, turbulent HC regimes. A such regime would be of particular importance for geophysical applications such as the overturning circulation. In particular, the present study further points to the idea that a deep circulation cannot be sustained in the limit of large Rayleigh numbers and that other effects must be taken into account to provide the necessary physics to model the meridional overturning circulation.
The companion paper gathers the results from this study together with an experimental study at a large Prandtl number. In particular, a regime diagram is provided that highlights all known limiting regimes of HC.
Acknowledgements
The authors thank Dr B. White, Professor D. Lohse and Dr O. Shiskina and three anonymous reviewers for insightful comments and encouragement.
Funding
This work was supported by the National Science Foundation (grant nos. OCE–1155558 and OCE-1736989).
Declaration of interests
The authors report no conflict of interest.
Appendix. Turbulent BL analysis
Although the scaling analysis in § 5.2.2 shows good agreement with the simulations, we provide further information regarding the velocity profiles and the local scaling laws of the BL in the large-Rayleigh-number regime.
We first depict near-wall BL profiles in the unstable region. These mean-velocity profiles were calculated for the largest Rayleigh number at 
$Ra = 1.92 \times 10^{15}$ and a Prandtl number 
$Pr=10^{-2}$ which is the most turbulent cases from this study. Horizontal velocity profiles are shown as 
$U^+(z) = (U(z^+)-U(z=0))/u^*$ where 
$u^*=\sqrt {\nu (\partial U/\partial z)}$ is obtained from the maximum shear stress obtained from the maximum shear in the laminar region near the wall and where 
$z^+=\nu /u^{*}$. The profiles are compared with the laminar viscous sublayer solution 
$U^+=y^+$ and the turbulent profile 
$U^+=C_\kappa ^{-1}\log (y^+)+5$. The results show that the resolution of the grid in the BL region for the first mesh size is below 
$y^+=1$. This was verified for all cases with large Rayleigh numbers. At 
$y^+>100$, the flow starts to follow the logarithmic flow region but at 
$y^+>300$, the BL is modified by the return flow and no longer follows the profile of turbulent BLs with zero pressure gradient. There also exists a peak at 
$y^+\approx 25$ 25 % greater than the beginning of the log profile that begins at 
$y^+\approx 100$. This feature is reminiscent of turbulent BLs (Monkewitz, Chauhan & Nagib Reference Monkewitz, Chauhan and Nagib2007) which might be associated with the onset of small-scale mixing caused by the statically unstable buoyancy profile.
In addition, we report the evolution of the displacement thickness 
$\delta ^*(x)$ and the momentum thickness 
$\varTheta ^*(x)$ compensated with 
$x^{-4/5}$, which is also a common feature of the turbulent BL scaling laws (see Schlichting & Gersten (Reference Schlichting and Gersten2016), for a detailed review). Figure 11(b) shows a compensated plot of both compensated quantities for the velocity profiles in figure 11(b) where the scaling laws vary by approximately 
$10\,\%$ around the mean value. The shape factor 
$\mathcal {H}_{1,2} = \delta ^*/\varTheta ^*$ is shown in figure 11(c) for increasing Rayleigh numbers at 
$Pr=10^{-1}$. It is shown that when the BL profiles become turbulent, the shape factor decreases to a value of 
$\mathcal {H}_{1,2} \approx 1.3$, also in agreement with zero-pressure gradient and homogeneous turbulent BL theory (Schlichting & Gersten Reference Schlichting and Gersten2016).
(a) Streamwise evolution of the turbulent BL profile for 
$Ra=1.92\times 10^{15}$ and 
$Pr=10^{-2}$. (b) Compensated displacement and momentum thickness of the turbulent BL measurements are shown in (a). (c) BL shape factor 
$\mathcal {H}_{1,2}$ along the streamwise direction at 
$Pr=10^{-1}$ for different values of 
$Ra$. (d) PDF of the cosine of the angle between vorticity 
$\omega$ and the vorticity stretching vector 
$W$.
Further confirmation that the BL is indeed fully turbulent comes from the analysis of the statistics of geometric invariants of the velocity gradient tensor (Meneveau Reference Meneveau2011). One such invariant is 
$(\omega {\,\cdot\, } W) /\|\omega \|\| W \|$, where 
$W_i \equiv \omega _j S_{j i}$ is the vortex stretching vector. The PDF calculated from our experiments agrees very well with the PDF of grid turbulence at 
$Re_\lambda =70$ reported by Tsinober (Reference Tsinober1998). Interestingly, this PDF remains essentially Rayleigh and the Prandtl number independent for sufficiently large values of the Rayleigh number. All of the above results point to the fact that, despite the free-slip boundary condition used for the simulations, the flow displays turbulent BL features. The turbulent scaling analysis essentially allows for correcting the transition from laminar to this weakly turbulent regime. This conclusion is also in line with the classical turbulence scaling derived in § 6.
Open access
