1. Introduction
Subglacial lakes are water environments trapped between ice sheets and bedrocks (Siegert et al. Reference Siegert, Ellis-Evans, Tranter, Mayer, Petit, Salamatin and Priscu2001). Over 400 subglacial lakes have been identified in Antarctica (Wright & Siegert Reference Wright and Siegert2012) and approximately 50 have been detected in Greenland (Bowling et al. Reference Bowling, Livingstone, Sole and Chu2019). Antarctica has 250 subglacial lakes that are stable, i.e. with water trapped for millions of years and in complete isolation from the Earth's climate. The remaining 150 are hydrologically active, i.e. are connected through networks of subglacial channels and communicate via filling and discharge with the surrounding ocean (Smith et al. Reference Smith, Fricker, Joughin and Tulaczyk2009). Here, we focus on stable subglacial lakes, which are of considerable interest to astrobiology since they could host microorganisms that might have developed novel survival strategies relevant to the oceans of icy moons (Cockell et al. Reference Cockell, Bagshaw, Balme, Doran, Mckay, Miljkovic, Pearce, Siegert, Tranter, Voytek and Wadham2011).
Subglacial lakes are heated by the Earth's geothermal flux, hence they are prone to vertical convection and can experience dynamic conditions. The water circulation in isolated subglacial lakes can also be driven or be affected by horizontal temperature gradients along the ice–water interface when it is tilted, due to the pressure dependence of the freezing temperature (Wells & Wettlaufer Reference Wells and Wettlaufer2008). The slope of the ice–water interface is typically of the order of or smaller than 
$10^{-2}$ (Siegert Reference Siegert2005), although here we will assume for simplicity that the ice–water interface is flat. Salt concentration levels are expected to be of the order of 0.1 
$\%$ or less in most subglacial lakes, such that the water is typically fresh (Siegert et al. Reference Siegert, Ellis-Evans, Tranter, Mayer, Petit, Salamatin and Priscu2001). A hypersaline lake has, however, been recently identified in the Canadian Arctic (Rutishauser et al. Reference Rutishauser, Blankenship, Sharp, Skidmore, Greenbaum, Grima, Schroeder, Dowdeswell and Young2018), suggesting that high salt concentrations remain possible. Subglacial lakes differ from ice-covered lakes because they typically have a much thicker ice cover and because they do not experience radiative heating (Ulloa, Wüest & Bouffard Reference Ulloa, Wüest and Bouffard2018).
Subglacial lakes under a thick ice cover, i.e. such as Lake Vostok, which lies beneath 4 km of ice (Siegert et al. Reference Siegert, Ellis-Evans, Tranter, Mayer, Petit, Salamatin and Priscu2001), are known to be unstable to vertical convection because the thermal expansion coefficient of water, 
$\beta$, is always positive at high pressures. Subglacial lakes under less than approximately 3 km of ice, such as Lake CECs (Rivera et al. Reference Rivera, Uribe, Zamora and Oberreuter2015), may on the contrary be stable against vertical convection because 
$\beta <0$ at low temperatures and for pressures lower than 
$p_*\approx 2848$ dbar (Thoma et al. Reference Thoma, Grosfeld, Smith and Mayer2010). Couston & Siegert (Reference Couston and Siegert2021) recently proposed that the geothermal flux, which is of the order of 50 mW m
$^{-2}$, is large enough to trigger convection in most subglacial lakes despite the nonlinearity of the equation of state. Convection typically occurs when the geothermal flux 
$F$ forces a bottom temperature in the static state 
$\bar {T}_b>T_d$, with 
$T_d$ the temperature of maximum density, such that 
$\beta (\bar {T}_b)>0$, which is a condition met by most lakes deeper than a few metres.
The existence of a density maximum at temperature 
$T_d>T_f$ with 
$T_f$ the freezing temperature means that low-pressure subglacial lakes self-organize into a lower convective layer coupled to an overlaying stably stratified fluid region. This two-layer dynamics has been extensively studied at atmospheric pressure, in which case 
$T_d\approx 4\,^{\circ }$C, both numerically (Lecoanet et al. Reference Lecoanet, Le Bars, Burns, Vasil, Brown, Quataert and Oishi2015; Toppaladoddi & Wettlaufer Reference Toppaladoddi and Wettlaufer2018; Wang et al. Reference Wang, Zhou, Wan and Sun2019) and experimentally (Large & Andereck Reference Large and Andereck2014; Léard et al. Reference Léard, Favier, Le Gal and Le Bars2020). Here, using direct numerical simulations (DNS), we investigate the turbulent dynamics of freshwater environments for different ice overburden pressures, 
$p_i$, which enclose and include 
$p_*$. Thus, our results generalize the study of two-layer freshwater systems to arbitrary ice overburden pressure.
An important point is that we consider a fixed top freezing temperature and fixed bottom heat flux condition, such that our boundary conditions are different from the classical isothermal top and bottom boundary conditions considered in the canonical Rayleigh–Bénard problem as well as by most numerical studies of mixed convective and stably stratified fluids (Couston et al. Reference Couston, Lecoanet, Favier and Le Bars2017; Toppaladoddi & Wettlaufer Reference Toppaladoddi and Wettlaufer2018; Wang et al. Reference Wang, Zhou, Wan and Sun2019). Laboratory and numerical experiments have shown that convection driven by a bottom isothermal boundary is statistically equivalent to convection driven by a bottom fixed-flux boundary (assuming a top isothermal boundary in both cases), provided that the temperature-based Rayleigh number is the same in both experiments and the fluid is in the Oberbeck–Boussinesq regime, i.e. its properties are independent of flow velocity and temperature (Verzicco & Sreenivasan Reference Verzicco and Sreenivasan2008; Johnston & Doering Reference Johnston and Doering2009). We will show that the same is true for subglacial lakes, even though they are not in the Oberbeck–Boussinesq regime, provided that an effective Rayleigh number is considered. We will demonstrate that there exist two limiting behaviours of the dimensional variables with the input heat flux and water depth depending on whether the thermal expansion coefficient 
$\beta$ is quasi-constant or linearly varying with the temperature difference driving convection. Importantly, our results support the idea that the convective and stably stratified layer dynamics are decoupled at leading order, which is a hypothesis that was recently invoked in order to predict flow velocities in Antarctic subglacial lakes (Couston & Siegert Reference Couston and Siegert2021).
We organize the paper as follows. We present the equations and numerical experiments in § 2. We analyse the DNS results and present the theoretical predictions in § 3. We discuss the geophysical implications in § 4 and conclude in § 5.
2. Problem formulation
2.1. Governing equations in dimensional form
We consider a Cartesian coordinate system 
$(x,y,z)$ centred on the lake's bottom boundary with 
$\boldsymbol {e}_z$ the upward-pointing unit vector of the 
$z$ axis, i.e. opposite to gravity, and we denote 
$H$ the ice thickness and 
$h$ the lake water depth (cf. figure 1a). For computational expediency we restrict our attention to two-dimensional motions, i.e. we assume 
$y$ invariance and neglect rotation. Here, as in most liquids, compressibility effects are weak and density fluctuations with temperature and pressure are small compared to the reference density 
$\rho _0=999$ kg m
$^{-3}$. As a result, the evolution of the lake's velocity vector 
${\boldsymbol {u}}$ and temperature 
$T$ is well approximated by the Navier–Stokes equations in the Boussinesq approximation and the incompressible energy equation, i.e. such that
\begin{gather} \partial_t {\boldsymbol{u}} - \nu\nabla^2{\boldsymbol{u}} + \boldsymbol{\nabla} (p/\rho_0) ={-} ({\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}){\boldsymbol{u}} - (\rho/\rho_0) g\boldsymbol{e}_z, \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{u}} = 0, \end{gather}
\begin{gather}\partial_t T - \kappa \nabla^2 T ={-} ({\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}) T, \end{gather}
where 
$p$ is the pressure, 
$\rho$ is the density, 
$\partial _t$ denotes time derivative and 
$\boldsymbol {\nabla }$ is the gradient operator. The physical parameters in (2.1) are the kinematic viscosity 
$\nu$, the reference density 
$\rho _0$, the gravitational acceleration 
$g$ and the thermal diffusivity 
$\kappa$ (cf. table 1). For the boundary conditions, we consider
\begin{equation} {\boldsymbol{u}}(z=0)={\boldsymbol{u}}(z=h)=\boldsymbol{0}, \quad \partial_zT(z=0)={-}F/k, \quad T(z=h)=T_f(p_i), \end{equation}
i.e. we assume no slip, fixed heat flux 
$F$ on the bottom boundary with 
$k$ the thermal conductivity and we set the temperature at the top of the lake equal to the temperature of freezing, 
$T_f$, which varies with the ice overburden pressure 
$p_i$.
Figure 1. (a) Problem schematic. The green shading highlights the region of the water column that stays stably stratified when 
$p_i \leq p_*$. (b) Thermal expansion coefficient 
$\beta (T,p)$. The solid black (respectively red) line shows 
$T_f$ (respectively 
$T_d$) while the dashed line shows the 
$p_*$ isobar. The arrows highlight the ice overburden pressures considered in the simulations. (c) Density variations with depth at 
$t=0$ for each of the four simulation cases 
$\mathcal {S}_i^1 (i=0,1,2,3)$ of the first experiment corresponding to the different ice pressures 
$p_i$ shown by arrows in (b).
Table 1. Physical parameters and polynomial approximations for 
$T_f$, 
$T_d$, 
$\rho _1$ and 
$C$ with temperatures in 
$^{\circ }$C, pressures in dbar, densities in 
$\textrm {kg}\ \textrm {m}^{-3}$ and 
$C$ in 
$\textrm {kg}\ \textrm {m}^{-3}\ ^{\circ }\textrm {C}^{-2}$.
We approximate the equation of state for the density of freshwater as a function of the lake pressure 
$p\geq p_i$ and temperature 
$T\geq T_f(p_i)$ using the bivariate polynomial
\begin{equation} \rho(p,T)=\rho_0 + \rho_1(p) + C(p)[ T-T_d(p)]^2, \end{equation}
where 
$T_d$ is the temperature of maximum density, i.e. such that 
$(\partial \rho /\partial T)|_p(T=T_d)=0$. We obtain (quadratic) polynomial expressions for 
$\rho$ (through 
$\rho _1$ and 
$C$), 
$T_d$ and 
$T_f$ as functions of pressure by minimizing their 
$\ell^2$ relative error norm compared to the exact thermodynamic values 
$\rho ^e$, 
$T_d^e$ and 
$T_f^e$ (superscript 
$^e$ denoting exact values) computed using TEOS-10 (McDougall & Barker Reference McDougall and Barker2011). The polynomial approximations for 
$\rho$, 
$T_d$ and 
$T_f$ (provided in table 1) result in errors smaller than 0.1 g m
$^{-3}$ and 0.002
$\,^{\circ }$C for 
$p,p_i\in [0,10\ 000]$ dbar and 
$T\in [T_f,T_f+15\,^{\circ }$C]. Figure 1(b) shows the pressure dependence of 
$T_f$ (solid black line) and 
$T_d$ (red line). Both 
$T_f$ and 
$T_d$ decrease with increasing pressure, but 
$T_d > T_f$, i.e. such that the water is densest at a non-freezing temperature, only for 
$p < p_*= 2848.5$ dbar, which we call the critical ice overburden pressure (dashed blue line). The form of the equation of state (2.3) highlights that the density can be non-monotonic with temperature and exhibits a maximum at 
$T=T_d$ within the water column provided that 
$T_d(p)>T_f(p_i)$. The condition 
$T=T_d$ requires 
$p_i < p_*$ and is most likely to be satisfied near the top of subglacial lakes since 
$p\geq p_i$ increases with depth by hydrostasy and 
$T_d$ decreases with 
$p$. The thermal expansion coefficient is
\begin{equation} \beta ={-}\frac{1}{\rho_0}\left.{\frac{\partial \rho}{\partial T}}\right|_{p} ={-}\frac{2C(p)[T-T_d(p)]}{\rho_0}. \end{equation}
Figure 1(b) clearly shows that 
$\beta >0$ for all temperatures when 
$p>p_*$, while 
$\beta$ can change sign with temperature for 
$p < p_*$, i.e. when the temperature of maximum density exceeds the freezing temperature. The critical ice-cover thickness associated with 
$p_*$ is 
$H_*=10^4p_*/(\rho _ig)=3166$ m (
$p_*$ in dbar) assuming a mean ice density 
$\rho _i=917$ kg m
$^{-3}$.
The density 
$\rho$ and the temperature of maximum density 
$T_d$ are functions of the full pressure 
$p$ (cf. table 1). However, for simplicity, here we will substitute 
$\rho (p,T)$ and 
$T_d(p)$ with 
$\rho (p_i,T)$ and 
$T_d(p_i)$ in the governing equations, i.e. such that 
$\rho$ and 
$T_d$ depend on the ice overburden pressure only. This approximation is legitimate for lakes that are not too deep, i.e. such that hydrostatic pressure variations are weak and considering 
$\rho (p,T)\approx \rho (p_i,T)$ and 
$T_d(p)\approx T_d(p_i)$ does not impact significantly buoyancy effects. All lakes considered in this work are shallow, i.e. the water depth does not exceed 8 metres, such that the approximation is valid. In particular, a simulation of a lake with a maximum depth of 8 metres yields almost identical results whether we make or relax the assumption 
$\rho (p,T)\approx \rho (p_i,T)$ and 
$T_d(p)\approx T_d(p_i)$ (cf. the Appendix). Note that approximating 
$\rho (p,T)\approx \rho (p_i,T)$ implies approximating 
$\beta (p,T)\approx \beta (p_i,T)$ as well.
Our study of natural convection in subglacial lakes is fundamentally a study of non-Oberbeck–Boussinesq (NOB) effects in thermal convection due to a temperature-dependent thermal expansion coefficient (2.4). Previous works on NOB effects due to a temperature-dependent thermal expansion coefficient that can change sign include Couston et al. (Reference Couston, Lecoanet, Favier and Le Bars2017); Toppaladoddi & Wettlaufer (Reference Toppaladoddi and Wettlaufer2018) and Wang et al. (Reference Wang, Zhou, Wan and Sun2019). Toppaladoddi & Wettlaufer (Reference Toppaladoddi and Wettlaufer2018) and Wang et al. (Reference Wang, Zhou, Wan and Sun2019) considered the equation of state for water at constant atmospheric pressure, i.e. such that their range of 
$\beta <0$ was fixed, whereas it varies with pressure in our case (see, e.g. figure 1b). Couston et al. (Reference Couston, Lecoanet, Favier and Le Bars2017) used a piecewise-linear equation of state and a variable stiffness parameter, which allowed them to consider different ranges for 
$\beta <0$. Our work is different from Couston et al. (Reference Couston, Lecoanet, Favier and Le Bars2017) because (i) we consider the full equation of state for water rather than an artificial equation of state and (ii) the bottom boundary condition is fixed heat flux in our work rather than fixed temperature, which we will show is an important point when 
$\beta$ varies with 
$T$. The dependence of viscosity 
$\nu$ and thermal diffusivity 
$\kappa$ with temperature are two other well-known NOB effects that can lead to noticeable deviations of thermal convection, including a top–down asymmetry, from the classical Rayleigh–Bénard experiment (Ahlers et al. Reference Ahlers, Brown, Araujo, Funfschilling, Grossmann and Lohse2006; Sugiyama et al. Reference Sugiyama, Calzavarini, Grossmann and Lohse2009). Nevertheless, here, we take 
$\nu$ and 
$\kappa$ as constants since their relative variations do not exceed 50 % over the range of 
$(p,T)$ considered, i.e. 0 dbar
$< p < 10^4$ dbar and 
$-5\,^{\circ }\textrm {C} < T < 5\,^{\circ }\textrm {C}$ (Forst, Werner & Delgado Reference Forst, Werner and Delgado2000; Huber et al. Reference Huber, Perkins, Friend, Sengers, Assael, Metaxa, Miyagawa, Hellmann and Vogel2012) and are not expected to have an effect on the lakes’ dynamics as important as the variations of the thermal expansion coefficient. Future works might consider relaxing this assumption.
2.2. Governing equations in dimensionless form
We use the water depth 
$h$ as the characteristic length scale, the diffusive time 
$\tau _{\kappa }=h^2/\kappa$ as the time scale, the velocity 
$u_{\kappa }=h/\tau _{\kappa }$ as the velocity scale and the pressure 
$p_{\kappa }=\rho _0u_{\kappa }^2$ as the pressure scale in order to identify dimensionless control parameters and non-dimensionalize the governing equations, which we recall are (2.1), (2.2a–c) and (2.3) with 
$\rho (p_i,T)$ substituted for 
$\rho (p,T)$. The temperature difference between the lake's top and bottom boundaries in the turbulent regime is unknown. However, we can use 
${\varDelta } = Fh/k$ as the temperature scale, which is the temperature difference across the lake's depth in the diffusive base state, which we denote by overbars and is given by 
$\bar {{\boldsymbol {u}}}=\boldsymbol {0}$, 
$\bar {T}=T_f+{\varDelta }(1-z/h)$ and hydrostatic pressure 
$\bar {p}=p_i+\int _z^h \bar {\rho }gdz'$. We use 
$T_f(p_i)$ as reference temperature and 
$p_i+[\rho _0+\rho _1(p_i)]g(h-z)$ as pressure gauge, i.e. such that we remove the leading-order mean buoyancy and hydrostatic pressure terms, which balance each other, in the governing equations. The dimensionless variables, which we denote by tildes, are then related to the dimensional variables through
\begin{equation} \left.\begin{gathered} (x,z)=h(\tilde{x},\tilde{z}),\quad t=\tau_{\kappa}\tilde{t}, \quad u=u_{\kappa}\tilde{u}, \\ p = p_i+[\rho_0+\rho_1(p_i)]g(h-z) + p_{\tau}\tilde{p}, \quad T = T_f + {\varDelta} \tilde{T}.\end{gathered} \right\} \end{equation}
Substituting (2.5a–e) into (2.1) and (2.2a–c) combined with (2.3) (with 
$p_i$ replacing 
$p$ in the expression for 
$\rho$), yields a set of dimensionless equations and boundary conditions, which we write as
\begin{gather} \partial_{\tilde{t}} \tilde{{\boldsymbol{u}}} - Pr\tilde{\nabla}^2\tilde{{\boldsymbol{u}}} + \widetilde{\boldsymbol{\nabla}} \tilde{p} ={-} (\tilde{{\boldsymbol{u}}}\boldsymbol{\cdot}\widetilde{\boldsymbol{\nabla}})\tilde{{\boldsymbol{u}}} + Pr\overline{Ra}_F\frac{( 1+\bar{S})}{2}\left[ \tilde{T}- \frac{\bar{S}}{( 1+\bar{S})} \right]^2\boldsymbol{e}_z, \end{gather}
\begin{gather}\partial_{\tilde{t}} \tilde{T} - \tilde{\nabla}^2 \tilde{T} ={-} (\tilde{{\boldsymbol{u}}}\boldsymbol{\cdot}\widetilde{\boldsymbol{\nabla}}) \tilde{T}, \end{gather}
\begin{gather}\widetilde{\boldsymbol{\nabla}}\boldsymbol{\cdot} \tilde{{\boldsymbol{u}}} = 0, \end{gather}
\begin{gather}\tilde{{\boldsymbol{u}}}(\tilde{z}=0)=\tilde{{\boldsymbol{u}}}(\tilde{z}=1)=\boldsymbol{0}, \quad \partial_{\tilde{z}}\tilde{T}(\tilde{z}=0)={-}1, \quad \tilde{T}(\tilde{z}=1)=0, \end{gather}with
\begin{equation} Pr=\frac{\nu}{\kappa}, \quad \overline{Ra}_F=\frac{gh^4F\bar{\beta}_b}{k\nu\kappa}, \quad \bar{S}=\frac{T_d-T_f}{\bar{T}_b-T_d}, \end{equation}the control parameters, and where
\begin{gather} \bar{\beta}_b={-}2C(\bar{T}_b-T_d)/\rho_0, \end{gather}
\begin{gather}\bar{T}_b=T_f+{\varDelta}, \end{gather}are the bottom thermal expansion coefficient and the bottom temperature of the diffusive base state, respectively.
The control parameters (2.7a–c) are the (constant) Prandtl number 
$Pr=12.8$, the base-state flux Rayleigh number 
$\overline {Ra}_F$, which is based on the heat flux 
$F$ and the bottom thermal expansion coefficient of the diffusive base state 
$\bar {\beta }_b$, and the base-state stiffness number 
$\bar {S}$, which compares the thermal expansion coefficient at the top of the lake, i.e. at temperature 
$T_f$, to the thermal expansion coefficient at the bottom in the diffusive base state, i.e. at temperature 
$\bar {T}_b$. The base-state flux Rayleigh number 
$\overline {Ra}_F$ is positive and convection is possible if 
$\bar {\beta }_b>0$, i.e. if the bottom temperature of the diffusive base state exceeds the temperature of maximum density. This condition is satisfied provided that the heat flux exceeds a minimum heat flux for fixed 
$p_i$ and 
$h$, or, equivalently, the water depth exceeds a minimum water depth for fixed 
$p_i$ and 
$F$, which we call the threshold heat flux and threshold water depth, respectively, and define as
\begin{align} F_t = max \left\{ \frac{k[ T_d(p_i)-T_f(p_i) ]}{h},0 \right\} , \quad h_t = max \left\{ \frac{k[ T_d(p_i)-T_f(p_i) ]}{F}, 0 \right\} , \end{align}
i.e. such that 
$F_t>0$ and 
$h_t>0$ if 
$p_i < p_*$ and 
$F_t=h_t=0$ if 
$p_i\geq p_*$. When 
$Fh>k[T_d(p_i)-T_f(p_i)]$, 
$\overline {Ra}_F$ increases monotonically with 
$p_i$, 
$h$ and 
$F$. The base-state stiffness number 
$\bar {S}$ can take any value in 
$[-\infty ,+\infty ]$. Specifically, 
$\bar {S}\leq -1$ when the lake is fully stable, i.e. 
$F < F_t$; 
$-1 < \bar {S}\leq 0$ when the lake is fully unstable, i.e. 
$p_i>p_*$; and, 
$\bar {S}>0$ when the lake is partially convective, i.e. 
$T_f < T_d < \bar {T}_b$, and self-organizes into a stably stratified upper layer and a convective bottom layer. When 
$\bar {S}>0$, we might expect that larger stiffness parameters 
$\bar {S}$ correlate with stronger resistance of the top stable layer to overshooting convective motions (see, e.g. Couston et al. Reference Couston, Lecoanet, Favier and Le Bars2017). When the lake is fully unstable, i.e. 
$-1< \bar {S} \leq 0$, 
$\bar {S}\approx -1$ indicates that the base-state thermal expansion coefficient is almost depth invariant, while 
$\bar {S}\approx 0$ indicates that it strongly varies with depth.
2.3. Numerical experiments
Since we impose the heat flux rather than the temperature on the bottom boundary of the lake, we expect the mean bottom temperature 
$\langle T_b \rangle =\langle T(z=0) \rangle$, with 
$\langle \cdot \rangle$ denoting the horizontal and temporal average, to become smaller than 
$\bar {T}_b$ as convection sets in. This means that the control parameters 
$\overline {Ra}_F$ and 
$\bar {S}$, while fully prescribing the system, do not provide an effective measure of buoyancy forcing compared to dissipation or density stratification in the stable layer when the flow is turbulent and at statistical steady state. As a result, our goals are to:
(i) investigate the variations of the bottom temperature
$\langle T_b \rangle$ at statistical steady state with the problem parameters;(ii) define an effective Rayleigh number
$Ra_{eff}$ based on $\langle T_b \rangle$, which can be used to predict the characteristic Reynolds number $Re$ and Nusselt number $Nu$ of the convective layer;(iii) investigate the influence of the stable layer on the convective dynamics through the use of, e.g. an effective stiffness parameter
$S_{eff}$.
We note that many studies have investigated how the classical scalings of Rayleigh–Bénard convection (between isothermal plates) change when changing the boundary conditions (e.g. fixed heat flux) or considering NOB effects (Chillà & Schumacher Reference Chillà and Schumacher2012). However, previous studies considering the effect of a fixed-flux boundary condition have been limited to the Oberbeck–Boussinesq regime (Verzicco & Sreenivasan Reference Verzicco and Sreenivasan2008; Johnston & Doering Reference Johnston and Doering2009), while those exploring the effect of a density maximum have used isothermal boundaries (Couston et al. Reference Couston, Lecoanet, Favier and Le Bars2017; Toppaladoddi & Wettlaufer Reference Toppaladoddi and Wettlaufer2018; Wang et al. Reference Wang, Zhou, Wan and Sun2019; Léard et al. Reference Léard, Favier, Le Gal and Le Bars2020). Thus, we expect that the analysis presented in this paper can provide new fundamental results on non-classical Rayleigh–Bénard convection while being also useful to the study of subglacial lakes.
Although parametric fluid dynamics studies are usually optimally designed when sweeping through parameters in dimensionless space, here, we explore the effect of the parameters on the flow dynamics by sweeping through the physical space 
$(p_i,F,h)$ rather than through the control parameter space 
$(\overline {Ra}_F,\bar {S})$. The reason is that we are interested in the variations and predictions of the bottom temperature 
$\langle T_b \rangle$ and flow velocities in the convective layer in terms of a specific range of either geophysically relevant or laboratory relevant lake parameters 
$p_i$, 
$F$ and 
$h$, which are difficult to cover with an exploration in 
$(\overline {Ra}_F,\bar {S})$ space due to the nonlinear relationships between 
$\overline {Ra}_F$, 
$\bar {S}$ and 
$F$, 
$h$ and 
$p_i$. Thus, we conduct our investigation of subglacial lake dynamics by sweeping through lines of constant heat flux and lines of constant water depth in physical space 
$(F,h)$ while considering ice pressures both above and below the critical ice pressure 
$p_*$, which separates the fully convective regime from the partially convective one (figure 2).
Figure 2. Graphical illustration of our two sets of numerical experiments in physical space with (a) 
$h=0.5$ m fixed in the first experiment and (b) 
$F$ (in W m
$^{-2}$) fixed in the second (cf. details in table 2). Thinner lines correspond to a larger ice overburden pressure. (c) Corresponding coverage in dimensionless space 
$(\bar {S},\overline {Ra}_F)$. The red shading highlights the region of fully convective lakes.
We consider subglacial lakes under 4 different ice overburden pressures, i.e. 
$p_i=0$, 2387, 2848.5 and 3549 dbar, corresponding to ice thicknesses 
$H=0$ (infinitesimally small ice layer), 
$H=2653$ m (relevant for subglacial Lake CECs, cf. Rivera et al. Reference Rivera, Uribe, Zamora and Oberreuter2015), 
$H=H_*$ and 
$H=3945$ m (relevant for subglacial lake Vostok, cf. Siegert et al. Reference Siegert, Ellis-Evans, Tranter, Mayer, Petit, Salamatin and Priscu2001), respectively. For each ice pressure considered we investigate the temperature variations in the lake and the root-mean-square velocity as functions of the geothermal flux 
$F$ and water depth 
$h$. We use two sets of experiments. First we focus on the case 
$h=0.5$ m and increase 
$F$ in successive stages. We denote the corresponding simulation cases 
$\mathcal {S}_i^1$, with 
$i=0,1,2,3$ increasing as 
$p_i$ increases (four multi-stage simulations), i.e., such that e.g. 
$\mathcal{S}_0^1$ corresponds to the simulation with 
$h=0.5$ m and 
$p_i=0$ dbar and with F increasing in successive stages. We pick 
$h=0.5$ m, which is a relatively standard height for water containers, such that the first set of simulations may be compared with future laboratory experiments provided that water can be pressurized. Second, we fix 
$F$ and increase 
$h$ in stages. We denote these simulations 
$\mathcal {S}_i^2$, again with 
$i=0,1,2,3$ increasing as 
$p_i$ increases (cf. table 2). Each stage of a simulation (with both h and F fixed) lasts one diffusive thermal time such that the results, averaged over the second half of a stage, describe the system at statistical steady state. Figure 2 highlights the physical parameter space covered by the numerical simulations as well as the corresponding coverage in dimensionless space.
Table 2. Dimensional parameters for the two sets of numerical experiments, which consider four distinct ice overburden pressures 
$p_i$ (in dbar) each and either a broad range of geothermal fluxes (7th column) or a broad range of water depths (last column). Ice thickness 
$H$ and water depths 
$h$, 
$h_c$ and 
$h_t$ are in metres and fluxes 
$F$, 
$F_c$ and 
$F_t$ are in W m
$^{-2}$; s.n. means simulation name. Note that figure 2 provides a graphical illustration of the dimensional and dimensionless parameter spaces explored.
For each simulation we first compute the threshold heat flux 
$F_t$ if 
$h$ is fixed (first set of experiments) or the threshold water depth 
$h_t$ if 
$F$ is fixed (second set of experiments) using (2.10a,b). Then, we evaluate the critical heat flux 
$F_c$ (respectively critical water depth 
$h_c$), which is required for the destabilizing buoyancy force to overcome viscous dissipation and thermal diffusion in (2.6) for the first (respectively second) set of experiments. For the calculation of 
$F_c$ and 
$h_c$ we use the eigentools package (https://github.com/jsoishi/eigentools) in Python, which is based on the eigenvalue-solver capability of the open-source pseudo-spectral code Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). Necessarily, 
$F_c > F_t$ and 
$h_c>h_t$. We report 
$F_t$ and 
$F_c$, as well as 
$h_t$ and 
$h_c$, and the range of supercritical heat fluxes and water depths considered for each simulation case in table 2. At 
$t=0$, we initialize the system with no velocities and a conductive (linear) temperature profile superimposed with small-amplitude white noise. The corresponding mean density profiles are shown in figure 1(c) for the four simulations of the first experiment. For both 
$\mathcal {S}_0^1$ and 
$\mathcal {S}_1^1$, for which 
$p_i < p_*$, the density increases with height, such that the conductive state is unstable to convection, in a lower subregion of the water column but decreases with height, hence is stably stratified, above. For 
$\mathcal {S}_2^1$ and 
$\mathcal {S}_3^1$ the density always increases with height such that the full water column is unstable to convection. For 
$\mathcal {S}_2^1$, 
$p_i=p_*$ and 
$T_f=T_d$ such that 
$\beta =0$ at 
$z=h$, which is why 
$\partial _z\rho =0$ at the top boundary.
We solve (2.6) with the open-source pseudo-spectral code Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020). We assume that the 
$x$ direction is periodic and has dimensional length 
$L_x=4h$. We recall that we use no-slip boundary conditions, an isothermal top boundary and a fixed heat flux bottom boundary. The horizontally averaged dynamic pressure, i.e. in excess of the hydrostatic pressure, is set to 0 at the top boundary. We use a Fourier basis with 
$n_x=512$ modes in the 
$x$ direction and a Chebyshev basis with 
$n_z=256$ in the 
$z$ direction before dealiasing for the most turbulent simulations. For the least turbulent simulations with a stable layer, i.e. 
$\mathcal {S}_0^1$, 
$\mathcal {S}_1^1$, 
$\mathcal {S}_0^2$, 
$\mathcal {S}_1^2$ with 
$F\leq 4F_t$, 
$F\leq 20F_t$, 
$h\leq 0.8$ m and 
$h\leq 1.3$ m, respectively, we decrease 
$n_x$ to 256. For the least turbulent fully convective simulations, i.e. 
$\mathcal {S}_2^1$, 
$\mathcal {S}_3^1$, 
$\mathcal {S}_2^2$, 
$\mathcal {S}_3^2$ with 
$F\leq 1250F_c$, 
$F\leq 2\times 10^5F_c$, 
$h\leq 4$ m and 
$h\leq 2$ m, respectively, we decrease both 
$n_x$ to 256 and 
$n_z$ to 128. We use a second-order two-step Runge–Kutta method for time integration and a Courant-Friedrichs-Lewy constant (CFL) between 0.2 and 0.4, with the lower CFL used for the most turbulent simulations.
3. Results
3.1. General flow features
The flow dynamics in a subglacial lake with a dimensionless thermal expansion coefficient 
$\tilde {\beta }=\beta /\bar {\beta }_b$, which changes sign within the water column (provided that 
$\tilde {T}_d>\tilde {T}_f=0$), is qualitatively different from the flow dynamics in a subglacial lake with 
$\tilde {\beta }$ positive throughout. Note that tildes denote dimensionless variables and that we will normalize all thermal expansion coefficients by 
$\bar {\beta }_b$ (cf. (2.8)). Figures 3(a) and 3(b) show several snapshots of the dimensionless temperature field 
$\tilde {T}$ for simulation 
$\mathcal {S}_{1}^1$, for which 
$\tilde {\beta }$ changes sign inside the water column, and simulation 
$\mathcal {S}_{3}^1$, for which 
$\tilde {\beta }>0$ everywhere. In figure 3(a), the water column is only partially unstable to convection. Convective motions (associated with 
$\tilde {\beta }>0$), which are shown with the red-to-blue colour map, coexist with a stably stratified layer (where 
$\tilde {\beta }<0$), which is shown by the yellow-to-green colour map. As 
$F$ increases (from top to bottom), the system transitions from a stationary laminar state to a turbulent state, and, at the same time, the bottom convective layer grows while the top stable layer shrinks. In figure 3(b), 
$\tilde {\beta }>0$ everywhere, such that the water column is convecting over the full depth and the flow dynamics is qualitatively similar to classical Rayleigh–Bénard convection.
Figure 3. (a) Snapshots of dimensionless temperature 
$\tilde {T}$ for simulation 
$\mathcal {S}_1^1$ with 
$F$ increasing from top to bottom. Here, 
$\tilde {T}_d$ is the dimensionless temperature of maximum density whereas 
$\langle \tilde {T}_b\rangle _h$ is the dimensionless horizontally averaged (but time-dependent) bottom temperature. (b) Same as (a) but for 
$\mathcal {S}_3^1$. (c) Dimensionless time and horizontally averaged temperature profiles 
$\langle \tilde {T} \rangle$ with depth at different stages (i.e. different 
$F$) for 
$\mathcal {S}_1^1$. (d) Same as (c) but for 
$\mathcal {S}_3^1$. The line colours go from dark to light as 
$F$ increases from small to large values (lines shifting from right to left as shown by the black arrows). The thin black lines show the conductive profiles at 
$\tilde {t}=0$. Panels (e) and (f) show the time evolution of 
$\langle \tilde {T}_b\rangle _h$ for 
$\mathcal {S}_1^1$ and 
$\mathcal {S}_3^1$, respectively. The vertical dashed lines highlight the times 
$\tilde {t}=i$ (
$i=1,2,3\ldots$) when the control parameter (
$F$ or 
$h$) starts increasing (smoothly) and the simulation stage changes, with a new statistical steady state reached before 
$\tilde {t}=i+0.5$.
We show in figures 3(c) and 3(d) the vertical profiles of the time and horizontally averaged dimensionless temperature 
$\langle \tilde {T} \rangle$ for each value (stage) of 
$F$ considered in simulations 
$\mathcal {S}_{1}^1$ and 
$\mathcal {S}_{3}^1$. When 
$F < F_c$, i.e. 
$F$ is subcritical, there is no motion nor mixing such that the dimensionless temperature profile is fully conductive, i.e. 
$\tilde {T}=1-z/h$, as shown by the black solid lines. As 
$F>F_c$ is increased (dark to light colours; following the direction of the arrow), convective motions emerge, intensify and mix the lake's unstable bulk more and more efficiently. The increased mixing results in a decreasing temperature of the lake's bulk and a decreasing temperature of the bottom boundary. For simulation 
$\mathcal {S}_1^1$ (figure 3c), the temperature profile remains conductive, i.e. linearly decreasing, in the top stably stratified layer where convective motions are inhibited (because 
$\tilde {\beta }<0$). The well-mixed convective region is small compared to the stably stratified region in 
$\mathcal {S}_1^1$ initially, but the situation reverses as 
$F$ increases. For simulation 
$\mathcal {S}_3^1$ (figure 3d), convection occurs everywhere such that 
$\langle \tilde {T} \rangle$ has a top–down symmetry and the bulk temperature is approximately the average of the top and bottom temperatures. We display the time history of the horizontally averaged bottom temperature 
$\langle \tilde {T}_b \rangle _h = \langle \tilde {T}(z=0) \rangle _h$ of simulations 
$\mathcal {S}_1^1$ and 
$\mathcal {S}_3^1$ in figures 3(e) and 3(f), respectively. The bottom temperature decreases in smooth steps every time the heat flux (or, alternatively, the water depth for the second experiment) increases. There are 10 stages in simulation 
$\mathcal {S}_1^1$, each lasting one diffusive time, and 7 stages in simulation 
$\mathcal {S}_3^1$. All simulations with ice overburden pressure 
$p_i < p_*$ are qualitatively similar to simulation 
$\mathcal {S}_1^1$, whereas simulations with 
$p_i\geq p_*$ are qualitatively similar to simulation 
$\mathcal {S}_3^1$. Note that time averaging of all variables of interest is performed over the second half of each simulation stage.
3.2. Effective temperature difference and thermal expansion coefficient driving the convection
The mean temperature on the bottom boundary is a key output of the simulations since it gives the range of temperatures involved in convective motions and contributing to the heat transport. The effective temperature difference driving the convection, which we denote by 
$\varDelta _{eff}$, may be taken as the difference between the mean bottom temperature and the maximum of the temperature of maximum density and freezing temperature, i.e. in dimensionless form,
\begin{equation} \tilde{\varDelta}_{eff}=\varDelta_{eff}/\varDelta=\langle\tilde{T}_b\rangle-\tilde{T}_d(\tilde{T}_d>0), \end{equation}
as 
$\tilde {\beta }>0$ for 
$\langle \tilde {T}_b \rangle \geq \tilde {T}\geq \langle \tilde {T}_b \rangle -\tilde {\varDelta }_{eff}$. Note that the term 
$(\tilde {T}_d>0)$ in (3.1) is to be understood as a Heaviside function (in fact, all greater than or less than signs in between parentheses should be understood as Heaviside functions in this paper), i.e. such that it is 1 if 
$\tilde {T}_d>0$ and 0 otherwise. For simulations with 
$\tilde {T}_d<0$, i.e. which are fully convective, the dimensionless effective temperature difference is simply equal to the dimensionless mean bottom temperature. For simulations with 
$\tilde {T}_d>0$, there can be no convection in the temperature range 
$0<\tilde {T}<\tilde {T}_d$, such that the effective temperature difference is equal to the mean bottom temperature minus the temperature of maximum density.
Since the mean bottom temperature is the highest (on average) temperature in the lake, it sets not only the effective temperature difference driving the convection but also the maximum value of the thermal expansion coefficient, which we write in dimensionless form as 
$\tilde {\beta }_b=\tilde {\beta }(\langle \tilde {T}_b\rangle )={\beta }_b/\bar {\beta }_b$ with subscript 
$_b$ denoting bottom variables. We recall that 
$\tilde {\beta }$ is also a function of 
$p_i$. However, the ice overburden pressure is fixed for each simulation, such that its influence on 
$\tilde {\beta }$ is not shown for simplicity. The effective thermal expansion coefficient 
$\tilde {\beta }_{eff}$ can be taken as the average between the bottom (maximum) thermal expansion coefficient and the thermal expansion coefficient at the top of the convective layer, which is 
$\tilde {\beta }_{f}=\tilde {\beta }(\tilde {T}_f=0)>0$ if the lake is fully convective and 0 if the lake has a stable layer (since in this case the mean temperature at the top of the convective layer is 
$\tilde {T}_d$), viz.
\begin{equation} \tilde{\beta}_{eff} = \frac{\tilde{\beta}_b+\tilde{\beta}_f(\tilde{T}_d<0)}{2}. \end{equation} We show in figure 4(a–d) the evolutions of the dimensionless effective temperature difference 
$\tilde {\varDelta }_{eff}$ and thermal expansion coefficient 
$\tilde {\beta }_{eff}$ as we increase the heat flux 
$F$ (for the first experiment) or the water depth 
$h$ (for the second experiment) in the simulations. Figure 4(a) shows that 
$\tilde {\varDelta }_{eff}$ decreases monotonically with the normalized heat flux 
$(F-F_t)/(F_c-F_t)$ in all simulations (of the first experiment). Two asymptotic behaviours, highlighted by the solid lines, emerge at relatively large values of 
$(F-F_t)/(F_c-F_t)$. The asymptotic behaviour is the same for simulations 
$\mathcal {S}_i^1$ (
$i=0,1,2$) but is different for simulation 
$\mathcal {S}_3^1$. A similar result is obtained with the second experiment, as shown in figure 4(b), i.e. 
$\tilde {\varDelta }_{eff}$ decreases monotonically with the normalized water depth 
$(h-h_t)/(h_c-h_t)$ and displays two asymptotic behaviours, although the difference between the asymptotic behaviour for 
$\mathcal {S}_i^2$ (
$i=0,1,2$) and 
$\mathcal {S}_3^2$ is tenuous (which is expected, as we will demonstrate in § 3.3). The origin of the two different asymptotic behaviours for 
$\tilde {\varDelta }_{eff}$ can be related to the evolution of 
$\tilde {\beta }_{eff}$ with the normalized heat flux and water depth shown in figure 4(c) and 4(d), respectively. On the one hand, the effective thermal expansion coefficient 
$\tilde {\beta }_{eff}$ decreases monotonically and displays a common asymptotic behaviour with the normalized heat flux for simulations 
$\mathcal {S}_i^1$ with 
$i=0,1,2$ (figure 4c) and with the normalized water depth for simulations 
$\mathcal {S}_i^2$ with 
$i=0,1,2$ (figure 4d). On the other hand, 
$\tilde {\beta }_{eff}\approx 1$ for simulations 
$\mathcal {S}_3^1$ (figure 4c) and 
$\mathcal {S}_3^2$ (figure 4d), although it can be seen that 
$\tilde {\beta }_{eff}$ starts decreasing with the normalized heat flux at large values for simulation 
$\mathcal {S}_3^1$ (figure 4c).
Figure 4. (a,b) Dimensionless effective temperature difference 
$\tilde {\varDelta }_{eff}$ (i.e. driving the convection) as a function of (a) the normalized geothermal flux 
$(F-F_t)/(F_c-F_t)$ for the simulations of the first experiment and (b) the normalized water depth 
$(h-h_t)/(h_c-h_t)$ for the simulations of the second experiment (cf. legends and table 2). (c,d) Same as (a,b) but for the dimensionless effective thermal expansion coefficient 
$\tilde {\beta }_{eff}$. (e–h) Show the same variables as (a–d) but in dimensional form, i.e. with 
$\varDelta _{eff}=\varDelta \tilde {\varDelta }_{eff}$ and 
${\beta }_{eff}=\bar {\beta }_b\tilde {\beta }_{eff}$. The solid lines show scaling laws as discussed in § 3.3 and listed in table 4. Note that the symbol size and line thickness are inversely proportional to 
$p_i$, i.e. large (respectively small) symbols and thick (respectively thin) lines highlight results for small (respectively large) 
$p_i$.
In order to understand why 
$\tilde {\beta }_{eff}$ either stagnates or decreases with the normalized heat flux or water depth, it is useful to look at the dimensional effective temperature difference 
$\varDelta _{eff}=\varDelta \tilde {\varDelta }_{eff}$ and the dimensional effective thermal expansion coefficient 
$\beta _{eff}=\bar {\beta }_{b}\tilde {\beta }_{eff}$ shown in figure 4(e–h). Figures 4(e) and 4(f) show that 
$\varDelta _{eff}$ increases in all simulations. This happens because 
$\tilde {\varDelta }_{eff}$ decreases more slowly with increasing 
$F$ or 
$h$ (which increase mixing), than the temperature scale 
$\varDelta =Fh/k$ increases with 
$F$ or 
$h$. Since 
$\varDelta _{eff}$ increases in all simulations, the bottom temperature also increases, and so does the thermal expansion coefficient 
$\beta _{eff}$ (figures 4g and 4h). However, 
$\beta _{eff}=\bar {\beta }_b[\varDelta _{eff}+2(T_f-T_d)(T_d < T_f)]/[2(\bar {T}_b-T_d)]$ is an affine function of 
$\varDelta _{eff}$ (cf. (3.2)). Thus, while we might expect that the effective thermal expansion coefficient always scales asymptotically like 
$\beta _{eff}\sim \varDelta _{eff}$ as 
$\varDelta _{eff} \rightarrow \infty$, there exists a range of effective temperature differences, or bottom temperatures, i.e. 
$0<\langle T_b\rangle -T_f\ll T_f-T_d$, when 
${T}_d < T_f$ for which 
$\beta _{eff}$ increases negligibly and remains approximately constant. The range of bottom temperatures for which 
$\beta _{eff}$ can be considered constant shrinks to 0 for subglacial lakes with 
$p_i\leq p_*$ since 
${T}_d\geq T_f$ in this case, and it increases as 
$p_i>p_*$ increases. The next section (§ 3.3) presents a derivation of two sets of closed-form expressions for several variables of interest, including the bottom temperature, in terms of the problem parameters, based on whether we assume that 
$\beta _{eff}$ is constant or is linearly proportional to 
$\varDelta _{eff}$. Then, § 4 discusses which predictions are applicable to subglacial lakes in Antarctica.
3.3. Predictions in the limit of decoupled convective and stably stratified dynamics
In this section we set out to define an effective (output) Rayleigh number 
$Ra_{eff}$ based on the simulation results and explore its dependence on the control parameters. We then demonstrate that the Nusselt number, 
$Nu$, which estimates the contribution of fluid motions to the transport of heat relative to conduction alone, and the Reynolds number, 
$Re$, which compares fluid inertia to viscous dissipation, display asymptotic behaviours with 
$Ra_{eff}$ similar to those observed in classical Rayleigh–Bénard convection. This allows for the derivation of predictive expressions for all output variables of interest in terms of the control parameters.
We assume that the convective and stably stratified layers are dynamically decoupled, such that we can simply define the effective Rayleigh number as
\begin{equation} Ra_{eff} = \frac{g\varDelta_{eff}\beta_{eff} h_{eff}^3}{\nu\kappa}, \end{equation}
where 
$h_{eff}$ is the effective (decoupled) convective layer depth and 
$\varDelta _{eff}$ and 
$\beta _{eff}$ are the effective temperature difference and thermal expansion coefficient discussed in § 3.2 and whose dimensionless forms are given by (3.1) and (3.2). Neglecting the influence from the stable layer on the convection means that the effective convective layer depth is well approximated by the full depth minus the mean thickness of the top stably stratified layer, which is equal to 
$h(T_d-T_f)/\varDelta$ if 
$T_d>T_f$ and 0 otherwise. Thus,
\begin{equation} h_{eff} = h\left[ 1 -\frac{T_d-T_f}{\varDelta}( T_d>T_f )\right], \end{equation}which can be rewritten as
\begin{equation} h_{eff} = h[ 1 -\tilde{T}_d(\tilde{T}_d>0) ] = \frac{h}{1+\bar{S}(\bar{S}>0)}, \end{equation}
using either the dimensionless temperature of maximum density or the base-state stiffness parameter (cf. (2.7a–c)). Using (3.5) we can rewrite 
$Ra_{eff}$, i.e. (3.3), as
\begin{equation} Ra_{eff} = \overline{Ra}_F\frac{\tilde{\varDelta}_{eff}\tilde{\beta}_{eff}}{[ 1+\bar{S}(\bar{S}>0)]^3}, \end{equation}
where 
$\overline {Ra}_F$ is the base-state flux-based Rayleigh number (cf. (2.7a–c)). Equations (3.3) and (3.6) can be expected to represent accurately the effective Rayleigh number when there is (almost) no contribution from the stable layer to the convective dynamics, but to be inaccurate when the stable layer is entrained into and modifies the properties of the lower convective bulk. Couston et al. (Reference Couston, Lecoanet, Favier and Le Bars2017) investigated the dynamics of mixed convective and stably stratified fluids over a broad range of input stiffnesses and showed that the Nusselt number was indeed enhanced due to an entrainment heat flux in the limit of small (effective) stiffness, suggesting that the effective Rayleigh number should be adjusted when entrainment from the stable layer into the convective bulk is significant. We will show in § 3.4 that the entrainment heat flux is small in all our simulations, such that the assumption of decoupled convective and stably stratified dynamics, leading to (3.3) and (3.6) for 
$Ra_{eff}$, is accurate at leading order. Note that the base-state stiffness 
$\bar {S}$, which can be related to the effective convective layer depth (cf. (3.5)), is not the effective stiffness of the system due to the variability of the bottom temperature and thermal expansion coefficient. The effective stiffness of the simulations is estimated a posteriori and discussed in § 3.4.
Figures 5(a)–5(d) show the effective Rayleigh number 
$Ra_{eff}$ as a function of the control parameters 
$\overline {Ra}_F$ and 
$\bar {S}$ or normalized heat flux 
$(F-F_t)/(F_c-F_t)$ and water depth 
$(h-h_t)/(h_c-h_t)$. Figures 5(a) and 5(b) show 
$Ra_{eff}$ for all simulations, i.e. combining the results of the first and second experiments, in addition to the effective Rayleigh number for a series of Oberbeck–Boussinesq (OB) simulations, which we denote by 
$\mathcal {S}_{OB}$ and which have different 
$\overline {Ra}_{F}$ and fixed 
$Pr=12.8$ (shown by the black crosses). The definition of the effective Rayleigh number for 
$\mathcal {S}_{OB}$ is simply 
$Ra_{eff}=\tilde {\varDelta }_{eff}\overline {Ra}_F$ since 
$\bar {S}=0$ and all physical variables, including the thermal expansion coefficient, are assumed constants in the OB approximation. It can be seen in figure 5(a) that 
$Ra_{eff}$ in simulations 
$\mathcal {S}_3^2$ (red triangles) and 
$\mathcal {S}_3^1$ (red circles; for relatively low values of 
$\overline {Ra}_F$) follows the same trend as 
$Ra_{eff}$ in simulation 
$\mathcal {S}_{OB}$ (black crosses), whose asymptotic behaviour is shown by the black solid line. This suggests that subglacial lakes that are fully convective, and for which the effective thermal expansion coefficient is approximately independent of 
$\overline {Ra}_{F}$, i.e. constant, behave similarly to OB fluids. For simulations 
$\mathcal {S}_i^j$ (
$i=0,1,2$; 
$j=1,2$), 
$Ra_{eff}$ also increases with 
$\overline {Ra}_{F}$ but following a trend (whose asymptotic behaviour is shown by the green solid line) that is markedly different from the OB results, which is due to the variability of the effective thermal expansion coefficient with the control parameters when 
$\bar {S}\geq 0$ (or 
$p_i\leq p_*$). Figure 5(b) shows 
$Ra_{eff}$ as a function of 
$(1+\bar {S})$ instead of 
$\bar {S}$ because the base-state stiffness enters the definition of 
$Ra_{eff}$ as 
$1/(1+\bar {S})$ (cf. (3.5)). It should be noted that 
$\overline {Ra}_{F}$ and 
$\bar {S}$ are not independent variables in our two numerical experiments, since we decided to explore the dynamics of subglacial lakes by sweeping along lines of constant heat flux and water depth rather than lines of constant 
$\overline {Ra}_{F}$ and 
$\bar {S}$. Thus, the increase of 
$Ra_{eff}$ for simulation 
$\mathcal {S}_3^1$ with 
$\bar {S}$ (small red circles in figure 5b) is due to the associated increase of 
$\overline {Ra}_{F}$, not 
$\bar {S}$. Similarly, the decrease of 
$Ra_{eff}$ with increasing 
$\bar {S}$ for simulations 
$\mathcal {S}_i^j$ (
$i=0,1$; 
$j=1,2$) is due to the associated decrease of 
$\overline {Ra}_{F}$ (cf. figure 2). The existence of two distinct scaling behaviours for 
$Ra_{eff}$ with the problem parameters is again clear in figure 5(c) where 
$Ra_{eff}$ follows a common asymptotic behaviour for simulations 
$\mathcal {S}_i^1$ (
$i=0,1,2$), shown by the parallel blue, orange and green (relatively thick) solid lines, and one other for simulation 
$\mathcal {S}_3^1$, shown by the thin red solid line. The difference in asymptotic behaviours for 
$Ra_{eff}$ between simulations 
$\mathcal {S}_i^2$ (
$i=0,1,2$) and simulation 
$\mathcal {S}_3^2$ in figure 5(d) is tenuous. This is because the difference in scalings for 
$Ra_{eff}$ with the water depth is much weaker than with the heat flux, as we demonstrate at the end of this section.
Figure 5. (a) Effective Rayleigh number 
$Ra_{eff}$ as a function of 
$\overline {Ra}_F$ for all simulations, including the Oberbeck–Boussinesq (OB) simulations 
$\mathcal {S}_{OB}$ whose results are shown by black crosses. Note that, as in figure 4, large (respectively small) symbols highlight simulation results with small (respectively large) 
$p_i$. The thick black and thin green solid lines show the predictive expressions for 
$Ra_{eff}$ as a function of 
$\overline {Ra}_F$ for simulations 
$\mathcal {S}_{OB}$ (OB regime) and 
$\mathcal {S}_{2}^1$ (linear expansion coefficient, or LEC, regime), respectively (cf. table 4); (b) 
$Ra_{eff}$ as a function of 
$1+\bar {S}$. The vertical dotted line highlights 
$\bar {S}=0$; (c) 
$Ra_{eff}$ for the first experiment only as a function of the normalized heat flux 
$(F-F_t)/(F_c-F_t)$; (d) 
$Ra_{eff}$ for the second experiment only as a function of the normalized water depth 
$(h-h_t)/(h_c-h_t)$. The solid lines in (c) and (d) show the theoretical scalings for 
$Ra_{eff}$ with 
$F$ and 
$h$ (cf. table 4) in the OB regime (thin red lines) and in the LEC regime (all other lines) with 
$\alpha$ and 
$\gamma$ listed in table 3.
Now that we have explored the dependence of the effective Rayleigh number on the problem parameters, we turn our attention to the Nusselt number 
$Nu$ and the Reynolds number 
$Re$. We define the Nusselt number as the ratio of the full heat flux 
$F$ divided by the conductive heat flux based on the output temperature difference between the top and bottom boundaries of the convective bulk at statistical steady state, i.e.
\begin{equation} Nu = \frac{F}{k\dfrac{ \varDelta_{eff}}{ h_{eff}}} = \frac{\tilde{h}_{eff}}{\tilde{\varDelta}_{eff}} = \frac{1}{\tilde{\varDelta}_{eff}[ 1+\bar{S}(\bar{S}>0)]}, \end{equation}
with 
$\tilde {h}_{eff}=h_{eff}/h$ and the Reynolds number as
\begin{equation} Re = \frac{V_{rms}h_{eff}}{\nu}, \end{equation}
with 
$V_{rms}=\sqrt {h_{eff}^{-1}\int _0^{h_{eff}} \langle |\boldsymbol {u}|^2 \rangle \, \textrm {d} z}$ the root-mean-square (r.m.s.) velocity within the convective layer.
We show the Nusselt number 
$Nu$ as a function of 
$Ra_{eff}$ in figure 6(a) for all simulations of the first experiment (with results shown by circles), all simulations of the second experiment (shown by triangles) and for simulation 
$\mathcal {S}_{OB}$ (black crosses). There is a clear universal asymptotic scaling of 
$Nu$ with 
$Ra_{eff}$ in all cases, which is highlighted by the black solid line that shows the best-fit power law for the OB results. More precisely, the Nusselt number in all simulations can be expressed to a good approximation as a power law of the form
\begin{equation} Nu = a Ra_{eff}^{\alpha}, \end{equation}
where pre-factor 
$a$ and exponent 
$\alpha$ are reported in table 3 for each simulation along with the relative error, which is typically of the order of 1 % and always less than 8 %. The exponent 
$0.26<\alpha <0.29$ shows little variability across the simulations, as can be seen from the similar slopes of the simulation results for large 
$Ra_{eff}$ (cf. figure 6a). The pre-factor is slightly more variable, i.e. 
$0.16 < a < 0.31$, and is relatively large for simulations 
$\mathcal {S}_{i}^j (i=0,1; j=1,2)$, i.e. with a stable layer, as can be seen from the upward shift of the large blue and orange circles and triangles relative to the other smaller symbols and the solid black line. The discrepancy of 
$Nu$ between simulations that are fully convective and simulations with a top stably stratified layer is due to the fact that the effective temperature difference 
$\tilde {\varDelta }_{eff}$ as defined by (3.1) underestimates the range of temperatures involved in the convective heat transport, which is extended in the presence of a stable layer. Substituting 
$\tilde {\varDelta }_{eff}$ with 
$3/2\tilde {\varDelta }_{eff}$ in the definitions of 
$Ra_{eff}$ and 
$Nu$ helps correct the discrepancy and yields the results shown by the light-blue and light-orange markers in figure 6(a), which better overlap with the other simulation results (cf. details in § 3.4). Similar to the Nusselt number, the Reynolds number also follows an almost universal scaling with 
$Ra_{eff}$ (see figure 6b), such that it can be predicted to a good approximation using a power law of the form
\begin{equation} Re = c Ra_{eff}^{\gamma}, \end{equation}
with almost the same pre-factor 
$c$ and exponent 
$\gamma$ across all simulations (cf. table 3). Note that the standard deviation 
$\sigma _{Re}$ (in time) of the Reynolds number (respectively 
$\sigma _{Nu}$ for the Nusselt number) over the (second) half of the diffusive time scale of each simulation stage is always in the range 
$10^{-2}Re<\sigma _{Re}<10^{-1}Re$ (respectively 
$10^{-3}Nu<\sigma _{Nu}<3\times 10^{-2}Nu$), i.e. small, except for simulation 
$\mathcal {S}_3^1$ with the smallest heat flux, for which 
$\sigma _{Re}\approx 1.4Re$, because 
$Re$ increases slowly toward its final value 
$Re\approx 0.5$ in this case.
Figure 6. (a) Nusselt number 
$Nu$ as a function of 
$Ra_{eff}$ for all simulations (we use the same symbol colours and size chart as in figures 4 and 5). (b) Reynolds number 
$Re$ as a function of 
$Ra_{eff}$. The black solid lines in (a,b) show the power law fit (3.9) and (3.10) for 
$Nu$ and 
$Re$ for the results of simulation 
$\mathcal {S}_{OB}$ (cf. table 3). (c) Compensated Nusselt number as a function of 
$\overline {Ra}_F$. (d) Compensated Reynolds number as a function of 
$\overline {Ra}_F$. The black (respectively green) solid lines in (c,d) show the power law fits provided in table 4 for the results of simulation 
$\mathcal {S}_{OB}$ (respectively 
$\mathcal {S}_2^1$), i.e. in the OB regime (respectively LEC regime). The light-coloured markers in (a,b,c) show the results for simulations 
$\mathcal {S}_{i}^{j}$ (
$i=0,1$; 
$j=1,2$) with 
$Nu$ multiplied by 2/3 and 
$Ra_{eff}$ multiplied by 3/2 (see details in § 3.4).
Table 3. Best-fit coefficients for the power laws 
$Nu=a Ra_{eff}^{\alpha }$ and 
$Re=cRa_{eff}^{\gamma }$ for 
$Ra_{eff}>10^5$ for all simulations, including the OB simulation 
$\mathcal {S}_{OB}$. Here, s.n. means simulation name while Rel. err. denotes the maximum relative error between the simulation results and the predictive best-fit power law.
With the predictive power law (3.9) for 
$Nu$ in terms of 
$Ra_{eff}$ in hand, we can derive a predictive equation for the effective temperature difference 
$\tilde {\varDelta }_{eff}$ in terms of the problem parameters. Substituting the expressions (3.6) and (3.7) for 
$Ra_{eff}$ and 
$Nu$ into (3.9), we obtain
\begin{equation} \frac{1}{\tilde{\varDelta}_{eff}[ 1+\bar{S}(\bar{S}>0)] } = a \left\{ \overline{Ra}_F\frac{\tilde{\varDelta}_{eff}\tilde{\beta}_{eff}}{[ 1+\bar{S}(\bar{S}>0)]^3} \right\}^{\alpha}, \end{equation}
which is an equation between 
$\tilde {\varDelta }_{eff}$, 
$\tilde {\beta }_{eff}$ and the control parameters. The effective thermal expansion coefficient (3.2) can be expressed in terms of 
$\tilde {\varDelta }_{eff}$ and 
$\bar {S}$ as
\begin{equation} \tilde{\beta}_{eff} = \frac{\tilde{\varDelta}_{eff}}{2}( 1+\bar{S} ) - \bar{S}( \tilde{T}_d < 0 ), \end{equation}
with the second term on the right-hand side being the non-zero dimensionless thermal expansion coefficient at the top of subglacial lakes that are fully convective, i.e. for which 
$\tilde {T}_d<0$ or 
$-1\leq \bar {S}< 0$. Substituting (3.12) into (3.11) yields an algebraic equation for 
$\tilde {\varDelta }_{eff}$ in terms of the control parameters 
$\overline {Ra}_F$ and 
$\bar {S}$, viz.
\begin{equation} \tilde{\varDelta}_{eff}^{1+\alpha} \left[ \frac{\tilde{\varDelta}_{eff}}{2}( 1+\bar{S} ) - \bar{S}( \tilde{T}_d < 0 ) \right]^{\alpha} = \frac{[ 1+\bar{S}(\bar{S}>0)]^{3\alpha-1}}{a \overline{Ra}_F^{\alpha}}. \end{equation}
Equation (3.13) is valid as long as the power law (3.9) provides a good approximation to the simulation results but is nonlinear in 
$\tilde {\varDelta }_{eff}$. Based on the results of § 3.2, we know that the thermal expansion coefficient in subglacial lakes exhibits two limiting behaviours, i.e. 
$\beta _{eff}$ is approximately constant, which we refer to as the OB regime, or 
$\beta _{eff}$ is linearly proportional to the effective temperature difference, which we refer to as the LEC regime. (3.13) can be simplified in both regimes. In the OB regime, 
$\beta _{eff}=\bar {\beta }_b$, i.e. 
$\tilde {\beta }_{eff}=1$, and 
$-1\leq \bar {S}<0$ (or 
$\tilde {T}_d<0$), such that (3.13) reduces to
\begin{equation} \tilde{\varDelta}_{eff}^{1+\alpha} = \frac{1}{a\overline{Ra}_F^{\alpha}}. \end{equation}
In the LEC regime, which is obtained when 
$\bar {S}>0$ (
$\tilde {T}_d>0$) or 
$|\bar {S}|\ll \tilde {\varDelta }_{eff}\leq 1$, the second term on the right-hand side of (3.12) is negligible and (3.13) can be approximated as
\begin{equation} \tilde{\varDelta}_{eff}^{1+2\alpha} = \frac{2^{\alpha}( 1+\bar{S} )^{2\alpha-1}}{a\overline{Ra}_F^{\alpha}}. \end{equation}
Substituting the expressions (2.7a–c) for 
$\overline {Ra}_F$ and 
$\bar {S}$ in terms of the control parameters into (3.14) yields an asymptotic scaling in the OB regime, which assumes constant 
$\bar {\beta }_b$ (such that 
$\overline {Ra}_F\sim Fh^4$), for the dimensional effective temperature difference as
\begin{equation} \varDelta_{eff} = \varDelta \tilde{\varDelta}_{eff} \sim F^{{1}/({1+\alpha})} h^{({1-3\alpha})/({1+\alpha})}, \end{equation}
which we recall is related to the mean bottom temperature through equation (3.1). Substituting (2.7a–c) into (3.15) and using the approximation 
$\bar {S}\approx 1$ in the limit 
$(Fh)/k\gg T_d$ (large thermal driving) when 
$\bar {S}>0$, yields instead
\begin{equation} \varDelta_{eff} = \varDelta \tilde{\varDelta}_{eff} \sim F^{{1}/({2\alpha+1})}h^{({1-3\alpha})/({2\alpha+1})}, \end{equation}
in the LEC regime (for which 
$\overline {Ra}_F\sim F^2h^5$). From the closed-form expressions for 
$\tilde {\varDelta }_{eff}$ (3.14) and (3.15) it is possible to derive the expressions for all variables of interest, including 
$Ra_{eff}$, 
$Nu$ and 
$Re$ but also dimensional variables such as 
$\beta _{eff}$ and 
$V_{rms}$, in terms of the problem parameters. We summarize all key expressions for the variables of interest in terms of 
$\overline {Ra}_F$ and 
$\bar {S}$ as well as the asymptotic scaling laws with the water depth and heat flux in table 4. Figure 6(c–d) shows 
$Nu$ and 
$Re$ in terms of 
$\overline {Ra}_F$ compensated such that the influence of the base-state stiffness predicted in table 4 is removed.
Table 4. Key expressions in terms of the control parameters 
$\overline {Ra}_F$ and 
$\bar {S}$ (cf. (2.7a–c)) and asymptotic scaling laws in terms of the water depth 
$h$ and heat flux 
$F$ for the dimensionless and dimensional variables of interest to the study of turbulent convection in subglacial lakes. The starting point is the algebraic equation for 
$\tilde {\varDelta }_{eff}$ (3.13). The difference between the limiting OB regime and LEC regime arises from the independence or variability of the effective thermal expansion coefficient 
$\tilde {\beta }_{eff}$ with the effective temperature difference 
$\tilde {\varDelta }_{eff}$. The asymptotic scalings with 
$F$ and 
$h$ are valid in the limit 
$(Fh/k)/|T_d-T_f| \ll 1$ (
$\bar {S}\approx -1$) for the OB regime and 
$(Fh/k)/|T_d-T_f| \gg 1$ (
$\bar {S}\approx 0$) for the LEC regime.
The asymptotic scalings with 
$F$ and 
$h$ are valid in the limit 
$(Fh/k)/|T_d-T_f| \ll 1$, i.e. such that 
$\bar {S}\approx -1$ for the OB regime and 
$(Fh/k)/|T_d-T_f| \gg 1$, i.e. such that 
$\bar {S}\approx 0$, for the LEC regime. The latter is satisfied when 
$F \rightarrow \infty$ and 
$h \rightarrow \infty$, i.e. for increasingly large heat flux and water depth. The former limit for the OB regime is only valid for a finite range of heat fluxes and water depth as the thermal expansion coefficient will eventually increase substantially as 
$F$ and 
$h$ keep increasing. The results of simulation 
$\mathcal {S}_3^1$ show how 
$\tilde {\beta }_{eff}$ in figure 4(c) and 
$Re$ in figure 6(d) transitions from varying according to the OB regime for relatively small 
$F$ but according to the LEC regime for relatively large 
$F$.
The validity of the scaling laws provided in table 4 can be verified from the good overlap between the simulation results and the solid lines, which show the asymptotic scalings with 
$F$ and 
$h$, in figures 4, 5(c) and 5(d). The relatively thick solid lines that are blue, orange and green, display the asymptotic scalings in the LEC regime, whereas thin red and black solid lines display asymptotic scalings in the OB regime. The simulation results have provided not only the scaling exponents 
$\alpha$ and 
$\gamma$, but also the pre-factors 
$a$ and 
$c$ for the predictive power laws (3.9) and (3.10). Thus, we can predict the actual value of any variable of interest using the equations listed in table 4 rather than just the scalings with the problem parameters. The validity of the full expressions derived for 
$Ra_{eff}$, 
$Nu$ and 
$Re$ can be verified from the good overlap between the simulation results and the black and green solid lines, which show the power law fits, in figure 5(a) and figure 6(a–d).
Here, we find 
$\alpha \approx 2/7$, which was first proposed as a natural scaling exponent for 
$Nu$ based on phenomenological arguments (Castaing et al. Reference Castaing, Gunaratne, Heslot, Kadanoff, Libchaber, Thomae, Wu, Zaleski and Zanetti1989), and 
$\gamma \approx 3/5$ for all simulations (cf. table 3). Substituting 
$\alpha = 2/7$ and 
$\gamma = 3/5$ in the asymptotic expressions for the dimensional variables with 
$F$ and 
$h$ in table 4 yields in the OB regime
\begin{equation} h_{eff} \sim h, \quad \varDelta_{eff} \sim F^{{7}/{9}}h^{{1}/{9}}, \quad \beta_{eff} \sim 1, \quad V_{rms} \sim F^{{7}/{15}}h^{{13}/{15}}, \end{equation}and in the LEC regime
\begin{equation} h_{eff} \sim h, \quad \varDelta_{eff} \sim F^{{7}/{11}}h^{{1}/{11}}, \quad \beta_{eff} \sim F^{{7}/{11}}h^{{1}/{11}}, \quad V_{rms} \sim F^{{42}/{55}}h^{{10}/{11}}. \end{equation}
It can be seen that the difference in scalings is typically greater with 
$F$ than with 
$h$, i.e. for instance, the difference in scaling exponent for 
$\varDelta _{eff}$ with 
$F$ is 
$7/9-7/11\approx 0.14$ whereas it is only 
$1/9-1/11\approx 0.02$ with 
$h$. The discrepancy in the difference of scaling exponents is the reason the two asymptotic regimes are usually easier to identify in simulation results from the first experiment than from the second experiment (compare, e.g. figures 4(a,b) and 5(c,d)). We find that 
$V_{rms}$ has a scaling with 
$F$ and 
$h$ that is steeper for subglacial lakes in the LEC regime than in the OB regime. The steeper scaling for 
$V_{rms}$ in the LEC regime occurs because the effective Rayleigh number 
$Ra_{eff}$ increases more rapidly with 
$F$ and 
$h$ in the LEC regime, due to the combined increase of the effective thermal driving and thermal expansion coefficient, than in the OB regime for which the thermal expansion coefficient remains constant.
3.4. Effective stiffness and entrainment dynamics
The power laws (3.9) and (3.10) for 
$Nu$ and 
$Re$ with 
$Ra_{eff}$ are in excellent agreement with the numerical results using pre-factors and exponents that are almost the same across all simulations, including the OB simulation (cf. table 3). This suggests that our definition of 
$Ra_{eff}$, i.e. (3.6), which assumes decoupled convective and stably stratified dynamics, is a good proxy for the true effective Rayleigh number of all laboratory-scale subglacial lakes, including lakes that have a top stable layer and can experience penetrative convection.
In this section, we investigate in details the influence of the stable fluid layer on the convection in order to predict whether the assumption of decoupled convective and stably stratified dynamics can hold at significantly higher heat flux and water depth. We first estimate the effective stiffness 
$S_{eff}$ of the convective–stably stratified interface in the simulations as
\begin{equation} S_{eff} = \frac{T_d-T_f}{\langle T_b \rangle-T_d}, \end{equation}
which is the ratio of the temperature difference across the stable layer to the temperature difference across the convective layer. In a mixed convective and stably stratified subglacial lake, i.e. for which 
$p_i\leq p_*$, this ratio is equal to the ratio of the opposite of the thermal expansion coefficient at the top of the stable layer (minimum negative thermal expansion coefficient) to the thermal expansion coefficient at the bottom of the convective layer (maximum positive thermal expansion coefficient). Thus, the effective stiffness 
$S_{eff}$ in (3.20) is equivalent to the 
$\varLambda$ parameter in Toppaladoddi & Wettlaufer (Reference Toppaladoddi and Wettlaufer2018) and is similar to the input stiffness parameter 
$S_i$ in Couston et al. (Reference Couston, Lecoanet, Favier and Le Bars2017), although we recall that Couston et al. (Reference Couston, Lecoanet, Favier and Le Bars2017) used a thermal expansion coefficient that is piecewise constant rather than linearly varying with temperature. We show 
$S_{eff}$ as a function of 
$Ra_{eff}$ in figure 7(a) for all simulations with a stably stratified layer. The effective stiffness decreases monotonically with 
$Ra_{eff}$, which is expected since the numerator in (3.20) is constant while the denominator, which is the dimensional thermal driving 
$\varDelta _{eff}$, increases with 
$Ra_{eff}$ (cf. figure 4e,f). The decrease of 
$S_{eff}$ with 
$Ra_{eff}$ suggests that the stably stratified layer may modify the properties of the convective bulk significantly for large enough 
$F$ and 
$h$. Couston et al. (Reference Couston, Lecoanet, Favier and Le Bars2017) demonstrated that the two layers become increasingly coupled as the input stiffness decreases, indeed, with, for instance, the mean temperature of the well-mixed bulk dropping below the temperature of maximum density in their experiment with 
$Ra=8\times 10^7$ and 
$S_i=1$.
Figure 7. (a) Effective stiffness 
$S_{eff}$ (cf. (3.20)) as a function of 
$Ra_{eff}$ for simulations 
$\mathcal {S}_i^j$ (
$i=0,1$; 
$j=1,2$), i.e. with a top stably stratified layer. (b) Ratio between the bulk temperature 
$\tilde {T}_{bulk}$ (cf. (3.21)) in excess of 
$\tilde {T}_d(\tilde {T}_d>0)$ and the mean bottom temperature 
$\langle \tilde {T}_b \rangle$ in excess of 
$\tilde {T}_d(\tilde {T}_d>0)$ as a function of 
$Ra_{eff}$. The results are shown for all simulations of the first and second experiments and the dashed line highlights the 1/4 ordinate. (c) Effective stiffness 
$S_{eff}$ as a function of the expected dimensionless convective layer depth 
$\tilde {h}_{eff}$ (cf. (3.5)) for all simulations with a top stably stratified layer.
In order to understand whether the decrease of 
$S_{eff}$ with 
$Ra_{eff}$ implies that the stable layer of subglacial lakes modifies the convective bulk dynamics for large (geophysical) 
$F$ and 
$h$, as may be expected based on the results of Couston et al. (Reference Couston, Lecoanet, Favier and Le Bars2017), we first investigate the evolution of the temperature of the well-mixed bulk, which we define as
\begin{equation} \tilde{T}_{bulk} = \frac{1}{\tilde{h}_{eff}}\int_0^{\tilde{h}_{eff}} \tilde{T} \, \textrm{d} \tilde{z}, \end{equation}
with 
$Ra_{eff}$. If there is no stable layer (
$\tilde {T}_d<0$) or if the stable layer does not influence the convective dynamics, we expect 
$[\tilde {T}_{bulk}-\tilde {T}_d(\tilde {T}_d>0)]\approx [\langle \tilde {T}_b\rangle -\tilde {T}_d(\tilde {T}_d>0)]/2$, i.e. the bulk temperature is the average of the temperatures at the top and bottom of the convection zone. Figure 7(b) shows that the bulk temperature is the average (approximately) of the top and bottom temperatures for simulations without a stable layer (small green and red symbols). For simulations with a stable layer, however, we find that the bulk temperature minus 
$\widetilde{T}_d(\widetilde{T}_d>0)$ of the two-layer system is typically lower than 
$[\langle \tilde {T}_b\rangle -\tilde {T}_d(\tilde {T}_d>0)]/2$. Specifically, we find that it decreases initially with 
$Ra_{eff}$ and then reaches the plateau 
$[\tilde {T}_{bulk}-\tilde {T}_d(\tilde {T}_d>0)]\approx [\langle \tilde {T}_b\rangle -\tilde {T}_d(\tilde {T}_d>0)]/4$ (black dashed line). For simulation 
$\mathcal {S}_1^1$, the normalized bulk temperature appears to anomalously increase with 
$Ra_{eff}$ at large 
$Ra_{eff}$ (two rightmost orange circles). The increase of the normalized bulk temperature toward the 0.5 mark at large 
$Ra_{eff}$ for 
$\mathcal {S}_1^1$ is due to the shrinking of the top stable layer. In fact, 
$Ra_{eff}$ increases because 
$F$ increases in the simulations of the first experiment, which forces a thinning of the stratified top layer. As the stably stratified layer shrinks, possibly to a point where it is thinner than the thermal boundary layer, upward plumes attempting to penetrate into the stable layer rapidly feels the effect of the top wall and lose their inertia, which reduces entrainment of the stable fluid and lowering of the bulk temperature. The vanishing of the stable layer thickness can be seen in figure 7(c), which shows that the convective layer depth 
$\tilde {h}_{eff}\rightarrow 1$ as 
$S_{eff}$ decreases, or, equivalently, as 
$Ra_{eff}$ increases. The observation of a plateau for 
$[\tilde {T}_{bulk}-\tilde {T}_d(\tilde {T}_d>0)]/[\langle \tilde {T}_b\rangle -\tilde {T}_d(\tilde {T}_d>0)]$ at 0.25 with 
$Ra_{eff}$ (let alone possible increase toward the 0.5 mark when the stable layer becomes too thin) suggests that there is an influence of the stable layer on the convection, but that this influence does not increase with increasing 
$Ra_{eff}$ or decreasing 
$S_{eff}$. We further demonstrate in the next paragraph that the entrainment of stable fluid into the convection zone does not increase with 
$Ra_{eff}$, such that the stable layer's influence on the convection is indeed limited and negligible at leading order.
The lowering of the bulk temperature is caused by turbulent plumes that occasionally penetrate into the stable layer and entrain some of the top cold fluid into the convection zone. The significance of the entrained stable fluid to the system's overall dynamics can be estimated from its contribution to the convective heat flux (Couston et al. Reference Couston, Lecoanet, Favier and Le Bars2017), as the downward motion of cold fluid can contribute a positive heat flux similar to upward-moving warm plumes. We define the depth-dependent (dimensionless) heat flux due to convective plumes, which we denote by 
$q_{cp}$, and the depth-dependent heat flux due to the entrained stable fluid, which we denote by 
$q_{ce}$, as
\begin{gather} q_{cp} = \langle \tilde{w}( \tilde{T}-\tilde{T}_d ) ( \tilde{T}>\tilde{T}_d ) \rangle, \end{gather}
\begin{gather}q_{ce} = \langle \tilde{w}( \tilde{T}-\tilde{T}_d ) ( \tilde{T}\leq\tilde{T}_d ) \rangle, \end{gather}
respectively. We show in figure 8(a–d) the entrained heat flux 
$q_{ce}$ (solid lines), the plume-driven convective heat flux 
$q_{cp}$ (dashed lines) and the conductive heat flux 
$q_d=-\langle \partial _{\tilde {z}} \tilde {T} \rangle$ (dotted lines) as functions of the normalized water depth 
$z/h$ for all simulation stages including a stable layer. The plume-driven convective heat flux is significant over a depth that becomes larger with each simulation stage (later stages are shown with lighter colours), due to the increase of either 
$F$ (figure 8a,b) or 
$h$ (figure 8c,d). The top of the convective layer is approximately where the entrained heat flux is maximum and where the conductive heat flux increases rapidly with 
$z/h$. Note that 
$q_{ce}+q_{cp}+q_d\approx 1$ for all 
$z/h$ as the total heat flux is depth invariant at statistical steady state (the discrepancy after temporal and horizontal averaging is smaller than a few per cent). We remark that the conductive heat flux 
$q_d$ can be negative in the lower half of the convective layer (see, e.g. dotted lines to the left of the 0 abscissa below 
$z/h\approx 0.2$ in figure 8a,b) and that the entrainment heat flux 
$q_{ce}$ can be slightly negative at the top of the convection zone (see, e.g. solid lines going to the left of the 0 abscissa at 
$z/h\approx 0.2, 0.55, 0.7, 0.8$ in figure 8c). The reversals of the conductive heat flux in the convective region and of the entrainment heat flux at the bottom of the stable layer are real physical phenomena (i.e. they are neither a numerical nor a statistical artefact) and we note that they have been observed in past laboratory experiments (Townsend Reference Townsend1964; Adrian Reference Adrian1975).
Figure 8. (a–d) Convective plume-driven 
$q_{cp}$ (dashed lines), entrainment-driven 
$q_{ce}$ (solid lines) and conductive 
$q_d$ (thin dotted lines) heat fluxes as functions of water depth 
$z/h$ (
$y$ axis) for simulations 
$\mathcal {S}_0^1$, 
$\mathcal {S}_1^1$, 
$\mathcal {S}_0^2$ and 
$\mathcal {S}_1^2$, respectively. (e–h) Top of the bottom conductive layer 
$\tilde {z}=L_{db}$ (big bottom circles), plume-driven convective layer 
$\tilde {z}=L_{db}+L_{cp}$ (lower horizontal bars), entrainment-driven convective layer 
$\tilde {z}=L_{db}+L_{cp}+L_{ce}$ (upper horizontal bars) and stably stratified layer 
$\tilde {z}=L_{db}+L_{cp}+L_{ce}+L_{dt}$ (top squares) as functions of 
$Ra_{eff}$ for simulations 
$\mathcal {S}_0^1$, 
$\mathcal {S}_1^1$, 
$\mathcal {S}_0^2$ and 
$\mathcal {S}_1^2$, respectively. The crosses show the height 
$\tilde {z}_d$ of the mean 
$\tilde {T}_d$ isotherm. The lines’ colour and symbol colour become lighter as the heat flux or water depth increases.
We estimate the contribution of each heat flux component (i.e. 
$q_{cp}$, 
$q_{ce}$ and 
$q_d$) to the full heat transport by computing the ratio of each heat flux component (integrated over depth) to the volume-averaged heat flux, which is unity in dimensionless space. We separate the contribution of the diffusive heat flux below and above the convection zone by integrating 
$q_d$ either over 
$[0,\tilde {h}_{eff}]$ or over 
$[\tilde {h}_{eff},1]$. As in Couston et al. (Reference Couston, Lecoanet, Favier and Le Bars2017), we interpret the ratio of each depth-integrated heat flux component to the mean heat flux as an equivalent (dimensionless) layer thickness (wherein the full heat flux is assumed transported by that heat flux component), in order to visualize the contribution of each heat flux component to the full heat transport. Specifically, we denote by
\begin{align} L_{db}= \int_0^{\tilde{h}_{eff}}q_d \, \textrm{d} z, \quad L_{cp} = \int_0^{1}q_{cp} \, \textrm{d} z, \quad L_{ce} = \int_0^{1}q_{ce} \, \textrm{d} z, \quad L_{dt} = \int_{\tilde{h}_{eff}}^{1}q_d \, \textrm{d} z, \end{align}
the equivalent layer thicknesses of the diffusive bottom boundary layer, of the plume-driven convective layer, of the entrained layer and of the diffusive top stable layer, respectively, which we envision as stacked on top of each other, with 
$L_{db}+L_{cp}+L_{ce}+L_{dt}\approx 1$. Figure 8(e–h) shows the top of each layer, i.e. 
$\tilde {z}=L_{db}$ (bottom circles), 
$\tilde {z}=L_{db}+L_{cp}$ (lower horizontal bars), 
$\tilde {z}=L_{db}+L_{cp}+L_{ce}$ (upper horizontal bars) and 
$\tilde {z}=L_{db}+L_{cp}+L_{ce}+L_{dt}$ (top squares) as functions of 
$Ra_{eff}$. We also display the height 
$\tilde {z}_d$ of the mean 
$\tilde {T}_d$ isotherm for each simulation (crosses), which is almost always between the top of the plume-driven layer and the top of the entrained layer. Importantly, the thickness of the entrained layer, which is given by the vertical distance between the two horizontal bars showing 
$L_{db}+L_{cp}$ and 
$L_{db}+L_{cp}+L_{ce}$, does not vary significantly with 
$Ra_{eff}$. This can be clearly seen in figure 9(a), which shows 
$L_{ce}$ as a function of 
$Ra_{eff}$. 
$L_{ce}$ is typically less then 0.05 and does not noticeably increase with 
$Ra_{eff}$. While the thickness of the entrained layer remains approximately constant, the thickness of the plume-driven convective layer thickness increases with 
$Ra_{eff}$, such that the entrainment parameter, 
$\mathcal {E}$, which we define as
\begin{equation} \mathcal{E} = \frac{q_{ce}}{q_{cp}+q_{ce}}, \end{equation}
and which compares the entrainment heat flux to the full convective heat flux, is (almost) monotonically decreasing with 
$Ra_{eff}$ (figure 9b). Figure 9(c) shows the height 
$\tilde {z}_d$ of the 
$\tilde {T}_d$ isotherm, which may be considered as an output mean convective–stably stratified interface height (similar to the top of the entrained layer 
$\tilde {z}=L_{db}+L_{cp}+L_{ce}$), as a function of the expected convective layer depth 
$\tilde {h}_{eff}=1/(1+\bar {S}(\bar {S}>0))$. We have 
$\tilde {z}_d\approx \tilde {h}_{eff}$ for all simulations and 
$\tilde {z}_d\approx L_{db}+L_{cp}\approx L_{db}+L_{cp}+L_{ce}$ since 
$L_{ce}$ is small (cf. figure 8e–h), which means that the convective and stably stratified layers are sufficiently decoupled that 
$\tilde {h}_{eff}$ is an accurate estimate for the thickness of the well-mixed turbulent bulk. Note that 
$\tilde {h}_{eff}$, 
$\tilde {z}_d$ and 
$L_{db}+L_{cp}+L_{ce}$ are not always close to each other, as is the case here, but can be significantly different when there is strong entrainment (Couston et al. Reference Couston, Lecoanet, Favier and Le Bars2017).
Figure 9. (a) Thickness of the entrained layer 
$L_{ce}$ (cf. (3.24a–d)) as a function of 
$Ra_{eff}$. (b) Entrainment parameter 
$\mathcal {E}$ as a function of 
$Ra_{eff}$. (c) Dimensionless height 
$\tilde {z}_d$ of the 
$\tilde {T}_d$ isotherm as a function of the expected convective layer thickness 
$\tilde {h}_{eff}$. The solid line shows 
$\tilde {z}_d=\tilde {h}_{eff}$.
The conclusions of this section are that (i) the entrainment of stable fluid can modify the temperature of the well-mixed bulk in subglacial lakes but that (ii) the influence of the stably stratified layer on the convective dynamics does not increase with the effective Rayleigh number or heat flux and water depth. Thus, we predict that the convective and stably stratified layers can be considered decoupled at leading order and that classical scaling laws for 
$Nu$ and 
$Re$ apply to natural subglacial lakes with large water depths and heat fluxes, provided that the effective Rayleigh number in (3.3) is considered. Further, we remark that the (finite) lowering of the bulk temperature from 
$(\tilde {T}_{bulk}-\tilde {T}_{d})\approx (\langle \tilde {T}_{b}\rangle -\tilde {T}_{d})/2$ to 
$(\tilde {T}_{bulk}-\tilde {T}_{d})\approx (\langle \tilde {T}_{b}\rangle -\tilde {T}_{d})/4$ in the presence of a thick stable layer can be accounted for by changing the effective temperature difference driving the convection from 
$\tilde {\varDelta }_{eff}$ to 
$3\tilde {\varDelta }_{eff}/2$. Multiplying the Nusselt number by 2/3 and the effective Rayleigh number by 3/2 for the simulation results with a stable layer indeed provides an improved collapse with the simulation results without a stable layer (cf. light-coloured symbols in figure 6a–c). We note that the lowering of the normalized bulk temperature in the presence of a stable layer, although relatively constant over the range of simulations considered (cf. figure 7b), should be expected (in general) to vary with the problem parameters (as is the case in Couston et al. Reference Couston, Lecoanet, Favier and Le Bars2017), such that the 3/2 correction proposed for the effective temperature difference is unlikely to be universal. Investigating the functional form of the correction coefficient for the effective temperature difference with the problem parameters in the context of real subglacial lakes would be interesting but is beyond the scope of the present work and is not expected to significantly change the current leading-order predictions, which neglect the influence of the stable layer on the convection.
4. Geophysical discussion
The closed-form expressions presented in table 4 are applicable to subglacial lakes that are either in the OB or LEC regime. We recall that the criterion for the OB regime (cf. (3.12)) is 
$\tilde {\beta }_{eff}\approx 1$, i.e. 
$\bar {S}\approx -1$, which can be satisfied by fully convective lakes only. The criterion for the LEC regime is instead 
$\tilde {\beta }_{eff}\approx \tilde {\varDelta }_{eff}(1+\bar {S})/2$, i.e. 
$\tilde {\varDelta }_{eff} \gg 2\bar {S}(\tilde {T}_d<0)/(1+\bar {S})$, which is satisfied when 
$\bar {S}\geq 0$, i.e. all lakes with a top stable layer, or 
$\tilde {\varDelta }_{eff} \gg |\bar {S}| \approx 0$, which may be verified for some fully convective lakes. We show in figure 10(a) 
$(1+\bar {S})$ as a function of ice pressure and lake water depth. We consider a heat flux of 
$F=68$ mW m
$^{-2}$, which is expected to be the average geothermal flux across Antarctica (Martos et al. Reference Martos, Catalán, Jordan, Golynsky, Golynsky, Eagles and Vaughan2017). For 
$1+\bar {S}<0$, 
$\bar {T}_b < T_d$, such that the thermal expansion coefficient is negative throughout the full water column and the lake is stable. We have 
$1+\bar {S}<0$ where the water depth 
$h$ and ice overburden pressure 
$p_i$ are small (hashed region in the top left corner), so small that 
$F=68$ mW m
$^{-2}$ is smaller than the threshold heat flux 
$F_t(p_i,h)$ (cf. (2.10a,b)). We have 
$1+\bar {S} \geq 1$ for all lakes with 
$p_i\leq p_*$ (as expected) and 
$1+\bar {S} \geq 0.9$ for many fully convective lakes, including the three well-known subglacial lakes Ellsworth, Concordia and Vostok (cf. three filled circles below the 
$p_*$ isobar). This means that the LEC regime is applicable to all mixed convective and stably stratified lakes but also possibly to some of the fully convective lakes if 
$\tilde {\varDelta }_{eff}$ is not too small. We have 
$1+\bar {S}\leq 0.1$, which is the criterion for the OB regime, for a limited number of fully convective lakes only (bottom left corner below the 0.1 isoline of the diagram), because all lakes that are deep can experience relatively large temperature differences, such that the thermal expansion coefficient becomes variable.
Figure 10. (a) Colour map of 
$1+\bar {S}$, with 
$\bar {S}$ the base-state stiffness parameter, in 
$(p_i,h)$ space relevant to Antarctic subglacial lakes with heat flux 
$F=68$ mW m
$^{-2}$. The hashed region in the top left corner highlights subglacial lakes that are fully stable because 
$F=68$ mW m
$^{-2}$ is smaller than the threshold heat flux 
$F_t$. The horizontal dashed line highlights the 
$p_*$ isobar and corresponds also to the isocontour 
$1+\bar {S}=1$. The dotted blue lines are isocontours of 
$1+\bar {S}$ and the coloured circles (see legend) highlight the location of some of the well-known subglacial lakes in parameter space (Couston & Siegert Reference Couston and Siegert2021). (b) Same as (a) but for the dimensionless effective thermal driving 
$\tilde {\varDelta }_{eff}$. The dotted blue lines still show the isocontours of 
$1+\bar {S}$ whereas the green solid lines show isocontours of 
$\tilde {\varDelta }_{eff}$.
In order to investigate the applicability of the LEC regime to real subglacial lakes, we display in figure 10(b) the dimensionless effective temperature difference 
$\tilde {\varDelta }_{eff}$ predicted based on (3.15) as a function of 
$(p_i,h)$ and superimpose isocontours of 
$1+\bar {S}$ (same blue dotted lines as seen in figure 10a). We find that 
$\tilde {\varDelta }_{eff}$ is always smaller or much smaller than 
$|\bar {S}|$, i.e. such that the LEC regime is not applicable unless 
$\bar {S} \geq 0$, except for a narrow band of ice pressures 
$p_i$ close to 
$p_*$ (horizontal dashed line). For instance, the isoline 
$\tilde {\varDelta }_{eff}=10^{-2}$ (green solid line) is the rightmost boundary of the region where 
$\tilde {\varDelta }_{eff}\geq 10^{-2}$, and it can be seen that this region does not overlap with the region of 
$|\bar {S}|\leq 10^{-2}$, which is above and to the right of the 
$1+\bar {S}=0.99$ isocontour (blue dotted line), unless 
$\bar {S}\geq 0$ (or 
$p_i\leq p_*$). Thus, the 
$-\bar {S}(\tilde {T}_d<0)$ term in (3.12) is not negligible when 
$\tilde {T}_d<0$ (or 
$-1\leq \bar {S}<0$) and neither the OB regime nor the LEC regime is applicable to deep fully convective lakes.
Although neither limiting regime is rigorously applicable to deep fully convective lakes, we can still use the closed-form expressions in table 4 in order to provide some first-order estimates for the r.m.s. velocity and temperature difference driving the convection in both shallow subglacial Lake CECs, where the LEC regime applies, and deep subglacial lake Vostok. The most recent geophysical information on Lake CECs is that the heat flux should be around 
$F\approx 100$ mW m
$^{-2}$ and that the maximum water depth is 
$h\approx 300$ m (personal communication with N. Donoso from Centro de Estudios Científicos, Chile). Using 
$F=100$ mW m
$^{-2}$, 
$h=300$ m and using the approximate pre-factors and exponents 
$a= 0.25$, 
$c= 0.006$, 
$\alpha = 2/7$ and 
$\gamma = 3/5$ for the closed-form expressions in the LEC regime yields 
$V_{rms}\approx 2.4$ mm s
$^{-1}$ and 
$\varDelta _{eff} \approx 0.04\,^{\circ }$C, which corresponds to a bottom temperature 
$T_b = T_d + \varDelta _{eff} = -1.12\,^{\circ }$C, i.e. 
$0.74\,^{\circ }$C above freezing (
$T_f=-1.86\,^{\circ }$C). For lake Vostok, the geothermal flux is closer to 
$50$ mW m
$^{-2}$ and the maximum water depth is of the order of 1000 m (Siegert et al. Reference Siegert, Ellis-Evans, Tranter, Mayer, Petit, Salamatin and Priscu2001). Using 
$F=50$ mW m
$^{-2}$, 
$h=1000$ m and 
$a= 0.25$, 
$c= 0.006$, 
$\alpha = 2/7$ and 
$\gamma = 3/5$ yields in the LEC regime 
$V_{rms}\approx 4.2$ mm s
$^{-1}$ and 
$\varDelta _{eff} \approx 0.03\,^{\circ }$C, with 
$\varDelta _{eff}$ directly equal to 
$T_b - T_f$, i.e. the bottom temperature in excess of the freezing point (
$T_f=-2.83\,^{\circ }$C), since there is no stable layer in this case. Note that we obtain 
$\varDelta _{eff} \approx 0.005\,^{\circ }$C and 
$V_{rms}\approx 25$ mm s
$^{-1}$, i.e. a significantly larger r.m.s. velocity, using predictions in the OB regime for lake Vostok because the effective thermal expansion coefficient is strongly overestimated in the OB regime for deep fully convective lakes.
The velocities in the LEC regime are of the same order as the velocities 
$V_{rms}\approx 1$ mm s
$^{-1}$ for Lake CECs and 
$V_{rms}\approx 4$ mm s
$^{-1}$ for lake Vostok predicted by previous studies, including (Wüest & Carmack Reference Wüest and Carmack2000), who applied scaling laws of vertical convection in the rapidly rotating regime to lake Vostok, and, more recently, Couston & Siegert (Reference Couston and Siegert2021), who predicted the physical properties of most subglacial lakes based on three-dimensional non-rotating scaling laws. The rough agreement between the velocities predicted in lake Vostok by Wüest & Carmack (Reference Wüest and Carmack2000) and this work is far from trivial since the former assumes that the flow dynamics is dominated by rotation while the latter neglects rotation. The importance of rotation in subglacial lakes is an open question, which is beyond the scope of this manuscript. The agreement between the velocities predicted in lakes CECs and Vostok by Couston & Siegert (Reference Couston and Siegert2021) and this work is also better than could be expected as there is a significant difference in the scaling exponent 
$\gamma \approx 3/5$ calculated in this work based on two-dimensional simulations and the scaling exponent 
$\gamma \approx 1/2$ used in Couston & Siegert (Reference Couston and Siegert2021) and obtained based on three-dimensional simulations (King, Stellmach & Buffett Reference King, Stellmach and Buffett2013). The discrepancy in exponents is compensated by the difference in pre-factors for the case of lakes CECs and Vostok considered here but may result in inaccurate predictions of three-dimensional r.m.s. velocities in other conditions. Another difference between this work and Couston & Siegert (Reference Couston and Siegert2021) is that the latter study used conservative estimates for the effective thermal expansion coefficient, such that their predicted velocities can be expected to be smaller than those we estimate based on the closed-form expression for the thermal expansion coefficient derived in the LEC regime. We refer the reader to Couston & Siegert (Reference Couston and Siegert2021) for a detailed geophysical discussion of first-order predictions of physical properties in most Antarctic subglacial lakes, which is based on the assumption – which we validated in this work through novel DNS – that the convective and stably stratified dynamics are decoupled at leading order.
We have demonstrated that the penetration of convective motions tends to lower the temperature of the turbulent well-mixed bulk when the top stable layer is not too thin, i.e. much thicker than the thermal boundary layer, despite the fact that the convection and stably stratified layer can be considered dynamically decoupled at leading order. Couston & Siegert (Reference Couston and Siegert2021) showed that the top stable layer of Antarctic subglacial lakes grows from 0 m at 
$p_i=p_*$ to approximately 40 m thickness at 
$p_i=p_{atm}\approx 10^5$ Pa. The thickness of the thermal boundary layer can be estimated as 
$h_{eff}/(2Nu)$ and is less than 1 m for most Antarctic subglacial lakes (Couston & Siegert Reference Couston and Siegert2021). Thus, we expect that most mixed convective and stably stratified subglacial lakes have bulk temperature that can be lowered by entrainment, but note that the decrease of the bulk temperature in three-dimensional simulations may differ from the lowering we reported in § 3.4 and figure 7(b).
5. Concluding remarks
We have investigated the dynamics of laboratory-scale subglacial lakes via DNS and demonstrated that the Nusselt number 
$Nu$ and Reynolds number 
$Re$ follow similar scalings laws as in classical Rayleigh–Bénard convection provided that an effective Rayleigh number 
$Ra_{eff}$ is considered. We obtained pre-factors and exponents similar to those obtained for two-dimensional fully convective Rayleigh–Bénard simulations in the OB regime with 
$Pr=1$ (Johnston & Doering Reference Johnston and Doering2009; Sugiyama et al. Reference Sugiyama, Calzavarini, Grossmann and Lohse2009), which means that the effect of the Prandtl number is relatively weak. As in Johnston & Doering (Reference Johnston and Doering2009), we remark that it is possible to define an effective flux-based Rayleigh number, i.e. 
$Ra_{Feff} = ( gh_{eff}^4F\beta _{eff})/( k\nu \kappa )$, which is equal to 
$\overline {Ra}_F$ in the OB regime (i.e. when the thermal expansion coefficient is constant) and is related to 
$Ra_{eff}$ through 
$Ra_{Feff} = NuRa_{eff}$.
We have shown that dimensional variables, such as the effective temperature difference, or bottom temperature, and r.m.s. velocity, scale differently with the problem parameters depending on whether the effective thermal expansion coefficient 
$\beta _{eff}$ is constant or linearly proportional to the effective temperature difference 
$\varDelta _{eff}$. We have called the dynamical regimes associated with the two limiting behaviours of 
$\beta _{eff}$ the OB regime and the LEC regime and derived explicit expressions for all variables of interest in terms of the problem parameters in both regimes (table 4). When 
$\beta _{eff}$ is an affine function of temperature, with non-negligible 
$y$-intercept, it is possible to infer 
$\varDelta _{eff}$ from (3.13) (and deduce all other variables thereafter) but an explicit expression is not available. We remark that the expressions for, e.g. 
$\varDelta _{eff}$, with the problem parameters are discontinuous between the OB regime and the LEC regime but can be continued and connected through the use of (3.13) in regions of the parameter space where neither regime is accurate.
The key results of our work are:
(i) the definition of an accurate
$Ra_{eff}$ for subglacial lakes leading to classical scaling laws for $Nu$ and $Re$ with $Ra_{eff}$;(ii) the demonstration that the convective and stably stratified layer dynamics are decoupled at leading order;
(iii) the identification of the two limiting OB and LEC regimes;
(iv) the derivation of closed-form expressions for all variables of interest in terms of the problem parameters in the OB and LEC regimes.
The numerical predictions of physical variables in subglacial lakes is briefly discussed in § 4. We emphasize that while we expect that the expressions hold in both two and three dimensions, the pre-factors 
$a$, 
$c$ and exponents 
$\alpha$, 
$\gamma$ may vary between two-dimensional and three-dimensional simulations. Clearly, future predictions of physical variables in subglacial lakes should consider power laws obtained in three dimensions, although we remark that the power law for 
$Nu$ is almost the same between two-dimensional and three-dimensional simulations of both classical Rayleigh–Bénard convection (Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009) and water convection close to the density maximum (Wang et al. Reference Wang, Zhou, Wan and Sun2019).
We have demonstrated that the penetration of convective motions in the stratified layer decreases with 
$Ra_{eff}$, such that there should be limited entrainment in most two-layer subglacial lakes. This does not mean that future studies should discard the stratified layer altogether. For instance, it would be interesting to investigate if internal waves excited by convection (Couston et al. Reference Couston, Lecoanet, Favier and Le Bars2018) can melt the ice ceiling such that ice-trapped oxygen and nutrients become available in subglacial lakes with a thin ice cover. The independence of the bulk temperature or thickness of the entrained layer with increasing 
$Ra_{eff}$ or decreasing 
$S_{eff}$, which we observed for simulations with a thick stable layer, seems at odds with the results of Couston et al. (Reference Couston, Lecoanet, Favier and Le Bars2017). Couston et al. (Reference Couston, Lecoanet, Favier and Le Bars2017) demonstrated a monotonic lowering of the bulk temperature and increase of the entrained layer thickness with a decrease of their input stiffness 
$S_i$, which, intuitively, is what may be expected if the inverse of the stiffness does provide a measure of the degree of coupling between the convective and stably stratified layers. We expect that the linear dependence of the thermal expansion coefficient with the temperature in our numerical simulations may be responsible for the observed discrepancy between our results and the results of (Couston et al. Reference Couston, Lecoanet, Favier and Le Bars2017), who instead used a piecewise-constant thermal expansion coefficient. On the one hand, in our numerical simulations, plumes approaching the 
$\tilde {T}_d$ isotherm lose their buoyancy even before entering the stable layer, which limits penetration. On the other hand, the base of the stable layer always has a small negative thermal expansion coefficient, i.e. a small Brunt–Väisälä frequency, such that penetration is always possible, i.e. including at large stiffness 
$S_{eff}$. Thus, the stiffness parameter can be expected to have a weaker effect on entrainment when the thermal expansion coefficient changes sign smoothly rather than discontinuously. A study designed specifically to investigate the effect of the functional form of the thermal expansion coefficient with temperature on entrainment would be a valuable addition to the literature on penetrative convection but is beyond the scope of this work.
Future works could investigate the effect of planetary rotation, which is most important near the poles and may play a significant role in the turbulent dynamics and mean vertical temperature profiles obtained in Antarctic subglacial lakes (Wüest & Carmack Reference Wüest and Carmack2000), and low to moderate salt concentrations. Considering planetary rotation will require three-dimensional simulations, unless we assume decoupled rotating convective and stably stratified dynamics and use existing results on the effect of planetary rotation in classical Rayleigh–Bénard convection (Plumley & Julien Reference Plumley and Julien2019). An important limitation of the present work is the assumption of a flat ice–water interface, such that considering a tilted ice ceiling and investigating the combined dynamics of the resulting baroclinic horizontal flow with the vertical convection is essential. There is also a possibility that the proposed scalings may become inaccurate as the water depth becomes large enough that the effect of pressure variations within the water column on buoyancy can no longer be neglected.
We expect that this paper and future studies improving our understanding of the hydrodynamic conditions in subglacial lakes could help identify subglacial environments that are physically favourable for a biome and guide future observations and sampling of subglacial lake water.
Acknowledgements
I gratefully acknowledge fruitful discussions with B. Favier and T. Alboussiére as well as many constructive comments from four anonymous reviewers, which helped me significantly extend and improve the initial manuscript.
Funding
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement 793450. I acknowledge PRACE for awarding me access to Marconi at CINECA, Italy.
Declaration of interest
The author reports no conflict of interest.
Appendix. Pressure effects
All simulations discussed in the main text and listed in table 2 assume that the thermal expansion coefficient can be approximated as 
$\beta (p,T)\approx \beta (p_i,T)$, i.e. such that pressure variations within the water column are neglected in the expression for 
$\beta$. We demonstrate that this is a valid approximation at leading order by showing in figure 11 the time history and vertical profiles of several variables obtained in simulation 
$\mathcal {S}_{14}^2$, which is simulation 
$\mathcal {S}_{1}^2$ with 
$h=4$ m and in an additional simulation, denoted 
$\mathcal {S}_{14+}^2$, which is the same as 
$\mathcal {S}_{14}^2$ but with 
$\beta (p,T)\approx \beta (p_i+\rho _0g(h-z),T)$, i.e. such that it includes pressure variations due to hydrostasy in the expression of the thermal expansion coefficient. We select 
$\mathcal {S}_{14}^2$ as a point of comparison as the temperature of maximum density, 
$T_d$, which is the variable that is most sensitive to 
$p$ in the expression for 
$\beta$, is reached inside the water column (i.e. there is a top stable layer) and because it is one of the simulations with the largest water depth (4 metres), such that pressure variations due to hydrostasy are relatively large.
Figure 11. (a) Time history of the dimensionless effective temperature difference 
$\langle \tilde {T}_b\rangle _h-\tilde {T}_d(p_i)$ for simulation 
$\mathcal {S}_{14}^2$ (thick orange line), which is simulation 
$\mathcal {S}_{1}^2$ discussed in the main text with 
$h = 4$ m and constant 
$\tilde {T}_d$, and simulation 
$\mathcal {S}_{14+}^2$ (thin black line), which includes hydrostatic effects in 
$\tilde {T}_d$. (b) Same as (a) but for the Reynolds number. (c) Same as (a,b) but for the dynamic pressure r.m.s. 
$p'_{rms}$ normalized by 
$10^{-9}$ dbar. (d) Mean vertical profiles of the dimensionless temperature 
$\tilde {T}$ (solid lines) and of the dimensionless temperature of maximum density 
$\tilde {T}_d$ (dashed lines) for simulations 
$\mathcal {S}_{14}^2$ (thick orange lines) and 
$\mathcal {S}_{14+}^2$ (thin black lines).
Figures 11(a)–11(c) show that the dimensionless effective temperature difference, the Reynolds number and the dynamic pressure r.m.s. 
$p'_{rms} =[\langle ( p - \langle p \rangle _h )^2 \rangle _h]^{1/2}$, which we evaluate in (approximately) the middle of the convection zone at 
$z=0.5$, are similar in both cases. The small discrepancies between the two cases are that the dimensionless effective temperature difference is slightly smaller in 
$\mathcal {S}_{14+}^2$ than in 
$\mathcal {S}_{14}^2$ and that 
$Re$ and 
$p'_{rms}$ are slightly larger in 
$\mathcal {S}_{14+}^2$ than in 
$\mathcal {S}_{14}^2$. This suggests that 
$\mathcal {S}_{14+}^2$ is slightly more turbulent than 
$\mathcal {S}_{14}^2$, which is expected since the decrease of 
$T_d$ with decreasing 
$z$ in 
$\mathcal {S}_{14+}^2$ results in a slightly larger thermal expansion coefficient in 
$\mathcal {S}_{14+}^2$ than in 
$\mathcal {S}_{14}^2$ at depth. The variability of 
$T_d$ within the water column for 
$\mathcal {S}_{14+}^2$ can be seen in figure 11(d), which displays the mean vertical profiles of 
$\tilde {T}_d$ and 
$\tilde {T}$. It can be seen that the decrease of 
$\langle \tilde {T}_d\rangle$ with depth in 
$\mathcal {S}_{14+}^2$ is significantly smaller than the variability of 
$\langle \tilde {T}\rangle$, which is, ultimately, the reason why the neglect of pressure variations in the expression for the thermal expansion coefficient is valid, at least for the laboratory-scale subglacial lakes considered in this paper. For deeper (geophysical) subglacial lakes, it may be expected that the decrease of 
$\tilde {T}_d$ with depth results in slightly more turbulent conditions than could be predicted assuming 
$\beta (p,T)\approx \beta (p_i,T)$. The exact discrepancy due to the neglect or consideration of hydrostatic pressure variations in the expression for 
$\beta$ in geophysical lakes is beyond the scope of this work. We note that considering hydrostatic effects in the expression for 
$\beta$ does not incur any computational overhead, such that they should be included in future works. Figure 11(c) shows that the dynamic pressure r.m.s. is much smaller than 1 dbar (by 9 orders of magnitude), which is equal to the change of hydrostatic pressure over 1 metre. Thus, dynamic pressure variations can be safely neglected from the expression for 
$\beta$ in simulations of laboratory-scale subglacial lakes as well as in simulations of (most) geophysical subglacial lakes.
