Louis-Alexandre Couston
email: louis.couston@ens-lyon.fr

Joseph Nandaha

Benjamin Favier

Competition between Rayleigh-Bénard and horizontal convection

Keywords: geological and geophysical flows, turbulent convection, convection in cavities. 1. Introduction

We investigate the dynamics of a fluid layer subject to an imposed bottom heat flux and a top monotonically-increasing temperature profile driving horizontal convection. We use direct numerical simulations and consider a large range of flux-based Rayleigh numbers 10 6 𝑅𝑎 𝐹 10 9 and imposed top horizontal to bottom vertical heat flux ratios 0 Λ 1. The fluid domain is a closed two-dimensional box with aspect ratio 4 Γ 16 and we consider no-slip boundaries and adiabatic side walls. We demonstrate a regime transition from Rayleigh-Bénard convection (RB) to horizontal convection (HC) at Λ ≈ 10 -2 , which is independent of 𝑅𝑎 𝐹 and Γ. At small Λ, the flow is organized in multiple overturning cells with approximately unit aspect ratio, while at large Λ a single cell is obtained. The RB-relevant Nusselt number scaling with 𝑅𝑎 𝐹 and the HC-relevant Nusselt number scaling with the horizontal Rayleigh number 𝑅𝑎 𝐿 = 𝑅𝑎 𝐹 ΛΓ 4 are in good agreement with previous results from classical RB convection and HC studies in the limit Λ 10 -2 and Λ 10 -2 , respectively. We demonstrate that the system is multi-stable near the transition Λ ≈ 10 -2 , i.e. the exact number of cells not only depends on Λ but also on the system's history. Our results suggest that subglacial lakes, which motivated this study, are likely to be dominated by RB convection, unless the slope of the ice-water interface, which controls the horizontal temperature gradient via the pressure-dependence of the freezing point, is greater than unity.

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 et al. 2009), 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 2008). The behavior of RB convection and HC is typically examined through the scaling of the Reynolds number 𝑅𝑒, which is a proxy for fluid velocity, and Nusselt number 𝑁𝑢, which is a proxy for the efficiency of heat transport by convection, with the control parameters, including notably, the Rayleigh number 𝑅𝑎, 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 𝑅𝑒 and 𝑁𝑢 as functions of 𝑅𝑎, as well as the Prandtl number 𝑃𝑟, which compares viscosity to thermal diffusivity, for both RB convection (Grossmann & Lohse 2000;Ahlers et al. 2009) and HC (Mullarney et al. 2004;Shishkina et al. 2016). However, many open questions remain. For instance, recent studies on RB convection aim to address the existence of the ultimate regime (Zhu et al. 2018) or its large-scale organization (Pandey et al. 2018;Wang et al. 2020;Vieweg et al. 2021), while 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 2002;Scotti & White 2011;Gayen et al. 2014;Passaggia et al. 2017;Rocha et al. 2020b) and the effect of rotation (Vreugdenhil et al. 2019;Gayen & Griffiths 2022).

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. For instance, RB convection and HC may concomitantly control the fluid dynamics within subglacial lakes in Greenland and Antarctica, which impact the dynamics of ice sheets and most likely host extremophiles of interest to astrobiology (Cockell et al. 2011;Livingstone et al. 2022). Subglacial lakes, which are typically fresh except when active or close to grounding lines (Priscu et al. 2021), 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. 2010). The fluid dynamics of planetary oceans may also be affected by both RB convection and HC, as oceans receive solar radiations that vary with latitude and typically experience geothermal heating (Wang et al. 2016). 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 1999;Mound & Davies 2017).

Past studies of the dual Rayleigh-Bénard-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 parameterized subgrid-scale processes (Thoma et al. 2007(Thoma et al. , 2009) ) and a laboratory experiment of lake Vostok (Wells & Wettlaufer 2008). There is yet no consensus on the type of fluid motions expected in subglacial lakes: Thoma et al. (2007,2009) predict HC-driven large-scale circulations in lakes Vostok and Concordia affected by rotation; Wells & Wettlaufer (2008) found that the dynamics in lake Vostok is better described by rotating RB convection, i.e. dominated by multiple columnar vortices; and, Couston & Siegert (2021) predicts non-rotating RB convection in most subglacial lakes. In open ocean Competition between Rayleigh-Bénard and horizontal convection 3 research, several studies (Hofmann & Maqueda 2009) 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 idealized numerical simulations (Mullarney et al. 2006) and laboratory experiments (Wang et al. 2016), wherein the MOC is driven by a horizontallyvarying 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 Γ = 6, Mullarney et al. (2006) 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 Γ = 1, Wang et al. (2016) 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 characterize 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 𝜆 along the top boundary multiplied by the thermal conductivity 𝑘 (which yields a horizontal heat flux), to the bottom heat flux 𝐹, which we write as Λ = 𝑘𝜆/𝐹. To this end, we run a large number of two-dimensional numerical simulations with variable bottom and surface forcing (listed in table 2), 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 manuscript is organized as follows.

In section § 2 we introduce the dimensional and dimensionless governing equations as well as the numerical method. In section § 3 we show results highlighting the regime transition between RB-like convection at Λ 10 -2 and horizontal convection at Λ 10 -2 , and we demonstrate the existence of multiple flow states for the same set of problem parameters. Finally, we conclude in § 4.

Problem formulation

We consider a two-dimensional rectangular fluid domain with Cartesian coordinates (𝑥, 𝑧) centred on the bottom boundary; 𝒆 𝑧 is the upward-pointing unit vector of the 𝑧 axis, which is opposite to gravity (figure 1(a)). We denote by 𝐻 and 𝐿 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 𝒖, pressure 𝑝 and temperature 𝑇 are governed by the Navier-Stokes equations in the Boussinesq approximation, i.e.

𝜕 𝑡 𝒖 -𝜈∇ 2 𝒖 + ∇( 𝑝/𝜌 0 ) = -(𝒖 • ∇) 𝒖 -(𝜌 /𝜌 0 )𝑔e 𝑧 , (2.1a) ∇ • 𝒖 = 0, (2.1b) 𝜕 𝑡 𝑇 -𝜅∇ 2 𝑇 = -(𝒖 • ∇) 𝑇, (2.1c)
where 𝜌 0 is the reference fluid density and 𝜌 is the density anomaly; 𝑔 is surface gravity, 𝜕 𝑡 denotes time derivative and ∇ 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 𝜌 = -𝜌 0 𝛼(𝑇 -𝑇 0 ) with 𝑇 0 the reference temperature and with the thermal expansion coefficient 𝛼, along with kinematic viscosity 𝜈 and thermal diffusivity 𝜅, 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

𝒖(𝑧 = 0) = 𝒖(𝑧 = 𝐻) = 𝒖 (𝑥 = ±𝐿/2) = 0, (2.2a) 𝜕 𝑧 𝑇 (𝑧 = 0) = - 𝐹 𝑘 , 𝑇 (𝑧 = 𝐻) = 𝑇 𝑓 (𝑥) = 𝜆𝐿 2 sin 𝜋𝑥 𝐿 , 𝜕 𝑥 𝑇 (𝑥 = ±𝐿/2) = 0. (2.2b)
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 horizontal convection, 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 𝑓 to denote the top temperature 𝑇 𝑓 , 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.

In order to identify the minimum number of independent parameters and explore their effect on the fluid dynamics, we non-dimensionalize the governing equations (2.1) and boundary conditions (2.2). We use the water depth 𝐻 as characteristic length scale, the diffusive time 𝜏 𝜅 = 𝐻 2 /𝜅 as time scale, the velocity 𝑢 𝜅 = 𝐻/𝜏 𝜅 as velocity scale, the temperature difference due to geothermal heating Δ = 𝐹𝐻/𝑘 as temperature scale, and the pressure 𝑝 𝜅 = 𝜌 0 𝑢 2 𝜅 as pressure scale. Using 𝑇 0 and 𝑝 0 + 𝜌 0 𝑔(𝐻 -𝑧) as temperature and pressure gauges (𝑝 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

(𝑥, 𝑧) = 𝐻 ( 𝑥, 𝑧), 𝑡 = 𝜏 𝜅 𝑡, 𝑢 = 𝑢 𝜅 𝑢, 𝑝 = 𝑝 0 + 𝜌 0 𝑔(𝐻 -𝑧) + 𝑝 𝜏 𝑝, 𝑇 = 𝑇 0 + Δ 𝑇 .
(2.3) Substituting (2.3) into (2.1) and (2.2) yields a set of dimensionless equations and boundary

Focus on Fluids articles must not exceed this page length

Competition between Rayleigh-Bénard and horizontal convection 5 conditions, which, dropping tildes, read

𝜕 𝑡 𝒖 -𝑃𝑟∇ 2 𝒖 + ∇𝑝 = -(𝒖 • ∇) 𝒖 + 𝑃𝑟 𝑅𝑎 𝐹 𝑇e 𝑧 , (2.4a) ∇ • 𝒖 = 0, (2.4b) 𝜕 𝑡 𝑇 -∇ 2 𝑇 = -(𝒖 • ∇) 𝑇, (2.4c)
and [START_REF]An overlapping Schwarz method for spectral element solution of the incompressible Navier-Stokes equations[END_REF]Deville et al. 2002), which has been extensively used recently in thermal convection studies (Scheel et al. 2013;Léard et al. 2020). The governing equations are cast into weak form and discretised in space by the Galerkin approximation. The Cartesian domain is discretised using 𝑛 𝑧 elements in the vertical direction and Γ𝑛 𝑧 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 𝑙 𝑑 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/2𝑙 𝑑 . Convergence has been tested by gradually increasing the polynomial order for a fixed number of elements (Scheel et al. 2013). 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.

𝒖(𝑧 = 0) = 𝒖(𝑧 = 1) = 𝒖(𝑥 = ±Γ/2) = 0, (2.5a) 𝜕 𝑧 𝑇 (𝑧 = 0) = -1, 𝑇 (𝑧 = 1) = ΛΓ 2 sin 𝜋𝑥 Γ , 𝜕 𝑥 𝑇 (𝑥 = ±Γ/2) = 0, ( 2
The initial state is motionless and has uniform temperature distribution superimposed with low-amplitude background noise (except in section § 3.5 where we perform continuation). We explore the fluid dynamics in the (𝑅𝑎 𝐹 , Λ) parameter space as shown in figure 1(b), i.e. with 10 6 𝑅𝑎 𝐹 10 9 and 0 Λ 1. Simulations with Λ = 0 yield canonical Rayleigh-Bénard results; for completeness we also run simulations without geothermal flux, i.e. setting 𝜕 𝑧 𝑇 = 0 at 𝑧 = 0, which yield canonical horizontal convection results. For simplicity, we set 𝑃𝑟 = 1 in all simulations. Most results are shown for Γ = 8, though we also ran simulations with Γ = 4 and 16 to partially assess the effect of the aspect ratio; 4 Γ 16 spans enclosure aspect ratios that are commonly used in HC studies (either experimental or numerical), including Γ = 6.2 (Mullarney et al. 2004;Gayen et al. 2013Gayen et al. , 2014;;Tsai et al. 2020) and Γ = 10 (e.g. Shishkina & Wagner 2016). Note that different colors highlight different 𝑅𝑎 𝐹 in figure 1(b), whereas squares, circles and diamonds denote results obtained for Γ = 4, 8 and 16, respectively; stars show the results of pure (no geothermal flux) horizontal convection simulations. As we will show, Λ is a better indicator of the regime dynamics than 𝑅𝑎 𝐿 /𝑅𝑎 𝐹 . The range of 0 Λ 1 values is motivated by subglacial lakes, as explained in Appendix B.

Results

Flow regimes

We first show in figures 2 and 3 snapshots of the velocity and temperature fields at statistical steady state for 𝑅𝑎 𝐹 = 10 8 and Γ = 8 with Λ increasing from top to bottom. The top plot in figure 2 shows the velocity field obtained for Λ = 0, which is the canonical RB case. The flow is organized in 4 pairs of counter-rotating rolls, with characteristic length scale approximately equal to twice the domain height, as expected from linear stability analysis (Chandrasekhar 1961). As Λ 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 Λ = 10 -2 (3rd plot from the top in figure 2), only 3 pairs of counter-rotating cells co-exist, and the 3 counter-clockwise cells are elongated because the upper leftward flow branch they support is enhanced by the top boundary. As Λ increases above 10 -2 , the flow becomes dominated by horizontal convection, 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 Λ = 1 with (2nd 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.

The temperature field in figure 3 also highlights the existence of RB-like cells for small Λ, which merge at intermediate Λ values. For Λ = 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 Λ = 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 Λ 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.

Mean temperature and Reynolds number

In sections § 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 𝑋 𝑥 to denote 𝑥-averaged variables, 𝑋 to denote volume-averaged variables, and over-line 𝑋 to denote temporal averages at statistical steady state (typically from 𝑡 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.

Λ = 0 Λ = 10 -3 Λ = 10 -2 Λ = 10 -1 Λ = 1 Pure HC (no geothermal flux) with 𝑅𝑎 𝐿 ≈ 4 × 10 11
Figure 4 shows the temporal evolution of the volume-averaged temperature 𝑇 (top row) and Reynolds number 𝑅𝑒 (bottom row). In this study, we use the kinetic energy density to construct the Reynolds number, i.e. such that

𝑅𝑒 ≡ √ 𝑢 2 + 𝑤 2 , (3.1)
which is close to the Reynolds number based on the velocity root mean square (not shown).

L.-A. Couston, J. Nandaha, B. Favier

Λ = 0 Λ = 10 -3 Λ = 10 -2 Λ = 10 -1 Λ = 1
Pure HC (no geothermal flux) with 𝑅𝑎 𝐿 ≈ 4 × 10 11 The time-averaged Reynolds number is

𝑅𝑒 ≡ √ 𝑢 2 + 𝑤 2 . (3.2)
Both 𝑇 and 𝑅𝑒 display a sharp transient followed by a statistical steady state (small fluctuations around a mean value independent of time) at approximately 𝑡 = 0.5; thus, here we use (conservatively) 𝑡 = 1 as the initial time for time averaging of output variables representative of the statistical steady state. For small Λ (figure 4(a)), the mean temperature increases from 0 up to a small but positive value as a result of geothermal heating, which dominates over horizontal convection. shown in figure 4(c)). The effect of Λ (color intensity) and Γ (line thickness) on 𝑅𝑒, which is relatively weak for our set of simulations (although clearly visible for 𝑅𝑎 𝐹 = 10 7 and 𝑅𝑎 𝐹 = 10 8 in figure 4), is commented on in greater details when discussing figure 5(d).

We show in figure 5 the time-averaged mean temperature 𝑇 and Reynolds number 𝑅𝑒 at statistical steady state as functions of the problem parameters. Figure 5(a) shows that the mean temperature decreases quickly with increasing Λ 10 -2 because down-welling below the cold top boundary (𝑇 (𝑥 = -Γ/2) = -ΛΓ/2 shown by black markers) becomes sufficiently strong to lower the bulk temperature. The mean temperature also decreases with increasing 𝑅𝑎 𝐹 (blue markers above orange, green and red markers) because increasing geothermal heating makes mixing more efficient, thus lowering temperature differences, while increasing horizontal convection increases mixing from the cold region of the top boundary.

Figures 5(b) shows that the power law curve 𝑅𝑒 = 𝑐 𝑅𝐵 𝑅𝑎 𝑑 𝑅𝐵 𝐹 (black solid line) provides a good prediction for the Reynolds number as a function of 𝑅𝑎 𝐹 for most simulations. Here, pre-factor 𝑐 𝑅𝐵 and exponent 𝑑 𝑅𝐵 are obtained from best fit with the results for Λ = 0 and Γ = 8 (see table 1 and Appendix C for a list and discussion of all pre-factors and exponents mentioned in the paper). The dependence of 𝑅𝑒(𝑅𝑎 𝐹 ) with Λ (and Γ) is small, especially at low Λ < 10 -2 (dark colors), which means that a small horizontal temperature gradient (or change of aspect ratio) has limited effect on the intensity of RBdriven flows. Conversely, figure 5(c) demonstrates that there is a wide spread of 𝑅𝑒(𝑅𝑎 𝐿 ) between simulations, even for relatively large Λ ∼ 0.1 (light colors), which means that the circulation is almost always affected by the bottom heat flux (if not driven by it), even in the HC-dominated regime obtained for Λ > 𝑂 (10 -2 ) (as we will show). Accordingly, the power law curve 𝑅𝑒 = 𝑐 𝐻 𝐶 𝑅𝑎 𝑑 𝐻𝐶 𝐿 (black solid line) obtained from pure HC results (shown by the stars) predicts 𝑅𝑒 accurately for large Λ 10 -1 only.

We plot the Reynolds number compensated by the RB scaling in figure 5(d) in order to highlight the effect of Λ and Γ on 𝑅𝑒. The spread of 𝑅𝑒 with Γ and Λ is overall small (less

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 Rayleigh-Bénard or horizontal convection. The Nusselt number of RB convection is defined as

Rapids articles must not exceed this page length

𝑁𝑢 𝑅𝐵 ≡ 𝐹𝐻 𝑘 𝑇 dim (𝑧 = 0) 𝑥 -min(𝑇 dim ) = 1 𝑇 (𝑧 = 0) 𝑥 -min(𝑇) , (3.3)
where superscript dim denotes a dimensional variable and min(𝑇) = -ΛΓ/2 is the coldest temperature achieved on the top boundary. As is customary in classical RB studies, 𝑁𝑢 𝑅𝐵 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(𝑇) instead of the mean 𝑇 = 0 value as a gauge for the top temperature in the denominator in equation (3.3) to ensure 𝑁𝑢 𝑅𝐵 > 0, since 𝑇 (𝑧 = 0) 𝑥 becomes negative for large-enough Λ (figure 5(a)). We note that the use of min(𝑇) affects the definition of the heat flux of the diffusive state; however, in the RB regime, i.e. for Λ 10 -2 , min(𝑇) is always much smaller than the mean bottom temperature, such that 𝑁𝑢 𝑅𝐵 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 while 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 RB-HC and HC-dominated simulations with 𝑅𝑎 𝐹 = 10 8 : the RB results (Λ = 0) show an oscillatory pattern of heat flux linked to the underlying overturning cells, the intermediate Λ = 10 -2 value leads to skewed oscillations and the large Λ = 1 value yields a much-larger monotonic and anti-symmetric heat flux pattern with respect to the middle of the domain 𝑥 = 0. Figure 6

𝑁𝑢 𝜒 𝐻 𝐶 = 𝑇 (𝑧 = 1)𝜕 𝑧 𝑇 (𝑧 = 1) -𝑇 (𝑧 = 0)𝜕 𝑧 𝑇 (𝑧 = 0) 𝑥 𝑇 diff (𝑧 = 1)𝜕 𝑧 𝑇 diff (𝑧 = 1) -𝑇 diff (𝑧 = 0)𝜕 𝑧 𝑇 diff (𝑧 = 0) 𝑥 , (3.7)
with an analytical expression for the denominator given in equation (D 5). We would like to note that while 𝑁𝑢 𝜒 𝐻 𝐶 = 𝜒/𝜒 diff is thermodynamically compelling, it is very close to 𝑁𝑢 abs

𝐻 𝐶

and 𝑁𝑢 half 𝐻 𝐶 for pure HC simulations. For mixed RBH simulations, differences exist but are due at leading order to the diffusive normalization (see details in Appendix E).

Figure 7(a) shows 𝑁𝑢 𝑅𝐵 (equation (3.3)) as a function of 𝑅𝑎 𝐹 . Simulation results obtained for Λ 10 -2 all collapse very well on the power law curve shown by the black solid line, which was obtained from best fit for Λ = 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 𝑁𝑢 𝑅𝐵 number in figure 7 (grey area), which we will show marks the transition from RB to horizontal convection. Note that the decrease of 𝑁𝑢 𝑅𝐵 with Λ is primarily the result of an increase of the denominator in equation (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). 10 -2 , which is exactly the scaling of the diffusive normalization 1/𝜒 diff -1 ∝ Λ 2 once Taylor expanded in the small Λ 2 limit (cf. equation (D 5)). Second, figure 8(c) demonstrates that replacing 𝜒 diff with 𝜒 dim = 𝜋 2 Λ 2 /8, which is dimensionally-equivalent but discards geothermal heating and approximates tanh (𝜋/Γ) ≈ 𝜋/Γ (large Γ limit), in the definition of 𝑁𝑢 𝜒 𝐻 𝐶 (equation (3.7)) yields a perfect overlap of all simulation results obtained for large Λ 10 -2 with the power law 𝑎 𝐻 𝐶 𝑅𝑎 𝑏 𝐻𝐶 𝐿 . Ultimately, the spread of 𝑁𝑢 𝜒 𝐻 𝐶 with Λ and Γ at large Λ is due to the diffusive normalization 𝜒 diff (denominator in (3.7)), not 𝜒, because 𝜒 diff remains sensitive to aspect ratio (Γ) and flow topology (Λ) for a much wider range of parameters than 𝜒 (as is the case for other definitions of the Nusselt number, see Appendix E). This result is in agreement with previous studies (Sheard & King 2011) who found no Γ dependence using a flux-based definition of the Nusselt number normalized by Λ rather than the diffusive solution in the convection-dominated regime, which spans a large range of 𝑅𝑎 𝐿 encompassing our simulation parameters (see also Hossain et al. 2019, for the effect of Γ 1 on the transition from diffusion-to convection-dominated HC).

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 Λ ≈ 10 -2 marks the transition from RB convection to HC. Since the leftward flow near the top boundary is an emblematic feature of HC, we use the variations in -𝑢(𝑧 = 0.9) 𝑥 is of the same order as 𝑅𝑒 for relatively large Λ.

𝑥 of the horizontal velocity at 𝑧 = 0.9 as diagnostic. We first show the time-and 𝑥-averaged horizontal flow at 𝑧 = 0.9 normalized by 𝑅𝑒 as a function of Λ in figure 9. For Λ 10 -2 , -𝑢(𝑧 = 0.9) ≈ 𝑅𝑒 > 0, which indicates that the large-scale leftward current dominates the dynamics. For small Λ 10 -2 , including Λ = 0, -𝑢(𝑧 = 0.9) 𝑅𝑒 (either positive or negative) but is not always zero, because 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 Λ = 10 -2 due to the superposition of the RB and HC dynamics, which cannot be inferred from figure 9 but will be discussed in section § 3.5.

The calculation of the characteristic length scale ℓ of overturning motions from 𝑢(𝑧 = 0.9) is illustrated in figure 10 for a RB-like case (Λ = 10 -3 ; top row) and a HC-like case (Λ = 1; bottom row) with 𝑅𝑎 𝐹 = 10 8 and Γ = 8. Figure 10(a) shows the (𝑥, 𝑡)-Hovmöller diagram of 𝑢(𝑧 = 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 Competition between Rayleigh-Bénard and horizontal convection 15 time. The auto-correlation function of 𝑢 in 𝑥 is defined as

R 𝑢𝑢 (𝑥) ≡ ∫ Γ/2 -Γ/2 𝑢(𝑥 )𝑢(𝑥 + 𝑥) d𝑥 ∫ Γ/2 -Γ/2 𝑢(𝑥 )𝑢(𝑥 ) d𝑥 , (3.8)
where 𝑢(𝑥 + 𝑥 > Γ/2) is set to 0 and 𝑥 > 0 is the lag. The auto-correlation function evaluated at 𝑧 = 0.9 shows an oscillatory pattern (figure 10(b)), like 𝑢(𝑧 = 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 ℓ 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 Λ = 10 -3 , figure 10(c) shows that ℓ ≈ 1 (dashed solid line), which is approximately the value expected for an unconfined RB roll (rotating clockwise or anti-clockwise). Figures 10(d)-(f) show the same results as figure 10(a)-(c) but for Λ = 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 𝑥. Thus, the characteristic length scale equals the domain size, i.e. ℓ = Γ, 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 𝑢(𝑥 ) with 𝑢(𝑥 + 𝑥) in the denominator). However, such a definition creates large boundary effects (not shown), which make the calculation of ℓ 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.

The characteristic length scale ℓ obtained from the auto-correlation function of 𝑢(𝑧 = 0.9) is shown for all simulations as a function of Λ in figure 11. For Λ 10 -2 , ℓ ≈ 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 Λ values, bursts of large-scale currents sweeping away RB cells (further explained in section § 3.5). For Λ 10 -2 , ℓ = Γ, 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 horizontal convection.

Bi-stability and bursts near the transition

This section focuses on the transitional regime between the RB and HC dynamical regimes obtained for Λ 10 -2 and Λ 10 -2 , respectively. To this end, we present numerical results that use the method of discrete continuation. We fix 𝑅𝑎 𝐹 = 10 6 and Γ = 12, and we gradually vary the flux ratio Λ in the range [0, 0.03]. For each value of Λ, 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

𝑢 𝑧>0.5 = 2 Γ ∫ Γ/2 -Γ/2 ∫ 1 0.5
𝑢 d𝑥 d𝑧, (3.9)

and mean temperature difference between the right-and left-hand sides of the domain

Δ𝑇 𝐻 = 𝑇 𝑥>Γ/4 -𝑇 𝑥<-Γ/4 = 4 Γ ∫ Γ/2 Γ/4 ∫ 1 0 𝑇 d𝑥 d𝑧 - 4 Γ ∫ -Γ/4 -Γ/2 ∫ 1 0 𝑇 d𝑥 d𝑧, (3.10)
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 equation (3.9) performed on the lower half of the domain would yield exactly the opposite mean value. We start by gradually increasing Λ 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 Λ ≈ 0.016 (figure 12(a)). This transition corresponds to the first destabilisation of the Rayleigh-Bénard cells with approximately unit aspect ratio, which merge into horizontally elongated cells. For our particular choice of Γ = 12, the flow is initially organised in twelve cells and switches to six more elongated cells. A secondary transition between these six cells and only one extended cell is observed at Λ = 0.02. The system remains in this state for larger values of Λ. The two transitions at Λ = 0.016 and Λ = 0.02 are also observed in the mean horizontal temperature difference (figure 12(b)). The mean horizontal temperature difference Δ𝑇 𝐻 initially follows the temperature difference imposed along the top boundary, which is equal to ΛΓ (dashed line in figure 12(b)). As Λ increases, Δ𝑇 𝐻 slowly deviates downward until it recovers the imposed value ΛΓ just after the transition at Λ = 0.016. After the second transition, Δ𝑇 𝐻 clearly follows a different trend and increases less rapidly with Λ, which is indicative of the growing efficiency of HC in its ability to mix the imposed horizontal temperature difference.

To explore the possibility of multi-stability, we gradually decrease Λ from a given equilibrium state. We do so independently for each observed transition, i.e. starting first at Λ = 0.016 and then at Λ = 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 Λ, which is highlighted by the green shaded region in figure 12(a),(b). The bi-stability is characterized by the separation of the ascending and descending branches, which are each connected to distinct flow organizations: for 0.01 < Λ < 0.02, there always exist at least two flow states for the same Λ that have a different number of convective rolls; each roll having an aspect ratio that can vary between 1 and Γ. Thus, the actual number of convective cells in both horizontal and vertical directions crucially depends on the history of the system. The possibility to have different number of rolls for the same Λ 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 Λ = 0.014 (figure 12(c)) and Λ = 0.018 (figure 12(d)). For completeness, we also show in figure 12(e) an example of the final state of pure horizontal convection reached once Λ > 0.02. While 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, since Δ𝑇 𝐻 can either increase (orange line above blue line in figure 12(b)) or decrease (blue line above green line in figure 12(b)) with decreasing roll number. The former behavior 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, six rolls are more efficient at maintaining a large Δ𝑇 𝐻 than twelve rolls.

We would like to make a few notes regarding the results of figures 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 downwelling versus weak upwelling). 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 Λ < 0.01, we still observe bi-stability since the system does not recover exactly its initial state (characterizing the increasing branch) as we continue decreasing Λ. 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 Λ ≈ 0.016: the ascending trajectory bifurcates exactly where the trajectory descending from Λ = 0.02 falls back onto the descending trajectory initialised from Λ = 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 what can be deciphered by our study, especially as Γ 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 Γ 2 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 et al. 2022), which is beyond the scope of the present study.

We expect the hysteretic behaviour to persist for aspect ratios larger than Γ = 𝑂 (10), since this allows for even more distinct convective cell configurations (see Wang et al. (2020) 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 𝑅𝑎 𝐹 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, while 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 𝑅𝑎 𝐹 = 10 8 and Γ = 8. For Λ = 0.01, i.e. right at the lower bound of the bi-stable regime observed at 𝑅𝑎 𝐹 = 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 equation (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 Rayleigh-Bénard cells, and are reminiscent of the burst dynamics observed in RB convection between stress-free plates (Goluskin et al. 2014). The two bottom 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, while figure 13(c) shows a state featuring a burst of horizontal convection (note that a similar HC burst can be seen in the movie available as Supplementary Material). Figure 13(d) shows the time history of 𝑢(𝑧 = 0.9), which can help visualize 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 Λ 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.

We conclude our analysis of RBH dynamics by discussing and comparing the Probability Density Function (PDF) of the horizontal mean flow defined by equation (3.9) in different regimes (including bi-modality) near the transition. We consider fixed (relatively large) 𝑅𝑎 𝐹 = 10 8 , Γ = 8 and variable Λ close to the transition, as well as pure RB convection (Λ = 0) and pure HC (no geothermal flux) for comparison. We first show in figure 14(a skewed distribution with an elongated negative tail for Λ = 10 -2 . As Λ further increases, the PDFs become more and more 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 RB-HC regime, i.e. which is absent in pure RB or HC regimes. Figure 14(b) shows the PDF of the mean horizontal upper flow for three cases at 𝑅𝑎 𝐹 = 10 8 and Γ = 8: one with Λ = 0 (pure RB case), one with Λ = 0.05 but no geothermal flux (pure HC) and finally the RBH case with Λ = 0.05 and geothermal heating. We clearly see that both RB and HC cases produce mean-flow PDFs that are approximately Gaussian and centered around zero (dotted line) and a small negative value (dashed line), respectively. Conversely, the RBH case yields a PDF 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. 2006) 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 Λ, 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 Λ 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 𝑅𝑎 𝐹 and large Γ limit, to future studies.

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 Λ 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 Λ 3 × 10 -3 (purple shading in figures 7, 8, 9 and 11) to HC at large Λ 3 × 10 -2 (grey shading). Importantly, the transition occurs near Λ = 10 -2 independently of the aspect ratio Γ and Rayleigh number 𝑅𝑎 𝐹 , 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 equations (3.3) and (3.7)) are good indicators of the flow regime; they deviate quickly from the classical scaling laws of Competition between Rayleigh-Bénard and horizontal convection 21 RB convection and HC in the regime where they are not relevant, i.e. Λ 10 -2 for 𝑁𝑢 𝑅𝐵 and Λ 10 -2 for 𝑁𝑢 𝜒 𝐻 𝐶 . The characteristic length scale ℓ derived from the auto-correlation function ( § 3.4) is another excellent indicator as it is equal to 1 in the RB regime and Γ in the HC regime. The Reynolds number is not a good indicator for our set of simulations, however, because it is primarily controlled by 𝑅𝑎 𝐹 , the intensity of buoyancy driving, with small variations only due to Λ (flow topology) and Γ.

The fact that the transition does not depend on Γ (at least for Γ 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. Since 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 𝑅𝑎 𝐹 = 10 8 but are qualitatively representative of simulations using other 𝑅𝑎 𝐹 ), 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 1970;Krishnamurti & Howard 1981;Busse 1983). For instance, arrays of convective rolls are known to be unstable to the so-called shearing instability (Hughes & Proctor 1990;Rucklidge & Matthews 1996;Goluskin et al. 2014) 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 Λ = 0, because of our closed domain and no-slip boundary conditions, which are known to inhibit the instability for elongated domains (Fitzgerald & Farrell 2014;Van Der Poel et al. 2014). However, for Λ > 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 Λ ≈ 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 𝑅𝑎 𝐹 = 10 6 , we have found at least two stable branches for 0.01 Λ 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 organization of the rolls can differ substantially between cases. For larger 𝑅𝑎 𝐹 = 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 𝑅𝑎 𝐹 = 10 8 and Λ = 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 Λ = 10 -2 with 𝑅𝑎 𝐹 is rigorously valid only for 10 6 𝑅𝑎 𝐹 10 9 as considered in this study. Exploring the RBH dynamics at larger 𝑅𝑎 𝐹 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 𝑃𝑟 = 𝑂 (10), hence is much larger than 𝑃𝑟 = 1, but its effect on the transition may be limited. Indeed, increasing 𝑃𝑟 for 𝑃𝑟 1 has little effect on 𝑁𝑢 𝑅𝐵 and 𝑁𝑢 𝜒 𝐻 𝐶 , whereas 𝑅𝑒 decreases almost like 𝑃𝑟 -1 in both RB convection and HC, at least for moderate Rayleigh numbers (Shishkina & Wagner 2016;Li et al. 2021). Thus, we hypothesize that increasing 𝑃𝑟 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 𝑃𝑟 (Goluskin et al. 2014), suggesting that the upper bound of the RB regime (rightmost edge of purple shading in, i.e. figure 11) may be in fact sensitive to 𝑃𝑟.

The dual Rayleigh-Bénard-Horizontal 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. (2015)) and because the contribution of geothermal heating to the ocean abyssal stratification and circulation is now known to be significant (Mashayek et al. 2013;De Lavergne et al. 2016). For subglacial lakes, our results suggest that RB convection should dominate since Λ 10 -2 for realistic ice-water interface slopes (Appendix B). It is unclear whether Rayleigh-Bénard convection will transition toward HC at lower or higher Λ in three dimensions than in two dimensions. Preliminary threedimensional simulations suggest a transition around Λ ≈ 10 -2 again, which means that our results may be (at least qualitatively) informative for real systems in spite of the twodimensional limitation. Non-rectangular geometry and Earth's rotation, whose impact on subglacial lake dynamics remains unclear (Couston & Siegert 2021), are other effects not considered in this work that may be more favorable 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 2021;Olsthoorn et al. 2021), 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 flat waterbedrock 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 characterize past climates by analyzing sediment cores extracted from the lake bed (Hodgson et al. 2009). In lakes dominated by HC, vertical profiles should be ideally performed at both ends of the lake in order to probe both the downwelling and upwelling 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 downwelling region. Dating sediments from continuous cores (resulting from successive layer depositions) is typically easier than from non-continuous cores (Hodgson et al. 2009). Thus, initial dating could be based on cores extracted from below the downwelling region, henceforth enabling dating sediment cores extracted from below the upwelling 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 homogenizes 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 upwelling and downwelling regions.

to thank Adrien Villaret for his careful review and Andy Smith and Keith Makinson for useful discussions about subglacial lake exploration.

The 𝑅𝑒(𝑅𝑎 𝐹 ) scaling can be recast as a 𝑅𝑒(𝑅𝑎 Δ ) scaling, where 𝑅𝑎 Δ is the classical Rayleigh-Bénard parameter based on the temperature difference Δ between the top and bottom plates (averaged in 𝑥 and 𝑡 and computed a posteriori). To see this, we first remark that the different definitions of the problem parameters yield the relationship (in the limit Λ → 0) 𝑁𝑢 𝑅𝐵 = 𝑅𝑎 𝐹 /𝑅𝑎 Δ (Johnston & Doering 2009), such that 𝑅𝑎 𝐹 = 𝑎 1/(1-𝑏 𝑅𝐵 ) 𝑅𝐵 𝑅𝑎 𝑏 𝑅𝐵 /(1-𝑏 𝑅𝐵 ) Δ . This yields the scaling law 𝑁𝑢 = 0.316𝑅𝑎 0.236 Δ using the pre-factor and exponent reported in table 1. The obtained 𝑁𝑢(𝑅𝑎 Δ ) scaling law is in relatively good agreement with Grossmann & Lohse (2000)'s unifying theory: the exponent is slightly larger (and the pre-factor is slightly smaller) than that predicted by regime 𝐼 𝐼 𝑢 , i.e. 1/5, which is expected for our low 1.8 × 10 5 𝑅𝑎 Δ 4.8 × 10 7 numbers (see updated regime diagram in Stevens et al. 2013, and note that our scaling law happens to be in better agreement with the 1/4 exponent of regime 𝐼 𝑢 expected at larger 𝑃𝑟 > 1). We expect the difference between our scaling exponent and that of regime 𝐼 𝐼 𝑢 to be due to the two-dimensional flow assumption, which is known to affect 𝑁𝑢 𝑅𝐵 , especially at low-to-moderate 𝑃𝑟 (Van Der Poel et al. 2013) and the relative proximity of regime 𝐼𝑉 𝑢 at slightly larger 𝑅𝑎 Δ , which features a 1/3 exponent. We note that our exponent is smaller (and the pre-factor is larger )) as functions of 𝑅𝑎 𝐿 . Pure HC simulation results are indistinguishable from the scaling law (solid line) derived from the 𝑁𝑢 𝜒 𝐻 𝐶 data, which means that the heat transfer efficiency by horizontal convection 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 figures 15(b) and 15(c)). However, for Λ 10 -2 , these differences are due to the diffusive normalization (i.e. appearing at the denominator), whose dependence on geothermal heating remains significant for Λ 𝑂 (1) (unlike the numerator) and depends on the choice of Nusselt number definition. When renormalizing the numerator in the definition of each Nusselt number by the diffusive solution without geothermal heating (i.e. considering 𝐹 = 0), those differences (at large Λ) are removed, as can be seen from figures 15(d)-(f) in which all renormalized Nusselt numbers divided by the scaling law 𝑎 𝐻 𝐶 𝑅𝑎 𝑏 𝐻𝐶 𝐿 converge toward unity. Note that differences due to aspect ratio Γ are also due to the diffusive normalization, as explained at the end of § 3.3.

Figure 1 :

 1 Figure 1: (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), while 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 horizontal convection because of the pressure-dependence of the freezing temperature. (b) Phase diagram of this study. The circles show the location of the simulations in (Λ, 𝑅𝑎 𝐹 ) space for Γ = 8; for Γ = 4 (resp. Γ = 16) the markers are diamonds (squares) and are slightly shifted downward (upward) with respect to 𝑅𝑎 𝐹 = 10 8 (𝑅𝑎 𝐹 = 10 7 ) for better visibility. The stars show pure HC simulations run with Λ = 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 Rayleigh-Bénard convection and horizontal convection, respectively.

Figure 2 :

 2 Figure 2: Snapshots of the velocity magnitude√𝑢 2 + 𝑤 2 and velocity vector field (shown by arrows) at the end of the simulations with 𝑅𝑎 𝐹 = 10 8 and Λ increasing from top to bottom. The bottom figure shows the results obtained for pure HC, i.e. with Λ = 1 but with an insulating bottom boundary (no geothermal flux). The characteristic length scale of overturning motions is estimated from the horizontal flow at 𝑧 = 0.9 (see § 3.4), which is highlighted by the blue dashed lines.

Figure 3 :

 3 Figure 3: Snapshots of the temperature field at the end of the simulations with 𝑅𝑎 𝐹 = 10 8 and Λ increasing from top to bottom. The solid (resp. dashed) lines show positive (resp. negative) contours of the streamfunction, which is set to 0 at the bottom left corner. The bottom figure shows the results obtained without geothermal heating. Note that the color bars are not the same between plots.

Figure 4 :

 4 Figure 4: Time history of volume-averaged (top row) temperature and (bottom row) Reynolds number. Simulation results are split into two different subplots (left and right columns) for better visibility (as indicated by the plot titles). Different colors correspond to different Rayleigh numbers 𝑅𝑎 𝐹 , lighter colors indicate higher Λ, and line width increases with Γ. Time-averaged variables representative of the statistical steady state are integrated from 𝑡 = 1 (shown by the vertical dashed lines) onward.

  Figure 5: (a) Mean temperature at statistical steady state for all simulations as a function of the flux ratio Λ. The black symbols show the coldest temperature on the top boundary, i.e. 𝑇 = -ΛΓ/2. (b),(c) Reynolds number 𝑅𝑒 as a function of 𝑅𝑎 𝐹 and 𝑅𝑎 𝐿 , respectively, with scalings 𝑅𝑒 = 𝑐 𝑅𝐵 𝑅𝑎 𝑑 𝑅𝐵 𝐹 and 𝑅𝑒 = 𝑐 𝐻 𝐶 𝑅𝑎 𝑑 𝐻𝐶 𝐿 shown by the black solid lines. (d) Compensated Reynolds number as a function of Λ. The star symbols show results obtained for pure HC simulations.

  Figure 6: (a) Time-averaged heat flux along the top boundary for 𝑅𝑎 𝐹 = 10 8 , Γ = 8 and four different Λ values (see 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 Λ = 10 -2 .

  Figure 7: (a) RB-based Nusselt number 𝑁𝑢 𝑅𝐵 as a function of 𝑅𝑎 𝐹 with scaling law 𝑁𝑢 𝑅𝐵 = 𝑎 𝑅𝐵 𝑅𝑎 𝑏 𝑅𝐵 𝐹 shown by the black solid line. (b) Compensated RB-based Nusselt number as a function of Λ.

  Figure 8(a) shows 𝑁𝑢 𝜒 𝐻 𝐶 as a function of 𝑅𝑎 𝐿 . Simulation results obtained for Λ 10 -2 all tend asymptotically (for fixed Λ, or color intensity) toward power laws parallel (in loglog space) to the one shown by the black solid line (see e.g. the dashed line), which was obtained from best fit for Λ = 1, Γ = 8 and without geothermal flux (cf. perfect overlap with star symbols), and whose exponent is compatible with the 1/5 exponent predicted by Rossby (1965).

  Figure 8(a) shows that 𝑁𝑢 𝜒 𝐻 𝐶 can help distinguish simulations dominated by RB convection from simulations dominated by HC: at small Λ (dark symbols with typically small 𝑅𝑎 𝐿 ), 𝑁𝑢 𝜒 𝐻 𝐶 < 1, while at large Λ (light symbols with typically large 𝑅𝑎 𝐿 ), 𝑁𝑢 𝜒 𝐻 𝐶 1. The large spread of 𝑁𝑢 𝜒 𝐻 𝐶 with Λ and Γ for fixed 𝑅𝑎 𝐿 in the HC limit (i.e. at large Λ) is somewhat unexpected but can be explained: it is due to the diffusive normalization 𝜒 diff . First, figure 8(b) shows that dividing 𝑁𝑢 𝜒 𝐻 𝐶 by the pure HC scaling 𝑎 𝐻 𝐶 𝑅𝑎 𝑏 𝐻𝐶 𝐿 results in a steep scaling with Λ (+2 slope shown by the solid line) for Λ

Figure 8 :

 8 Figure 8: (a) HC-based Nusselt number 𝑁𝑢

Figure 9 :

 9 Figure 9: Temporally-and horizontally-averaged horizontal flow near the top boundary at 𝑧 = 0.9 normalized by the Reynolds number as a function of Λ.-𝑢(𝑧 = 0.9) 𝑥 is of the same order as 𝑅𝑒 for relatively large Λ.

  Figure 10: (a) Horizontal flow near the top boundary (𝑧 = 0.9) as a function of (𝑥, 𝑡) for 𝑅𝑎 𝐹 = 10 8 , Γ = 8 and Λ = 10 -3 . (b) Auto-correlation function R 𝑢𝑢 in 𝑥 of the horizontal flow shown in (a) as a function of time 𝑡. (c) Time average auto-correlation R 𝑢𝑢 as a function of lag in 𝑥. 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 Λ = 1. For a monotonically decreasing auto-correlation function, the characteristic length scale equals the domain length Γ.

Figure 11 :

 11 Figure 11: Characteristic length scale ℓ at statistical steady state as a function of Λ for all simulations. The dashed lines show the aspect ratios (or domain lengths), i.e. Γ = 4, 8 and 16. At large Λ, ℓ = Γ, which is an indication of horizontal convection dominating the dynamics, while at small Λ, ℓ ≈ 1 (dotted line) indicates multiple overturning cells, as in classical RB convection.

  L.-A.Couston, J. Nandaha, B. Favier (a) (

Figure 12 :

 12 Figure 12: Bifurcation diagrams based on averaged velocity and temperature. (a) Mean horizontal flow averaged over time and the upper half of the domain (equation (3.9)) as a function of the heat flux ratio Λ. The parameters are 𝑅𝑎 𝐹 = 10 6 and Γ = 12. Round symbols correspond to gradually increasing Λ while cross symbols correspond to gradually decreasing Λ. 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-(𝑥 > Γ/4) and left-hand (𝑥 < -Γ/4) sides of the domain (equation (3.10)). (c) Snapshots in (𝑥, 𝑧) of the temperature field with streamlines shown as black contours for Λ = 0.014 when Λ is going up (top panel) and down (bottom panel). (d) Same as (c) but for Λ = 0.018. (e) Same as (c) but for Λ = 0.028 (branch with Λ going up only). (f) Horizontal flow at 𝑧 = 0.9 as a function of 𝑥 and 𝑡 (initial time set to 0 arbitrarily) for the simulation right after the first transition of the increasing-Λ (blue) branch on panels 12(a),(b) (Λ = 0.016; highlighted by triangles). Note that merging events, e.g. at 𝑥 ≈ 2 and 𝑡 ≈ 7, appear discontinuous because they occur on short time scales.

  Figure 13: (a) Mean horizontal flow averaged over the upper half of the domain as a function of time for 𝑅𝑎 𝐹 = 10 8 , Γ = 8 and Λ = 10 -2 . (b), (c) Snapshots of the temperature field at times 𝑡 = 3.1 (vertical dashed line in the top panel) and 𝑡 = 4 (vertical dotted line), respectively. (d) Spatio-temporal plot of 𝑢(𝑧 = 0.9) for the same simulation.

Figure 14 :

 14 Figure 14: (a) Probability density function (PDF) at statistical steady state of 𝑢 𝑧>0.5 for 𝑅𝑎 𝐹 = 10 8 , Γ = 8 with Λ = 0.007, 0.01, 0.02, 0.05. (b) Same as (a) but for pure RB convection (Λ = 0; dotted line), pure HC (Λ = 0.05 but no geothermal flux; dashed line), and RBH convection (Λ = 0.05; solid line).

  Figures 15(a)-(c) show 𝑁𝑢 abs 𝐻 𝐶 , 𝑁𝑢 half 𝐻 𝐶 and 𝑁𝑢 𝜒 𝐻 𝐶 (defined in equations (3.4)-(3.7)) as functions of 𝑅𝑎 𝐿 . Pure HC simulation results are indistinguishable from the scaling law (solid line) derived from the 𝑁𝑢

Figure 15 :

 15 Figure 15: The top row shows the three Nusselt numbers defined in equations (3.4)-(3.7) as functions of 𝑅𝑎 𝐿 with the solid black line showing the scaling law derived from the 𝑁𝑢 𝜒 𝐻 𝐶 data obtained for pure HC simulations. The bottom row shows the same Nusselt numbers divided by the scaling law and renormalized by the diffusive solution derived without geothermal heating (𝐹 = 0); e.g. 𝑁𝑢 𝜒 𝐻 𝐶 (𝐹 = 0) = 𝑁𝑢 𝜒 𝐻 𝐶 × 𝜒 diff /𝜒 diff (𝐹 = 0).

Table 1 :

 1 Pre-factor and exponent of all power laws referred to in the paper and shown as black solid lines in figures 5(b),(c), 7(a) and 8(a). The power laws for 𝑅𝑒(𝑅𝑎 𝐹 ) and 𝑁𝑢(𝑅𝑎 𝐹 ) are based on four simulations with Λ = 0 and Γ = 8, whereas the power laws for 𝑅𝑒(𝑅𝑎 𝐿 ) and 𝑁𝑢(𝑅𝑎 𝐿 ) are based on three simulations without geothermal flux and Γ = 8 (see table 2 in Appendix A).

		power law	pre-factor	exponent
	𝑅𝑒 = 𝑐 𝑅𝐵 𝑅𝑎 𝑑 𝑅𝐵 𝐹	2.326 × 10 -1	0.441
	𝑅𝑒 = 𝑐 𝐻 𝐶 𝑅𝑎 𝑑 𝐻𝐶 𝐿	1.839 × 10 -2	0.396
	𝑁𝑢 𝑅𝐵 = 𝑎 𝑅𝐵 𝑅𝑎 𝑏 𝑅𝐵 𝐹	3.938 × 10 -1	0.191
	𝑁𝑢	𝜒 𝐻 𝐶 = 𝑎 𝐻 𝐶 𝑅𝑎 𝑏 𝐻𝐶 𝐿	3.181 × 10 -1	0.226

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

Funding. This research received funding from LabEx LIO (ANR-10-LABX-0066), Université de Lyon.

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 sections § 3.1-3.4 are provided in table 2. The details of the simulations used for the multi-stability analysis are provided in § 3.5 and summarized on the last line of table 2.

Appendix B. Parameter range relevant to subglacial lakes

In this Appendix we motivate the choice of small Λ considered in the paper. We are particularly interested in subglacial lakes, which lie beneath several kilometers 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. 2022). 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 icewater interface because it depends on local pressure, which is larger where the ice is thicker (Thoma et al. 2010). A linear approximation for the freezing temperature of freshwater with no dissolved air as a function of local ice pressure 𝑝 𝑖 (in decibar) or ice thickness ℎ 𝑖 = 10 4 𝑝 𝑖 /(𝜌 𝑖 𝑔) (in meter), with 𝜌 𝑖 = 917 kg/m 3 the ice density and 𝑔 = 9.81 m/s 2 the surface gravity, is 𝑇 𝑓 = 6.67 × 10 -2 -8.21 × 10 -4 𝑝 𝑖 = 6.67 × 10 -2 -9.12 × 10 -4 ℎ 𝑖 , (B 1)

as obtained from best-fit with the exact freezing temperature calculated with the python package gsw (McDougall & Barker 2011) for 0 < 𝑝 𝑖 < 5000 dbar (ℎ 𝑖 ranging from 0 to 5558 m). Equation (B 1) shows that 𝑇 𝑓 decreases with 𝑝 𝑖 or ℎ 𝑖 (𝑇 𝑓 ≈ -4 • C when 𝑝 𝑖 = 5000 dbar) and yields a horizontal temperature gradient equal to

for a subglacial lake with an ice-water interface slope 𝑑ℎ 𝑖 /𝑑𝑥 = -𝛾, i.e. thicker for smaller 𝑥, as assumed in this study. We find that 𝛾 = 0.016, 0.003, 0.03, 0.002, 0.003, for the five well-documented subglacial lakes reported in Couston & Siegert (2021), i.e. CECs, SPL, Ellsworth, Vostok and Concordia, respectively. This range of 𝛾 values yields approximately 3 × 10 -5 < Λ = 𝑘𝜆/𝐹 < 3 × 10 -4 , which lies within the range 0 Λ 1 explored in this study, with 𝑘 = 0.56 W/(m • C) the thermal conductivity of water and 𝐹 = 50 mW/m 2 an average value for Earth's geothermal flux. Thus, for subglacial lakes, Λ 10 -2 , unless the slope is of order unity (in which case Λ 10 -2 ), which is most likely unrealistic. For completeness, the flux-based Rayleigh number of a subglacial lake below thick ice with depth 𝐻 = 100 m, 𝐹 = 50 mW/m 2 , 𝜈 = 1.7 × 10 -6 m 2 /s, 𝜅 = 1.33 × 10 -7 m 2 /s, 𝛼 = 10 -4 • C -1 (typical values inferred from Couston (2021)) is 𝑅𝑎 𝐹 ≈ 3.9 × 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).

yes 10 6 0 ↔ 0.03 0 ↔ 6 × 10 8 10 7 7 DNS Table 2: Physical and numerical parameters of the simulations, with a dash denoting same value as in the row above for readability. Γ is the aspect ratio, 𝑅𝑎 𝐹 is the flux-based Rayleigh number, Λ is the heat flux ratio, 𝑅𝑎 𝐿 is the horizontal Rayleigh number, 𝑛 𝑧 is the number of elements in the vertical direction, 𝑙 𝑑 is the polynomial order and 𝑑𝑡 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 (2001) from 𝑡 = 0 to 𝑡 = 1 and then restarted with numerical parameters as indicated in the table and integrated from 𝑡 = 1 to 𝑡 = 1.2 via DNS. The last line shows the parameters of the simulations discussed in section 3.5, which focus on bi-stability.

Appendix D. Dissipation of temperature variance and diffusive solution

Recently, Rocha et al. (2020b) showed that a compelling definition of the Nusselt number of horizontal convection should involve the dissipation of temperature (or buoyancy 𝑏 = -𝑇) variance, 𝜒 = |∇𝑇 | 2 , because it is related to the vertically-averaged horizontal heat flux set up by the imposed heterogeneous temperature profile at (here) 𝑧 = 1. They further showed that 𝜒 = 𝑇 (𝑧 = 1)𝜕 𝑧 𝑇 (𝑧 = 1) 𝑥 , 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 |∇𝑇 | 2 . Here we derive a similar result for the case with geothermal heating. Multiplying (2.1c) by 𝑇 and rearranging yields

The first term becomes zero when time averaging at statistical steady state, while the third term becomes zero when performing a volume average because of the no-slip condition on L.-A. Couston, J. Nandaha, B. Favier the walls. Thus,

where the last equality is obtained after integrating in 𝑧 and enforcing no-flux conditions on the side walls. Our expression for 𝜒 differs from that of Rocha et al. (2020b) because of the second term on the right-hand-side of equation (D 2), which is non-zero here because our derivation takes into account the heat flux (𝜕 𝑧 𝑇 = -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 equation (2.5), the diffusive solution without geothermal heating (𝜕 𝑧 𝑇 = 0 at 𝑧 = 0) is

while, with geothermal heating (𝜕 𝑧 𝑇 = -1 at 𝑧 = 0), we simply add a linearly-varying vertical profile, i.e.

The denominator in equation (3.7) based on equation (D 4) then reads