L.-A Couston
email: couston@irphe.univ-mrs.fr.

D Lecoanet

B Favier

M Le Bars

DYNAMICS OF MIXED CONVECTIVE-STABLY-STRATIFIED FLUIDS

We study the dynamical regimes of a density-stratified fluid confined between isothermal noslip top and bottom boundaries (at temperatures Tt and T b ) via direct numerical simulation. The thermal expansion coefficient of the fluid is temperature dependent and chosen such that the fluid density is maximum at the inversion temperature T b > T i > Tt. Thus, the lower layer of the fluid is convectively unstable while the upper layer is stably stratified. We show that the characteristics of the convection change significantly depending on the degree of stratification of the stable layer. For strong stable stratification, the convection zone coincides with the fraction of the fluid that is convectively unstable (i.e. where T > T i ), and convective motions consist of rising and sinking plumes of large density anomaly, as is the case in canonical Rayleigh-Bénard convection; internal gravity waves are generated by turbulent fluctuations in the convective layer and propagate in the upper layer. For weak stable stratification, we demonstrate that a large fraction of the stable fluid (i.e. with temperature T < T i ) is instead destabilized and entrained by buoyant plumes emitted from the bottom boundary. The convection thus mixes cold patches of low density-anomaly fluid with hot upward plumes, and the end result is that the T i isotherm sinks within the bottom boundary layer and that the convection is entrainment-dominated. We provide a phenomenological description of the transition between the regimes of plume-dominated and entrainment-dominated convection through analysis of the differences in the heat transfer mechanisms, kinetic energy density spectra, and probability density functions for different stratification strengths. Importantly, we find that the effect of the stable layer on the convection decreases only weakly with increasing stratification strength, meaning that the dynamics of the stable layer and convection should be studied self-consistently in a wide range of applications.

Introduction

Buoyancy-driven convection can result from a number of physical mechanisms leading to a change of the fluid density, including evaporation, heating, and sedimentation. Convection occurs naturally in the interior of stars and planets, as well as in oceans and atmospheres, but is also sometimes used to improve heat transfer properties in buildings and industrial processes. The transport of mass, momentum and heat are all important aspects of convective flows.

Because density variations are often small compared to the mean fluid density, buoyancy-driven convection studies typically neglect density effects other than in the buoyancy force, an approximation known as the Boussinesq assumption. The Oberbeck-Boussinesq assumption (OB), which additionally assumes constant physical properties (such as viscosity, conductivity, and thermal expansion), is also often considered and has resulted in much of our fundamental understanding of convective flows, spear-headed by the canonical problem of (thermally-driven) Rayleigh-Bénard convection (RBC) [START_REF] Ahlers | Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection[END_REF].

While the Boussinesq assumption is applicable in a number of physically-interesting fluid systems, some dynamics of realistic convective fluids are not attained under the Oberbeck-Boussinesq approximation. In thermal convection, for instance, symmetric top and bottom boundary conditions necessarily result in a top/down mid-plane symmetry of the flow under OB assumption. Helium convection experiments have demonstrated asymmetric fluid motions with respect to the mid-plane due to temperature-dependent viscosity [START_REF] Urban | Effect of boundary layers asymmetry on heat transfer efficiency in turbulent Rayleigh-Benard convection at very high Rayleigh numbers[END_REF]; an effect discarded within the OB framework. Non-Oberbeck-Boussinesq (NOB) effects are of current research interest because a number of fluids (e.g. water, helium, ethane) have temperature-dependent physical properties, in particular close to their critical point (such as in glycerol, see e.g. [START_REF] Horn | On non-Oberbeck Boussinesq effects in three-dimensional Rayleigh Bénard convection in glycerol[END_REF]. Of all NOB effects, those due to temperature-dependent expansion coefficients may be most striking (in ethane, see e.g. [START_REF] Ahlers | Non-Oberbeck-Boussinesq effects in turbulent thermal convection in ethane close to the critical point[END_REF], and at the same time can be readily observed in nature. This is especially true because the expansion coefficient of water varies significantly close to its density maximum at temperature T 0 = 4 • C, a temperature ubiquitous in oceans and lakes which can thus routinely experience NOB effects [START_REF] Nycander | The nonlinear equation of state of sea water and the global water mass distribution[END_REF][START_REF] Vettoretti | Thermohaline instability and the formation of glacial North Atlantic super polynyas at the onset of Dansgaard-Oeschger warming events[END_REF]. The density maximum for water results from the quadratic equation of state (EoS) ρ ∝ ρ 0 [1 -β(T -T 0 ) 2 ], which is a simplified yet realistic EoS for water at atmospheric pressure [START_REF] Roquet | Defining a Simplified Yet Realistic Equation of State for Seawater[END_REF]. Because of the quadratic temperature nonlinearity, the water with a positive vertical temperature gradient reverses from buoyancy decreasing with height for T ≤ T 0 = 4 • C to buoyancy increasing with height for T ≥ T 0 = 4 • C. This results in a layer subject to convective instability beneath a stably stratified layer. The existence of an adjacent stable stratification can have a significant influence on the convective dynamic, most notably allowing for penetrative convection, as demonstrated theoretically via stability analysis decades ago (e.g. [START_REF] Veronis | Penetrative convection[END_REF]. More recently, laboratory experiments (Le [START_REF] Bars | Experimental study of internal wave generation by convection in water[END_REF] and numerical simulations [START_REF] Lecoanet | Numerical simulations of internal wave generation by convection in water[END_REF] have improved our understanding of such penetrative convection following the seminal study by [START_REF] Townsend | Natural convection in water over an ice surface[END_REF] and have shown a clear picture of plumes interacting with a stable layer and generating internal waves, an effect inaccessible within the OB framework.

NOB convection with a buoyancy reversal is not limited to water but also occurs in a number of alcoholglycol-water mixtures used as analogues of gas dissolution in porous media [START_REF] Backhaus | Convective instability and mass transport of diffusion layers in a Hele-Shaw geometry[END_REF][START_REF] Hidalgo | Scaling of convective mixing in porous media[END_REF][START_REF] Hewitt | Ultimate regime of high Rayleigh number convection in a porous medium[END_REF], glass-forming liquids such as BeF 2 [START_REF] Hemmati | Interpretation of the molten BeF2 viscosity anomaly in terms of a high temperature density maximum, and other waterlike features[END_REF], and in the troposphere where an inversion layer forms atop the cloud-top mixing layer as a result of evaporation [START_REF] Mellado | The evaporatively driven cloud-top mixing layer[END_REF]. The heterogeneous convective/wave dynamic is also relevant to many astrophysical and geophysical settings: In stars such as our Sun, convection only occurs in the relatively cool outer shell, and couples with an underlying stable region where gravity waves propagate [START_REF] Brummell | Penetration and Overshooting in Turbulent Compressible Convection[END_REF][START_REF] Rogers | Penetrative Convection within the Anelastic Approximation[END_REF][START_REF] Alvan | Coupling between internal waves and shear-induced turbulence in stellar radiation zones: the critical layers[END_REF][START_REF] Pinçon | Astrophysics Generation of internal gravity waves by penetrative convection[END_REF]; In the Earth, observations within the convective liquid outer core point to a possible stratified layer near both the core-mantle boundary and the inner-core boundary [START_REF] Hirose | Composition and State of the Core[END_REF][START_REF] Buffett | Geomagnetic fluctuations reveal stable stratification at the top of the Earth's core[END_REF]. Buoyancy reversal in astrophysical and geophysical systems is due to changes from super-to sub-adiabatic mean temperature profiles, related to changes of physical properties with temperature and pressure.

In addition to the effects of the Rayleigh and Prandtl numbers on thermal convection, buoyancy reversal dynamics are strongly affected by an additional physical parameter, the relative strength (S) of the stable stratification compared to the destabilizing buoyancy difference. A number of questions related to buoyancy reversal effects in thermal convection are related to the effect of varying S. Such questions, addressed in this paper, include: for which S is it safe to assume that the convection is not affected by the dynamic of the stable layer? How is the generation of internal waves affected by S? How does the height of the convection zone change with S? Is the heat transfer sensitive to S?

The idea of possibly decoupling the convection zone from the stable fluid layer is of particular interest because the convective and internal gravity wave dynamics can be on vastly different time and length scales (Le [START_REF] Bars | Experimental study of internal wave generation by convection in water[END_REF][START_REF] Lecoanet | Numerical simulations of internal wave generation by convection in water[END_REF]. In numerical simulations it would be advantageous to run a turbulent convection-only simulation, and then use the output flow statistics from this simulation as a source term for generating internal waves in a separate run. This idea was put forward theoretically by [START_REF] Goldreich | Wave generation by turbulent convection[END_REF]: assuming a Kolmogorov-type spectrum representative of Reynolds stresses due to sweeping plumes below the unstable/stable interface, they bypassed the simulation of the convection and were able to directly analyze in details the generated internal wavefield. Of course this approach requires a careful investigation of whether a Kolmogorov spectrum is appropriate and whether Reynolds stress is indeed the main generating mechanism. Ongoing work and recent results on water convection seem to support that waves are mainly generated by Reynolds stresses, at least for some (high) values of S [START_REF] Lecoanet | Numerical simulations of internal wave generation by convection in water[END_REF]. Another simplified approach consists in deriving one-dimensional (along the vertical) models with turbulence diffusivity parameterization, as done in, for instance, [START_REF] Wunsch | Stochastic simulations of buoyancy-reversal experiments[END_REF].

Here we address via direct numerical simulations some of the questions pertaining to the effect of the stable stratification strength on buoyancy reversal convection. Because our goal is to obtain a clear picture of the various convective dynamics for a wide range of stratification parameters, we focus in this paper on two-dimensional simulations. Reducing the system to two dimensions allows us to run a large number of simulations for several diffusive times, which is necessary to reach a statistically steady state. We describe the model setup in §2, and in §3 we analyze the onset and transition to steady state. The qualitative features of simulations with high, moderate, and small stratification strengths are then described in §4. A detailed analysis is presented in §5 and concluding remarks are offered in §6.

Problem formulation

We consider the Navier-Stokes equations under Boussinesq approximation and in two dimensions. We take x (z) to be the horizontal (vertical upward) direction, with x and ẑ the unit vectors. The dimensionless equations are

∂ t u + u • ∇u = -∇p + Pr∇ 2 u -PrRaρẑ, (2.1a) ∂ t T + u • ∇T = ∇ 2 T, (2.1b) ∇ • u = 0. (2.1c)
Length and time scales are non-dimensionalized using the full vertical extent of the domain and the thermal diffusion time 2 /κ, where κ is the constant thermal diffusivity of the fluid. T is the temperature, u the velocity vector, and ρ is the density anomaly compared to the reference density ρ 0 . We assume isothermal and no-slip boundary conditions on the top and bottom horizontal plates and periodicity of all fluid variables in the horizontal x direction. The dimensionless parameter Pr = ν/κ is the Prandtl number, with ν the constant fluid viscosity. Ra is the Rayleigh number, defined as

Ra = α * c g(T * b -T * i ) 3 νκ , (2.2)
where g is the constant gravitational acceleration, α * c is the (dimensional) constant thermal expansion coefficient for temperatures T * > T * i , and T * b -T * i > 0 is the (dimensional) temperature scale obtained from the difference between the temperature on the bottom plate (T * b ) and inversion temperature (T * i defined below). Using T * i as the absolute reference temperature then gives the dimensionless temperature T in (2.1) in terms of the dimensional one, T * , as

T = (T * -T * i )/(T * b -T * i ), implying T b = 1 and T i = 0.
As in classical RBC, we consider an equation of state for the density anomaly of the form ρ = -αT . However, we deviate from the OB paradigm as we allow the dimensionless thermal expansion coefficient α(T ) to vary with temperature, as is the case for e.g. water. Specifically, as can be seen in figure 1, we consider ρ = -α(T )T = -T, T ≥ 0, ST, T < 0, (2.3) such that α changes sign at T = 0, and with S > 0 a free parameter of our model. The change of sign of α results in the density anomaly being non-monotonic and maximum at T = 0, which we thus refer to as the inversion temperature (at the dimensional temperature T * i ). The form of α is the simplest model that can reproduce a non-monotonic EoS similar to water (shown in figure 1), and that allows varying degrees of stratification of the top stable layer. While S 1 is obtained when considering a stronger stratification than is the case for water, i.e. a steeper effect of varying temperature on the density, S 1 implies a weaker stratification. It can be noted that a smoothed piecewise profile for α may have worked equally well as the discontinuous profile considered in equation (2.3), though it would have introduced additional parameters in the problem. Numerically, equation (2.3) did not pose any problem, although a relatively high vertical resolution was necessary, because even though α is discontinuous, relevant terms appearing in the equations (i.e. T and its first-order and second-order derivatives as well as the buoyancy αT ) are all continuous. The variations of ρ with T are shown in figure 1 for S = 1, along with the stiffest S = 2 8 (almost horizontal line) and most flexible S = 2 -8 (almost vertical line) cases considered in this paper. For S = 2 8 , note that the stratification in density is so strong in the stable layer that the density anomaly cannot be shown in figure 1 even for a unitary increment in temperature; at the top of the domain where we will consider T t = -20, we will have ρ = -5120.

The system described by equations (2.1) and (2.3) is completely defined by Pr, Ra, S and the dimensionless top temperature T t (recall T b = 1). We will call S the stiffness parameter as it is related to the buoyancy resistance of the upper fluid to overshooting convective plumes. Defining the buoyancy frequency as N = √ -PrRa∂ z ρ (based on our notations), we can see that S is also related to the ratio of the buoyancy frequency in the stable and unstable fluid through Equations of state for the density anomaly ρ as a function of temperature T (T decreasing upward) for three different stiffness parameters S (solid lines). For comparison with water, the quadratic equation of state ρ = -T 2 is also shown (dashed line); note that the S parameter that best represents water depends on the temperature range considered.

S = (N 2 /∂ z T )| zs /(N 2 /∂ z T )| zc ,
temperature profile T = 1 -z(1 -T t ). In this state, the fluid is convectively unstable below z = 1/(1 -T t ), and stably stratified above, and the interface separating the unstable and stable regions corresponds to the T = 0 (inversion temperature) isotherm. We will also consider the neutral buoyancy level Z NB , the height at which fluid parcels rising adiabatically and without inertia from the bottom boundary would reach equilibrium.

In the next sections we solve equations (2.1) with ρ given by (2.3) using the pseudo-spectral code Dedalus1 (see Burns et al., (unpublished, and further details in appendix). We first briefly discuss the transient evolution ( §3), and then we explore the dynamical regimes and statistical properties of the coupled convectivestably stratified system at thermal equilibrium ( §4-5) for different stiffness parameters, S = 2 i with i ∈ [-8, 8] an integer, and three different reference Rayleigh numbers, i.e. Ra = 8 × 10 6 , 8 × 10 7 , 8 × 10 8 (details of all 51 simulations are provided in the appendix). For simplicity, the Prandtl number is fixed to Pr = 1, and the top temperature is fixed to T t = -20 in all simulations. The choice of T t = -20 is rather arbitrary but serves the purpose of having stable and convective layers of similar height for the prototypal case shown in figure 3b (Ra = 8 × 10 7 , S = 1), thus minimizing confinement effects on the convection and on wave propagation in the stably-stratified region. The box width relative to the height is fixed to 2, such that the effective aspect ratio for the convection is approximately 4, which limits horizontal confinement effects.

Transient evolution and the height of the convection zone

We first show in figure 2 the transient evolution of the temperature and vorticity fields between t/τ c = 2.1 and t/τ c = 50.5 starting from the conductive state T = 1 -z(1 -T t ) with white noise at t = 0 for physical parameters Ra = 8 × 10 7 and S = 1 (we recall that we set Pr = 1 and T t = -20 for all simulations in this paper); τ c = 1/f c = 2π/ √ RaPr is the reference turnover time scale for buoyancy forces, which will be used throughout (along with the turnover frequency f c ) to normalize time variables and vorticity (note that it can be obtained from equations (2.1) assuming small perturbations, upon substitution of ρ by (2.3), and assuming a unitary temperature gradient). After roughly two turnover times (figures 2a,c), the convective instability develops close to the bottom boundary layer; the convective cells then expand over time upward into the overlaying stable layer (figures 2b,d). The convective region expands because the convective flux in the convection zone is larger than the diffuse flux in the stable zone. When the convection expands, the convective flux stays about constant, but the diffusive heat flux in the stable region increases because the temperature gradient increases. The depth of the convection zone reaches a constant mean value when the heat flux through the convection zone is on average equal to the heat flux through the stable layer, which is known as the state of thermal equilibrium.

The characteristics of the convection vary with the height of the convective region. Determining the height of the convective region at thermal equilibrium is thus a major goal of this paper. The Rayleigh number defined by (2.2) is a reference Rayleigh number, but is not a strict analogue of the Rayleigh number in classical RBC because the height of the convective region is an output of the simulation, not an input. Indeed, in our definition of the Rayleigh number we use the temperature difference T * b -T * i and thermal expansion coefficient α * c related to the convection, but choose the total domain depth as length scale. The inversion depth /(1 -T t ) of the conductive state could be used alternatively, but from figure 2 and previous studies on convection in water close to its density maximum, we know that the convective layer depth grows with time from the static state, and that it reaches a mean height at dynamical equilibrium that can substantially vary with the problem parameters (see e.g. Large and Andereck , 2014, and figure 2). As a result, neither nor /(1 -T t ) can be expected to provide an accurate measure of the convective height. Providing a prediction for the convective layer depth including coupling between the two layers is the topic of section §5.6.

In a previous analysis of convection in water close to its density maximum, [START_REF] Moore | Nonlinear penetrative convection[END_REF] showed that a prediction of the (dimensionless) convective height h could be made under several assumptions. Assuming no dynamical coupling between the stratified upper layer and the convection, the method equates the diffusive heat flux in the upper layer (approximated as -(T t -T i )/(1 -h)) with the heat flux in the convection zone of the form CRa β ef f (T b -T i )/h, where Ra ef f is the Rayleigh number based on the effective convective depth. The form of the convective heat flux is suggested by classical RBC studies which show that for fixed Pr the Nusselt number, Nu, can be approximated as Nu = CRa β with C and β two constants. Solving for the convective depth h with Ra ef f = Rah 3 then yields (with T i = 0 and T b = 1)

1 -h = -T t h 1-3β CRa β , (3.1)
which, under the assumption that 1 -3β 1 (an acceptable approximation since 0.28 < β < 0.31 for laminar thermal boundary layers, cf. [START_REF] Grossmann | Scaling in thermal convection: A unifying theory[END_REF]; [START_REF] Ahlers | Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection[END_REF]), reduces to

1 -h ≈ -T t CRa β . (3.2)
The approximate expression (3.2) for 1 -h predicts that the convection height increases with the reference Rayleigh number Ra but decreases with T t . Both effects are expected since as T t decreases, the stably stratified layer grows and the unstable layer correspondingly shrinks in the conductive state. The prediction, however, is independent of the stratification parameter S, as the dynamical coupling between the stably stratified and convective regions are neglected. We will show that this simplification is appropriate only in the large stiffness regime (S 1), and that the convection height cannot be in general inferred from this model.

Dynamics at thermal equilibrium for different stiffnesses

We now turn our attention to the system dynamics at thermal equilibrium, i.e. once statistical steady-state and depth-invariant heat flux are achieved. Unlike in figure 2, in order to avoid the thermal equilibration time for the stable layer, we start our simulations from a temperature profile close to what is expected at thermal equibrium (see appendix).

We show in figure 3a-c three snapshots obtained at thermal equilibrium for (a) high S = 2 8 , (b) moderate S = 1, and (c) small S = 2 -8 stiffness. We select snapshots for Ra = 8 × 10 7 , since this intermediate Rayleigh number case will be the focus of the detailed analysis in section §5 (Pr = 1, T t = -20). In figures 3a,b, the density anomaly ρ is shown from the bottom up to the instantaneous neutral buoyancy level (i.e. z = Z NB (x), where ρ(x, Z NB ) = -1, labeled as the interface in figures 3), while the vorticity is shown above it. The idea of showing both density and vorticity such as in figure 3 is not new [START_REF] Lecoanet | Numerical simulations of internal wave generation by convection in water[END_REF], and allows to combine in a single plot the information on buoyancy effects in the convection zone and waves in the stable layer. In figure 3c, the density anomaly is shown everywhere because plumes can in theory rise all the way to the top of the domain; the stratification is indeed so weak that the lightest fluid is within the bottom boundary layer. Figure 3a, obtained for a high stiffness parameter, shows two distinct regions above and below Z NB , which coincides with the T i isotherm. The bottom layer consists of RB convective cells made of up-down symmetric plumes rising between the bottom plate and the neutral buoyancy level. The top layer is stably stratified and supports waves. The interface is almost perfectly straight and we expect that the waves are excited by Reynolds stress forcing [START_REF] Lecoanet | Numerical simulations of internal wave generation by convection in water[END_REF]. Global internal gravity modes are present, with approximately three wavelengths fitting along the vertical direction.

Figure 3b, obtained for stiffness S = 1, exhibits the same features as figure 3a. Plumes are emitted from the bottom plate and the interface, while internal waves are observed above Z NB . The aspect ratio of the convective cells is roughly 1, as in figure 3a, but because the interface lies higher, only 2 cells are obtained, compared to 3 for S = 2 8 . Compared to the high-stiffness case, however, the instantaneous interface is much more distorted as a result of upward-going plumes impinging on the relatively weak stable stratification. Rising and sinking plumes do not have the same symmetry either: down-going heavy plumes look like filaments, entraining chunks of lighter fluid with them. These observations are also reflected in the difference of location of the T i isotherm and neutral buoyancy level, which can be seen as a first approximation of the overshooting length scale of strong penetrative convection. Note that this overshooting length scale is defined from the mixing of stable fluid into the convective region (sinking of T i isotherm within the convection zone), and is therefore not exactly equivalent to the overshooting length scale based on the mixing of convective fluid into the stable region considered in penetrative convection studies on top of a background state. The more chaotic dynamical transition between the convective and stably-stratified regions clearly results in a less organized internal gravity wavefield than in the high-stiffness case of figure 3a.

Figure 3c, obtained for the smallest stiffness parameter S = 2 -8 considered, shows a very different behavior than the previous two cases. In this simulation, the bulk of the fluid is everywhere heavier than near the bottom plate. A stark difference also comes from the fact that there are only upward-going plumes. This is because the T i isotherm is within the thermal bottom boundary layer, so convection only occurs near z = 0. The sinking of the T i isotherm is an important result that can be explained from the fact that the top stable layer, which is initially quiescent in the simulations, is so heavy for small S that it can squeeze the buoyancy-driven convective layer down as soon as it is destabilized by overshooting convective plumes. Because the stable layer is only marginally stratified for small S, the entire stable fluid is put into motion and slumps down when the bottom part is destabilized by the plumes. This explains why entrained fluid motions extend almost everywhere in the fluid. Note that the cold entrained fluid moves seemingly passively, i.e. slower than the plumes, but because of its inertia and large temperature anomaly still dominates dynamically almost everywhere and brings the bulk density anomaly close to zero. Looking at the temperature field at the same simulation time (cf. figure 4) shows that hot rising plumes still impact high parts of the fluid as a result of their inertia. Thus, plumes play an essential role in the dynamics by maintaining the entrained fluid in motion, but not necessarily in the heat transfer. As they rise further up, plumes broaden (thereby carving the large-scale entrained fluid visible in figure 3c) and cool down to the point that they can contribute negatively to the (positive) heat transfer. Nonetheless, as will be demonstrated later, convection is efficient from z = 0 up to z = 0.9, because the large-scale entrained structures of low density anomaly (in part visible from contour lines of iso-density in figure 3c) are efficient at transporting cold fluid down. In the Above Z NB , the fluid is stably stratified, and we plot the vorticity. In (a) convective cells with a top/down symmetry are seen in the lower layer, generating global wave modes in the upper stably-stratified layer. In (b), the convection zone is larger but significant coupling between the two layers is obtained as shown by the strong variations of Z NB with x. In (c) the system changes significantly with rising plumes close to the bottom interacting with large-scale structures (shown by contour lines) higher in the fluid (note that Z NB is outside of the domain such that we show the density field everywhere in this case). We also show the mean temperature and anomaly density profiles by the dashed and solid lines as later discussed and reported in figure 6. Movies of figures 2 and 3 can be found online at sites.google.com/site/fludyco in the outreach section.

weakly-stratified limit, the overshooting length scale is of the same order as the entire fluid depth. Note that the top of the fluid layer is stably-stratified, although weakly.

Statistical analysis and regime properties

The simulation snapshots of figure 3 show that the stiffness parameter S plays a key role in the coupling between the convective and stably-stratified regions. The high-stiffness regime shows strong top-down symmetry in the convection zone, suggesting negligible feedback from the waves on the convection. This symmetry is no longer present for moderate stiffness S ∼ 1, and the dynamics completely change for low stiffness.

We now quantitatively study the effect of the stiffness S on the interaction between the convective and stably stratified regions and on the heat transfer. The results are obtained for simulations in thermal equilibrium, i.e. such that the statistics are temporally converged and the total heat flux is depth-invariant. Reaching thermal equilibrium requires running simulations over several dissipation times, which requires long run times. Use of judicious initial conditions allows a more rapid convergence to thermal equilibrium (see appendix for more details). We will make extensive use of the x average operator • x = dx/2, (2 is the box width); the volume average operator • = dzdx/2; and the time average operator, • . Time averages are typically performed over tens to thousands turnover times depending on S (see appendix for details). We first present in figure 5 simulation results obtained for 3 different Rayleigh numbers (Ra = 8 × 10 6 , 8 × 10 7 , 8 × 10 8 ) and 17 different stiffness parameters (S from 2 8 to 2 -8 ; see table 1 in appendix). In figure 5a we show the bulk temperature, i.e. T bulk = z dz T x , averaged in z over the convection zone (i.e. where the convective heat transfer q c accounts for at least 95% of the total heat flux q). For all three Rayleigh numbers we can see that the bulk temperature increases with the stiffness parameter. In the limit of large stiffness, T bulk → 0.5, which is the temperature expected from a convection-only simulation with top/bottom temperatures of T t = 0 and T b = 1. Thus the large-stiffness limit of the bulk temperature tends to the classical RBC bulk temperature even though the EoS is nonlinear. This is consistent with figure 3a, where the density anomaly in the bulk is roughly ρ = -0.5, corresponding to T bulk = 0.5. In the smallstiffness limit, the bulk temperature decreases below the T i = 0 isotherm, suggesting that significant mixing of low-temperature fluid occurs within the convective region. The overlap of the results for all Ra further indicates that this mixing is relatively independent of the degree of turbulence. As long as the lower layer is convectively unstable and plumes are emitted, they can destabilize the upper stable layer entirely, because it is only marginally stratified. For high stiffness, higher Ra tends to increase the entrainment/mixing and lower the bulk temperature. The data for all Ra are well approximated in the range S ∈ [2 -8 , 2 8 ] by a best polynomial fit law in terms of ln(S), i.e.

T bulk (S) = -0.0075 ln(S) 2 + 0.11 ln(S) + 0.1, (5. 1) shown as the dash-dot line in figure 5a. This empirical law will be used in §5.6 to provide an estimate for the depth of the convection zone.

We define the z-dependent convective and diffusive heat fluxes to be q c = wT x , (5.2a)

q d = -∂ z T x (5.2b)
and we recall that the total heat flux q = q c +q d can be shown to be independent of z when averaged over long time periods (recall that time and therefore w are normalized by the thermal time, and that Pr = 1). Figure 5b then shows the total heat transfer Q = qdz as a function of the stiffness parameter (z-averaging removes the small 2% maximum relative discrepancy obtained using the infinity norm in our simulations, showing the good statistical convergence). Because we work with dimensionless variables, we note that Q is normalized by a diffusive heat flux, and thus may be interpreted as a Nusselt number. We refrain from interpreting Q as a Nusselt number, however, because the Nusselt number definition requires a choice of length scale (see section §5.5). The heat transfer Q increases with decreasing stiffness S for all three Rayleigh numbers, which is in agreement with the observation that the diffusive heat flux in the stable layer should increase as a result of the growing convection zone (see figure 3). As expected, we also observe that Q increases with Ra for all S. In the high-stiffness limit, the slope |dQ/dS| → 0, which suggests that the classical RBC regime might be recovered for large S. This further implies that the depth of the convection zone h becomes independent of S since the heat flux in the stable layer is dominated by diffusion, so Q ≈ -T t /(1 -h). However, how large S should be to recover RBC clearly depends on Ra, since for the highest Ra simulations, Q still changes appreciably for S = 256. For small S, on the other hand, Q increases continuously with decreasing S. Thus, different mechanisms occur when varying S, as could be inferred from figures 3 and 4. 5.1. Vertical profiles of the mean temperature and density. In figure 6, the mean temperature T x and density ρ x profiles are shown as a function of z. Because the convection is steady for Ra = 8 × 10 6 at high S, and vertical confinement effects may affect the higher Ra = 8 × 10 8 simulations (due to the higher convective flux resulting in a larger h, and to overshoot), from here onward we choose to focus on the case Ra = 8 × 10 7 (obtaining converged statistics would also require a much longer time for Ra = 8 × 10 8 than for Ra = 8 × 10 7 ). The trends regarding the effect of S are, nevertheless, similar for all three Rayleigh numbers considered. We use colors ranging from blue to red to represent increasing stiffness. In figure 6a, we plot the temperature between -4 and 1 because the temperature decreases linearly with z for T x ≤ -4 (recall T t = -20). The highest-stiffness case (dark red curve) shows a symmetric temperature profile with z for T x ∈ [0, 1], i.e. between the bottom and inversion temperature, similar to what is obtained from a purely convection simulation (see unfilled-square symbols). The mean temperature in the well-mixed region is roughly 0.5, as observed in figure 5a. As the stiffness decreases, the mean temperature profiles become less symmetric; eventually the bulk temperature drops below 0. It is relatively surprising that for S ≤ 2 -2 the T i = 0 isotherm is found within the thermal bottom boundary layer, so the fluid is on average stablystratified from the outside of the bottom boundary layer all the way to the top of the domain. Nonetheless, it is in agreement with the observations made on figure 3c,i.e. that entrained cold fluid motions squeeze the buoyancy-driven convection zone close to the bottom boundary layer.

The stable stratification of the entire fluid domain, except the bottom boundary layer, for low stiffness S can be clearly seen in the (blueish) ρ x profiles of figure 6b. For S ≤ 2 -2 the density increases rapidly from -1 at z = 0 to 0, and then becomes roughly uniform or slowly decreasing with z (as observed in figure 3c). For high stiffness, on the other hand, large positive density gradients are obtained in two distinct regions separated by the bulk region, such that density increases over a relatively large vertical extent. It should be noted that because density anomaly gradients vary strongly where the dynamics change from stable to unstable, the buoyancy frequency N varies strongly there too. However, since we chose a constant background density and an EoS for which α is constant in the stable layer, N is mostly constant deep in the stable region. Had we chosen a non-uniform background density or an EoS for which N would be variable deep in the stable region on length scales large in comparison to the vertical wavelength of the wave, then wave propagation would follow WKB ray theory [START_REF] Dosser | Anelastic internal wave packet evolution and stability[END_REF]. Unstable (positive) density gradients are only found close to the bottom boundary for low stiffness S because buoyancy-driven convection is limited to a thin bottom layer when the stratification is weak. Nonetheless, it is clear from figure 6b that the vertical extent over which the density is approximately uniform always increases with decreasing stiffness S, such that there might be a mixed region of non-negligible extent for all values of S considered. We will confirm that convective motions take place over a volume fraction that grows monotonously with decreasing stiffness in the next section. 5.2. Heat transfer and entrainment. In order to explain the increase of Q with decreasing S from figure 5b, as well as the-first increase, then decrease-of the layer depth over which buoyancy-driven convection occurs (figure 6), we split the convective heat transfer q c into two contributions, i.e. (5.3b) where H is the Heaviside step function (recall T i = 0). The heat transfer q cp gives the contribution from convective motions occurring as a result of buoyancy effects, and is equivalent to the convective heat flux in classical Rayleigh-Bénard convection. q ce measures the contribution to the heat flux from fluid that is entrained (hence the superscript e), i.e. from fluid parcels expected to be in a stable region (T < T i ) and that yet have a net negative velocity. We recall that the sum q cp + q ce + q d = q = Q is constant for all z. . Vertical profiles of the heat transfer due to (a) convective plumes q cp , (b) entrainment q ce , and (c) diffusion q d for Ra = 8 × 10 7 . Different curves correspond to different stiffness S. For high stiffness, we obtain a symmetry in the heat transfer in the convection zone due to plumes (q cp ) in (a), suggesting classical Rayleigh-Bénard convection. The contribution from entrainment to the total heat transfer can be seen in (b) to be dominant in the low-stiffness cases and negligible in the high-stiffness regime. The total heat transfer q = q cp + q ce + q d is constant along the z axis up to 2% maximum relative discrepancy. As in figure 6, the unfilled-square symbols show results from a convection-only simulation with a vertical height comparable to the convection depth of the highest-stiffness case.

q cp = wT H [+(T -T i )] x , (5.3a) q ce = wT H [-(T -T i )] x ,
From figure 7a, we find for the high-stiffness cases (S ∈ [2 2 , 2 8 ]) a symmetric profile for q cp in the convection zone, with a rapid increase close to the bottom boundary and a rapid decrease higher up, where density gradients change sign (see figure 6b). This is a typical profile of the total convective heat flux in classical RBC (see unfilled-square symbols). As S decreases from 2 8 to 2 2 , however, q cp increases in parallel with increasing convective layer depth (which can be taken approximately as the height at which q cp decreases quickly). This is qualitatively different from RBC, which exhibits a weak dependence of q cp on h. For S ≤ 2 1 , q cp loses its symmetry and decreases as the stiffness decreases (gray to dark blue). For low stiffness, the buoyancy-driven heat flux reaches a maximum close to the bottom boundary layer, and decays rapidly with z away from it.

The dominant contribution to convective heat transfers does not come from buoyancy effects for low stiffness S, but from fluid motions in an overall stably-stratified environment, as can be seen from figure 7b. The heat transfer q ce is almost uniformly 0 for high-stiffness cases, but is large for low stiffness S; both the maximum of q ce and the vertical extent increase with decreasing stiffness. As S decreases from 2 8 to approximately 1, q ce becomes non-negligible close to the density inversion height as a result of plumes overshooting in the stable layer and inducing some entrainment of the stable fluid (T < T i ) (see e.g. [START_REF] Brummell | Penetration and Overshooting in Turbulent Compressible Convection[END_REF]. For low stiffness simulations, the entrained fluid has a density roughly equal to the density of the fluid just outside of the thermal boundary layer (see figure 6b), so the entrained fluid has a strong influence on the rising thermal plumes. Buoyancy-driven effects become confined closer and closer to the bottom boundary as S decreases, but at the same time plumes can entrain more of the stably-stratified fluid since the stable buoyancy frequency decreases. This competition results in higher heat transfers for smaller stiffness parameters (as seen in 5b), and a smaller diffusive upper layer close to the top boundary of the fluid (figure 7c). It should be noted that the efficiency of heat transfers by the entrained fluid depends on the temperature anomaly carried by the stable fluid. Here, the top temperature is T t = -20, which means that the entrained fluid can have temperature anomalies in the range [-20, 0]. This range is broader than the typical temperature anomaly of buoyant plumes of 0.5 to 1 (see e.g. figure 10), hence partly explaining the relative high values of q ce compared to q cp in spite of the relatively slow velocities of the entrained fluid (see section 5.4).

Figure 8. Volume fractions occupied by the bottom thermal boundary layer, the plume zone, the entrainment zone, and the stably-stratified top layer as functions of the stiffness parameter S and for Ra = 8 × 10 7 (see figure 15 in appendix for Ra = 8 × 10 8 ). The layer thicknesses are computed based on the profiles shown in figure 7 and according to equations (5.4). For each stiffness parameter, the position where the average temperature in the simulations is equal to the inversion temperature is shown by a star. It can be clearly seen that the inversion temperature isotherm is on average close to or within the bottom thermal boundary layer for S ≤ 2 -2 . The double-headed arrow on the left-hand side shows the convection layer thickness h prediction based on equation (3.1). The filled diamonds show the empirical prediction of h based on equation (5.9) presented in section 5.6.

The respective contributions from plumes, entrained fluid and diffusive effects to the heat transfer for varying stiffness are summarized in figure 8 for Ra = 8 × 10 7 (and in figure 15 in appendix for Ra = 8 × 10 8 ). Each region represents the volume fraction where either q cp , q ce , or q d , dominates over the others, and each region thickness is computed based on the corresponding heat flux variables as

L cp = 1 0 q cp Q dz, L ce = 1 0 q ce Q dz, L db = z * 0 q d Q dz z * , L ds = 1 z * q d Q dz 1 -z * . (5.4)
The thicknesses of the two diffusive regions, i.e. of the thermal bottom boundary layer and the upper stable layer, denoted by L db and L ds , are computed from the same diffusive heat flux q d , but we separate the integral values over the two zones at z = z * with z * taken in the middle of the mixed region (i.e. defined as {z, q d (z) ≤ 5%Q}), in order to avoid an overlap. The choice of z * is rather arbitrary, but results in an appropriate estimate for the depths of the bottom boundary layer and upper stable layer. Importantly, we have L cp + L ce + L db + L ds = 1, which is the total domain height.

Figure 8 clearly shows that plumes dominate convective heat transfers for high stiffness (left-hand-side of figure 8). For low stiffness (right-hand-side), however, it is the entrained fluid that dominates over most of the mixed region, and plume-driven heat fluxes are confined near the bottom plate. The vertical extent of the mixed region increases monotonically with decreasing S, and the stratified upper layer correspondingly shrinks, which is in agreement with increasing heat transfer (cf. figure 5). The white stars show the position of the T i isotherm in our simulations, and it can be seen that this isotherm sinks close to or within the thermal bottom boundary layer for S ≤ 2 -2 . The fact that the change of the T i position is much more sudden than the change of the plume region relative to the entrainment region indicates that plumes carrying heat upward are still somewhat effective despite the low bulk temperature; this is in agreement with figure 4 in which plumes that are still relatively hot can be seen as high as z ∼ 0.6. In the high-stiffness regime it is expected that the convective layer thickness h is well approximated by equation (3.1). This is verified in figure 8 as the predicted h (shown by the double-headed arrow), including the plume region as well as two symmetric thermal boundary layers, is in good agreement with the highest-stiffness S = 2 8 results (note that we use β = 0.27 and C = 1/5.5 in equation (3.1) based on convection-only simulation results that we ran with Dedalus and obtained for Rayleigh numbers in the range [10 6 , 10 9 ]). The fact that the convection height for S = 2 7 is different than for S = 2 8 , however, suggests that the stiffness must be quite large (> O(100)) to be in the asymptotic high-stiffness regime. 5.3. Kinetic energy density spectra. The regimes of entrainment-and plume-dominated heat transfer exhibit different types of flows and dynamical signatures. Figure 9 shows the temporal kinetic energy density spectrum averaged in x ( x operator), i.e. 2πf K x = πf û2 (f, x, z)+ ŵ2 (f, x, z) x , for the three simulations shown in figure 3 (• denotes temporal Fourier transform, f the frequency).

The kinetic energy density spectrum for the high-stiffness case (figure 9a) shows a convective layer with a mid-plane symmetry, similar to classical RBC. As expected, the energy peaks at frequency f /f c ∼ 1, where f c = √ RaPr/(2π) is the reference turnover frequency. The kinetic energy density decreases rapidly above the mean neutral buoyancy level (z ∼ 0.33, see figure 3a) in the stably stratified region. Diffusive and viscous damping is most rapid for low-frequency internal waves (such as those at the convection frequency), so only relatively high-frequency waves carry energy all the way to the top of the domain. The kinetic energy density is small above the buoyancy frequency, shown as the black solid line. Interestingly, the arced dashed lines in figure 9a, which highlight low-energy (z, f /f c ) paths, indicate the presence of nodes of standing high-frequency internal waves.

Figure 9b shows the kinetic energy density spectrum for the moderate-stiffness case. As for figure 9a, the energy density peaks at f /f c ∼ 1 in the convective region, but also at f /f c ∼ 0.04. A second peak at lower frequency can also be seen for the high-stiffness case (figure 9a), though less intense, that is due to a largescale circulation, which is stable on very long timescales. For the moderate-stiffness case, the second peak may be due to the large-scale circulation, but is also possibly related to the coupling of the convective and wave dynamics; it may be due to the slow motion downward of low-density-anomaly stable fluid entrained by the penetrating plumes. The buoyancy frequency (black solid line) in the moderate-stiffness case is relatively close to the reference turnover frequency f c , such that the kinetic energy density of the internal-wave field that peaks at f /f c ∼ O(1) quickly decreases with frequency. Compared to the high-stiffness case, internal waves are yet more strongly generated as a result of the penetrating plumes, complementing the generation by Reynolds stress due to weaker plumes sweeping beneath the neutral buoyancy interface (which is the main generation mechanism for stronger stratification, see [START_REF] Lecoanet | Numerical simulations of internal wave generation by convection in water[END_REF]). We note that in the case of water cooled from below (T * b = 0 • C) and heated from above (T * t = 25 • C), the kinetic energy density spectrum obtained shows similarities with both figure 9a and9b. This is because the stiffness parameter approximating the quadratic EoS for water (with

T * i = 4 • C) is then S ≈ (T * t -T * i )/(T * i -T * b ) = 5.
25. The kinetic energy density spectrum for the low-stiffness case (S = 2 -8 ) is shown in figure 9c. Fluid motions in the upper stable layer (z ≥ 0.8) have a lot of energy at frequencies larger than N , showing that even though the stratification is stable, the motions are not wave-like. This is in part due to the fact that the stratification is weak, such that large convective motions from the entrained fluid are not strongly affected by the restoring buoyancy forces. The maximum of kinetic energy density spectrum in the bulk is approximately two orders of magnitude smaller for S = 2 -8 compared to the high-stiffness case. This is because buoyancy-driven fast plumes are confined near the bottom boundary (the energy peak is close to 1 near z = 0), while slower larger-scale structures of lower frequencies dominate the fluid bulk (as already suggested from the large energy peak at ∼ 0.1 in figure 9b); this can also be seen from the decrease in the frequency of the kinetic energy density maximum from ∼ 1 close to the bottom boundary layer to ∼ 0.1 at z ∼ 0.2. 5.4. Probability density functions. Probability density functions (PDFs) of flow variables describe important aspects of the flow characteristics. In experiments, PDFs of the flow variables are accessible using thermometers and PIV techniques, and geophysicists can infer flow statistics for e.g. Earth's deep interior from measurements of the magnetic field. Here we show that different convective regimes exhibit different flow statistics.

We show in figure 10 the PDFs of the temperature for all our simulations with Ra = 8 × 10 7 . All PDFs are obtained based on temperatures interpolated on a uniform grid within the entire domain, and have been normalized such that the integral value P(T )dT = 1. As the stiffness increases and becomes large, the peak of the PDF tends toward T ≈ 0.5 and becomes more and more symmetric about its peak value, as is the case for convection-only simulations (see PDF shown by the unfilled-square symbols in figure 10). Moreover, the PDFs for the high-stiffness cases appear mostly exponential, in agreement with earlier studies [START_REF] Siggia | High Rayleigh Number Convection[END_REF][START_REF] Kerr | Rayleigh number scaling in numerical convection[END_REF]. As the stiffness S decreases, however, the peak shifts toward negative temperatures and the PDF becomes skewed toward the lower temperature values. The decrease for low stiffness of the peak temperature is in agreement with the low average temperature reported in figure 6, and the asymmetry tells us that the large-scale convecting structures have densities only marginally smaller than 0 (the maximum density, see figure 6b). This asymmetry is a feature of the entrainment regime, and is not due to sampling a subset of the total convective domain (see e.g. [START_REF] Emran | Fine-scale statistics of temperature and its derivatives in convective turbulence[END_REF]. Asymmetric temperature distributions have already been observed in laboratory experiments of convection in water (see Le [START_REF] Bars | Experimental study of internal wave generation by convection in water[END_REF], in which case the asymmetry is most likely related to entrainment and mixing of the stable fluid in the convection zone. Note that the oscillations of the temperature distributions for high-stiffness cases are most likely due to statistics that are not yet fully converged in the stable layer.

The PDFs of the vertical velocity w obtained for different stiffnesses are shown in figure 11a (for clarity, we do not include results from all simulations). The PDFs are for vertical velocities interpolated on a uniform grid within the convection zone only, and normalized such that P(w)dw = 1. As in classical RBC (cf. unfilled-square symbols), the PDF of w for high stiffness is symmetric with respect to positive/negative values. The PDF for low stiffness is, on the other hand, skewed towards the positive vertical velocities w > 0, and a similar albeit weaker asymmetry is also observed for S = 2 0 . The asymmetry for S = 2 -8 is toward high positive vertical velocities because plumes, which are much faster than the entrained fluid, are only ascending. That the velocity distributions for the more flexible cases are more peaked around small velocities than for the stiffer cases suggests that patches of the stable fluid that are entrained within the convection zone have small velocities. In an attempt to further characterize the dynamics of the plume-and entrainment-dominated convection, we show in figure 11b the PDFs of vertical velocity w after normalization with respect to their mean and standard deviation. For the high-stiffness cases, we see that the slow speeds (in absolute magnitude) follow an exponential distribution while the distribution at higher speeds decays like a Gaussian (cf. dashed line). For the low-stiffness cases the result is the opposite: slow speeds seem to follow a Gaussian distribution (cf. dotted line), while high speeds follow an asymmetric exponential distribution. The statistics for weak stratification suggest a mixed dynamics. Large-scale structures that are relatively slow (but with a Reynolds number ∼ 100 still relatively large based on a length scale of ∼ 1 and typical rms velocity of ∼ 100) have a velocity distribution similar to isotropic random fluctuations, possibly due to the fact that the motions of the entrained fluid are not driven by buoyancy and thus are isotropic. Fast fluid As S is increased, the peak of P(T ) tends to 0.5 (dotted line), the bulk temperature expected for classical RBC. As S is decreased, the peak of P(T ) drops below the inversion temperature (dashed line). The profile of P(T ) is symmetric about the peak for hight stiffness (black solid line shows exponential fit), but asymmetric for low stiffness. Small local maxima of P(T ) in the range T ∈ [0.5, 1] are expected to be representative of the plumes' temperature. As in figure 6, the unfilled-square symbols show results from a convection-only simulation with a vertical height comparable to the convection depth of the highest-stiffness case. Note that P(T ) for the convection-only results is normalized such that the integral value is equal to the integral value of P(T ) for S = 2 8 over the range

T ∈ [0, 1].
motions (i.e. plumes rising in the dense fluid bulk), on the other hand, are intermittent and strongly skewed towards w > 0. For the stiffest cases, fluid motions in the convection zone are expected to be turbulent, which is in agreement with the Gaussian distribution for the highest speeds. That the PDFs for the highstiffness cases (and also the convection-only simulations) are not purely Gaussian for slow speeds may be explained from the fact that we have compiled statistics for w not only in the bulk but also partly on the edge of the boundary layers, which are not as turbulent as the bulk.

Figures 12a-c show the joint PDFs of w and T for S = 2 8 , 2 0 , 2 -8 . For large stiffness (figure 12a), the joint PDF shows a symmetry for T ∈ [0, 1] around zero velocity and T = 0.5, because the heat transfer is due to both rising and sinking plumes. In figure 12b, the peak convective heat transfer occurs for values of T closer to 0, and the larger values of the joint PDF for T > 0 and w > 0 highlight the increased importance of hot rising plumes compared to cool descending ones. In the weak-stratification regime (figure 12c), low-velocity sinking fluid parcels are primarily responsible for the convective heat flux, in addition to the slightly-hot fast rising plumes (seen in figure 3c). The asymmetry of positive/negative vertical velocities is evident and again demonstrates that buoyancy-driven plumes only go upward. 5.5. Heat transfer scaling. We now discuss the Nusselt-Rayleigh scaling obtained from our simulations, which is an essential element of the study of convective heat transfers. The Nusselt number Nu is traditionally defined as the total heat transfer Q, normalized by the purely conductive heat flux, and thus Nu measures the heat transfer enhancement due to convective motions. In our case, the purely conductive heat flux across the domain boundaries is

-(T t -1), such that Nu = -Q T t -1 (5.5)
In classical RBC, it is generally assumed that one can write Nu(Ra, Pr) = CRa β Pr γ , with Ra and Pr the two (separable) problem parameters. Here, we have two additional parameters (S and T t ), but as we restrict Pr = 1 and T t = -20, we only consider Nu ≡ Nu(Ra, S).

Figure 13 shows Nu/Ra 0.28 as a function of S for all three reference Rayleigh numbers. The exponent 0.28 for Ra is obtained from a best fit power law regression of Nu for our simulations with Ra = 8×10 7 , 8×10 8 for For low stiffness (blue) the distribution looks Gaussian (dotted line) for slow speeds with asymmetric exponential tails. As in figure 6, the unfilled-square symbols show results from a convection-only simulation with a vertical height comparable to the convection depth of the highest-stiffness case. each S, and is within the range obtained in classical RBC [START_REF] Ahlers | Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection[END_REF]. We do not expect the RBC scaling to work equally well for all S, and the r 2 -value of this first fit is therefore small (i.e. ∼ 0.52) due to the relatively large dispersion of the exponents in the range [0.26, 0.31] obtained for different S. Nevertheless, there is then a relatively good collapse of the data for different S with exponent -0.086 (with final r 2 -value of 0.99). That the exponent is negative is expected since as the stiffness decreases we showed that the heat transfer increases. Note, however, that this exponent has no theoretical grounding yet. For the smaller Ra = 8 × 10 6 , Nu/Ra 0.28 is above the collapsed data, especially for the high-stiffness cases, because the convection pattern is steady. This is consistent with classical RBC studies which showed that the Nusselt number has a steeper scaling with Ra for smaller Ra, close to the onset of convection. (1 -T t ) (5.6) suggested by figure 13 can be useful in predicting the size of the convection zone in mixed convective-stablystratified systems. However, the Nusselt numbers obtained are small (we find Nu ∈ [1,7]), compared to those demonstrated in classical RBC, where Nusselt numbers have been reported in the range Nu ∈ [10,30] for Pr = 1 and Rayleigh numbers in 8 × 10 6 -8 × 10 8 (see e.g. [START_REF] Grossmann | Scaling in thermal convection: A unifying theory[END_REF]. In order to reconcile the expected enhancement of the heat transfer due to convective motions in our simulations in the high-stiffness limit with RBC results, we consider the rescaled Nusselt number (subscript C standing for convection)

Nu C = Q 1/h , (5.7)
where h = (2L db + L cp + L ce ) is a measure of the depth of the convection zone based on (5.4). Nu C thus corresponds to the heat flux normalized by the diffusive heat flux through the depth of the convective layer only.

In figure 14 we plot Nu C as a function of the effective Rayleigh number Ra ef f = Rah 3 (also characterizing this convective layer only). We find that Nu C ≈ CRa 0.28 ef f for each fixed S where C is a constant that depends on S. For the high-stiffness cases (S = 2 8 ) we find that Nu C is in the range 3-30 for Ra ef f ∈ [2×10 4 , 2×10 8 ], in relatively good agreement with classical RBC results (shown by the filled squares). For fixed reference Ra number, d log Nu C /d log(Ra ef f ) appears to approach 0.28 as S → ∞, similar to classical RBC, which suggests that the coupling between the two layers diminishes. For small stiffness, however, d log Nu C /d log(Ra ef f ) increases up to 1 for Ra = 8 × 10 7 and S = 2 -8 , indicating that the weakly stratified fluid significantly enhances the convective heat transfer. 5.6. Extension of Moore's model to include the effect of S. With the help of (5.6) we can now extend Moore's approximate model [START_REF] Moore | Nonlinear penetrative convection[END_REF] for the prediction of the convective layer depth. This is done by equating the diffusive heat flux in the stable region with the convective heat flux, i.e.

-[T t -T bulk (S)]

1 -h ≈ -Nu(S)(T t -1), (5.8) with consideration of the effect of S on Nu(S) and T bulk (S) (approximately the temperature at the top of the convective layer) through equations (5.6) and (5.1). We obtain h ≈ 1 -[T bulk (S) -T t ] 0.26Ra 0.28 S -0.086 , (5.9) and the prediction (5.9) (with T bulk substituted by (5.1)) is shown in figure 8 by the filled diamonds. It should be noted that in the high-stiffness limit Nu, T bulk and h seem to become independent of S, so (5.9) cannot be valid as S → ∞. That some variations of Nu are observed for S ∈ [2 6 , 2 8 ] shows that we are not yet in the high-stiffness regime even with S = 2 8 , and suggests that RBC is quantitatively recovered only when S O(100). The classical RBC scaling of heat transfer is nevertheless recovered for the highest stiffness considered, when choosing an appropriate effective Rayleigh number (see figure 14).

Conclusions

We have presented two-dimensional direct numerical simulation results of thermal convection with buoyancy reversal, and demonstrated that the convective and internal-wave dynamics can be strongly coupled depending on the stable stratification strength. Our results are expected to hold qualitatively for arbitrary profiles of the thermal expansion coefficient with temperature, despite only considering a simple piecewise linear EoS. Similarities with convection in water, for which the EoS is quadratic, have been demonstrated by comparing our results with previous numerical and laboratory experiments [START_REF] Lecoanet | Numerical simulations of internal wave generation by convection in water[END_REF]; Le [START_REF] Bars | Experimental study of internal wave generation by convection in water[END_REF].

The different dynamical regimes obtained for different S are best summarized by figure 8 (see also figure 15 for higher Rayleigh simulations), and the two most important findings are that (i) for weak stratification (S ≤ O(10 -1 )), entrained fluid motions that are maintained by rising plumes close to the bottom boundary dominate the convection zone, and (ii) classical RB convection is recovered only for S 100 since otherwise the influence of the stable layer on the mean temperature profiles and flow statistics is significant. For the smallest stiffnesses considered, the entire stable layer slumps down and becomes entrained by the buoyant plumes, such that the convective dynamics is qualitatively very different from what is expected from classical RBC or even penetrative convection. As S increases, however, entrained fluid motions vanish and the convective dynamics becomes again dominated by buoyant plumes; as expected, the overshooting length scale and mixing of the stable fluid decreases with increasing S. For large stratification (i.e. S ≥ O(100)), the convection is relatively similar to RBC, the convective layer depth can be estimated assuming no feedback from the stably-stratified fluid (see equation (3.2)), and flow statistics are similar to convection-only statistics.

In the high-stiffness limit, the convection is not affected much by the waves it generates, such that classical RBC statistics might be used as forcing in simulations of internal waves within the stable zone. Whether the theoretical approach in [START_REF] Lecoanet | Internal gravity wave excitation by turbulent convection[END_REF] or numerical method based on Reynolds stress [START_REF] Lecoanet | Numerical simulations of internal wave generation by convection in water[END_REF] is valid is, however, beyond the scope of this manuscript. In fact, as S decreases to S ≤ O(10), Reynolds stress may not be the only mechanism exciting internal waves, so it may be necessary to consider other mechanisms such as the mechanical oscillator effect (see e.g. Ansong and Sutherland , 2010, for a discussion of the mechanical oscillator effect in relation with the generation of atmospheric internal waves by storms).

We have limited our study to fixing the top temperature to T t = -20, which results in approximately equalsize convective and stable layers for many parameter values. The effects of changing T t are qualitatively predictable and are not expected to change the importance of the convective/wave coupling. Imagine a domain with larger T t (e.g., T t = -10). If the convection height is the same as for the case T t = -20, the temperature gradient in the stable layer is larger (i.e. less negative), such that the heat flux is smaller in the stable than in the convective layer. To obtain a thermal equilibrium the convection height h must thus be different. In particular, h should decrease as T t < 0 decreases (in agreement with equation (3.2)), because the diffusive heat flux of the upper layer increases with h (i.e. Q d ≈ -T t /(1 -h)) while the convective heat flux of the lower layer stays the same (at least in the high-stiffness case). Preliminary results for S = 1 have confirmed this claim. For T t = -20 we have h ≈ 0.60 and Q ≈ 48, whereas for T t = -40 we find h ≈ 0.27 and Q ≈ 55, and for T t = -10, we find h ≈ 0.89 and Q ≈ 44. Note that the increase of Q with decreasing T t suggests that the effect of T t on h is not quite as dramatic as indicated by (3.2), at least for S = 1.

The effect of the Prandlt number, which we did not explore (we set Pr = 1), is likely to be more complex than T t and will be reported elsewhere. Nonetheless, we might expect that decreasing Pr will result in stronger convective/wave coupling. Plumes overshooting from the convective region are expected to penetrate further into the stratified layer as a result of rapid thermal diffusion processes in the low Péclet limit. Significant overshooting could potentially result in turbulent motions within the stably-stratified layer, a topic of interest to astrophysical systems that would be worth investigating [START_REF] Brummell | The effect of the Prandtl number on penetrative convection[END_REF][START_REF] Garaud | The stability of stratified spatially periodic shear flows at low Péclet number[END_REF].

We have based our analysis on the stiffness parameter S, which is equal to the ratio of the thermal expansion coefficients, or N 2 /(dT /dz), in the stable and convective regions. In real systems, we can only measure the end states of nonlinear processes, which means that N 2 and dT /dz are expected to be both approximately zero in the well-mixed (convective) zone and hence difficult to measure. Extracting S from real systems might thus require knowing the thermal expansion coefficient within the convection based on the fluid composition directly, while measuring both N 2 and dT /dz in the stably-stratified zone. A different stiffness parameter, S H , has been used in previous astrophysical studies of compressible convection [START_REF] Hurlburt | Penetration below a convective zone[END_REF]. Anelastic simulations have a background temperature profile T (z), and S H is defined as the ratio of N 2 /(dT /dz) in the two regions, multiplied by the ratio of T (z) in different parts of the domain. For domains smaller than a typical pressure scale height (the Boussinesq regime), T (z) is almost constant, and thus S ≈ S H . Massive stars have stable layers with typical stiffness S H ∼ 10 5 , and therefore are good examples of coupled systems in which the stable layer has a large stiffness [START_REF] Shiode | The observational signatures of convectively excited gravity modes in main-sequence stars[END_REF]. Knowing whether the convection is similar to RBC or include some types of entrained motions would yet still require a self-consistent study, because the Rayleigh number in stars is so large that the required stiffness for the entrainment to be negligible should be very large too. In the case of Earth's liquid outer core, the stratification strength of the stable layer-if it exists-is still uncertain (possibly not too different from the Coriolis frequency, [START_REF] Vidal | Quasi-geostrophic modes in the Earth's fluid core with an outer stably stratified layer[END_REF]) and the Rayleigh number is relatively large [START_REF] Gubbins | The rayleigh number for convection in the earths core[END_REF]. Therefore, self-consistent studies might be required in the context of planetary interiors in order to estimate the importance of the stable layer's feedback on the convection. It is possible that convective motions in the oceans and atmosphere [START_REF] Marshall | Open-ocean convection: Observations, theory, and models[END_REF]; [START_REF] Parodi | A theory for buoyancy and velocity scales in deep moist convection[END_REF] may sometimes fall within the class of low-to-moderate stiffness convective regimes (i.e. including some entrainment), depending on the properties of the stably-stratified fluid [START_REF] Grise | A global survey of static stability in the stratosphere and upper troposphere[END_REF]; [START_REF] King | Buoyancy frequency profiles and internal semidiurnal tide turning depths in the oceans[END_REF].

The plume-and entrainment-dominated regimes have been shown to have unique dynamical signatures that can be of interest to astrophysicists, meteorologists, oceanographers, and geophysicists. Understanding how the flow signatures relate to the mean state of a fluid is of significant importance since measurements of, e.g., the stable stratification in planets and stars are typically difficult, and the present study can provide information on the dynamics of a system with a buoyancy reversal. We have neglected a number of physical effects that can have a significant impact on astrophysical and geophysical fluids. In the case of stars and planetary interiors, the mean background convective state does not have a constant density profile, but a constant entropy profile, and compressibility and magnetic effects can be significant. Non-uniform background density profiles can lead to upward buoyancy plumes stronger than downward ones (and viceversa), and strong variations of the buoyancy frequency even deep in the stable region, resulting in internal waves of varying amplitudes [START_REF] Dosser | Anelastic internal wave packet evolution and stability[END_REF]. Double diffusive effects may also play a major role in some cases. In the ocean, for instance, density is a function of both temperature and salt such that temperature gradients can be stabilizing while salt stratification is destabilizing, and vice-versa. Our generic model of convective-stably-stratified fluids provides a basis for future works aimed at considering such physical effects, but also spherical geometry and rotation, along with buoyancy reversal. The potential emergence of a mean flow in the stable layer similar to the Quasi-Biennial Oscillation of the Earth's equatorial stratosphere [START_REF] Baldwin | The quasi-biennial oscillation[END_REF] is also a topic of significant importance which will be explored in the near future. the stably-stratified fluid layer above needs ∼ O(1) thermal time to reach equilibrium. Setting up the initial temperature field close to thermal equilibrium in the stable region thus reduces the simulation time significantly by allowing us to obtain converged flow statistics rapidly, i.e. without having to wait several thermal times for the entire system to be at equilibrium. We note that the initial state of the system is of much less importance for a convective system without buoyancy reversal since in this case the transient initial phase consists of an exponential growth of the most unstable modes until rapid saturation. We would like to note that the cadence at which data must be output is controlled by the shortest of the convective time scale τ c = 2π/ √ RaPr and the internal-wave time scale τ w = τ c / √ S. In order to construct flow statistics that capture the fastest dynamics here we typically use min(τ c , τ w )/10 as the time step between two data outputs. It should be noted that when the output cadence is high (S large), statistics are at equilibrium with relatively short time averages. When the output cadence is long, however, longer time averages are necessary to ensure statistical steady-state. This implies that about the same number of data outputs are required to obtain converged statistics for all parameters S.

Ra

Ra

S n x × n z CFL dt 8 × 10 7 2 -8 512 × 256 0.7 O(7 × 10 -7 ) 8 × 10 7 2 -7 512 × 256 0.7 O(7 × 10 -7 ) 8 × 10 7 2 -6 512 × 256 0.7 O(7 × 10 -7 ) 8 × 10 7 2 -5 512 × 256 0.7 O(7 × 10 -7 ) 8 × 10 7 2 -4 512 × 256 0.7 O(7 × 10 -7 ) 8 × 10 1. (Cont.) List of all 51 simulations discussed in this paper. Ra is the reference Rayleigh number (see equation (2.2)). Simulations are run with a stiffness parameter S, and with n x Fourier and n z Chebyshev basis functions in the horizontal and vertical directions. Timesteps were taken using the Courant-Friedrichs-Lewy (CFL) condition, with a safety factor given in the table. dt is the typical time step, and each simulation was run for at least O(1) thermal time (about dt -1 timesteps). The aspect ratio of all simulations is 2, Pr = 1, and T t = -20.

Figure 1 .

 1 Figure1. Equations of state for the density anomaly ρ as a function of temperature T (T decreasing upward) for three different stiffness parameters S (solid lines). For comparison with water, the quadratic equation of state ρ = -T 2 is also shown (dashed line); note that the S parameter that best represents water depends on the temperature range considered.

Figure 2 .

 2 Figure 2. Transient evolution of (a)-(b) the temperature field, and (c)-(d) the normalized vorticity field ωτ c at times (left) t/τ c = 2.1 and (right) t/τ c = 50.5. The simulation is started from the static equilibrium of linear temperature profile, and the physical parameters are Pr = 1, Ra = 8 × 10 7 , S = 1, and T t = -20. τ c is the reference turnover time scale (see text).

Figure 3 .

 3 Figure 3. See caption on next figure.

Figure 3 .

 3 Figure 3. (Cont.) Simulation snapshots for Ra = 8×10 7 and stiffness (a) S = 2 8 , (b) S = 1, (c) S = 2 -8 representatives of the system's dynamic at thermal equilibrium, i.e. obtained after several thousands convective turnover times. The density anomaly is shown in the convection zone, below the neutral buoyancy height Z NB (x) (i.e. where ρ(x, Z NB ) = -1).Above Z NB , the fluid is stably stratified, and we plot the vorticity. In (a) convective cells with a top/down symmetry are seen in the lower layer, generating global wave modes in the upper stably-stratified layer. In (b), the convection zone is larger but significant coupling between the two layers is obtained as shown by the strong variations of Z NB with x. In (c) the system changes significantly with rising plumes close to the bottom interacting with large-scale structures (shown by contour lines) higher in the fluid (note that Z NB is outside of the domain such that we show the density field everywhere in this case). We also show the mean temperature and anomaly density profiles by the dashed and solid lines as later discussed and reported in figure6. Movies of figures 2 and 3 can be found online at sites.google.com/site/fludyco in the outreach section.

Figure 4 .Figure 5 .

 45 Figure 4. Snapshot of the temperature field corresponding to the density anomaly plot of figure 3c.

Figure 6 .

 6 Figure 6. (a) Mean temperature T x and (b) mean density anomaly ρ x profiles as functions of z for Ra = 8 × 10 7 and for different stiffness parameters S. For strong stratification, the density anomaly in the convection zone is constant within a bulk bounded by two thin layers where the density quickly increases, reminiscent of classical RBC. For weak stratification, the fluid density only increases within the bottom boundary layer and is almost uniform or weakly decreasing above. As in figure 5, the vertical dashed line in (a) highlights the inversion temperature T i = 0, and the dotted line indicates the bulk temperature expected from classical Rayleigh-Bénard convection. The unfilled-square symbols show results from a convection-only simulation with rigid no-slip top/bottom boundaries and a vertical height h = 0.32, which is an estimate of the depth of the convection zone for S = 2 8 (obtained from h = (2L db + L cp + L ce ), see (5.4)).

  Figure7. Vertical profiles of the heat transfer due to (a) convective plumes q cp , (b) entrainment q ce , and (c) diffusion q d for Ra = 8 × 10 7 . Different curves correspond to different stiffness S. For high stiffness, we obtain a symmetry in the heat transfer in the convection zone due to plumes (q cp ) in (a), suggesting classical Rayleigh-Bénard convection. The contribution from entrainment to the total heat transfer can be seen in (b) to be dominant in the low-stiffness cases and negligible in the high-stiffness regime. The total heat transfer q = q cp + q ce + q d is constant along the z axis up to 2% maximum relative discrepancy. As in figure6, the unfilled-square symbols show results from a convection-only simulation with a vertical height comparable to the convection depth of the highest-stiffness case.

Figure 9 .

 9 Figure 9. Spectrogram of the kinetic energy density spectrum 0.5f û2 + ŵ2 x as a function of frequency and z for each of the three simulations shown in figure 3 (Ra = 8 × 10 7 ; S = 2 8 , 2 0 , 2 -8 from left to right). The frequencies are normalized by the reference turnover frequency f c = √ RaPr/2π. The black solid (white dashed) line represents the mean value of the real (imaginary) part of the normalized buoyancy frequency, i.e. N/f c with N = √ -RaPr∂ z ρ x /2π based on our notations. The arced dashed lines in (a) highlight global internal-wave modes.

Figure 10 .

 10 Figure10. Probability density function of temperature P(T ) for different stiffness parameters S and Ra = 8 × 10 7 . As S is increased, the peak of P(T ) tends to 0.5 (dotted line), the bulk temperature expected for classical RBC. As S is decreased, the peak of P(T ) drops below the inversion temperature (dashed line). The profile of P(T ) is symmetric about the peak for hight stiffness (black solid line shows exponential fit), but asymmetric for low stiffness. Small local maxima of P(T ) in the range T ∈ [0.5, 1] are expected to be representative of the plumes' temperature. As in figure6, the unfilled-square symbols show results from a convection-only simulation with a vertical height comparable to the convection depth of the highest-stiffness case. Note that P(T ) for the convection-only results is normalized such that the integral value is equal to the integral value of P(T ) for S = 2 8 over the range T ∈ [0, 1].

Figure 11 .

 11 Figure 11. (a) Probability density function P of velocity w from the convective region for different stiffness parameters S and for Ra = 8 × 10 7 . Note that the convective region encompasses both the entrainment zone and the plume zone (cf. figure 8). (b) Same as (a) but normalized by the standard deviation σ and mean µ of the distributions. For high stiffness (red solid line), the distribution looks exponential for slow speeds (dashed line).For low stiffness (blue) the distribution looks Gaussian (dotted line) for slow speeds with asymmetric exponential tails. As in figure6, the unfilled-square symbols show results from a convection-only simulation with a vertical height comparable to the convection depth of the highest-stiffness case.

Figure 12 .

 12 Figure 12. Joint probability density function of temperature and vertical velocity for (a) S = 2 8 , (b) S = 2 0 , (c) S = 2 -8 and Ra = 8 × 10 7 . Note that the oscillations in (c) are most likely due to statistics that are not yet fully converged in the stable layer for the highest-stiffness case.

Figure 13 .

 13 Figure13. Variations of Nu/Ra 0.28 with S where Nu = -Q/(T t -1) for three reference Rayleigh numbers. The dashed line shows the best fit power law regression of Nu in terms of S for Ra = 8 × 10 7 , 8 × 10 8 .

Figure 14 .

 14 Figure 14. Rayleigh-Nusselt scaling based on Nu C = Qh and effective Rayleigh number Ra ef f = Rah 3 where h is an estimate of the convective layer depth (see equation (5.7) and text for details). The circles, triangles and stars correspond to reference Rayleigh number Ra = 8 × 10 6 , 8 × 10 7 , 8 × 10 8 , respectively. The dash-dot line has slope 0.28, while the dashed line has slope 1. Results of five convection-only simulations are shown by filled squares; note that the lowest Ra ef f convection-only simulation is stationary in time; also, the scaling exponent for the convection-only results is 0.28 (discarding the lowest Ra ef f case).

  7 2 -3 512 × 256 0.5 O(3 × 10 -7 ) 8 × 10 7 2 -2 512 × 256 0.5 O(3 × 10 -7 ) 8 × 10 7 2 -1 512 × 256 0.5 O(3 × 10 -7 ) 8 × 10 7 2 0 512 × 256 0.5 O(3 × 10 -7 ) 8 × 10 7 2 1 512 × 256 0.5 O(3 × 10 -7 ) 8 × 10 7 2 2 512 × 256 0.5 O(3 × 10 -7 ) 8 × 10 7 2 3 512 × 256 0.5 O(3 × 10 -7 ) 8 × 10 7 2 4 512 × 256 0.5 O(3 × 10 -7 ) 8 × 10 7 2 5 512 × 256 0.35 O(2 × 10 -7 ) 8 × 10 7 2 6 512 × 256 0.35 O(2 × 10 -7 ) 8 × 10 7 2 7 512 × 256 0.35 O(2 × 10 -7 ) 8 × 10 7 2 8 512 × 256 0.35 O(2 × 10 -7 )

Table 1 .

 1 S n x × n z CFL dt 8 × 10 6 2 -8 256 × 256 0.7 O(3 × 10 -6 ) 8 × 10 6 2 -7 256 × 256 0.7 O(3 × 10 -6 ) 8 × 10 6 2 -6 256 × 256 0.7 O(3 × 10 -6 ) 8 × 10 6 2 -5 256 × 256 0.7 O(3 × 10 -6 )8 × 10 6 2 -4 256 × 256 0.7 O(3 × 10 -6 ) 8 × 10 6 2 -3 256 × 256 0.7 O(3 × 10 -6 ) 8 × 10 6 2 -2 256 × 256 0.7 O(3 × 10 -6 ) 8 × 10 6 2 -1 256 × 256 0.7 O(3 × 10 -6 ) 8 × 10 6 2 0 256 × 256 0.5 O(10 -6 ) 8 × 10 6 2 1 256 × 256 0.5 O(10 -6 ) 8 × 10 6 2 2 256 × 256 0.5 O(10 -6 ) 8 × 10 6 2 3 256 × 256 0.5 O(10 -6 ) 8 × 10 6 2 4 256 × 256 0.5 O(10 -6 ) 8 × 10 6 2 5 256 × 256 0.5 O(10 -6 ) 8 × 10 6 2 6 256 × 256 0.5 O(10 -6 ) 8 × 10 6 2 7 256 × 256 0.5 O(10 -6 ) 8 × 10 6 2 8 256 × 256 0.5 O(10 -6 ) See caption below.

Table 1 .

 1 (Cont.) See caption below. × 10 8 2 6 1024 × 512 0.35 O(2 × 10 -8 ) 8 × 10 8 2 7 1024 × 512 0.35 O(2 × 10 -8 ) 8 × 10 8 2 8 1024 × 512 0.35 O(2 × 10 -8 )

	Ra	S	n x × n z	CFL	dt
	8 × 10 8 2 -8 512 × 256	0.5	O(10 -7 )
	8 × 10 8 2 -7 512 × 256	0.5	O(10 -7 )
	8 × 10 8 2 -6 512 × 256	0.5	O(10 -7 )
	8 × 10 8 2 -5 512 × 256	0.5	O(10 -7 )
	8 × 10 8 2 -4 512 × 256	0.5	O(10 -7 )
	8 × 10 8 2 -3 512 × 256	0.5	O(10 -7 )
	8 × 10 8 2 -2 512 × 256	0.5	O(10 -7 )
	8 × 10 8 2 -1 512 × 256	0.5	O(10 -7 )
	8 × 10 8 2 0	512 × 256	0.5	O(10 -7 )
	8 × 10 8 2 1	512 × 512	0.5 O(4 × 10 -8 )
	8 × 10 8 2 2	512 × 512	0.5 O(4 × 10 -8 )
	8 × 10 8 2 3	512 × 512	0.5 O(4 × 10 -8 )
	8 × 10 8 2 4	512 × 512	0.5 O(4 × 10 -8 )
	8 × 10 8 2 5	512 × 512	0.5 O(4 × 10 -8 )
	8 Table

Dedalus is available at http://dedalus-project.org.

Acknowledgement This research is funded by the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program through grant agreement No 681835-FLUDYCO-ERC-2015-CoG. D.L. is supported by a PCTS fellowship, and a Lyman Spitzer Jr. fellowship. This project has been carried out in the framework of the Labex MEC (ANR-10-LABX-0092) and of the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the program "Investissements dAvenir" supervised by the Agence Nationale de la Recherche. This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program "Investissements d'Avenir" supervised by the Agence Nationale de la Recherche, and to the HPC ressources of IDRIS (clusters Turing and Ada) under allocation A0020407543. Additional computing resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. This project was supported in part by NASA under TCAN grant number NNX14AB53G, and by the National Science Foundation under Grant No. NSF PHY-1125915.

Appendix A: Numerical method and simulation details

We use the open-source pseudo-spectral code Dedalus Burns et al. ((unpublished) to solve equations (2.1). The spectral decomposition consists of Chebyshev polynomial functions in ẑ direction and complex exponential functions (Fourier series) in the periodic x direction. We note n z and n x the number of modes in the ẑ and x directions, respectively. The use of Chebyshev polynomials along the vertical axis results in increased resolution near the top and bottom domain boundaries, such that an accurate treatment of the thermal and viscous boundary layers is obtained with relatively few spectral modes. A 2-step implicit/explicit Runge-Kutta scheme is used for time integration [START_REF] Ascher | Implicit-explicit Runge-Kutta methods for timedependent partial differential equations[END_REF]. We list all simulation runs whose results are presented in this paper in table 1 with relevant numerical and physical parameters.

In order to reach thermal equilibrium relatively rapidly, we initiate all numerical simulations at t = 0 with low-amplitude noise added to a background zero-velocity field (u = 0 ¯) and a background temperature field of the form (except those in section 3, which are started from the conductive state):

, z ≤ z conv , (6.1a)

In equation (6.1a) z conv is an estimate of the convective layer depth and T bulk the bulk temperature both obtained from preliminary simulations at low resolution, which allow us to initialize the problem close to thermal equilibrium. An accurate estimate for z conv and T bulk results in an accurate estimate for the heat flux -(T t -T bulk )/(1 -z conv ), and hence the temperature field in the stably-stratified layer, which is the fluid region that is longest to reach equilibrium. While the convection reaches a dynamical equilibrium relatively quickly (i.e. in ∼ O(100) turnover times, which corresponds to O(0.01) thermal time for Ra = 10 8 ),