1. Introduction
An important feature of freshwater lakes is that they have a nonlinear equation of state (EOS). This nonlinear EOS is nearly quadratic with temperature near the temperature of maximum density, which is above the freezing temperature of the water (e.g. 
$\tilde {T}_{md} \approx 3.98\,^\circ \textrm {C}$ for distilled water at atmospheric pressure). The significance of this nonlinearity for lakes has been recognised for over a century (Whipple Reference Whipple1898). As a result of the density maximum, cooling the surface of a water body that has a mean internal temperature below 
$\tilde {T}_{md}$ will stabilise the water column, resulting in the characteristic reverse temperature stratification found during the winter months in temperate lakes (Farmer Reference Farmer1975). Stratification with temperatures on opposite sides of 
$\tilde {T}_{md}$ will lead to cabbeling (nonlinear mixing), which has important implications for convection (Carmack Reference Carmack1979; Farmer & Carmack Reference Farmer and Carmack1981; Couston et al. Reference Couston, Lecoanet, Favier and Le Bars2017, Reference Couston, Lecoanet, Favier and Le Bars2018). Convection of this type is also relevant in other fields such as Arctic melt ponds (Kim et al. Reference Kim, Moon, Wells, Wilkinson, Langton, Hwang, Granskog and Rees Jones2018) and geologically sequestered carbon dioxide (Hewitt, Neufeld & Lister Reference Hewitt, Neufeld and Lister2013). This paper aims to understand how thermal convection is altered near 
$\tilde {T}_{md}$.
The theoretical studies of such a system date back to Veronis (Reference Veronis1963). He considered the linear stability of a fluid layer with fixed temperatures at the top and bottom boundaries; temperatures fixed on either side of 
$\tilde {T}_{md}$. Shortly thereafter, Townsend (Reference Townsend1964) performed a complementary experimental study, again with fixed temperatures at the top and bottom. Both studies showed that the nonlinear EOS results in a stable upper layer above a convectively unstable lower layer, which will preferentially mix the lower-layer temperature stratification (similar to figure 1(a) except with a bottom thermal boundary layer). Several studies have shown that convection can generate internal waves in the stable upper layer, which are particularly relevant for astrophysical applications (Bars et al. Reference Bars, Lecoanet, Perrard, Ribeiro, Rodet, Aurnou and Gal2015; Lecoanet et al. Reference Lecoanet, Le Bars, Burns, Vasil, Brown, Quataert and Oishi2015). Recently, Toppaladoddi & Wettlaufer (Reference Toppaladoddi and Wettlaufer2017) and Wang et al. (Reference Wang, Zhou, Wan and Sun2019) built on this previous work and performed two-dimensional and three-dimensional simulations of convection with a nonlinear EOS. Both studies highlighted that in addition to the traditional dimensionless parameters (the Rayleigh number and the Prandtl number) considered for this flow set-up, the nonlinear EOS introduces an additional independent parameter. The additional non-dimensional parameter quantifies the ratio of the temperature variation across the stable and unstable layers. We will show that a similar ratio is important here. In each of these studies, the top and bottom temperatures are fixed, and the results focus on the statistical steady state.
(a) Representative mean temperature (blue line) and density (magenta dashed-dot line) profiles with a nonlinear EOS. Note the presence of an upper stable thermal boundary layer (blue). The model's piecewise-linear profiles (4.2) are also included as a solid black line. (b) A diagram of the numerical domain. (c) Comparative mean temperature and density profiles for a linear EOS.
Most lakes do not reach a steady state, but warm and cool throughout the year. In addition, the dominant heat loss in these freshwater systems is through the water surface with a relatively insulated bottom. Motivated by these considerations, we study a box of warm fluid (
$\tilde {T}>\tilde {T}_{md}$) cooled from above (
$\tilde {T}<\tilde {T}_{md}$). Here, unlike the previous studies, the lower boundary is insulating, and the temperature stratification is transient. As is typical of these theoretical studies (Veronis Reference Veronis1963; Townsend Reference Townsend1964; Olsthoorn, Tedford & Lawrence Reference Olsthoorn, Tedford and Lawrence2019), we consider an EOS that is quadratic with temperature, which approximates the full EOS of Chen & Millero (Reference Chen and Millero1986) for temperatures below 
$10\,^\circ \textrm {C}$. In this paper, we refer to a cooling of the surface, though the problem is symmetric for a box of cold water that is heated from above, at least for a quadratic EOS assumed here. We want to understand how convection and the resultant heat transport is changed in the presence of this nonlinear density relationship, near 
$\tilde {T}_{md}$. In particular, we address the three following questions.
(i) Does the nonlinear EOS affect the vertical transport of heat within the water?
(ii) Which parameters determine the heat flux out of the water surface?
(iii) Can we predict the vertical heat transport and kinetic energy produced by the turbulent convection?
We begin, in § 2, with a discussion of the numerical methods used in this paper. Section 3 is an analysis of the three-dimensional direct numerical simulations, highlighting that the vertical transport of heat is significantly different for a nonlinear EOS than that predicted for a linear EOS. In § 4 we discuss the relevant parameters that control the heat flux and present a predictive model for this system. We show that the model agrees well with the data from the numerical simulations in § 5. Finally, we conclude in § 6.
2. Numerical set-up
We consider a body of water that is insulated from below and cooled from above. We restrict our analysis to freshwater, where the temperature of maximum density 
$\tilde {T}_{md}$ is above its freezing temperature. Figure 1(b) is a schematic of the numeric domain of interest. We ignore the effects of salinity, pressure and higher-order terms in the EOS such that the density (
$\tilde {\rho }$) is given as
\begin{equation} \tilde{\rho}= \tilde{\rho}_0 - C_T\left(\tilde{T} - \tilde{T}_{md} \right)^2. \end{equation}
Here, 
$C_T$ is a constant and 
$\tilde {\rho }_0$ is the density at 
$\tilde {T} = \tilde {T}_{md}$.
We are interested in parameterising the convective heat flux through the top boundary. As such, a natural time scale for this set-up is the diffusive timescale of heat 
$\tau _{\kappa }$ over the domain height 
$H$. That is,
\begin{equation} \tau_{\kappa} = \frac{H^2}{\kappa}, \end{equation}
where 
$\kappa$ is the diffusivity of heat. We then non-dimensionalise the spatial coordinates (
$\tilde {\boldsymbol {x}}$), the fluid velocity (
$\tilde {\boldsymbol {u}}$), time (
$\tilde t$) and temperature (
$T$) as
\begin{equation} {\boldsymbol{x}} = \frac{\tilde{{\boldsymbol{x}}}}{H} , \quad {\boldsymbol{u}} = \frac{\tilde{{\boldsymbol{u} }} H}{\kappa}, \quad t = \frac{\tilde{ t}}{\tau_{\kappa}}, \quad T = \frac{\tilde{T} - \tilde{T}_{md} }{\Delta \tilde{T}}. \end{equation}
Boldface variables denote vector quantities and 
$\Delta \tilde {T} = \tilde {T}_{md} - \tilde {T}(\tilde {z}=H)$. As an aside, we would like to highlight that, throughout the paper, we have dropped factors of depth 
$H =1$. Although we believe that this is a reasonable simplification, we highlight this convention to avoid any confusion for the reader.
The equations of motion for this flow are the Navier–Stokes equations under the Boussinesq approximation. These equations are written
\begin{gather} \left( \frac{\partial }{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \right)\boldsymbol{u} ={-} \boldsymbol{\nabla} P + {Ra}_0 \,{Pr}\,T^2 \hat{\boldsymbol{k}} + {Pr} \nabla^2 \boldsymbol{u}, \end{gather}
\begin{gather}\left( \frac{\partial }{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) T = \nabla^2 T, \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0. \end{gather}
We have defined the Rayleigh number (
${Ra}_0$) and Prandtl number (Pr) as
\begin{equation} {Ra}_0 = \frac{g C_T \Delta \tilde{T}^2}{\rho_0}\frac{H^3}{\kappa \nu}, \quad {Pr} = \frac{\nu}{\kappa}, \end{equation}
where 
$\nu$ is the molecular viscosity of water and 
$g$ is gravitational acceleration. The value of the Prandtl number for freshwater varies from 
$Pr\approx 7$ at 
$\tilde {T} = 20 \,^\circ \textrm {C}$ to 
$Pr \approx 13.4$ at 
$\tilde {T} = 0\,^\circ \textrm {C}$, largely owing to variations in 
$\nu$. Incorporating the functional dependence of 
$\nu$ on 
$T$ is outside the scope of this paper. Hay & Papalexandris (Reference Hay and Papalexandris2019) performed convective simulations with an evolving 
$Pr$, in a different context, that do not show significant changes to the vertical heat flux from the constant Pr case. Future work will discuss the role of the changing viscosity at low temperatures. In this paper, we assume a constant 
$\nu$ and prescribe 
$Pr=9$.
We enforce a fixed temperature on the surface and an insulating bottom condition. That is
\begin{equation} T\left|_{z = 1} ={-}1, \quad \frac{\partial T}{\partial z}\right|_{z=0} = 0, \end{equation}along with no-slip top and bottom velocity boundary conditions.
As we show, the well-mixed bottom water temperature 
$T_B$ is an important parameter in this system. As illustrated in figure 1(a), although the system is convectively unstable, the temperature profile below the top boundary layer is nearly uniform. In the numerical simulations, we determine 
$T_B$ by fitting the horizontally averaged temperature profile to an erf function. That is,
\begin{equation} \bar{T} \approx \left( 1 + T_B\right) \hbox{erf} \left( \frac{z-1}{\eta}\right) - 1, \end{equation}
where 
$\eta \approx \delta _{BL}$ (defined later) is a fit parameter and 
$\bar {\cdot } = (1/{\hbox {Lx Ly}})\int \cdot \, \textrm {d} x\,\hbox {d}y$. Over time, 
$T_B$ will decrease, whereas 
${Ra}_0$, by definition, will remain fixed.
In most cases, we selected the initial fluid temperature such that there was an initial symmetry between the bottom water temperature (
$T_B=1$) and the top boundary condition about 
$T=0$. For reference, for a water depth of 
$H=0.05 \ \text {m}$ with 
$\tilde {T}(\tilde {z}=H) =0\,^\circ \text {C}$, the initial water temperature would be 
$\tilde {T}\approx 8\,^\circ \textrm {C}$ and 
${Ra}_0\approx 8\times 10^{5}$. Thus, the simulations presented in this paper (see table 1) are on the scale of potential laboratory experiments. In the lowest 
${Ra}_0$ case, we selected 
$T_B (t=0) = 2$ so that the initial instability grew significantly quicker than the background diffusion.
Parameters for each numerical simulation.
Although the top boundary temperature will remain fixed (
$T(x,y,z=1,t)=-1$), 
$T_B$ will cool throughout the numerical simulations. While 
$T_B>0$, there exist two sub-domains to the mean temperature stratification: an upper stable stratification of depth 
$\delta _{St}$ where 
$\bar {T}<0$ and a lower hydrostatically-unstable stratification where 
$\bar {T}>0$. Figure 1(a) is a plot of a representative temperature profile within the fluid domain. The total top boundary-layer thickness between the well-mixed uniform temperature and the upper boundary is 
$\delta _{BL}>\delta _{St}$ (see figure 1a).
Before continuing, we highlight that for a linear EOS, 
$\tilde {T}_{md}$ cannot be internal to the fluid domain, and an upper stable layer cannot form. As such, for a linear EOS, we would expect a temperature profile to resemble that in figure 1(c).
2.1. Numerical implementation
We performed direct numerical simulations using SPINS (Subich, Lamb & Stastna Reference Subich, Lamb and Stastna2013). SPINS solves the Navier–Stokes equations under the Boussinesq approximation using pseudospectral spatial derivatives and a third-order time-stepping scheme. As the top boundary layer controls the dynamics of the initial instability and the subsequent convection, we implement a Chebyshev grid in the vertical, which clusters grid points at the domain boundaries. We assume periodic horizontal boundary conditions, implemented with fast Fourier transforms.
We performed five numerical simulations with a nonlinear EOS at different Rayleigh numbers. We performed one additional simulation with a linear EOS for comparison. The details of these numerical simulations are provided in table 1. In all six cases, the Rayleigh number was large enough for the system to become unstable (see Appendix C for the linear stability analysis). We initially perturb the three velocity components with a random perturbation sampled from a Normal distribution scaled by 
$10^{-2}$. The numerical resolution (
$\textrm {Nx} \times \textrm {Ny} \times \textrm {Nz}$) was selected, such that 
$\max ({\Delta x}/{\eta _B})<3$, where we compute the Batchelor scale 
$\eta _B = (\bar {\varepsilon })^{-{1}/4} \,Pr^{-{1}/2}$, for viscous dissipation rate 
$\varepsilon$ (see (4.8a,b)) and horizontal grid spacing 
$\varDelta x$. The vertical grid is clustered towards the boundaries and 
$\max ({\Delta z}/{\eta _B})<\max ({\Delta x}/{\eta _B})$, in all cases. This resolution criterion is similar to that employed in Hay & Papalexandris (Reference Hay and Papalexandris2019), Kaminski & Smyth (Reference Kaminski and Smyth2019), Olsthoorn et al. (Reference Olsthoorn, Tedford and Lawrence2019) and Smyth & Moum (Reference Smyth and Moum2000). The resolution sufficiency is highlighted in Appendix D.
3. Simulation results
In this section, we describe both the qualitative and quantitative dynamics of the numerical simulations as they relate to the transport of heat within the fluid domain. In particular, we show that the convection is self-similar in 
$T_B$ for the range of Rayleigh numbers considered. Our discussion is focused on a single representative case: case 3, 
$Ra=9.0\times 10^5$. The results are similar for the other simulations, except where otherwise noted.
Figure 2 contains snapshots of the temperature field for case 3. The left column of figure 2 contains plots of the temperature field at different times 
$t=\{2.50\times 10^{-3}, 5.00\times 10^{-3}, 7.50\times 10^{-3}, 1.00\times 10^{-2}, 1.25\times 10^{-2}, 5.00\times 10^{-2}\}$. These plots were made with VisIt's (Childs et al. Reference Childs2012) volume plot option that varies the transparency of the temperature field according to its value. The downwelling plumes are visible in the 
$X$–
$Z$ slices (figure 2b,e,h,k,n,q). Initially, the temperature stratification is linearly stable, and the temperature profile simply diffuses. At 
$t=0.0025$ (figure 2a–c), the temperature stratification matches that predicted by pure diffusion. Once the top thermal boundary layer has grown sufficiently, the system is unstable and small perturbations to the velocity and temperature field will grow. These near-modal perturbations are visible at 
$t=5.00\times 10^{-3}$ (figure 2d–f). Once the perturbations grow large enough, they begin to merge, forming dense plumes that transport cold fluid to the bottom of the domain. The columnar plumes are highly variable in both space and time, and the resultant convection will mix the bottom fluid. Eventually (figure 2p–r), the bottom water temperature decreases sufficiently such that the vertical heat flux rapidly decreases.
Snapshots of the temperature field for case 3: 
${Ra}_0=9.0\times 10^5, Pr=9$. The volume plots (a,d,g,j,m,p) encapsulate the three-dimensional structure of the flow field. The middle figure panels (b,e,h,k,n,q) contour plots of vertical (
$x-z$) temperature slices at 
$y=-1$. Similarly, the right figure panels (c,f,i,l,o,r) are contour plots of horizontal (
$x-y$) temperature slices at 
$z=0.9$. Snapshots are given at (a–c) 
$t=2.50\times 10^{-3}$, (d–f) 
$t=5.00\times 10^{-3}$, (g–i) 
$t=7.50\times 10^{-3}$, (j–l) 
$t=1.00\times 10^{-2}$, (m–o) 
$t=1.25\times 10^{-2}$ and (p–r) 
$t=5.00\times 10^{-2}$. Note the jump in output times highlighted by the horizontal dashed line.
Slices of the spanwise (
$X$–
$Y$) structure of the temperature field at 
$z = 0.9$ are presented in figure 2(c,f,i,l,o,r). Although a complete analysis of this structure is beyond the scope of this paper, we highlight it here as it demonstrates that the spacing between the downwelling plumes is increasing over time. The horizontal length scale of the three-dimensional flow structure increases as 
$T_B\to 0$, as viscous dissipation preferentially diffuses small-scale motions. Eventually, the spacing between the plumes reaches the size of the domain, which limits the run time of these simulations. We have performed several simulations at different domain sizes and have determined that the present results are not domain-size dependent.
We highlight the mean structure of the temperature stratification in figure 3(a–f), which are plots of individual horizontally averaged temperature profiles at times 
$t=\{2.50\times 10^{-3}, 5.00\times 10^{-3}, 7.50\times 10^{-3}, 1.00\times 10^{-2}, 1.25\times 10^{-2}, 5.00\times 10^{-2}\}$ (identical to figure 2). The evolution of the temperature stratification is further highlighted in the contours of horizontally averaged temperature over time (figure 3g). These contours show that the boundary-layer depth increases over time and that the bottom water temperature decreases. Note that there is a weak thermal gradient near the bottom of the domain owing to the solid boundary. This thermal gradient becomes weaker with increasing Ra, owing to the increased energy within the system.
The horizontally averaged temperature profiles are plotted at the output time of the snapshots of figure 2: (a) 
$t=2.50\times 10^{-3}$, (b) 
$t=5.00\times 10^{-3}$, (c) 
$t=7.50\times 10^{-3}$, (d) 
$t=1.00\times 10^{-2}$, (e) 
$t=1.25\times 10^{-2}$ and (f) 
$t=5.00\times 10^{-2}$. Contours of the horizontally averaged temperature 
$(\bar {T})$ over time are also plotted (g). The vertical dashed lines in panel (g) represent the time of the vertical temperature profiles. Data are shown for case 3: 
${Ra}_0=9.0\times 10^5$.
3.1. Vertical heat transport
Now that we have presented the essential features of the temperature evolution, we proceed to quantify the heat loss resulting from the thermal convection. The rate of heat loss of the water in the domain is
\begin{equation} \frac{\mathrm{d}}{\mathrm{d} t} \int_0^1 \int_A T \,\hbox{d} A \,\hbox{d} z ={-}F A, \end{equation}
where 
$F$ is the average outward heat flux at the domain surface, and 
$A=\hbox {Lx Ly}$ is the area of the water surface. The vertical transport of heat is approximately constant over the upper conductive boundary layer.
The vertical heat flux within the convective layer is significantly greater than the diffusive flux. We define 
$\sigma$ as the ratio of average heat flux to the diffusive heat flux over the convective domain. That is, we write
\begin{equation} \sigma \approx \frac{F/\left( 1 - \delta_{St}\right)}{\left(T_B - T_0\right)/\left( 1 - \delta_{St}\right)} = \frac{F}{T_B - T_0}, \end{equation}
where 
$T_0$ is the temperature of maximum density within the domain. As we show, for a quadratic EOS, 
$\sigma$ is constant for a significant portion of the simulation time, which indicates that the temperature decays exponentially. For a linear EOS, 
$T_0$ will be equal to the top boundary temperature (
$T_0=-1$). In this case, 
$\sigma$ is identical to the Nusselt number (Nu), which is the ratio of the average vertical heat flux to the diffusive heat flux across the entire domain. For the nonlinear EOS discussed here, where the temperature of maximum density is interior to the fluid domain, 
$T_0 = 0$.
As an aside, it is worth noting that 
$F$ is a function of the temperature difference 
$T_B - T_0$. As the water becomes uniform, 
$( T_B - T_0) \to 0$, the vertical flux 
$F\to 0$ such that 
$\sigma$ is bounded.
As a simple model, we first consider the scenario where the water column is uniformly mixed to some temperature 
$T_U$. If 
$\sigma$ were a constant 
$\sigma _0$, then
\begin{equation} \frac{\mathrm{d}}{\mathrm{d} t} T_U A ={-}\sigma_0 \cdot \left(T_U - T_0\right) A \implies T_U = B \exp \left( -\sigma_0 t\right) + T_0, \end{equation}
where 
$B$ is a constant of integration. That is, 
$T_U$ decays exponentially in time. In principle, there is no a priori reason to suspect that 
$\sigma$ will be constant and, in general, it is not. However, as we show, for the initial convective evolution, 
$\sigma$ is constant, and the temperature will decay exponentially.
Before continuing, we return to the first of the main questions we are trying to answer in this paper. Does the nonlinear EOS affect the vertical transport of heat out of the domain? We performed a single numerical simulation with a linear EOS with 
${Ra}_0 = 9.0 \times 10^5$ and 
$Pr=9$. The boundary and initial conditions remain unchanged. For a linear EOS, there are only two free parameters: a Prandtl number (Pr) and a time-dependent Rayleigh number. The relevant Rayleigh number for a linear EOS includes the density difference across the domain 
$\Delta \rho$ as
\begin{equation} {Ra}_{Lin} = \frac{g \Delta \rho}{\rho_0} \frac{H^3}{\kappa \nu} = {Ra}_0 \left(1 + T_B\right). \end{equation}
Here, we have related the constant Rayleigh number (
${Ra}_0$) defined in this paper with a more traditionally defined time-varying Rayleigh number (
${Ra}_{Lin}$), used with a linear EOS. The equivalent effective Rayleigh number (
${Ra}_{eff}$) for the quadratic EOS is given by
\begin{equation} {Ra}_{eff} = {Ra}_0 T_B^2. \end{equation}
We will return to this in our discussion of 
$\sigma$ later.
For a linear EOS, it is empirically determined that
\begin{equation} {Nu}_{Lin} \propto {Ra}_{Lin}^n. \end{equation}
While significant controversy persists over the exact value of 
$n$ (e.g. see Plumley & Julien Reference Plumley and Julien2019), we will specify 
$n=0.28$ as it best fits the data discussed.
Figure 4 is a comparison plot of (a) temperature and (b) vertical heat flux as a function of time, with a comparison between a linear and quadratic EOS. The scaling (3.6) is included in panel (b). Here, and for the rest of the paper, we plot in grey the initial transition period before reaching a quasi-equilibrium state (where the vertical buoyancy flux approximately balances viscous dissipation). We observe that the curvature in the EOS fundamentally changes the temperature evolution of the system over time. The surface heat flux is significantly greater for the linear EOS, resulting in the bottom water temperature 
$T_B$ decaying much faster for a linear EOS. For the quadratic EOS, the presence of a temperature of maximum density also restricts the convective mixing such that 
$T_B\to 0$ (note that molecular diffusion will eventually reduce 
$T_B\to -1$ as 
$t\to \infty$). If we return to the dimensional example of a 0.05 m deep container with a surface temperature at 
$0\,^\circ \textrm {C}$, then by 
$t=0.05$ (
${\approx }15$ min) there is a 
$1.5\,^\circ \textrm {C}$ difference between the internal temperatures of the linear and quadratic EOS; nearly a 40 % increase in the heat loss.
A comparison plot of (a) 
$T_B$ and (b) 
$F$ versus time between convection with a linear or a quadratic EOS. The temperature of maximum density is included as a horizontal dashed line in (a). The scaling (3.6) is included as a dashed line in panel (b). In both cases, 
${Ra}_0=9.0\times 10^5$. We plot in grey the initial (diffusive) transition period.
The Nusselt number dependency in (3.6) is inadequate to describe the temperature evolution for a quadratic EOS. In the next section, we derive a model for the vertical heat flux (
$F$) and turbulent kinetic energy (
$\hbox {TKE}$) density of convection with a nonlinear EOS. We show that the convection is fundamentally dependent on three independent parameters, as opposed to the two needed with a linear EOS.
4. Scaling laws
We begin with a Reynolds decomposition of the temperature and velocity field into a mean temperature profile and fluctuations from it, where
\begin{equation} \boldsymbol{u} = \boldsymbol{0} + \boldsymbol{u}', \quad T = \bar{T} + T'. \end{equation}
In the numerical simulations, we take 
$\overline{\cdot}$ as the horizontal average through the domain. For this scaling analysis, we simplify 
$\bar {T}$ as a piecewise linear profile
\begin{equation} \bar{T} = \begin{cases} -1 - \dfrac{1 + T_B}{\delta_{BL}} \left(z - 1 \right) & z> 1 - \delta_{BL},\\ T_B & z\leqslant 1- \delta_{BL}, \end{cases}\end{equation}
as in figure 1(a). Note that 
$\delta _{BL}$ is the total transition thickness with 
$\delta _{BL}>\delta _{St}$.
4.1. Boundary-layer thickness 
$\delta _{St}$
The diffusion of heat through the top boundary is
\begin{equation} F ={-} \left.\frac{\partial \bar{T}}{\partial z} \right|_{z=1} = \frac{1 + T_B}{\delta_{BL}} = \frac{1}{\delta_{St}}. \end{equation}
That is, in this non-dimensionalisation, the boundary-layer thickness determines solely the outward heat flux. As a reminder, 
$\delta _{St}$ is the top stable boundary-layer thickness between the top boundary and where 
$\bar {T}=0$ (the temperature of maximum density).
Substituting the Reynolds decomposition (4.1a,b) into the temperature evolution equation (2.5), we derive the evolution equation for the mean temperature profile,
\begin{equation} \frac{\partial \bar{T}}{\partial t} ={-}\frac{\partial }{\partial z} \left( -\frac{\partial \bar{T}}{\partial z} + \overline{T'w'} \right). \end{equation}Balancing the heat fluxes within the domain and top boundary layer, we determine
\begin{equation} -\underbrace{\frac{\partial T_B}{\partial t} \left(1-\left(1 + T_B\right)\delta_{St}\right)}_{\hbox{Interior Cooling}} + \underbrace{\frac{\left(1 + T_B\right)^2}{2} \frac{\partial \delta_{St}}{\partial t} }_{\hbox{Increased}\ \delta_{St}} = \underbrace{\frac{1}{\delta_{St}}. }_{\hbox{Outward Heat Flux}} \end{equation}
Appendix A provides a derivation of (4.5). We can see from (4.5) that a fraction of the heat loss is derived from the cooling of the interior fluid. The remainder of the heat loss is given by an increase in the boundary-layer thickness 
$\delta _{St}$. The surface heat loss is controlled by the boundary-layer thickness 
$\delta _{St}$.
Building on § 3.1, the average temperature within the domain is 
$T_B$ such that
\begin{equation} \frac{\textrm{d} T_B}{\textrm{d}t} ={-}\sigma T_B. \end{equation}
As convection greatly enhances the vertical flux of temperature, we find that 
$\sigma \gg 1$. We continue a discussion of 
$\sigma$ later. To leading order in 
$O(1/\sigma )$, the dominant balance between the interior cooling and the outward heat flux determines that
\begin{equation} \delta_{St} \sim \frac{1}{\sigma T_B}, \quad F = \frac{1}{\delta_{St}} \sim \sigma T_B, \quad \sigma\gg 1. \end{equation}4.2. TKE
Similar to (4.4), the volume-integrated TKE density (
$\hbox {TKE} = \frac {1}{2}\boldsymbol {u}^\prime \boldsymbol {\cdot } \boldsymbol {u}^\prime$) evolution equation is written
\begin{equation} \frac{\mathrm{d}\langle \hbox{TKE}\rangle_V }{\mathrm{d} t} ={-} {Ra}_0 \,Pr \langle w^\prime \rho^\prime \rangle_V - \varepsilon, \quad \varepsilon =Pr \langle \boldsymbol{\nabla} \boldsymbol{u}^\prime : \boldsymbol{\nabla} \boldsymbol{u}^\prime \rangle_V . \end{equation}
The 
$(:)$ operator is the double dot product. Unlike the case of a linear EOS, the vertical buoyancy flux (
$\overline {w^\prime \rho ^\prime }$) is not directly proportional to the vertical temperature flux (
$\overline {w^\prime T^\prime }$). For a quadratic EOS, the vertical buoyancy flux is the sum of two components,
\begin{equation} \overline{w^\prime \rho^\prime} ={-} 2 \overline{w^\prime T^\prime} \, \bar{T} - \overline{w^\prime T^\prime T^\prime}. \end{equation}
Along with the mixing coefficient (
$\varGamma$, the ratio of the vertical buoyancy flux to viscous dissipation), we define a buoyancy flux ratio (
$\lambda$) to quantify the relative contribution of each buoyancy flux term. These are defined as
\begin{equation} \varGamma = \frac{ -{Ra}_0 \,Pr \langle \overline{w^\prime \rho^\prime}\rangle_V }{\varepsilon}, \quad \lambda = \frac{-\langle \overline{w^\prime T^\prime T^\prime}\rangle_V }{\langle 2 \overline{w^\prime T^\prime} \ \bar{T} \rangle_V }. \end{equation}Assuming self-similarity of TKE, we can show that (see Appendix B for details)
\begin{equation} \langle \hbox{TKE}\rangle_V \sim \tfrac 12 \left( \varGamma^{{-}1} - 1 \right) \left( 1 - \varLambda \right) {Ra}_0 \,Pr T_B^2, \quad \sigma \gg 1, \end{equation}
where 
$\varGamma$ is assumed constant and
\begin{equation} \varLambda = \frac{2}{T_B^2} \int_{T_B} T' \lambda(T') \,\textrm{d}T'. \end{equation}This model for the kinetic energy is consistent with the empirical Reynolds number scaling found in Wang et al. (Reference Wang, Zhou, Wan and Sun2019).
Which parameters control the heat flux out of the water surface? The scaling laws (4.7a–c) and (4.11a,b) depend on three undetermined coefficients: the decay rate 
$\sigma$, the buoyancy flux ratio 
$\lambda$ and the mixing coefficient 
$\varGamma$. These three coefficients are, as we show, functions of the three parameters 
${Ra}_0$, 
$Pr$ and 
$T_B$. In the next section, we derive the functional form of 
$\sigma$, 
$\lambda$ and 
$\varGamma$, and show that the models provided here agree well with the simulated quantities.
5. Model comparison
As discussed previously, 
$T_B$ initially decays exponentially with time. Figure 5(a) is a plot of 
$T_B$ as a function of time on a semi-log axis for all five of the nonlinear EOS simulations listed in table 1. We compute the instantaneous decay rate as
\begin{equation} \sigma = \frac{-1}{T_B} \frac{\mathrm{d}T_B}{\mathrm{d} t} , \end{equation}
and we observe (figure 5b) that, after an initial transient, 
$\sigma$ is nearly constant over a range of 
$T_B$, with 
$\sigma \gg 1$ (
$T_B\gtrsim 0.37$). During this period we perform a linear regression to determine the 
${Ra}_0$-number dependence of the flow. Figure 5(c) is a plot of 
$\sigma$, scaled by 
${Ra}_0^{0.28}$, which collapses the data for all of the simulations with a nonlinear EOS. For lower values of 
$T_B$ (
$T_B\lesssim 0.37$), 
$\sigma$ decreases with 
$T_B$. Once the system has achieved a quasi-steady state, we fit piecewise curves to the different cases and approximate
\begin{equation} \sigma = 0.97 {Ra}_0^{0.28} \begin{cases} 1, & T_B \geqslant 0.37\\ \left(\dfrac{T_B}{0.37}\right)^{0.63}, & T_B < 0.37 \end{cases}. \end{equation}
We find a regime change at 
$T_B\approx 0.37$, between a constant exponential decay (
$T_B\gtrsim 0.37$), and when 
$\sigma$ rapidly decreases (
$T_B\lesssim 0.37$). We note that the model presented in § 4 is applicable only for 
$\sigma \gg 1.$ As 
$T_B\to 0$, 
$\sigma$ decreases and higher-order corrections are increasingly significant. It is also worth noting that case 1 (
${Ra}_0=9.0\times 10^4$) has the lowest value of 
$\sigma$ and 
${Ra}_0$, and does not collapse as well as the other cases.
(a) Plot of the bottom water temperature 
$T_B$ with time for all numerical runs with a nonlinear EOS. (b) Plot of the instantaneous decay rate (
$\sigma$) as a function of 
$T_B$. (c) Plot of the scaled decay rates (
$\sigma$), which collapse on to a single curve. We include the scaling (5.2) as a dotted black line in panel (c). We plot in grey the initial transition period before reaching a quasi-equilibrium state. Note that the 
$x$-axis has been reversed in panels (b,c), with 
$T_B$ decreasing towards the right such that it reads in the same direction as time.
In connection with the linear EOS, we might expect that 
$\sigma$ should scale with an effective Rayleigh number (
${Ra}_{eff}$) (as was highlighted in Anders et al. Reference Anders, Vasil, Brown and Korre2020). This is partially correct. To justify this statement, we approximate equation (5.2) as
\begin{equation} \frac{\sigma}{\sigma_0} \approx \min \left\{ {Ra}_{max}^{0.28}, {Ra}_{eff}^{0.28}\right\} , \quad {Ra}_{max} = {Ra}_{eff}(T_B = T_{max}), \end{equation}
where 
$\sigma _0 = 0.97 T_{max}^{-0.63}$ and 
$T_{max} = 0.37$. That is, the temperature decay rate 
$\sigma$ does increase with 
${Ra}_{eff}$ below some threshold value. For large enough 
${Ra}_{eff}$, the convection is sufficiently turbulent that 
$\sigma$ remains constant. Thus, the convection is self-limiting. We believe this results from the increased diffusion at the top boundary owing to the stable thermal layer that is not present in a linear EOS. Importantly, as highlighted in Toppaladoddi & Wettlaufer (Reference Toppaladoddi and Wettlaufer2017) and Wang et al. (Reference Wang, Zhou, Wan and Sun2019), a third parameter must be defined to characterise the convection. We define this third parameter as 
$T_B$.
The mixing coefficient 
$\varGamma$ and buoyancy flux ratio 
$\lambda$, defined in (4.10a,b), are also functions of the non-dimensional parameters. Figure 6 is a plot of (a) 
$\lambda$ and (b) 
$\varGamma$ as a function of 
$T_B$. We find that 
$\lambda$ increases as 
$T_B$ decreases. Fitting 
$\lambda$ as a function of 
$T_B$ with piecewise-functions, we estimate that while in quasi-steady state,
\begin{equation} \lambda = 0.24 \begin{cases} T_B^{{-}1/2}, & T_B \geqslant 0.35,\\ 1, & T_B < 0.35 . \end{cases} \end{equation}
As the data are noisy, our estimated power-law dependencies are approximate approximations, preferentially weighting the fit values of the higher 
${Ra}_0$ cases. We anticipate future work to illuminate a theoretical prediction for the appropriate scaling laws. In addition, it is important to note that case 1 (
${Ra}_0=9.0\times 10^{4}$) does not follow the trend of the other cases. For case 1, the low Rayleigh number results in a significant diffusive contribution to the total heat flux and is in a weakly unstable regime. As in (5.2), we find that there is a regime change that occurs at 
$T_B\approx 0.37$. Similarly, we find that
\begin{equation} \varGamma\approx 0.96 \end{equation}
for all 
${Ra}_0$, though large fluctuations are present. This value of 
$\varGamma$ indicates that the rate of viscous dissipation is nearly equal to the vertical buoyancy flux. In steady state, by definition, 
$\varGamma =1$ which is consistent with our observation that the cooling box is nearly in steady state.
Plots of (a) the buoyancy flux ratio 
$\lambda$ and (b) the mixing coefficient 
$\varGamma$ as a function of 
$T_B$. We include the fit of (5.4) as a dashed line in (a). We further include a dashed line at 
$\varGamma = 0.96$ in (d). We plot in grey the initial transition period before reaching a quasi-equilibrium state.
5.1. Model agreement
Can we predict the vertical heat transport and kinetic energy produced by the turbulent convection? In § 4, we derived scaling laws for 
$\delta _{St}$, 
$F$ and 
$\hbox {TKE}$ as a function of 
${Ra}_0$, 
$Pr$ and 
$T_B$. Figure 7 is a comparison plot between the model equations (4.7a–c) and (4.11a,b) and the numerically computed (a) 
$\delta _{St}$, (b) 
$F$ and (c) 
$\hbox {TKE}$. The model agrees well with the data from the direct numerical simulations. We have scaled the 
$y$-axes by the appropriate powers of 
${Ra}_0$ and 
$Pr$. Note that the scaling coefficient 
$0.02$ in panel (c) is equivalent to 
$\varGamma =0.96$. The parameters 
$\sigma$ and 
$\varLambda$ are defined using the fit equations (5.2) and (5.4), respectively.
We further find that as the bottom water 
$T_B\to 0$, the simulated values of 
$\delta$ and 
$F$ appear to diverge from the model. We recall the model equations (4.7a–c) are valid for 
$\sigma \gg 1$. As the decay rate 
$\sigma$ increases with 
${Ra}_0$, the larger values of 
${Ra}_0$ show better agreement between the model and the simulated data. Similarly, case 4 (
${Ra}_0 = 9\times 10^4$) exhibits the largest deviation from the model fit as it also is the lowest 
${Ra}_0$ simulation. As the bottom water temperature decreases (
$T_B\to 0$) 
$\sigma$ decreases rapidly. The first-order approximation (4.7a–c) is not expected to perform well in that limit. Higher-order corrections can account for this discrepancy, but that is outside the scope of this paper.
The flow transition that occurs at 
$T_B\approx 0.37$ (see (5.2) and (5.4)) results in a ‘kink’ in the modelled predictions, as seen in figure 7. We do not yet have a prediction for this limiting temperature value of 
$T_{max} = 0.37$. However, this transition is important in the evolution of 
$\delta _{St}$, 
$F$ and, to a lesser extent, 
$\hbox {TKE}$.
6. Conclusions
We considered a box of warm fluid, cooled from the surface through a fixed-temperature boundary condition. In the box, the density is quadratic with temperature. The top boundary temperature and the initial domain temperature are selected to be on opposite sides of the temperature of maximum density, leading to the generation of convection and the formation of an upper stable layer and a lower convectively unstable layer. As the convection mixes the lower-layer fluid, its near homogeneous temperature decreases (
$T_B\to 0$).
We developed a model for our system and scaling laws for 
$\delta _{St}$, 
$F$ and 
$\hbox {TKE}$ as follows,
\begin{gather} \delta_{St} \sim \frac{1}{\sigma T_B}, \quad F = \frac{1}{\delta_{St}} \sim \sigma T_B, \quad \sigma \gg 1, \end{gather}
\begin{gather} \langle \hbox{TKE}\rangle_V \sim \tfrac{1}{2}\left( \varGamma^{{-}1} - 1 \right) \left(1 - \varLambda \right) {Ra}_0 \,Pr \,T_B^2. \end{gather}
Analysing the numerical simulations, we determined the 
${Ra}_0$ and 
$T_B$ dependence of the decay 
$\sigma$ and the integrated-buoyancy flux ratio 
$\varLambda$. We showed that the mixing coefficient 
$\varGamma$ was effectively constant over the whole range of parameters considered here. Once determined, the model equations agreed well with the numerically computed 
$\delta _{St}$, 
$F$ and 
$\hbox {TKE}$.
The main goal of this paper was to understand how convection is changed when the EOS is nonlinear. We have shown that:
(i) the surface heat flux is dramatically different, in both magnitude and parameter dependence, between a nonlinear and linear EOS;
(ii) in addition to the Rayleigh and Prandtl number, convection with a quadratic EOS depends on a third non-dimensional parameter, the bottom water temperature
$T_B$;(iii) our model (6.1a–c) and (6.2) accurately predicts the heat flux (
$F$), boundary-layer thickness ($\delta _{St}$) and $\hbox {TKE}$ based on these three parameters.
This work forms a crucial step towards understanding how a nonlinear EOS modifies convection. The quadratic EOS limits the vertical heat flux and the kinetic energy, compared with a linear EOS, and depends on an additional non-dimensional parameter 
$T_B$. We are in the process of constructing a cold convection facility, capable of convectively cooling the surface of a fresh body of water, to run complementary laboratory experiments. This work provides a framework, including the essential parameters and model considerations, to understand much more complicated systems and, subsequently, freshwater systems in the environment.
Acknowledgements
We want to thank A. Wells, H. Ulloa and L.-A. Couston for their feedback on this work.
Funding
This work was funded in part by the Natural Sciences and Engineering Research Council of Canada, the Isaak Killam Trust and Syncrude Canada Ltd.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of the scaling law for $\delta _{St}$
First, assuming a piecewise-linear profile for 
$\bar {T}$, we can establish that
\begin{equation} \delta_{BL} = \left(1+T_B\right) \delta_{St} .\end{equation}The rate of change of the total temperature within the domain is then
\begin{align} \frac{\mathrm{d}}{\mathrm{d} t} \int_0^1 \bar{T}\,\hbox{d} z&= \frac{\mathrm{d}}{\mathrm{d} t} \left( \int_0^{1 - \delta_{BL}} T_B\hbox{d} z + \int_{1-\delta_{BL}}^1\left({-}1 - \frac{1 }{\delta_{St}} \left(z - 1 \right) \right) \hbox{d} z\right) \end{align}
\begin{align} &= \frac{\mathrm{d}T_B}{\mathrm{d} t} \left(1-\delta_{BL}\right) - \frac{1}{2} \frac{\delta_{BL}^2}{\delta_{St}^2} \frac{\partial \delta_{St}}{\partial t} . \end{align}Noting that the top temperature gradient prescribes the rate of change of heat within the domain, and including (A1), we arrive at
\begin{equation} \frac{\mathrm{d}T_B}{\mathrm{d} t} \left(1-\left(1 + T_B\right)\delta_{St}\right) - \frac{\left(1 + T_B\right)^2}{2} \frac{\partial \delta_{St}}{\partial t} ={-}\frac{1}{\delta_{St}}. \end{equation}If we further evaluate
\begin{equation} \frac{\textrm{d} T_B}{\textrm{d}t} \sim{-}\sigma T_B, \quad \delta_{St}\ll1, \quad \sigma\gg1, \end{equation}then to leading order
\begin{equation} \delta_{St} \sim \frac{1}{\sigma T_B}, \quad \sigma\gg1. \end{equation}Appendix B. Derivation of the scaling law for $\hbox {TKE}$
We first recall the Reynolds decomposition
\begin{equation} \boldsymbol{u} = \boldsymbol{0} + \boldsymbol{u}', \quad T = \bar{T} + T'. \end{equation}The density flux is then written out as
\begin{equation} \rho ={-} T^2 ={-} \left( \bar{T}^2 + 2 \bar{T} T' + T^\prime T^\prime\right) \implies \overline{w^\prime \rho^\prime} ={-} 2 \overline{w^\prime T^\prime} \ \bar{T} - \overline{w^\prime T^\prime T^\prime}. \end{equation}
For 
$z<1-\delta _{BL}$, the temperature stratification is well mixed, and thus
\begin{align} \left.\begin{gathered} \overline{T'w'} ={-}\frac{\mathrm{d}T_B}{\mathrm{d} t} z = \sigma T_B z, \quad z<1-\delta_{BL} \implies \langle 2 \overline{w^\prime T^\prime} \ \bar{T}\rangle_V \approx \int_0^{1} 2 \sigma T_B^2 z \,\hbox{d} z = \sigma T_B^2, \\ \delta_{BL}\ll1. \end{gathered}\right\} \end{align}
The approximations can be formalised by performing the full integrals and expanding the 
$\delta _{BL}$ as a perturbation series in 
${1}/{\sigma }$. Therefore, the TKE equation reduces to
\begin{equation} \frac{\mathrm{d}\langle \hbox{TKE}\rangle_V }{\mathrm{d} t} ={-}\left(\varGamma^{{-}1} - 1\right) \left( 1 - \lambda\right) \sigma T_B^2. \end{equation}The final step in deriving equation (4.11a,b) is to use self-similarity to determine
\begin{equation} \frac{\mathrm{d}\langle \hbox{TKE} ({Ra}_0,Pr,T_B)\rangle_V }{\mathrm{d} t} = \frac{\mathrm{d}\langle \hbox{TKE} \rangle_V }{\mathrm{d} T_B} \frac{\mathrm{d}T_B}{\mathrm{d} t} . \end{equation}
Substituting for 
$\varGamma$ and 
$\varLambda$ and integrating with respect to 
$T_B$, results in
\begin{equation} \langle \hbox{TKE}\rangle_V \sim \tfrac {1}{2}\left( \varGamma^{{-}1} - 1 \right) \left( 1 - \varLambda \right) {Ra}_0 \,Pr\, T_B^2, \quad \sigma \gg 1, \end{equation}
where 
$\varGamma$ is assumed constant and
\begin{equation} \varLambda = \frac{2}{T_B^2} \int_{T_B} T' \lambda(T') \,\textrm{d}T'. \end{equation}Appendix C. Linear stability analysis
The initial evolution of the temperature profiles described in § 3 is diffusive. For a deep box, the diffusive temperature solution is
\begin{equation} \bar{T} ={-}1 - \left(1+T_B\right) \hbox{erf} \left( \frac{z}{\delta}\right), \quad \delta = \sqrt{4 t},\end{equation}
where 
$T_B$ is fixed for the purposes of this stability analysis. Here, as we have a deep box, we redefine 
$z=0$ as the top boundary with the domain of interest below. In this analysis, we need to include diffusion of this background profile in the linear stability. We follow the approach of Nijjer, Hewitt & Neufeld (Reference Nijjer, Hewitt and Neufeld2018) and define the similarity variable
\begin{equation} \xi = \frac{z}{\sqrt{t}}, \end{equation}such that the background density profile is
\begin{equation} \bar{T} ={-}1 - \left(1+T_B\right) \hbox{erf} \left( \frac{\xi}{2}\right). \end{equation}
We want to know the growth rate of infinitesimal modal perturbations to the diffusive background state. The linear vertical velocity (
$w_\epsilon$) and temperature (
$T_\epsilon$) perturbations are assumed to have the form
\begin{equation} \begin{bmatrix} w_\epsilon \\ T_\epsilon \end{bmatrix} = \begin{bmatrix} \hat{w} \\ \hat{T} \end{bmatrix}\tau(t) \exp \left[ \textrm{i} \boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{x}\right]. \end{equation}Following the approach of Drazin & Reid (Reference Drazin and Reid2004), we linearise the equations of motion (2.4)–(2.6), which will result in a linear eigenvalue problem of the form
\begin{equation} \boldsymbol{\underline{A}} \begin{bmatrix} \hat{w}\\ \hat{T} \end{bmatrix} = \lambda \boldsymbol{\underline{B}} \begin{bmatrix} \hat{w}\\ \hat{T} \end{bmatrix}, \end{equation}
where 
$\boldsymbol {\underline {A}},\boldsymbol {\underline {B}}$ are matrix operators and 
$\lambda = ({1}/{\tau })({\textrm {d} \tau }/{\textrm {d}t} |_{t=t0})$ is the growth rate of the perturbations at some time 
$t0$. The solutions to these eigenvalue equations are a function of five parameters: 
$t0,k$, Ra, Pr and 
$T_B$. We solve these equations using an in-house built eigenvalue solver using Chebyshev differentiation matrices. We found that 50 grid points were sufficient to determine the growth rate of the system. We impose the same top boundary conditions as in the numerical simulations. Spurious solutions that did not satisfy the appropriate boundary conditions were removed.
Figure 8(a) is a plot of the growth rate (
$\lambda$) of the linear instability as a function of wavenumber 
$k$ and 
$t0$ for 
${Ra}=10^6$, 
$Pr=9$ and 
$T_B = 1$. We first note that there exists a minimum 
$t0_{min}$, below which the system is linearly stable. For 
$t>t0_{min}$, the system is unstable for a finite range of wavenumbers and a peak 
$\lambda$ at 
$t\gg t0_{min}$. Figure 8(b) is a plot of the maximum 
$\lambda$ over all wavenumbers, for each time 
$t0$ at different 
$T_B$. The maximum 
$\lambda$ decreases with 
$T_B$ until the system becomes stable at finite 
$T_B$.
(a) Contour plot of the growth rate (
$\lambda$) of the linear instability as a function of 
$t0$ and wavenumber (
$k$). (b) Plot of the maximum growth rate for all 
$k$ as a function of 
$t0$ for decreasing 
$T_B$. Both (a,b) are computed for 
${Ra}=10^6$ and 
$Pr=9$. (c) Plot of the minimum time to instability 
$t0_{min}$ below which the system is stable as a function of 
$T_B$. We include a vertical line at 
$T_B=0.35$ to show a similar regime change to the nonlinear dynamics. (d) Plot of the maximum growth rate for all 
$k$ and 
$\delta$ as a function of 
$T_B$. (e) Plot of the minimum 
$T_B$ below which the system is linearly stable versus 
${Ra}$.
By computing 
$t0_{min}$, we determine the earliest time at which the diffusive system becomes linearly unstable. Figure 8(c) is a plot of 
$t0_{min}$ as a function of 
$T_B$ for different Ra. Fitting the data, we find that
\begin{equation} t_0 = {Ra}^{{-2}/{3}}\begin{cases} 4 T_B^{{-}3/2}, & T_B<0.35\\ 3 T_B^{{-5}/{4}}, & T_B\geqslant0.35 \end{cases}\end{equation}
This fit is included as dashed/dotted lines in figure 8(c). As with the nonlinear simulations, we observe that there exists a critical regime change that occurs around 
$T_B\approx 0.35$. Further, there is a minimum 
$T_B$ [
$T_{B,min}$], where the system is linearly stable for all time and wavenumbers. Figure 8(d) is a plot of the maximum growth rate (
$\lambda _{max}$) for all 
$k$ and 
$t0$ as a function of 
$T_B$. For 
$T_B$ close to 1, the growth rate of the system follows
\begin{equation} \lambda_{max} \approx 0.4 {Ra}^{2/3} T_B^{4/3} = 0.4 {Ra}_{eff}^{1/3}. \end{equation}
Again, this fit is included as dashed lines in figure 8(d). As 
$T_B\to T_{B,min}$, 
$\lambda _{max}$ diverges from the fit (C7), decreasing rapidly to zero.
For a linear EOS, we know that there is a minimum Rayleigh number below which the system is stable. Similarly, we have found a minimum condition for instability with a quadratic EOS that depends on both 
${Ra}$ and 
$T_B$. Figure 8(e) is a plot of 
$T_{B,min}$ for different values of 
${Ra}$. Note that this panel has been rotated so that 
$T_{B,min}$ can be easily compared with panels (c,d). We estimate
\begin{equation} T_{B,min}\approx 0.5 {Ra}^{{-}0.44}. \end{equation}
Note that in terms of 
${Ra}_{eff}$, this suggests that the convection is stable where 
${Ra}_{eff}\lesssim \frac {1}{4}$. There exists a minimum stability condition below which the system is linearly stable for these diffusive temperature profiles for all time.
For 
${Ra}_{eff}$ larger than this minimum stability criterion, the analysis highlights that there still exists a 
$t0_{min}$ (or equivalently a minimum interface thickness), below which the system remains linearly stable. For 
$t>t0_{min}$, perturbations about the base state will grow and result in convective mixing of the temperature field. We can compare this minimum time with the numerical simulations.
Figure 9(a) is a plot of the volume-integrated 
$\hbox {TKE}$ density. As with the linear stability, the initial temperature stratification is stable, such that 
$\hbox {TKE}$ initially decays (see figure 9(a) inset). Once the top boundary layer is sufficiently deep, flow instability results in a large peak in the kinetic energy, which quickly decays back to an equilibrium value. From there, the 
$\hbox {TKE}$ slowly decays. We define the time when the kinetic energy reaches a minimum (vertical dashed line in figure 9(a), inset) as 
$t_{min}$. Note that this is only an approximation for the time when the system becomes linearly stable. Figure 9(b) is a plot of 
$t_{min}$ as a function of the Rayleigh number. A fit of the data to the expected power-law suggests that, for the numerical simulations,
\begin{equation} t_{min} \approx 4.4 {Ra}_0^{{-}2/3}, \end{equation}
which is close to that expected from linear theory from (C6) with 
$T_B=1$.
(a) Plot of the kinetic energy as a function of time for case 3 (
${Ra}_0=9.0\times 10^5$). The inset presents a subset of the data on a log-linear axis to show the initial TKE growth. The vertical dashed line indicates the time of minimum kinetic energy (
$t_{min}$). (b) Plot of 
$t_{min}$ as a function of Rayleigh number (cases 1–5). The scaling law (C9) is included as a dashed black line in (b).
Appendix D. Resolution and domain-size verification
D.1. Domain dependency
As mentioned in the text, we verified that the present results are not dependent on the box size of the numerical simulations. Figure 10 is similar to figure 5, with data from two different sized domains. We show that when we double the numerical domain, the temperature decay rate is nearly identical, despite the different initial random noise.
Similar to figure 5, we demonstrate that the decay rate of change is similar for different sized domains.
D.2. Resolution sufficiency
To demonstrate that the numerical simulations are indeed resolved, we have run an additional low-resolution simulation. We reran case 3 (
${Ra}_0 = 9\times 10^5$), but with a grid resolution of 
$\textrm {Nx} \times \textrm {Ny} \times \textrm {Nz} = 128 \times 128 \times 128$, or half the resolution in each direction. Figure 11 is a comparison plot between the full-resolution case 3 and the low-resolution run described here. We observe that the rate of change of temperature is essentially unchanged with one-eighth the number of grid points.
Similar to figure 5, we demonstrate that the simulation resolution is sufficient to resolve the rate of change of temperature.
