1. Introduction
Buoyancy-driven flows are ubiquitous in nature and industrial processes. Geothermal heating generates plumes in Earth's mantle that drive plate tectonics (Jaupart, Labrosse & Mareschal Reference Jaupart, Labrosse and Mareschal2007), variations of solar irradiance with latitude drive large-scale winds in Earth's atmosphere (Trenberth & Caron Reference Trenberth and Caron2001), brine rejection from sea ice increases the density of near-surface water masses that contribute to the large-scale ocean circulation (Jacobs Reference Jacobs2004; Abernathey et al. Reference Abernathey, Cerovecki, Holland, Newsom, Mazloff and Talley2016) and metal assembly or coating through laser-induced heat deposition involve multi-component phase changes and fluid convection that can affect material microstructure (Gan et al. Reference Gan, Yu, He and Li2017).
Rayleigh–Bénard (RB) convection and horizontal convection (HC) are two canonical buoyancy-driven flow configurations that have attracted significant interest. RB convection considers fluid motions between horizontal plates held at different temperatures, such that the diffusive base state is gravitationally unstable (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009), whereas HC considers an inhomogeneous heat flux or temperature distribution along a single horizontal boundary, which is baroclinically unstable and drive a vigorous boundary-layer flow (Hughes & Griffiths Reference Hughes and Griffiths2008). The behaviour of RB convection and HC is typically examined through the scaling of the Reynolds number 
$Re$, which is a proxy for fluid velocity, and Nusselt number 
$Nu$, which is a proxy for the efficiency of heat transport by convection, with the control parameters, including notably, the Rayleigh number 
$Ra$, which measures the available potential energy relative to dissipation processes. Research over the past few decades has seen the development of detailed phase diagrams for the flow regimes and scalings of 
$Re$ and 
$Nu$ as functions of 
$Ra$, as well as the Prandtl number 
$Pr$, which compares viscosity to thermal diffusivity, for both RB convection (Grossmann & Lohse Reference Grossmann and Lohse2000; Ahlers et al. Reference Ahlers, Grossmann and Lohse2009) and HC (Mullarney, Griffiths & Hughes Reference Mullarney, Griffiths and Hughes2004; Shishkina, Grossmann & Lohse Reference Shishkina, Grossmann and Lohse2016). However, many open questions remain. For instance, recent studies on RB convection aim to address the existence of the ultimate regime (Zhu et al. Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018) or its large-scale organisation (Pandey, Scheel & Schumacher Reference Pandey, Scheel and Schumacher2018; Wang et al. Reference Wang, Verzicco, Lohse and Shishkina2020; Vieweg, Scheel & Schumacher Reference Vieweg, Scheel and Schumacher2021), whereas active topics of research in HC include the emergence and properties of turbulence, which is spatially heterogeneous, the definition of a thermodynamically compelling Nusselt number, which is not as straightforward as in RB convection (Paparella & Young Reference Paparella and Young2002; Scotti & White Reference Scotti and White2011; Gayen, Griffiths & Hughes Reference Gayen, Griffiths and Hughes2014; Passaggia, Scotti & White Reference Passaggia, Scotti and White2017; Rocha et al. Reference Rocha, Constantinou, Llewellyn Smith and Young2020b) and the effect of rotation (Vreugdenhil, Gayen & Griffiths Reference Vreugdenhil, Gayen and Griffiths2019; Gayen & Griffiths Reference Gayen and Griffiths2022).
The fluid dynamics resulting from the superposition of RB convection with HC has received surprisingly limited attention, in spite of being of fundamental interest and having potentially important applications in the environment. Notably, RB convection and HC may concomitantly control the fluid dynamics within subglacial lakes in Greenland and Antarctica, which affect the dynamics of ice sheets and most likely host extremophiles of interest to astrobiology (Cockell et al. Reference Cockell2011; Livingstone et al. Reference Livingstone2022). Subglacial lakes, which are typically fresh except when active or close to grounding lines (Priscu et al. Reference Priscu2021), are exposed to geothermal heating as well as horizontal temperature gradients along the ice–water interface when the ice thickness above is spatially variable. Thicker ice produces larger pressure at the ice–water interface and thus lower interface temperature, because the freezing (or fusion) temperature of water decreases with pressure (Thoma et al. Reference Thoma, Grosfeld, Smith and Mayer2010). As recently shown (Ulloa et al. Reference Ulloa, Ramón, Doda, Wüest and Bouffard2022), sloped surface waterbodies can similarly experience both RB convection and HC when cooling, as heat losses from the surface vary with water depth. The fluid dynamics of planetary oceans may also be affected by both RB convection and HC, as oceans receive solar radiation that varies with latitude and typically experience geothermal heating (Wang et al. Reference Wang, Huang, Zhou and Xia2016). Another field where thermal convection might be affected by boundary inhomogeneities in a direction transverse to gravity is the Earth's liquid outer core. Heterogeneous heat fluxes along the core–mantle boundary, which are due to large-scale convective patterns within the mantle, can sustain large-scale azimuthal flows and affect heat fluxes across the fluid layer (Sumita & Olson Reference Sumita and Olson1999; Mound & Davies Reference Mound and Davies2017).
Past studies of the dual RB–horizontal (RBH) configuration include a handful of simulations and experiments relevant to subglacial lakes and open oceans. RBH dynamics underly a series of realistic numerical simulations of subglacial lakes that used a large-scale ocean code with parameterised subgrid-scale processes (Thoma, Grosfeld & Mayer Reference Thoma, Grosfeld and Mayer2007; Thoma et al. Reference Thoma, Grosfeld, Filina and Mayer2009) and a laboratory experiment of lake Vostok (Wells & Wettlaufer Reference Wells and Wettlaufer2008). There is yet no consensus on the type of fluid motions expected in subglacial lakes: Thoma et al. (Reference Thoma, Grosfeld and Mayer2007, Reference Thoma, Grosfeld, Filina and Mayer2009) predicted HC-driven large-scale circulations in lakes Vostok and Concordia affected by rotation; Wells & Wettlaufer (Reference Wells and Wettlaufer2008) found that the dynamics in lake Vostok is better described by rotating RB convection, i.e. dominated by multiple columnar vortices; and Couston & Siegert (Reference Couston and Siegert2021) predicted non-rotating RB convection in most subglacial lakes. In open ocean research, several studies (Hofmann & Maqueda Reference Hofmann and Maqueda2009) have shown using general circulation models (GCMs) that the geothermal flux affects the global ocean circulation, which is otherwise primarily driven by winds and heat fluxes at the air-sea interface. The impact of geothermal heating on the meridional overturning circulation (MOC) of the Atlantic Ocean has also been investigated through idealised numerical simulations (Mullarney, Griffiths & Hughes Reference Mullarney, Griffiths and Hughes2006) and laboratory experiments (Wang et al. Reference Wang, Huang, Zhou and Xia2016), wherein the MOC is driven by a horizontally varying source of buoyancy or buoyancy flux. Both studies focused on dynamical regimes dominated by HC but demonstrated significant effects of bottom heating: in a box with aspect ratio 
$\varGamma =6$, Mullarney et al. (Reference Mullarney, Griffiths and Hughes2006) found that the volume flux driven by HC is 145 % higher with a bottom heat flux equal to just 10 % the amount of heat extracted from the top boundary, while in a box with aspect ratio 
$\varGamma =1$, Wang et al. (Reference Wang, Huang, Zhou and Xia2016) found a 260 % increase of the volume flux with bottom heating equal to 6.8 % the amount of heat input through half of the top boundary.
In this paper, our goal is to identify and characterise the transition from RB convection to HC as a function of the control parameters, including, most importantly, the ratio of the imposed horizontal temperature gradient 
$\lambda$ along the top boundary multiplied by the thermal conductivity 
$k$ (which yields a horizontal heat flux), to the bottom heat flux 
$F$, which we write as 
$\varLambda =k\lambda /F$. To this end, we run a large number of two-dimensional numerical simulations with variable bottom and surface forcing (listed in table 2 in Appendix A), and diagnose the resulting dynamics through the Reynolds and Nusselt numbers, as well as the characteristic length scale of overturning motions (see §§ 3.1–3.4). We also run a set of simulations over very long time scales (tens of diffusive time scales) near the transition, in order to demonstrate that RBH convection is multi-stable for some parameters (§ 3.5). The article is organised as follows. In § 2 we introduce the dimensional and dimensionless governing equations as well as the numerical method. In § 3 we show results highlighting the regime transition between RB-like convection at 
$\varLambda \ll 10^{-2}$ and HC at 
$\varLambda \gg 10^{-2}$, and we demonstrate the existence of multiple flow states for the same set of problem parameters. Finally, we conclude in § 4.
2. Problem formulation
We consider a two-dimensional rectangular fluid domain with Cartesian coordinates 
$(x,z)$ centred on the bottom boundary; 
$\boldsymbol {e}_z$ is the upward-pointing unit vector of the 
$z$-axis, which is opposite to gravity (figure 1a). We denote by 
$H$ and 
$L$ the fluid depth and domain length. We consider a pure fluid (similar to fresh water), i.e. without dissolved salts. The evolution of the fluid velocity 
${\boldsymbol {u}}$, pressure 
$p$ and temperature 
$T$ are governed by the Navier–Stokes equations in the Boussinesq approximation, i.e.
\begin{gather}
\partial_t {\boldsymbol{u}} - \nu{\nabla^2}{\boldsymbol{u}}
+ \boldsymbol{\nabla} (p/\rho_0) ={-}
\left({\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}\right){\boldsymbol{u}}
- (\rho'/\rho_0) g\boldsymbol{e}_z,
\end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{u}} = 0, \end{gather}
\begin{gather}\partial_t T - \kappa {\nabla^2} T ={-} \left({\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}\right) T, \end{gather}
where 
$\rho _0$ is the reference fluid density and 
$\rho '$ is the density anomaly; 
$g$ is gravitational acceleration, 
$\partial _t$ denotes time derivative and 
$\boldsymbol {\nabla }$ is the gradient operator. We consider constant thermodynamic and transport parameters, which is known as the Oberbeck approximation, such that the equation of state is simply 
$\rho '=-\rho _0\alpha (T-T_0)$ with 
$T_0$ the reference temperature and with the thermal expansion coefficient 
$\alpha$, along with kinematic viscosity 
$\nu$ and thermal diffusivity 
$\kappa$, taken constant. All boundaries are no slip. We impose a uniform heat flux on the bottom boundary, a variable temperature profile along the top plate and adiabatic side walls, such that the boundary conditions read
\begin{gather}
{\boldsymbol{u}}(z=0)={\boldsymbol{u}}(z=H)={\boldsymbol{u}}\left(
x={\pm} L/2\right)=\boldsymbol{0},
\end{gather}
\begin{gather}\partial_zT(z=0)=-\dfrac{F}{k},
\quad T(z=H)=T_f(x)=T_0+\dfrac{\lambda
L}{2}\sin\left(\dfrac{{\rm \pi} x}{L}\right), \quad
\partial_xT\left( x=\pm L/2\right)=0.
\end{gather}
It may be noted that the sinusoidal profile of temperature imposed on the top boundary differs from the step-wise or linear profile used in previous studies of HC, the latter being consistent at leading order with a constant-tilt ice–water interface. We chose a sinusoidal profile because it is consistent with the adiabatic vertical walls on the sides. We use subscript 
$f$ to denote the top temperature 
$T_f$, because the top of the water column must be at the freezing temperature in a subglacial lake. For simplicity, we consider a flat horizontal top boundary. Note that we often refer to the heat flux imposed on the bottom boundary as the geothermal flux as our study is motivated by subglacial lakes.
(a) Problem schematic. A rectangular fluid domain is subject to bottom heating and a horizontal temperature gradient along the top boundary. The former generates warm rising plumes contributing to multiple overturning cells (thin arrows), whereas the latter drives a single cell with intense down-welling near and below the cold end of the top boundary (thick arrow). In a subglacial lake, a tilted ice–water interface (shown by the dotted line), due to variable ice thickness above the lake water, would drive HC because of the pressure-dependence of the freezing temperature. (b) Phase diagram of this study. The circles show the location of the simulations in 
$(\varLambda,Ra_F)$ space for 
$\varGamma =8$; for 
$\varGamma =4$ (respectively, 
$\varGamma =16$) the markers are diamonds (squares) and are slightly shifted downward (upward) with respect to 
$Ra_F=10^8$ (
$Ra_F=10^7$) for better visibility. The stars show pure HC simulations run with 
$\varLambda =1$ (with markers shifted to the right for better visibility). The purple and grey shadings highlight regions of the parameter space wherein the regime dynamics is similar to RB convection and HC, respectively.
In order to identify the minimum number of independent parameters and explore their effect on the fluid dynamics, we non-dimensionalise the governing equations (2.1) and boundary conditions (2.2). We use the water depth 
$H$ as the characteristic length scale, the diffusive time 
$\tau _{\kappa }=H^2/\kappa$ as the time scale, the velocity 
$u_{\kappa }=H/\tau _{\kappa }$ as the velocity scale, the temperature difference due to geothermal heating 
$\Delta = FH/k$ as the temperature scale and the pressure 
$p_{\kappa }=\rho _0u_{\kappa }^2$ as the pressure scale. Using 
$T_0$ and 
$p_0+\rho _0g(H-z)$ as temperature and pressure gauges (
$p_0$ is the reference pressure), such that we remove the leading-order mean buoyancy and hydrostatic pressure terms that balance each other, we then define dimensionless variables (denoted by tildes) as
\begin{align} (x,z)= H(\tilde{x},\tilde{z}), \quad t=\tau_{\kappa}\tilde{t}, \quad u= u_{\kappa}\tilde{u}, \quad p = p_0+\rho_0 g(H-z) + p_{\tau}\tilde{p}, \quad T = T_0 + {\Delta} \tilde{T}. \end{align}Substituting (2.3a–e) into (2.1) and (2.2) yields a set of dimensionless equations and boundary conditions, which, dropping tildes for clarity, read
\begin{gather} \partial_{t} {\boldsymbol{u}} - Pr\nabla^2{\boldsymbol{u}} + \boldsymbol{\nabla} p ={-} \left({\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}\right){\boldsymbol{u}} + PrRa_FT\boldsymbol{e}_z, \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{u}} = 0, \end{gather}
\begin{gather}\partial_{t} T - \nabla^2 T ={-} \left({\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}\right) T, \end{gather}and
\begin{gather}{\boldsymbol{u}}(z=0)={\boldsymbol{u}}(z=1)={\boldsymbol{u}}(x={\pm} \varGamma/2)=\boldsymbol{0}, \end{gather}
\begin{gather}\partial_zT(z=0)={-}1, \quad T(z=1)=\dfrac{\varLambda \varGamma}{2}\sin\left(\dfrac{{\rm \pi} x}{\varGamma}\right), \quad \partial_xT(x={\pm} \varGamma/2)=0, \end{gather}wherein
\begin{equation} Pr\equiv\dfrac{\nu}{\kappa}, \quad Ra_F\equiv\dfrac{\alpha g H^4 F}{k\nu\kappa}, \quad \varLambda\equiv\dfrac{\lambda k}{F}, \quad \varGamma\equiv\dfrac{L}{H}, \end{equation}
are the control parameters with 
$Pr$ the Prandtl number, 
$Ra_F$ the flux-based Rayleigh number, 
$\varLambda$ the flux ratio and 
$\varGamma$ the aspect ratio. HC studies typically use a horizontal Rayleigh number as a proxy for the strength of the buoyancy anomaly along the top (or bottom, as the case may be) boundary, which is defined with our notation as
\begin{equation} Ra_L\equiv Ra_F\varLambda\varGamma^4. \end{equation}
In this study, classical RB simulations are run with 
$\varLambda =0$, whereas classical HC simulations are run with 
$\varLambda > 0$ (typically 
$\varLambda =1$) and 
$Ra_F> 0$ (such that 
$Ra_L> 0$) but with an insulating bottom boundary, i.e. such that the temperature condition (2.5) for 
$z=0$ is replaced with 
$\partial _z T(z=0)=0$. Note that although 
$Ra_L$ is a commonly employed Rayleigh number in HC studies, a temperature-based Rayleigh number, typically written as 
$Ra_{\varDelta }$, is usually preferred over the flux-based Rayleigh number 
$Ra_F$ in RB convection studies (primarily because RB convection studies typically consider fixed-temperature boundary conditions). Here we most often use 
$Ra_F$ as a control parameter to describe RB results, but use 
$Ra_{\varDelta }$, which is related to 
$Ra_F$ through 
$Ra_{\varDelta }=Ra_F/Nu_{RB}$ (Johnston & Doering Reference Johnston and Doering2009) with 
$Nu_{RB}$ the RB-specific Nusselt number (defined in § 3.3), when appropriate (e.g. in Appendix C).
We solve (2.4) with the open-source spectral-element code Nek5000 (Fischer Reference Fischer1997; Deville, Fischer & Mund Reference Deville, Fischer and Mund2002), which has been used extensively recently in thermal convection studies (Scheel, Emran & Schumacher Reference Scheel, Emran and Schumacher2013; Léard et al. Reference Léard, Favier, Le Gal and Le Bars2020). The governing equations are cast into weak form and discretised in space by the Galerkin approximation. The Cartesian domain is discretised using 
$n_z$ elements in the vertical direction and 
$\varGamma n_z$ elements in the horizontal direction. Elements have been refined close to all boundaries to properly resolve viscous and thermal boundary layers. The velocity is discretised within each element using Lagrange polynomial interpolants based on tensor-product arrays of Gauss–Lobatto–Legendre quadrature points. The polynomial order 
$l_d$ of the expansion basis on each element varies between 
$7$ and 
$11$ in this study. We use the 3/2 rule for dealiasing, i.e. with extended dealiased polynomial order 
$3/2l_d$. Convergence has been tested by gradually increasing the polynomial order for a fixed number of elements (Scheel et al. Reference Scheel, Emran and Schumacher2013). The nonlinear terms are treated explicitly by a second-order extrapolation scheme whereas the viscous terms are treated implicitly by a second-order backward differentiation scheme. The list of simulations performed with corresponding physical and numerical parameters is provided in table 2 in Appendix A.
The initial state is motionless and has uniform temperature distribution superimposed with low-amplitude background noise (except in § 3.5 where we perform continuation). We explore the fluid dynamics in the 
$(Ra_F,\varLambda )$ parameter space as shown in figure 1(b), i.e. with 
$10^6 \leq Ra_F \leq 10^9$ and 
$0 \leq \varLambda \leq 1$. Simulations with 
$\varLambda =0$ yield canonical RB results; for completeness we also run simulations without geothermal flux, i.e. setting 
$\partial _z T =0$ at 
$z=0$, which yield canonical HC results. For simplicity, we set 
$Pr=1$ in all simulations. Most results are shown for 
$\varGamma =8$, though we also ran simulations with 
$\varGamma =4$ and 16 to partially assess the effect of the aspect ratio; 
$4\leq \varGamma \leq 16$ spans enclosure aspect ratios that are commonly used in HC studies (either experimental or numerical), including 
$\varGamma = 6.2$ (Mullarney et al. Reference Mullarney, Griffiths and Hughes2004; Gayen et al. Reference Gayen, Griffiths, Hughes and Saenz2013, Reference Gayen, Griffiths and Hughes2014; Tsai et al. Reference Tsai, Hussam, King and Sheard2020) and 
$\varGamma = 10$ (e.g. Shishkina & Wagner Reference Shishkina and Wagner2016). Note that different colours highlight different 
$Ra_F$ in figure 1(b), whereas squares, circles and diamonds denote results obtained for 
$\varGamma =4,$ 8 and 16, respectively; stars show the results of pure (no geothermal flux) HC simulations. As we will show, 
$\varLambda$ is a better indicator of the regime dynamics than 
$Ra_L/Ra_F$. The range of 
$0\leq \varLambda \leq 1$ values is motivated by subglacial lakes, as explained in Appendix B.
3. Results
3.1. Flow regimes
We first show in figures 2 and 3 snapshots of the velocity and temperature fields at statistical steady state for 
$Ra_F=10^8$ and 
$\varGamma =8$ with 
$\varLambda$ increasing from top to bottom. The top plot in figure 2 shows the velocity field obtained for 
$\varLambda =0$, which is the canonical RB case. The flow is organised into four pairs of counter-rotating rolls, with characteristic length scale approximately equal to twice the domain height, as expected from linear stability analysis (Chandrasekhar Reference Chandrasekhar1961). As 
$\varLambda$ increases, the overturning cells become distorted because of the buoyancy anomaly on the top boundary, which triggers preferentially leftward flows in the upper half of the domain. For 
$\varLambda =10^{-2}$ (third plot from the top in figure 2), only three pairs of counter-rotating cells co-exist, and the three counter-clockwise cells are elongated because the upper leftward flow branch they support is enhanced by the top boundary. As 
$\varLambda$ increases above 
$10^{-2}$, the flow becomes dominated by HC, which includes intense down-welling below the cold end (left) of the top boundary that trigger vigorous local overturning, and a large-scale counter-clockwise current, which inhibits the growth of RB-like cells in the rest of the domain. The velocity fields obtained for 
$\varLambda =1$ with (second plot from bottom) and without (bottom plot) bottom heat flux are similar, although the former displays stronger meanders of the large-scale flow near the bottom boundary.
Snapshots of the velocity magnitude 
$\sqrt {u^2+w^2}$ and velocity vector field (shown by arrows) at the end of the simulations with 
$Ra_F=10^8$ and 
$\varLambda$ increasing from (a) to (e): (a) 
$\varLambda =0$; (b) 
$\varLambda =10^{-3}$; (c) 
$\varLambda =10^{-2}$; (d) 
$\varLambda =10^{-1}$; (e) 
$\varLambda =1$. Panel ( f) shows the results obtained for pure HC (no geothermal flux), i.e. with 
$\varLambda =1$ but with an insulating bottom boundary (no geothermal flux), such that 
$Ra_L\approx 4\times 10^{11}$. The characteristic length scale of overturning motions is estimated from the horizontal flow at 
$z=0.9$ (see § 3.4), which is highlighted by the blue dashed lines.
Snapshots of the temperature field at the end of the simulations with 
$Ra_F=10^8$ and 
$\varLambda$ increasing from (a) to (e): (a) 
$\varLambda =0$; (b) 
$\varLambda =10^{-3}$; (c) 
$\varLambda =10^{-2}$; (d) 
$\varLambda =10^{-1}$; (e) 
$\varLambda =1$. The solid (respectively, dashed) lines show positive (respectively, negative) contours of the streamfunction, which is set to zero at the bottom left corner. Panel ( f) shows the results obtained with pure HC (no geothermal flux) with 
$\varLambda=1$, such that 
$Ra_L\approx 4\times 10^{11}$. Note that the colour bars are not the same between plots.
The temperature field in figure 3 also highlights the existence of RB-like cells for small 
$\varLambda$, which merge at intermediate 
$\varLambda$ values. For 
$\varLambda =1$ (2 bottom plots of figure 3), the bulk temperature becomes significantly smaller than 0, which is the average temperature of the top boundary, because mixing occurs preferentially on the left of the domain where the top fluid is cold and the down-welling is intense; a warm layer develops near the top right-hand side of the domain but does not mix with the bulk as it is locally stably stratified. It is noteworthy to remark that the bulk temperature is cooler on the left-hand side than on the right-hand side of the domain at intermediate 
$\varLambda =10^{-2}$ (3rd plot from the top). This horizontal temperature gradient of the bulk, which is maintained across multiple cells, builds up until it becomes so great that a single counter-rotating large-scale flow takes over. However, the large-scale flow is short lived because 
$\varLambda$ is small, such that another cycle of RB cells emerges until the horizontal temperature gradient becomes too large again. This subtle bursting dynamics points toward the possible existence of hysteresis, which we investigate in § 3.5.
3.2. Mean temperature and Reynolds number
In §§ 3.2 and 3.3 we first show that all simulations reach a statistical steady state. Then we explore the scaling trends of the Reynolds and Nusselt numbers with the problem parameters, and we demonstrate that the latter can be used to distinguish RB- from HC-dominated simulations. We use 
$\langle X \rangle _x$ to denote 
$x$-averaged variables, 
$\langle X \rangle$ to denote volume-averaged variables and overline 
$\bar {X}$ to denote temporal averages at statistical steady state (typically from 
$t\geq 1$ onward). Whenever relevant, we show the standard deviation due to temporal fluctuations of averaged quantities with vertical error bars. Note, however, that the standard deviation is always small, such that the error bars are often smaller than the marker size and thus barely visible.
Figure 4 shows the temporal evolution of the volume-averaged temperature 
$\langle T \rangle$ (top row) and Reynolds number 
$\widehat {Re}$ (bottom row). In this study, we use the kinetic energy density to construct the Reynolds number, i.e. such that
\begin{equation} \widehat{Re} \equiv \sqrt{\langle u^2 + w^2 \rangle}, \end{equation}which is close to the Reynolds number based on the velocity root mean square (not shown). The time-averaged Reynolds number is
\begin{equation} Re \equiv \sqrt{\overline{\langle u^2 + w^2 \rangle}}. \end{equation}
Both 
$\langle T \rangle$ and 
$\widehat {Re}$ display a sharp transient followed by a statistical steady state (small fluctuations around a mean value independent of time) at approximately 
$t=0.5$; thus, here we use (conservatively) 
$t=1$ as the initial time for time averaging of output variables representative of the statistical steady state. For small 
$\varLambda$ (figure 4a), the mean temperature increases from zero up to a small but positive value as a result of geothermal heating, which dominates over HC. For large 
$\varLambda$ (figure 4b), down-welling beneath the cold end of the top boundary cools down the fluid more efficiently than geothermal heating can warm it, such that the bulk temperature becomes negative. The Reynolds number exhibits stronger fluctuations in time than the mean temperature, which are on the overturning time scale. As expected, 
$\widehat {Re}$ increases with 
$Ra_F$, as can be seen from the stacks of lines shown by blue (
$Ra_F=10^6$), orange (
$Ra_F=10^7$), green (
$Ra_F=10^8$) and red (
$Ra_F=10^9$) colours that are successively on top of each other (results shown in figure 4(d) would be above those shown in figure 4c). The effect of 
$\varLambda$ (colour intensity) and 
$\varGamma$ (line thickness) on 
$\widehat {Re}$, which is relatively weak for our set of simulations (although clearly visible for 
$Ra_F=10^7$ and 
$Ra_F=10^8$ in figure 4), is commented on in greater detail when discussing figure 5(d).
Time history of volume-averaged (a,b) temperature and (c,d) Reynolds number. Simulation results are split into two different subplots (a,c and b,d) for better visibility (as indicated by the plot titles). Different colours correspond to different Rayleigh numbers 
$Ra_F$, lighter colours indicate higher 
$\varLambda$ and line width increases with 
$\varGamma$. Time-averaged variables representative of the statistical steady state are integrated from 
$t=1$ (shown by the vertical dashed lines) onward.
(a) Mean temperature at statistical steady state for all simulations as a function of the flux ratio 
$\varLambda$. The black symbols show the coldest temperature on the top boundary, i.e. 
$T=-\varLambda \varGamma /2$. (b,c) Reynolds number 
$Re$ as a function of 
$Ra_F$ and 
$Ra_L$, respectively, with scalings 
$Re=c_{RB}Ra_F^{d_{RB}}$ and 
$Re=c_{HC}Ra_L^{d_{HC}}$ shown by the black solid lines. (d) Compensated Reynolds number as a function of 
$\varLambda$. The star symbols show results obtained for pure HC simulations.
We show in figure 5 the time-averaged mean temperature 
$\overline {\langle T \rangle }$ and Reynolds number 
$Re$ at statistical steady state as functions of the problem parameters. Figure 5(a) shows that the mean temperature decreases quickly with increasing 
$\varLambda \geq 10^{-2}$ because down-welling below the cold top boundary (
$T(x=-\varGamma /2)=-\varLambda \varGamma /2$ shown by black markers) becomes sufficiently strong to lower the bulk temperature. The mean temperature also decreases with increasing 
$Ra_F$ (blue markers above orange, green and red markers) because increasing geothermal heating makes mixing more efficient, thus lowering temperature differences, whereas increasing HC increases mixing from the cold region of the top boundary.
Figure 5(b) shows that the power law curve 
$Re=c_{RB}Ra_F^{d_{RB}}$ (black solid line) provides a good prediction for the Reynolds number as a function of 
$Ra_F$ for most simulations. Here, pre-factor 
$c_{RB}$ and exponent 
$d_{RB}$ are obtained from best fit with the results for 
$\varLambda =0$ and 
$\varGamma =8$ (see table 1 and Appendix C for a list and discussion of all pre-factors and exponents mentioned in the article). The dependence of 
$Re(Ra_F)$ with 
$\varLambda$ (and 
$\varGamma$) is small, especially at low 
$\varLambda <10^{-2}$ (dark colours), which means that a small horizontal temperature gradient (or change of aspect ratio) has limited effect on the intensity of RB-driven flows. Conversely, figure 5(c) demonstrates that there is a wide spread of 
$Re(Ra_L)$ between simulations, even for relatively large 
$\varLambda \sim 0.1$ (light colours), which means that the circulation is almost always affected by the bottom heat flux, even in the HC-dominated regime obtained for 
$\varLambda >O(10^{-2})$ (as we will show). Accordingly, the power law curve 
$Re= c_{HC}Ra_L^{d_{HC}}$ (black solid line) obtained from pure HC results (shown by the stars) predicts 
$Re$ accurately for large 
$\varLambda \geq 10^{-1}$ only.
Pre-factor and exponent of all power laws referred to in the paper and shown as black solid lines in figures 5(b), 5(c), 7(a) and 8(a). The power laws for 
$Re(Ra_F)$ and 
$Nu(Ra_F)$ are based on four simulations with 
$\varLambda =0$ and 
$\varGamma =8$, whereas the power laws for 
$Re(Ra_L)$ and 
$Nu(Ra_L)$ are based on three simulations without geothermal flux and 
$\varGamma =8$ (see table 2 in Appendix A).
We plot the Reynolds number compensated by the RB scaling in figure 5(d) in order to highlight the effect of 
$\varLambda$ and 
$\varGamma$ on 
$Re$. The spread of 
$Re$ with 
$\varGamma$ and 
$\varLambda$ is overall small (less than 40 % increase) except for 
$Ra_F=10^7$ (up to 100 % increase). At low 
$\varLambda$, the increase of 
$Re$ with 
$\varGamma$ is modest (approximately 10 % or less), in agreement with previous two-dimensional RB studies that have shown a monotonic increase of 
$Re$ with 
$\varGamma$ in two-dimensional bounded domains with a 
$\sim$20 %-increase plateau reached at 
$\varGamma \approx 10$ (Van Der Poel et al. Reference Van Der Poel, Stevens, Sugiyama and Lohse2012). At large 
$\varLambda$, the increase of 
$Re$ with 
$\varGamma$ becomes more important (up to 50 %) as HC, which is controlled by 
$Ra_L\propto \varGamma ^4$, starts dominating. An increase of 
$Re$ with 
$\varLambda$ may be expected at large 
$\varLambda$, i.e. once HC dominates, because 
$Ra_L\propto \varLambda$. Here, the increase of 
$Re$ with 
$\varLambda$ is significant only for the set of simulations with 
$Ra_F=10^7$ and 
$\varGamma =16$, i.e. which attain the largest 
$Ra_L$ (figure 5c). This suggests that our simulations transition from RB to HC (around 
$\varLambda \approx 10^{-2}$, as is shown) primarily through a change of flow structure, not intensity (
$Re$ almost fixed). The exact scaling of 
$Re$ with 
$\varLambda$ and 
$\varGamma$ in the HC regime is beyond the scope of our study given the limited amount of data in the large 
$\varLambda$ and 
$\varGamma$ limit. Note, however, that the enhancement of 
$Re$ by geothermal heating in the HC regime seen in figure 5(d) (compare star symbols and circles at large 
$\varLambda$) is compatible with previous studies (Mullarney et al. Reference Mullarney, Griffiths and Hughes2006; Wang et al. Reference Wang, Huang, Zhou and Xia2016). Both studies indeed reported a 
$O(1)$ increase due to geothermal heating in the volume flux (measured through the maximum of the streamfunction), which may be expected to scale linearly with 
$Re$, for 
$\varLambda Nu^{\chi }_{HC}\sim O(10)$, which is equivalent to 
$\varLambda \sim O(0.1)$–
$O(1)$ as the Nusselt number of HC (defined in § 3.3) 
$Nu^{\chi }_{HC}\sim O(10)$–
$O(100)$ for 
$Ra_L\sim O(10^{9})$–
$O(10^{12})$ (Mullarney et al. Reference Mullarney, Griffiths and Hughes2004).
3.3. Heat transfer efficiency via Nusselt numbers
In order to assess the efficiency of heat transfer, we define two distinct Nusselt numbers designed to measure heat transfer due to either RB convection or HC. The Nusselt number of RB convection is defined as
\begin{equation} Nu_{RB} \equiv \dfrac{FH}{k\left[\overline{\langle T^{{dim}}(z=0) \rangle_x}-\min(T^{{dim}})\right]} = \dfrac{1}{\overline{\langle T(z=0) \rangle_x}-\min(T)}, \end{equation}
where superscript 
${dim}$ denotes a dimensional variable and 
$\min (T)=-\varLambda \varGamma /2$ is the coldest temperature achieved on the top boundary. As is customary in classical RB studies, 
$Nu_{RB}$ compares the effective vertical heat flux that comes out of the system to the vertical heat flux that would be obtained from conduction only due to the temperature difference between the top and bottom boundaries. Here we use 
$\min (T)$ instead of the mean 
$T=0$ value as a gauge for the top temperature in the denominator in (3.3) to ensure 
$Nu_{RB}>0$, because 
$\overline {\langle T(z=0) \rangle _x}$ becomes negative for large-enough 
$\varLambda$ (figure 5a). We note that the use of 
$\min (T)$ affects the definition of the heat flux of the diffusive state; however, in the RB regime, i.e. for 
$\varLambda \ll 10^{-2}$, 
$\min (T)$ is always much smaller than the mean bottom temperature, such that 
$Nu_{RB}$ does converge toward the commonly defined Nusselt number of RB convection.
Unlike RB convection, HC sets up a large-scale circulation, which produces large asymmetry between the left-hand and right-hand sides of the fluid domain: heat extraction occurs below the cold (left) end of the top boundary whereas heat is replenished on the warm (right) end via a slow return flow and a thick boundary layer. This asymmetry can be seen in figure 6(a), which shows the time-averaged heat flux on the top boundary for RB-dominated, mixed RBH convection and HC-dominated simulations with 
$Ra_F=10^8$: the RB results (
$\varLambda =0$) show an oscillatory pattern of heat flux linked to the underlying overturning cells, the intermediate 
$\varLambda =10^{-2}$ value leads to skewed oscillations and the large 
$\varLambda =1$ value yields a much-larger monotonic and anti-symmetric heat flux pattern with respect to the middle of the domain 
$x=0$. Figure 6(b) shows similar patterns for the temperature on the bottom boundary. The temperature increase from left to right imposed on the top boundary is imprinted on the bottom boundary quite clearly for large 
$\varLambda$ but also for 
$\varLambda$ as small as 
$10^{-2}$.
(a) Time-averaged heat flux along the top boundary for 
$Ra_F=10^8$, 
$\varGamma =8$ and four different 
$\varLambda$ values (see the legend). (b) Same as (a) but for the temperature on the bottom boundary with the horizontal mean removed. The thin black solid line shows the temperature on the top boundary for 
$\varLambda =10^{-2}$.
Different Nusselt numbers have been used in the past to measure heat transfer by HC, such as the Nusselt number based on the intensity in absolute value of heat extraction and deposition at the top boundary (Rossby Reference Rossby1965, Reference Rossby1998), i.e.
\begin{equation} Nu^{{abs}}_{HC} \equiv \dfrac{\overline{\langle \left|\partial_z T (z=1)\right| \rangle_x}}{ \langle \left|\partial_z T_{{diff}}(z=1)\right| \rangle_x}, \end{equation}and the Nusselt number based on the heat flux above the cold (or warm) half of the boundary (Sheard & King Reference Sheard and King2011), i.e.
\begin{equation} Nu^{{half}}_{HC} \equiv \dfrac{\overline{\langle \partial_z T (z=1) \rangle_{x<0}} }{ \langle \partial_z T_{{diff}}(z=1) \rangle_{x<0} }, \end{equation}
where subscript 
$_{{diff}}$ means of the diffusive base state. Although 
$Nu^{{abs}}_{HC}$ and 
$Nu^{{half}}_{HC}$ are intuitive, they are not based on fundamental variables of the energy budget. Thus, following the recommendation by Rocha et al. (Reference Rocha, Constantinou, Llewellyn Smith and Young2020b) we use another definition of the Nusselt number, i.e.
\begin{equation} Nu^{\chi}_{HC} \equiv \dfrac{\chi}{\chi_{{diff}}}, \end{equation}
where 
$\chi =\overline {\langle |\boldsymbol {\nabla } T|^2 \rangle }$ is the dissipation of temperature variance, which is related to the horizontal heat transport and entropy production. We show in Appendix D that 
$\chi$ is proportional to the correlation of the heat flux with the temperature on the bounding horizontal plates (as already shown by Rocha et al. (Reference Rocha, Constantinou, Llewellyn Smith and Young2020b) without geothermal heating), such that the Nusselt number can be calculated from line integrals as
\begin{equation} Nu^{\chi}_{HC} = \dfrac{\overline{\langle T(z=1)\partial_z T(z=1) - T(z=0)\partial_z T(z=0) \rangle_x}}{\overline{\langle T_{{diff}}(z=1)\partial_z T_{{diff}}(z=1) - T_{{diff}}(z=0)\partial_z T_{{diff}}(z=0) \rangle_x}}, \end{equation}
with an analytical expression for the denominator given in (D5). We would like to note that although 
$Nu^{\chi }_{HC}=\chi /\chi _{{diff}}$ is thermodynamically compelling, it is very close to 
$Nu^{{abs}}_{HC}$ and 
$Nu^{{half}}_{HC}$ for pure HC simulations. For mixed RBH simulations, differences exist but are due at leading order to the diffusive normalisation (see the details in Appendix E).
Figure 7(a) shows 
$Nu_{RB}$ (3.3) as a function of 
$Ra_F$. Simulation results obtained for 
$\varLambda \leq 10^{-2}$ all collapse very well on the power law curve shown by the black solid line, which was obtained from best fit for 
$\varLambda =0$ (see table 1) and whose exponent is compatible with the classical scaling law of RB convection for moderate Rayleigh numbers (details in Appendix C). The plot of the compensated 
$Nu_{RB}$ number in figure 7(b) highlights the deviation of the RB-based Nusselt number from the RB scaling with increasing 
$\varLambda$. In particular, the difference between 
$Nu_{RB}$ and 
$a_{RB}Ra_F^{b_{RB}}$ becomes of order 1 when 
$\varLambda \geq 10^{-2}$ (grey area), which thus marks the transition from RB convection to HC. Note that the decrease of 
$Nu_{RB}$ with 
$\varLambda$ is primarily the result of an increase of the denominator in (3.3), which is due to the fact that the coldest temperature on the top boundary increases more quickly than the mean bottom temperature (in absolute value).
(a) RB-based Nusselt number 
$Nu_{RB}$ as a function of 
$Ra_F$ with scaling law 
$Nu_{RB}=a_{RB}Ra_F^{b_{RB}}$ shown by the black solid line. (b) Compensated RB-based Nusselt number as a function of 
$\varLambda$.
Figure 8(a) shows 
$Nu^{\chi }_{HC}$ as a function of 
$Ra_L$. Simulation results obtained for 
$\varLambda \geq 10^{-2}$ all tend asymptotically (for fixed 
$\varLambda$, or colour intensity) toward power laws parallel (in log–log space) to that shown by the black solid line (see, e.g., the dashed line), which was obtained from best fit for 
$\varLambda =1$, 
$\varGamma =8$ and without geothermal flux (cf. perfect overlap with star symbols); the exponent 
$b_{HC}$ is compatible with the 1/5 exponent predicted by Rossby (Reference Rossby1965). Figure 8(a) shows that 
$Nu^{\chi }_{HC}$ can help distinguish simulations dominated by RB convection from simulations dominated by HC: at small 
$\varLambda$ (dark symbols with typically small 
$Ra_L$), 
$Nu^{\chi }_{HC} < 1$, whereas at large 
$\varLambda$ (light symbols with typically large 
$Ra_L$), 
$Nu^{\chi }_{HC}\gg 1$. The large spread of 
$Nu^{\chi }_{HC}$ with 
$\varLambda$ and 
$\varGamma$ for fixed 
$Ra_L$ in the HC limit (i.e. at large 
$\varLambda$) is somewhat unexpected but can be explained: it is due to the diffusive normalisation 
$\chi _{{diff}}$. First, figure 8(b) shows that dividing 
$Nu^{\chi }_{HC}$ by the pure HC scaling 
$a_{HC}Ra_L^{b_{HC}}$ results in a steep scaling with 
$\varLambda$ (
$+2$ slope shown by the solid line) for 
$\varLambda \gg 10^{-2}$, which is exactly the scaling of the diffusive normalisation 
$1/\chi _{{diff}} - 1 \propto \varLambda ^2$ once Taylor expanded in the small 
$\varLambda ^2$ limit (cf. (D5)). Second, figure 8(c) demonstrates that replacing 
$\chi _{{diff}}$ with 
$\chi _{{dim}}={\rm \pi}^2\varLambda ^2/8$, which is dimensionally equivalent but discards geothermal heating and approximates 
$\tanh \left({\rm \pi}/\varGamma \right)\approx {\rm \pi}/\varGamma$ (large 
$\varGamma$ limit), in the definition of 
$Nu_{HC}^{\chi }$ (3.7) yields a perfect overlap of all simulation results obtained for large 
$\varLambda \gg 10^{-2}$ with the power law 
$a_{HC}Ra_L^{b_{HC}}$. All this shows that the spread of 
$Nu_{HC}^{\chi }$ with 
$\varLambda$ and 
$\varGamma$ at large 
$\varLambda$ is due to the diffusive normalisation 
$\chi _{{diff}}$ (denominator in (3.7)), not 
$\chi$, because 
$\chi _{{diff}}$ remains sensitive to aspect ratio (
$\varGamma$) and flow topology (
$\varLambda$) for a much wider range of parameters than 
$\chi$ (as is the case for other definitions of the Nusselt number, see Appendix E). This result is in agreement with previous studies (Sheard & King Reference Sheard and King2011) who found no 
$\varGamma$ dependence using a flux-based definition of the Nusselt number normalised by 
$\varLambda$ rather than the diffusive solution in the convection-dominated regime, which spans a large range of 
$Ra_L$ encompassing our simulation parameters (see also Hossain, Vo & Sheard Reference Hossain, Vo and Sheard2019, for the effect of 
$\varGamma \gg 1$ on the transition from diffusion- to convection-dominated HC).
(a) HC-based Nusselt number 
$Nu^{\chi }_{HC}$ as a function of 
$Ra_L$ with scaling law shown by the black solid line (derived from pure HC simulation results). Note that the dashed line shows a similar scaling albeit with a different pre-factor. (b) Compensated Nusselt number as a function of 
$\varLambda$. The solid line has a 
$+2$ slope and shows the scaling of 
$\chi _{{diff}}$ for 
$\varLambda \ll 1$. (c) Same as (b) but for a Nusselt number with denominator 
$\chi _{{diff}}$ replaced by 
$\chi _{{dim}}={\rm \pi}^2\varLambda ^2/8$, which discards the effect of geothermal heating and aspect ratio.
3.4. Characteristic length scale from auto-correlation
In this section we estimate the characteristic length scale of the overturning cells in order to further demonstrate that 
$\varLambda \approx 10^{-2}$ marks the transition from RB convection to HC. As the leftward flow near the top boundary is an emblematic feature of HC, we use the variations in 
$x$ of the horizontal velocity at 
$z=0.9$ as diagnostic. We first show the time- and 
$x$-averaged horizontal flow at 
$z=0.9$ normalised by 
$Re$ as a function of 
$\varLambda$ in figure 9. For 
$\varLambda \geq 10^{-2}$, 
$\overline {\langle -u(z=0.9) \rangle } \approx Re >0$, which indicates that the large-scale leftward current dominates the dynamics. For small 
$\varLambda \ll 10^{-2}$, including 
$\varLambda =0$, 
$\overline {\langle -u(z=0.9) \rangle } \ll Re$ as expected, but is not always zero. Indeed a closed domain with moderate aspect ratio can deform RB-like overturning cells durably, such that a mean flow exists in the upper half of the domain, which is balanced by an equivalently strong return flow in the bottom half. Note that the mean horizontal flow exhibits complex fluctuations in time near the transition 
$\varLambda =10^{-2}$ due to the superposition of the RB and HC dynamics, which cannot be inferred from figure 9 but is discussed in § 3.5.
Temporally and horizontally averaged horizontal flow near the top boundary at 
$z=0.9$ normalised by the Reynolds number as a function of 
$\varLambda$. Here 
$\overline {\langle -u(z=0.9) \rangle _x}$ is of the same order as 
$Re$ for relatively large 
$\varLambda$.
The calculation of the characteristic length scale 
$\ell$ of overturning motions from 
$u(z=0.9)$ is illustrated in figure 10 for a RB-like case (
$\varLambda =10^{-3}$; top row) and a HC-like case (
$\varLambda =1$; bottom row) with 
$Ra_F=10^{8}$ and 
$\varGamma =8$. Figure 10(a) shows the 
$(x,t)$-Hovmöller diagram of 
$u(z=0.9)$ for the RB-like simulation at statistical steady state. Four pairs of counter-rotating rolls can be identified, which have approximately the same width and slightly meander in time. The auto-correlation function of 
$u$ in 
$x$ is defined as
\begin{equation} \mathcal{R}_{uu}(x) \equiv \dfrac{\int_{-\varGamma/2}^{\varGamma/2} u(x') u(x'+x) \,\mathrm{d} x'}{\int_{-\varGamma/2}^{\varGamma/2} u(x') u(x')\,\mathrm{d} x'}, \end{equation}
where 
$u(x'+x>\varGamma /2)$ is set to 0 and 
$x>0$ is the lag. The auto-correlation function evaluated at 
$z=0.9$ shows an oscillatory pattern (figure 10b), like 
$u(z=0.9)$, which is damped as the lag increases due to the presence of fluctuations and the scarcity of data for large lag. We estimate the characteristic length scale 
$\ell$ of the overturning cells from the first minimum of the time-averaged auto-correlation function, which is always the largest in absolute value for all our simulations. For 
$\varLambda =10^{-3}$, figure 10(c) shows that 
$\ell \approx 1$ (dashed solid line), which is approximately the value expected for an unconfined RB roll (rotating clockwise or anti-clockwise). Figure 10(d–f) show the same results as figure 10(a–c) but for 
$\varLambda =1$. The horizontal flow near the top boundary is now approximately negative everywhere, such that the auto-correlation function monotonically decreases with the lag in 
$x$. Thus, the characteristic length scale equals the domain size, i.e. 
$\ell =\varGamma$, as shown by the vertical dashed line in figure 10(f). We note that the auto-correlation function (3.8) may be defined differently to account for the scarcity of data at large lag (e.g. replacing 
$u(x')$ with 
$u(x'+x)$ in the denominator). However, such a definition creates large boundary effects (not shown), which make the calculation of 
$\ell$ more complicated (because of large variations for large lag), although ultimately unchanged since the locations of most extrema are not modified for small-to-moderate lags.
(a) Horizontal flow near the top boundary (
$z=0.9$) as a function of 
$(x,t)$ for 
$Ra_F=10^8$, 
$\varGamma =8$ and 
$\varLambda =10^{-3}$. (b) Auto-correlation function 
$\mathcal {R}_{uu}$ in 
$x$ of the horizontal flow shown in (a) as a function of time 
$t$. (c) Time average auto-correlation 
$\overline {\mathcal {R}_{uu}}$ as a function of lag in 
$x$. The dashed line shows the first minimum (trough), which we identify as the characteristic length scale of the horizontal flow. (d–f) Same as (a–c) but for 
$\varLambda =1$. For a monotonically decreasing auto-correlation function, the characteristic length scale equals the domain length 
$\varGamma$.
The characteristic length scale 
$\ell$ obtained from the auto-correlation function of 
$u(z=0.9)$ is shown for all simulations as a function of 
$\varLambda$ in figure 11. For 
$\varLambda \ll 10^{-2}$, 
$\ell \approx 1$ (dotted line), as expected for classical RB simulations, although there is a small spread between 0.85 and 1.4, which is a result of bounded-box effects or, at non-zero 
$\varLambda$ values, bursts of large-scale currents sweeping away RB cells (further explained in § 3.5). For 
$\varLambda \gg 10^{-2}$, 
$\ell =\varGamma$, which is here equal to either 4, 8 or 16 (levels shown by dashed lines in figure 11) and an indication that the dynamics is dominated by HC.
Characteristic length scale 
$\ell$ at statistical steady state as a function of 
$\varLambda$ for all simulations. The dashed lines show the aspect ratios (or domain lengths), i.e. 
$\varGamma =4$, 8 and 
$16$. At large 
$\varLambda$, 
$\ell =\varGamma$, which is an indication of HC dominating the dynamics, whereas at small 
$\varLambda$, 
$\ell \approx 1$ (dotted line) indicates multiple overturning cells, as in classical RB convection.
3.5. Bi-stability and bursts near the transition
This section focuses on the transitional regime between the RB and HC regimes obtained for 
$\varLambda \ll 10^{-2}$ and 
$\varLambda \gg 10^{-2}$, respectively. To this end, we present numerical results that use the method of discrete continuation. We fix 
$Ra_F=10^6$ and 
$\varGamma =12$, and we gradually vary the flux ratio 
$\varLambda$ in the range 
$[0,0.03]$. For each value of 
$\varLambda$, we integrate the system for 
$20$ diffusive timescales to ensure that the system has reached a statistical quasi-steady state. This integration time is increased up to 
$100$ diffusive timescales for cases close to bifurcation points, hence limiting this study to only moderately large Rayleigh number and aspect ratio. In order to track bifurcations between different states, we consider the following averaged horizontal flow
\begin{equation} \overline{\langle u\rangle}_{z>0.5}=\frac{2}{\varGamma}\int_{-\varGamma/2}^{\varGamma/2}\int_{0.5}^1 \bar{u} \,\mathrm{d} x \,\mathrm{d}z, \end{equation}and mean temperature difference between the right- and left-hand sides of the domain
\begin{equation} \Delta T_H=\overline{\langle T\rangle}_{x>\varGamma/4}-\overline{\langle T\rangle}_{x<{-}\varGamma/4}=\frac{4}{\varGamma}\int_{\varGamma/4}^{\varGamma/2}\int_{0}^1 \bar{T} \,\mathrm{d} x \,\mathrm{d}z-\frac{4}{\varGamma}\int_{-\varGamma/2}^{-\varGamma/4}\int_{0}^1 \bar{T} \,\mathrm{d} x \,\mathrm{d}z, \end{equation}where time averaging is only performed once the statistical steady state has been reached. Due to the temperature profile imposed on the top boundary (2.5), we expect the time-averaged mean flow (3.9) to be generally negative and the mean temperature difference (3.10) to be positive. Note that by mass conservation, a similar integration as (3.9) performed on the lower half of the domain would yield exactly the opposite value.
We start by gradually increasing 
$\varLambda$ from 0 to 
$0.03$. The results are shown in figure 12(a,b) as round symbols. We recall that each point corresponds to at least 
$20$ diffusive timescales and up to 
$100$ timescales for cases on each side of a given transition. The mean horizontal flow in the upper half of the fluid domain gradually increases in absolute value until a first transition occurs at 
$\varLambda \approx 0.016$ (figure 12a). This transition corresponds to the first destabilisation of the RB cells with approximately unit aspect ratio, which merge into horizontally elongated cells. For our particular choice of 
$\varGamma =12$, the flow is initially organised into 12 cells and switches to six more elongated cells. A secondary transition between these six cells and only one extended cell is observed at 
$\varLambda =0.02$. The system remains in this state for larger values of 
$\varLambda$. The two transitions at 
$\varLambda =0.016$ and 
$\varLambda =0.02$ are also observed in the mean horizontal temperature difference (figure 12b). The mean horizontal temperature difference 
$\Delta T_H$ initially follows the temperature difference imposed along the top boundary, which is equal to 
$\varLambda \varGamma$ (dashed line in figure 12b). As 
$\varLambda$ increases, 
$\Delta T_H$ slowly deviates downward until it recovers the imposed value 
$\varLambda \varGamma$ just after the transition at 
$\varLambda =0.016$. After the second transition, 
$\Delta T_H$ clearly follows a different trend and increases less rapidly with 
$\varLambda$, which is indicative of the growing efficiency of HC in its ability to mix the imposed horizontal temperature difference.
Bifurcation diagrams based on averaged velocity and temperature. (a) Mean horizontal flow averaged over time and the upper half of the domain (3.9) as a function of the heat flux ratio 
$\varLambda$. The parameters are 
$Ra_F=10^6$ and 
$\varGamma =12$. Round symbols correspond to gradually increasing 
$\varLambda$ while cross symbols correspond to gradually decreasing 
$\varLambda$. Each symbol corresponds to a simulation lasting from 
$20$ and up to 
$100$ (close to bifurcations) diffusive timescales. The shaded area shows the region of bi-stability. (b) Same as (a) but for the depth-averaged difference of temperature between the right- (
$x>\varGamma /4$) and left-hand (
$x<-\varGamma /4$) sides of the domain (3.10). (c) Snapshots in 
$(x,z)$ of the temperature field with streamlines shown as black contours for 
$\varLambda =0.014$ when 
$\varLambda$ is going up (top panel) and down (bottom panel). (d) Same as (c) but for 
$\varLambda =0.018$. (e) Same as (c) but for 
$\varLambda =0.028$ (branch with 
$\varLambda$ going up only). ( f) Horizontal flow at 
$z=0.9$ as a function of 
$x$ and 
$t$ (initial time set to 0 arbitrarily) for the simulation right after the first transition of the increasing-
$\varLambda$ (blue) branch on panels 12(a,b) (
$\varLambda =0.016$; highlighted by triangles). Note that merging events, e.g. at 
$x\approx 2$ and 
$t \approx 7$, appear discontinuous because they occur on short time scales.
To explore the possibility of multi-stability, we gradually decrease 
$\varLambda$ from a given equilibrium state. We do so independently for each observed transition, i.e. starting first at 
$\varLambda =0.016$ and then at 
$\varLambda =0.02$. We show the corresponding descending branches in figure 12(a,b) using cross symbols. We clearly observe hysteretic behaviour over a large range of heat flux ratio 
$\varLambda$, which is highlighted by the green shaded region in figure 12(a,b). The bi-stability is characterised by the separation of the ascending and descending branches, which are each connected to distinct flow organisations: for 
$0.01<\varLambda <0.02$, there always exist at least two flow states for the same 
$\varLambda$ that have a different number of convective rolls; each roll having an aspect ratio that can vary between 1 and 
$\varGamma$. Thus, the actual number of convective cells crucially depends on the history of the system. The possibility to have different number of rolls for the same 
$\varLambda$ can be seen from the lower-left panels on figure 12 where we show the temperature field and streamlines for two bi-stable states at 
$\varLambda =0.014$ (figure 12c) and 
$\varLambda =0.018$ (figure 12d). For completeness, we also show in figure 12(e) an example of the final state of pure HC reached once 
$\varLambda >0.02$. Although the mean flow is always stronger for descending than for ascending branches, the number of rolls has a non-trivial effect on the mean horizontal temperature difference, because 
$\Delta T_H$ can either increase (orange line above blue line in figure 12b) or decrease (blue line above green line in figure 12b) with decreasing roll number. The former behaviour may be explained from the persistence of HC dynamics along the descending branch (green line), which efficiently mixes the horizontal temperature gradient, whereas the latter suggests, counter-intuitively, that, for our choice of parameters, 6 rolls are more efficient at maintaining a large 
$\Delta T_H$ than 12 rolls.
We would like to make a few notes regarding the results of figure 12(a–e). First, the merging of rolls does not percolate from either the left- or right-hand side of the domain, where HC drives distinct dynamics (intense down-welling versus weak up-welling). Rather, it occurs in the bulk through nucleation, as shown, e.g., in figure 12( f) wherein the first merging occurs between the second and third pair of rolls from the right boundary. Second, for low values of 
$\varLambda <0.01$, we still observe bi-stability because the system does not recover exactly its initial state (characterising the increasing branch) as we continue decreasing 
$\varLambda$. However, the two solutions only differ by the spatial organisation of the convective rolls, not their numbers. Third, we do not observe tri-stability close to 
$\varLambda \approx 0.016$: the ascending trajectory bifurcates exactly where the trajectory descending from 
$\varLambda =0.02$ falls back onto the descending trajectory initialised from 
$\varLambda =0.016$. This does not mean that tri-stability cannot be obtained in RBH convection systems. In fact, we expect the bifurcation diagram to be richer than can be deciphered by our study, especially as 
$\varGamma$ increases. Finally, even though we have integrated these bi-stable states for up to 
$100$ diffusive timescales based on the height of the fluid domain, we cannot exclude the possibility that rare spontaneous transitions can occur between them (partly because the diffusive timescale based on the horizontal size of the domain scales with 
$\varGamma ^2\gg 1$, hence is much longer). One way to more firmly prove that bi-stable states exist would require a detailed measurement of the transition time between different states, as was recently done in thin-layer turbulence (de Wit, van Kan & Alexakis Reference de Wit, van Kan and Alexakis2022), which is beyond the scope of the present study.
We expect the hysteretic behaviour to persist for aspect ratios larger than 
$\varGamma =O(10)$, because this allows for even more distinct convective cell configurations (see Wang et al. Reference Wang, Verzicco, Lohse and Shishkina2020 for classical RB convection). However, it is unclear whether such dynamics will persist at large Rayleigh numbers. The analysis of bi-stable dynamics at large 
$Ra_F$ is obviously very intensive numerically as it requires integrating a fully turbulent system for hundreds of diffusive timescales. Nevertheless, we explore a particular case close to a transition to show that, although we do not observe bi-stability, we do observe bi-modality in certain range of parameters. For example, let us consider the much more turbulent case 
$Ra_F=10^8$ and 
$\varGamma =8$. For 
$\varLambda =0.01$, i.e. right at the lower bound of the bi-stable regime observed for 
$Ra_F=10^6$, we observe spontaneous transitions between states with large and weak horizontal mean flows, which can be seen from the time history of the horizontal mean flow (defined in (3.9)) in figure 13(a). The bi-modality is clear: the horizontal mean flow displays abrupt (negative) HC-like peaks before rapidly relaxing to a state of weaker RB-like mean flow. The mean-flow bursts correspond to abrupt intrusions of warm fluid accumulating below the heated part of the right corner, which can temporarily disrupt the RB cells, and are reminiscent of the burst dynamics observed in RB convection between stress-free plates (Goluskin et al. Reference Goluskin, Johnston, Flierl and Spiegel2014). The two bottom left panels of figure 13 show snapshots of the temperature field of a single simulations exhibiting bi-modality: figure 13(b) displays convective rolls reminiscent of classical RB convection, whereas figure 13(c) shows a state featuring a burst of HC (note that a similar HC burst can be seen in Supplementary Movie available at https://doi.org/10.1017/jfm.2022.613). Figure 13(d) shows the time history of 
$u(z=0.9)$, which can help visualise the horizontal development of the bursts: the bursts always start from the right boundary, then the system relaxes from the left boundary. For this particular value of 
$\varLambda$ and duration of the simulation, the convection cells always recover their RB-like configuration and HC remains intermittent. Whether the system remains in this bi-modal state or eventually converge towards one of the two attractors at long times is an open question.
(a) Mean horizontal flow averaged over the upper half of the domain as a function of time for 
$Ra_F=10^8$, 
$\varGamma =8$ and 
$\varLambda =10^{-2}$. (b,c) Snapshots of the temperature field at times 
$t=3.1$ (vertical dashed line in the top panel) and 
$t=4$ (vertical dotted line), respectively. (d) Spatio-temporal plot of 
$u(z=0.9)$ for the same simulation (slightly earlier times).
We conclude our analysis of RBH dynamics by discussing and comparing the probability density function (p.d.f.) of the horizontal mean flow defined by (3.9) in different regimes (including bi-modality) near the transition. We consider fixed (relatively large) 
$Ra_F=10^8$, 
$\varGamma =8$ and variable 
$\varLambda$ close to the transition, as well as pure RB convection (
$\varLambda =0$) and pure HC (no geothermal flux) for comparison. We first show in figure 14(a) the mean flow p.d.f. at statistical steady state for simulations close to the transition. We clearly see an abrupt transition between a Gaussian distribution for 
$\varLambda =7\times 10^{-3}$ and a skewed distribution with an elongated negative tail for 
$\varLambda =10^{-2}$. As 
$\varLambda$ further increases, the p.d.f.s become increasingly skewed toward large negative values, but maintain a large spread confirming the intermittent nature of the mean flow dynamics. Interestingly, the burst dynamics is a specific property of the mixed RBH convection regime, i.e. which is absent in pure RB or HC regime. Figure 14(b) shows the p.d.f. of the mean horizontal upper flow for three cases at 
$Ra_F=10^8$ and 
$\varGamma =8$: one with 
$\varLambda =0$ (pure RB case), one with 
$\varLambda =0.05$ but no geothermal flux (pure HC) and finally the RBH case with 
$\varLambda =0.05$ and geothermal heating. We clearly see that both RB and HC cases produce mean-flow p.d.f.s that are approximately Gaussian and centered around zero (dotted line) or a small negative value (dashed line), respectively. Conversely, the RBH case yields a p.d.f. with a large spread, which is consistent with rare but intense mean-flow bursts and thus a characteristic feature of RBH dynamics. The emergence of strong mean flows (persistent or bursty) in RBH convection is consistent with earlier results (Mullarney et al. Reference Mullarney, Griffiths and Hughes2006) and confirms the underlying tendency of RB convection to drive large-scale horizontal flows, which are here enhanced by the imposed horizontal temperature gradient. One might expect that the bursts eventually disappear at large 
$\varLambda$, i.e. once the mean flow driven by the upper horizontal temperature gradient becomes dominant compared to the maximum value observed during each burst event. The underlying 
$\varLambda$ threshold would correspond to the lower bound of the HC-dominated regime, highlighted by grey shadings in figures 7, 8, 9 and 11. We leave the detailed investigation of this threshold, which will require exploration of the burst dynamics in the large 
$Ra_F$ and large 
$\varGamma$ limit, to future studies.
(a) Probability density function (p.d.f.) at statistical steady state of 
${\langle u \rangle }_{z>0.5}$ for 
$Ra_F=10^8$, 
$\varGamma =8$ with 
$\varLambda =0.007$, 0.01, 0.02 and 0.05. (b) Same as (a) but for pure RB convection (
$\varLambda =0$; dotted line), pure HC (
$\varLambda =0.05$ but no geothermal flux; dashed line) and RBH convection (
$\varLambda =0.05$; solid line).
4. Discussion and conclusions
We have investigated the dual RBH convection problem via direct numerical simulations in order to identify the transition from RB convection to HC in parameter space. The 
$\varLambda$ parameter, which is the ratio of the top horizontal heat flux divided by the bottom heat flux, clearly controls the system dynamics: RBH transitions from RB convection at small 
$\varLambda \leq 3\times 10^{-3}$ (purple shading in figures 7, 8, 9 and 11) to HC at large 
$\varLambda \geq 3\times 10^{-2}$ (grey shading). Importantly, the transition occurs near 
$\varLambda =10^{-2}$ independently of the aspect ratio 
$\varGamma$ and Rayleigh number 
$Ra_F$, which serves as a proxy for energy input from both the bottom and top boundary.
The Nusselt numbers of RB convection and HC (defined in (3.3) and (3.7)) are good indicators of the flow regime; they deviate quickly from the classical scaling laws of RB convection and HC in the regime where they are not relevant, i.e. 
$\varLambda \gg 10^{-2}$ for 
$Nu_{RB}$ and 
$\varLambda \ll 10^{-2}$ for 
$Nu^{\chi }_{HC}$. The characteristic length scale 
$\ell$ derived from the auto-correlation function (§ 3.4) is another excellent indicator as it is equal to 1 in the RB regime and 
$\varGamma$ in the HC regime. The Reynolds number is not a good indicator for our set of simulations, however, because it is primarily controlled by 
$Ra_{F}$, the intensity of buoyancy driving, with small variations only due to 
$\varLambda$ (flow topology) and 
$\varGamma$.
The fact that the transition does not depend on 
$\varGamma$ (at least for 
$\varGamma \geq 4$, as considered in this study) means that the competition mechanism is mostly local and does not depend on the horizontal extent nor the number of convective rolls present. As we observe intense mean flows near the transition, i.e. more intense than those observed in pure RB or HC regimes at similar parameter values (the results of figure 14(b) are for 
$Ra_F=10^8$ but are qualitatively representative of simulations using other 
$Ra_F$), we suspect that the Reynolds stresses originating from the convective rolls are enhanced by the imposed horizontal temperature gradient (by favouring counter-clockwise rolls), leading to efficient mean flow amplification. Research on horizontal mean flows coupling with RB convection cells in periodic domains has a long history (Thompson Reference Thompson1970; Krishnamurti & Howard Reference Krishnamurti and Howard1981; Busse Reference Busse1983). For instance, arrays of convective rolls are known to be unstable to the so-called shearing instability (Hughes & Proctor Reference Hughes and Proctor1990; Rucklidge & Matthews Reference Rucklidge and Matthews1996; Goluskin et al. Reference Goluskin, Johnston, Flierl and Spiegel2014) when both horizontal plates are stress free, whereby Reynolds stresses positively couple with the roll tilt induced by the emerging mean flow. We do not expect to directly observe this shear instability, especially for 
$\varLambda =0$, because of our closed domain and no-slip boundary conditions, which are known to inhibit the instability for elongated domains (Fitzgerald & Farrell Reference Fitzgerald and Farrell2014; Van Der Poel et al. Reference Van Der Poel, Ostilla-Mónico, Verzicco and Lohse2014). However, for 
$\varLambda >0$, we suspect that the imposed symmetry breaking in the horizontal direction activates a mechanism reminiscent of the shearing instability, which might explain why a weak horizontal temperature gradient can easily disrupt the convective rolls. The subtle interplay between the horizontal mean flow triggered by the imposed temperature gradient along the top boundary and the convective rolls most likely underpins the origin of the transition from RB convection to mixed RBH convection, as well as the complex bi-stable and bi-modal dynamics observed at 
$\varLambda \approx 10^{-2}$. Therefore, it would be of interest to revisit the secondary instability of convection rolls in the presence of an imposed horizontal temperature gradient driving a mean flow using, e.g., a weakly nonlinear theory in future work.
We have shown that the system is multi-stable near the transition. For 
$Ra_F=10^6$, we have found at least two stable branches for 
$0.01\leq \varLambda \leq 0.02$. Each branch corresponds to a different number of overturning cells. Interestingly, we found that multiple flow states can also exist for a fixed number of rolls, as the spatial organisation of the rolls can differ substantially between cases. For larger 
$Ra_F=10^8$, the bi-stability seems to disappear, at least for simulations lasting only a few diffusive time scales, possibly because turbulent fluctuations force transitions between the two states so efficiently that their basin of attraction overlap. For 
$Ra_F=10^8$ and 
$\varLambda =10^{-2}$ the dynamics is better described as bi-modal, i.e. akin to RB convection with bursts of HC.
Obviously, the independence of the transition observed at 
$\varLambda =10^{-2}$ with 
$Ra_F$ is rigorously valid only for 
$10^6 \leq Ra_F\leq 10^9$ as considered in this study. Exploring the RBH dynamics at larger 
$Ra_F$ and in three dimensions is thus required to extend the applicability of our results to geophysical fluids, such as Earth's atmosphere and subglacial lakes. The Prandtl number for the latter is 
$Pr=O(10)$, hence is much larger than 
$Pr=1$, but its effect on the transition may be limited. Indeed, increasing 
$Pr$ for 
$Pr \geq 1$ has little effect on 
$Nu_{RB}$ and 
$Nu^{\chi }_{HC}$, whereas 
$Re$ decreases almost like 
$Pr^{-1}$ in both RB convection and HC, at least for moderate Rayleigh numbers (Shishkina & Wagner Reference Shishkina and Wagner2016; Li et al. Reference Li, He, Tian, Hao, Huang, Li, He, Tian, Hao and Huang2021). Thus, we hypothesise that increasing 
$Pr$ may slow down fluid motions without changing the type of convection. That being said, we note that the burst dynamics in RB convection between stress-free plates has been observed to disappear for large enough 
$Pr$ (Goluskin et al. Reference Goluskin, Johnston, Flierl and Spiegel2014), suggesting that the upper bound of the RB regime (rightmost edge of purple shading in, e.g., figure 11) may be in fact sensitive to 
$Pr$.
The dual RBH convection problem is likely to receive increased attention in the coming years because subglacial lakes will soon be explored and monitored (including, notably, lake CECs; see Rivera et al. Reference Rivera, Uribe, Zamora and Oberreuter2015) and because the contribution of geothermal heating to the ocean abyssal stratification and circulation is now known to be significant (Mashayek et al. Reference Mashayek, Ferrari, Vettoretti and Peltier2013; De Lavergne et al. Reference De Lavergne, Madec, Le Sommer, Nurser and Naveira Garabato2016). For subglacial lakes, our results suggest that RB convection should dominate since 
$\varLambda \ll 10^{-2}$ for realistic ice–water interface slopes (Appendix B). It is unclear whether RB convection will transition toward HC at lower or higher 
$\varLambda$ in three dimensions than in two dimensions. Preliminary three-dimensional simulations suggest a transition around 
$\varLambda \approx 10^{-2}$ again, which means that our results may be (at least qualitatively) informative for real systems in spite of the two-dimensional limitation. Non-rectangular geometry and Earth's rotation, whose effect on subglacial lake dynamics remains unclear (Couston & Siegert Reference Couston and Siegert2021), are other effects not considered in this work that may be more favourable to HC, hence should be investigated in future studies. The nonlinear equation of state of freshwater is also an important physical ingredient that can completely modify the fluid dynamics of subglacial lakes (Couston Reference Couston2021; Olsthoorn, Tedford & Lawrence Reference Olsthoorn, Tedford and Lawrence2021), which would be interesting to consider in RBH convection.
Despite the numerous approximations of our work, we provide in closing some thoughts on the implications of the fluid dynamics regime (considering both RB convection and HC) for the future exploration of subglacial lakes. For simplicity, our discussion assumes a flat water–bedrock interface, constant and positive thermal expansion coefficient and neglects the effect of rotation. In addition to measuring flow velocities, temperature and salinity along vertical profiles, future explorations of subglacial lakes will most likely investigate populations of suspended particulates (including microorganisms) through direct sampling of the water column and analysis of accreted ice (which may host particulates due to freezing of the lake water onto the ice ceiling), and characterise past climates by analyzing sediment cores extracted from the lake bed (Hodgson et al. Reference Hodgson, Roberts, Bentley, Carmichael, Smith, Verleyen, Vyverman, Geissler, Leng and Sanderson2009). In lakes dominated by HC, vertical profiles should be ideally performed at both ends of the lake in order to probe both the down-welling and up-welling regions. The former may be affected by subglacial water discharge from upstream (if any) whereas the latter is best to search for suspended particulates, whether in the water column or in the ice above. Sediments brought in the lake via upstream subglacial discharge will likely accumulate below the down-welling region. Dating sediments from continuous cores (resulting from successive layer depositions) is typically easier than from non-continuous cores (Hodgson et al. Reference Hodgson, Roberts, Bentley, Carmichael, Smith, Verleyen, Vyverman, Geissler, Leng and Sanderson2009). Thus, initial dating could be based on cores extracted from below the down-welling region, henceforth enabling dating sediment cores extracted from below the up-welling region that are more ancient. In lakes dominated by RB convection, the lack of upstream/downstream asymmetry combined with the possible migration of RB cells means that the drilling strategy can be based on constraints other than the water circulation, which homogenises lake conditions. In this case, four drilling sites separated (horizontally) by a quarter water depth along the main direction of the lake would most likely provide good coverage of up-welling and down-welling regions.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2022.613.
Acknowledgements
We gratefully acknowledge support from the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for the computing resources. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high-performance computing resources. This work was performed using HPC/AI resources from GENCI-IDRIS/TGCC (Grant 2021-A0120407543). The authors would like to thank Adrien Villaret for his careful review and Andy Smith and Keith Makinson for useful discussions about subglacial lake exploration.
Funding
The authors are grateful to the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) Lyon for its financial support within the program “Investissements d'Avenir” of the French government operated by the National Research Agency (ANR).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Simulations details
The physical and numerical parameters of all simulations that led to the results discussed in §§ 3.1–3.4 are provided in table 2. The details of the simulations used for the multi-stability analysis are discussed in § 3.5 and summarised on the last line of table 2.
Physical and numerical parameters of the simulations, with a dash denoting same value as in the row above for readability. Here 
$\varGamma$ is the aspect ratio, 
$Ra_F$ is the flux-based Rayleigh number, 
$\varLambda$ is the heat flux ratio, 
$Ra_L$ is the horizontal Rayleigh number, 
$n_z$ is the number of elements in the vertical direction, 
$l_d$ is the polynomial order and 
$dt$ is the typical time step at statistical steady state. Most simulations include a bottom heat flux, also referred to as geothermal heating (Geo. Flux), and are run via direct numerical simulation (DNS) from start to finish; some of the most demanding simulations are run via large-eddy simulation (LES) using the filtering approach described in Fischer & Mullen (Reference Fischer and Mullen2001) from 
$t=0$ to 
$t=1$ and then restarted with numerical parameters as indicated in the table and integrated from 
$t=1$ to 
$t=1.2$ via DNS. The last line shows the parameters of the simulations discussed in § 3.5, which focus on bi-stability.
Appendix B. Parameter range relevant to subglacial lakes
In this appendix, we motivate the choice of small 
$\varLambda$ considered in the paper. We are particularly interested in subglacial lakes, which lie beneath several kilometres of ice in Antarctica and Greenland and are subject to geothermal heating and horizontal temperature gradients along their bottom and top boundaries, respectively (Livingstone et al. Reference Livingstone2022). The latter arises when variable ice thickness above the lake water results in a tilted ice–water interface by hydrostatic equilibrium. The temperature of freezing varies along a tilted ice–water interface because it depends on local pressure, which is larger where the ice is thicker (Thoma et al. Reference Thoma, Grosfeld, Smith and Mayer2010). A linear approximation for the freezing temperature of freshwater with no dissolved air as a function of local ice pressure 
$p_i$ (in decibars) or ice thickness 
$h_i=10^4p_i/(\rho _ig)$ (in metres), with 
$\rho _i=917\,{\rm kg}\,{\rm m}^{-3}$ the ice density and 
$g=9.81\,{\rm m}\,{\rm s}^{-2}$ the gravity, is
\begin{equation} T_f = 6.67\times 10^{{-}2}-8.21\times 10^{{-}4} p_i = 6.67\times 10^{{-}2}-9.12\times 10^{{-}4} h_i, \end{equation}
as obtained from best-fit with the exact freezing temperature calculated with the Python package gsw (McDougall & Barker Reference McDougall and Barker2011) for 
$0< p_i<5000$ dbar (
$h_i$ ranging from 0 to 5558 m). Equation (B1) shows that 
$T_f$ decreases with 
$p_i$ or 
$h_i$ (
$T_f\approx -4 \, ^{\circ }{\rm C}$ when 
$p_i=5000$ dbar) and yields a horizontal temperature gradient equal to
\begin{equation} \lambda = \dfrac{{\rm d}T_f}{{\rm d} x} = 9.12\times 10^{{-}4} \gamma, \end{equation}
for a subglacial lake with an ice–water interface slope 
${\rm d}h_i/{{\rm d} x} = -\gamma$, i.e. thicker for smaller 
$x$, as assumed in this study. We find that 
$\gamma =0.016$, 0.003, 0.03, 0.002 and 0.003, for the five well-documented subglacial lakes reported in Couston & Siegert (Reference Couston and Siegert2021), i.e. CECs, SPL, Ellsworth, Vostok and Concordia, respectively. This range of 
$\gamma$ values yields approximately 
$3\times 10^{-5}<\varLambda =k\lambda /F<3\times 10^{-4}$, which lies within the range 
$0\leq \varLambda \leq 1$ explored in this study, with 
$k=0.56\, {\rm W}\,({\rm m}\,^{\circ }{\rm C})^{-1}$ the thermal conductivity of water and 
$F=50\,{\rm mW}\,{\rm m}^{-2}$ an average value for Earth's geothermal flux. Thus, for subglacial lakes, 
$\varLambda \ll 10^{-2}$, unless the slope is of order unity (in which case 
$\varLambda \geq 10^{-2}$), which is most likely unrealistic. For completeness, the flux-based Rayleigh number of a subglacial lake below thick ice with water depth 
$H=100$ m, 
$F=50\,{\rm mW}\,{\rm m}^{-2}$, 
$\nu =1.7\times 10^{-6}\,{\rm m}^2\,{\rm s}^{-1}$, 
$\kappa =1.33\times 10^{-7}\,{\rm m}^2\,{\rm s}^{-1}$, 
$\alpha =10^{-4}\,^{\circ }{\rm C}^{-1}$ (typical values inferred from Couston Reference Couston2021) is 
$Ra_F\approx 3.9 \times 10^{16}$, which cannot be reached in numerical simulations, even in two dimensions, with present-day computational resources.
Appendix C. Scaling laws in the RB and HC regimes
In this section we comment on the scaling laws inferred from our simulations in light of previously published results (see table 1).
The 
$Re(Ra_F)$ scaling can be recast as a 
$Re(Ra_{\varDelta })$ scaling, where 
$Ra_{\varDelta }$ is the classical RB parameter based on the temperature difference 
$\varDelta$ between the top and bottom plates (averaged in 
$x$ and 
$t$ and computed a posteriori). To see this, we first remark that the different definitions of the problem parameters yield the relationship (in the limit 
$\varLambda \rightarrow 0$) 
$Nu_{RB}=Ra_F/Ra_{\varDelta }$ (Johnston & Doering Reference Johnston and Doering2009), such that 
$Ra_F=a_{RB}^{1/(1-b_{RB})}Ra_{\varDelta }^{b_{RB}/(1-b_{RB})}$. This yields the scaling law 
$Nu=0.316Ra_{\varDelta }^{0.236}$ using the pre-factor and exponent reported in table 1. The obtained 
$Nu(Ra_{\varDelta })$ scaling law is in relatively good agreement with Grossmann & Lohse (Reference Grossmann and Lohse2000)'s unifying theory: the exponent is slightly larger (and the pre-factor is slightly smaller) than that predicted by regime 
$II_u$, i.e. 1/5, which is expected for our low 
$1.8\times 10^5\leq Ra_{\varDelta }\leq 4.8\times 10^7$ numbers (see updated regime diagram in Stevens et al. Reference Stevens, van der Poel, Grossmann and Lohse2013, and note that our scaling law happens to be in better agreement with the 1/4 exponent of regime 
$I_u$ expected at larger 
$Pr>1$). We expect the difference between our scaling exponent and that of regime 
$II_u$ to be due to the two-dimensional flow assumption, which is known to affect 
$Nu_{RB}$, especially at low-to-moderate 
$Pr$ (Van Der Poel, Stevens & Lohse Reference Van Der Poel, Stevens and Lohse2013) and the relative proximity of regime 
$IV_u$ at slightly larger 
$Ra_{\varDelta }$, which features a 1/3 exponent. We note that our exponent is smaller (and the pre-factor is larger) than other two-dimensional flux-based RB studies (Johnston & Doering Reference Johnston and Doering2009; Couston Reference Couston2021) (compare 0.236 with 
$2/7\approx 0.286$) because we considered lower Rayleigh numbers. Replacing 
$Ra_F$ with 
$Ra_{\varDelta }$ in 
$Re(Ra_F)$ yields 
$Re=0.140Ra_{\varDelta }^{0.545}$, which is in good agreement with previously reported results provided that the 
$\sim Pr^{-0.82}$ dependence (Li et al. Reference Li, He, Tian, Hao, Huang, Li, He, Tian, Hao and Huang2021) and 
$\varGamma$ dependence (Van Der Poel et al. Reference Van Der Poel, Stevens, Sugiyama and Lohse2012) of 
$Re$ are taken into account.
Our HC scalings are in relatively good agreement with the theoretical predictions 
$Nu^{\chi }_{HC}\sim Ra_L^{0.2}$ and 
$Re\sim Ra_L^{0.4}$ based on laminar conditions, with 
$Re$ usually based on the peak velocity within the HC boundary layer (Rossby Reference Rossby1965; Hughes & Griffiths Reference Hughes and Griffiths2008; Sheard & King Reference Sheard and King2011) instead of the mean kinetic energy density. The slightly higher exponent for the Nusselt number 
$a_{HC}=0.224>0.2$ is likely due to the unsteadiness of the horizontal flow given our relatively large 
$Ra_L$ (Sheard & King Reference Sheard and King2011), while remaining substantially smaller than the upper bound 1/3, which is the exponent of the ultimate regime (Rocha et al. Reference Rocha, Bossy, Llewellyn Smith and Young2020a). Again, we note that the small dependence of 
$Nu^{\chi }_{HC}$ with 
$\varGamma$, which has not been observed previously (Sheard & King Reference Sheard and King2011), comes from the normalisation, which is based on the diffusive solution rather than on a purely dimensional argument (see equation (3.7)).
Appendix D. Dissipation of temperature variance and diffusive solution
Recently, Rocha et al. (Reference Rocha, Constantinou, Llewellyn Smith and Young2020b) showed that a compelling definition of the Nusselt number of HC should involve the dissipation of temperature (or buoyancy 
$b=-T$) variance, 
$\chi = \overline {\langle |\boldsymbol {\nabla } T|^2 \rangle }$, because it is related to the vertically averaged horizontal heat flux set up by the imposed heterogeneous temperature profile at (here) 
$z=1$. They further showed that 
$\chi =\overline {\langle T(z=1)\partial _z T(z=1) \rangle _x}$, such that it can be evaluated through a line integral, which is more tractable in laboratory experiments and numerical simulations than a volume integral as required by 
$\overline {\langle |\boldsymbol {\nabla } T|^2 \rangle }$. Here we derive a similar result for the case with geothermal heating. Multiplying (2.1c) by 
$T$ and rearranging yields
\begin{equation} \partial_t \left( \dfrac{T^2}{2} \right) + |\boldsymbol{\nabla} T|^2 ={-} \boldsymbol{\nabla}\boldsymbol{\cdot}\left( \frac{{\boldsymbol{u}} T^2}{2} \right) + \boldsymbol{\nabla}\boldsymbol{\cdot}\left( T\boldsymbol{\nabla} T \right). \end{equation}The first term becomes zero when time averaging at statistical steady state, whereas the third term becomes zero when performing a volume average because of the no-slip condition on the walls. Thus,
\begin{equation} \chi=\overline{\langle|\boldsymbol{\nabla} T|^2\rangle} = \overline{\langle\boldsymbol{\nabla}\boldsymbol{\cdot}\left( T\boldsymbol{\nabla} T \right)\rangle} = \overline{\langle T(z=1)\partial_z T(z=1) - T(z=0)\partial_z T(z=0) \rangle_x} > 0, \end{equation}
where the last equality is obtained after integrating in 
$z$ and enforcing no-flux conditions on the side walls. Our expression for 
$\chi$ differs from that of Rocha et al. (Reference Rocha, Constantinou, Llewellyn Smith and Young2020b) because of the second term on the right-hand side of (D2), which is non-zero here because our derivation takes into account the heat flux (
$\partial _z T=-1$) on the bottom boundary.
The diffusive solution in HC, which we derive in terms of dimensionless variables in this section, can be readily obtained using the method of separation of variables. For the surface temperature profile of (2.5), the diffusive solution without geothermal heating (
$\partial _z T=0$ at 
$z=0$) is
\begin{equation} T_{{diff}} = \dfrac{\varLambda\varGamma \sin \left( \dfrac{{\rm \pi} x}{\varGamma} \right) \cosh \left( \dfrac{{\rm \pi} z}{\varGamma} \right) }{2 \cosh \left( \dfrac{\rm \pi}{\varGamma} \right)} , \end{equation}
whereas, with geothermal heating (
$\partial _z T=-1$ at 
$z=0$), we simply add a linearly varying vertical profile, i.e.
\begin{equation} T_{{diff}} = \dfrac{\varLambda\varGamma \sin \left( \dfrac{{\rm \pi} x}{\varGamma} \right) \cosh \left( \dfrac{{\rm \pi} z}{\varGamma} \right) }{2 \cosh \left( \dfrac{\rm \pi}{\varGamma} \right)} + 1-z. \end{equation}The denominator in (3.7) based on (D4) then reads
\begin{align}
\chi_{{diff}} = \langle \partial_z T_{{diff}}(z=1)
T_{{diff}}(z=1) - \partial_z T_{{diff}}(z=0)
T_{{diff}}(z=0)\rangle_x \!=\!
\dfrac{{\rm \pi}\varLambda^2\varGamma}{8}\tanh\left(\!
\dfrac{\rm \pi}{\varGamma}\!\right) +1,
\end{align}
which tends asymptotically to 
${\rm \pi}^2\varLambda ^2/8+1$ for large 
$\varGamma$ with the 
$+1$ term being inherited from geothermal heating (which is thus discarded for the calculation of 
$Nu^{\chi }_{HC}$ in pure HC simulations). Importantly, geothermal heating has a negligible effect on 
$\chi _{{diff}}$ only if 
$\varLambda \gg 1$, since in this case the 
$+1$ term can be discarded.
Appendix E. Nusselt numbers of HC
Figure 15(a–c) show 
$Nu^{{abs}}_{HC}$, 
$Nu^{{half}}_{HC}$ and 
$Nu^{\chi }_{HC}$ (defined in (3.4)–(3.7)) as functions of 
$Ra_L$. Pure HC simulation results are indistinguishable from the scaling law (solid line) derived from the 
$Nu^{\chi }_{HC}$ data, which means that the heat transfer efficiency by HC is the same for all three definitions of the Nusselt number in the HC-dominating limit. Differences exist between RBH simulation results (compare, e.g., light blue, orange and green circles in figure 15b,c). However, for 
$\varLambda \gg 10^{-2}$, these differences are due to the diffusive normalisation (i.e. appearing at the denominator), whose dependence on geothermal heating remains significant for 
$\varLambda \leq O(1)$ (unlike the numerator) and depends on the choice of Nusselt number definition. When renormalising the numerator in the definition of each Nusselt number by the diffusive solution without geothermal heating (i.e. considering 
$F=0$), those differences (at large 
$\varLambda$) are removed, as can be seen from figure 15(d–f) in which all renormalised Nusselt numbers divided by the scaling law 
$a_{HC}Ra_L^{b_{HC}}$ converge toward unity. Note that differences due to aspect ratio 
$\varGamma$ are also due to the diffusive normalisation, as explained at the end of § 3.3.
(a–c) Three Nusselt numbers defined in (3.4)–(3.7) as functions of 
$Ra_L$ with the solid black line showing the scaling law derived from the 
$Nu^{\chi }_{HC}$ data obtained for pure HC simulations. (d–f) Same Nusselt numbers divided by the scaling law and renormalised by the diffusive solution derived without geothermal heating (
$F=0$); e.g. 
$Nu^{\chi }_{HC}(F=0)=Nu^{\chi }_{HC}\times \chi _{{diff}}/\chi _{{diff}}(F=0)$.
Open access
