# Accounting for surface temperature variations in Rayleigh-Bénard convection

Jason Olsthoorn  Department of Civil Engineering, Queen’s University,
Kingston, Ontario, Canada, K7L 3N6<EMAIL_ADDRESS>
###### Abstract

Turbulent Rayleigh-Bénard convection is often modelled with a constant surface
temperature. However, the surface temperature of many geophysical systems,
such as lakes, is coupled to the atmospheric forcing. In this paper, we
account for this dynamic surface temperature through an additional parameter
$\beta$. Using an appropriately defined dynamical Rayleigh number
$\hbox{Ra}_{\Delta}$, we recover many of the results from the standard
Rayleigh-Bénard model. We hope that this work will simplify the application of
Rayleigh-Bénard theory in geophysical contexts, such as lakes.

##  I Introduction

In its usual configuration, Rayleigh-Bénard convection results from
sufficiently heating the bottom, and cooling the top, of a fluid. In their
experimental work [4], Bénard observed that the surface water deflections from
this convection formed hexagonal cells. In an attempt to model these cells,
Lord Rayleigh [14] then made the simplification that “the fluid is supposed to
be bounded by two infinite fixed planes […], where also the temperatures are
maintained constant.” Lord Rayleigh noted that this was a simplification from
“Bénard’s experiments, where, indeed, the [temperature] conditions are
different at the two boundaries.” In the years since that publication,
Rayleigh-Bénard convection has proven to be a complex dynamical system and
remains a topic of intense academic interest in fields ranging from
astrophysics to oceanography [7].

In many geophysical systems, such as lakes and oceans, the water surface can
be heated or cooled from the surrounding environment. For example, the surface
temperature of a lake depends on the atmospheric temperature, surface
radiation, and evaporation. During Autumn, the lake surface cools, which
drives convection that will transport heat from within the lake to the water
surface. That is, the surface temperature depends on both the atmospheric
forcing and the convection; it is not a constant as is typically considered in
Rayleigh-Bénard convection. However, we will show that this coupling can be
encapsulated by the Biot number, $\beta$, and that, with appropriately defined
dynamic parameters, the results of this modified setup are similar to those
found in the standard Rayleigh-Bénard model. Thus, we hope that this work will
simplify the application of Rayleigh-Bénard theory in geophysical contexts,
such as lakes.

Modifying the surface boundary conditions will change the linear stability of
the Rayleigh-Bénard model, which has been previously discussed in Sparrow et
al. [17] and subsequently in Foster [8]. The linear stability of the system is
determined by the Rayleigh number ($\hbox{Ra}_{0}$), which characterizes the
ratio of advective to diffusive transport. Sparrow et al. [17] demonstrated
that the critical Rayleigh number (RaC), the minimum Rayleigh number for
instability, monotonically increases as the upper boundary condition changed
from an insulating condition ($\beta=0$, Ra${}_{C}\approx 720$) to an
isothermal condition ($\beta\to\infty$, , Ra${}_{C}\approx 1800$). Similarly
changing the upper velocity boundary condition from free-slip to no-slip
resulted in a near uniform decrease in RaC ( a decrease of $\approx 100$, see
figure 1a in Sparrow et al. [17]). However, while the boundary conditions
affect RaC, Foster [8] demonstrated that the precise value of the thermal
boundary condition has a weak effect on the time to instability ($t_{0}$) for
the semi-infinite system. In this paper, we will focus on the nonlinear
behaviour of the system, after it becomes convectively unstable, both before
and at thermal equilibrium.

Once started, convection enhances the effective thermal conductivity between
the two boundaries in the Rayleigh-Bénard system, which increases with
$\hbox{Ra}_{0}$. As a result, Chillà et al. [5] argue that, at high
$\hbox{Ra}_{0}$, the finite thermal conductivity of the bounding plates used
in laboratory experiments will be unable to maintain a fixed temperature. That
is, the true boundary temperature in laboratory experiments depends upon the
convection once the effective thermal conductivity of the fluid increases
beyond a certain value. Subsequently, Verzicco [18] and Brown et al. [1]
demonstrated that an empirical correction factor that depends on the
conductivity ratio between the plates and the fluid can account for the
relative decrease in laboratory-measured heat transport at high
$\hbox{Ra}_{0}$. Interestingly, Wittenberg [21] argue that where the finite
conductivity of the plates is significant (as is the case when
$\hbox{Ra}_{0}\to\infty$), the appropriate $\hbox{Ra}_{0}$ is defined by the
temperature difference across the entire system, including the conductive
plates. We will define $\hbox{Ra}_{0}$ in a similar manner below. While the
methods and application of this paper are different from these laboratory
setups, we will similarly argue that the ratio $\beta/\hbox{Ra}_{0}^{p}$, ($p$
is a constant defined below) determines if the surface temperature is
effectively “fixed” at a constant temperature, or if the convection
significantly modifies the surface temperature.

More recently, Clarté et al. [6] performed a set of three-dimensional
convective simulations of a semi-insulated rotating convective shell with a
range of Biot numbers $\beta$. They argue that the dynamic surface boundary
condition can be replaced by a fixed flux condition when $\beta$ is
sufficiently small, or replaced by a fixed temperature condition when $\beta$
is sufficiently large. In this paper, we will extend their work by arguing
that the ratio $\beta/\hbox{Ra}_{0}^{p}$ is the important parameter that
determines the surface boundary value.

We present a model of convection between two fixed planes with a dynamic
surface boundary condition for temperature. This boundary condition will
assume that the far-field atmospheric conditions are fixed, with a dynamic
thermal boundary-layer near the water surface. We want to answer the following
key questions:

  1. 1.

What is the heat transfer rate at the water surface?

  2. 2.

What is the equilibrium surface water temperature?

  3. 3.

How quickly does the system reach equilibrium?

  4. 4.

How vigorous is the convection?

In the following discussion, we will describe the problem setup (§II), and the
theoretical framework to answer the questions above (§III). The results of the
2D numerical simulations will be discussed in §IV, before concluding in §V.

##  II Problem Setup

In analogy to a surface cooled lake, we model a body of water that is cooled
by exposure to the atmosphere through the top boundary. In this study, we will
fix the bottom temperature of the water to a known value
$\tilde{T}_{\tilde{z}=0}$. Further, we consider the case where the water is
shallow compared to its width, which we will model as periodic. This is a very
similar setup to the classic Rayleigh-Bénard problem with one major
modification: the surface temperature is coupled to the flux of heat through
the surface.

Figure 1: A contour plot of the simulated value of nondimensional temperature
$\Theta$ for $\hbox{Ra}_{0}=10^{7},\beta=8$ at time t=1500. Superimposed is a
cartoon of the mean temperature profile ($\overline{\Theta}$) between the
fixed tank bottom at $\Theta=1$ and the air temperature above.

Figure 1 is a schematic of the non-dimensionalized (see below) mean
temperature profile $\left(\overline{\Theta}\right)$ in this configuration.
The temperature is fixed at the bottom of the domain, which is warmer than the
atmosphere above. Due to the resultant convection, the temperature at the
centre of the water domain is well mixed, with temperature boundary layers at
the top and bottom boundary. The surface water temperature depends on how much
heat is being transported upwards by the convection.

###  II.1 Temperature Boundary Condition

The rate of change of the dimensional mean water temperature $\tilde{T}\
[\hbox{K}]$ in a volume $V\ [\hbox{m}^{3}]$ is determined by the heat flux
$\mathbf{q}\ [\hbox{W m}^{-2}]$ through the boundaries $A\ [\hbox{m}^{2}]$.
That is,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\tilde{t}}\int_{V}\rho_{0}c_{p}\tilde{T}\
\mathrm{d}V$ $\displaystyle=\int_{A}-\mathbf{q}\cdot\hat{\mathbf{n}}\
\mathrm{d}A^{\prime}.$ (1)

Here, $\hat{\mathbf{n}}$ is the outward normal vector, $\rho_{0}\ [\hbox{kg
m}^{-3}]$ is the reference density of water, and $c_{p}\ [\hbox{J
kg}^{-1}\hbox{K}^{-1}]$ is the specific heat capacity of water. As mentioned
above, we fix the lower boundary temperature ($\tilde{T}_{\tilde{z}=0}$). The
domain is horizontally periodic. The heat flux at the surface depends on the
driving environmental process.

In many natural systems, such as lakes and oceans, the surface heat loss is
given by the sum of sensible (conductive), radiative, and latent heat fluxes.
In this paper, we omit any additional heat inputs, such as precipitation, and
we will ignore short-wave radiation as it is not strictly a boundary effect.
The surface heat flux is then computed as the sum of net longwave radiation,
sensible heat flux, and evaporation, written:

$\displaystyle\mathbf{q}_{Surf}\cdot\mathbf{n}=-\epsilon_{w}\sigma_{SB}\tilde{T}_{\tilde{z}=\tilde{H}}^{4}+\epsilon\sigma_{SB}\tilde{T}_{Air}^{4}-K_{Air}\left(\tilde{T}_{\tilde{z}=\tilde{H}}-\tilde{T}_{Air}\right)+q_{L}.$
(2)

Here, $\sigma_{SB}\ [\hbox{W m}^{-2}\hbox{K}^{-4}]$ is the Stefan-Boltzmann
constant and $\tilde{T}_{Air}$ is the atmospheric temperature. The effective
air heat transfer coefficient $\left(K_{Air}\ [\hbox{W
m}^{-2}\hbox{K}^{-1}]\right)$, and the emissivity of the air ($\epsilon$ [-],
including reflection from the water surface) and water $\left(\epsilon_{w}\
[-]\right)$ are model parameters and need to be computed. However, these
parameters have been studied extensively in the literature (see Imboden and
Wüest [10] and elsewhere). Finally, the heat loss due to evaporation is
denoted as $q_{L}\ [\hbox{W m}^{-2}]$.

Hitchen and Wells [9] showed that in the absence of evaporation ($q_{L}=0$),
the surface heat flux ($\mathbf{q}_{Surf}$) can be linearized to the
following:

$\displaystyle\mathbf{q}_{Surf}\cdot\mathbf{n}=\gamma\left(\tilde{T}_{\tilde{z}=\tilde{H}}-\tilde{T}_{R}\right),$
(3)

where, in our notation,

$\displaystyle\gamma=\left(K_{Air}+4\epsilon_{w}\sigma_{SB}\tilde{T}_{Air}^{3}\right),\qquad\tilde{T}_{R}=\tilde{T}_{Air}+\frac{\sigma_{SB}\left(\epsilon-\epsilon_{w}\right)\tilde{T}_{Air}^{4}}{\left(K_{Air}+4\epsilon_{w}\sigma_{SB}\tilde{T}_{Air}^{3}\right)}.$
(4)

Here, $\gamma\ [\hbox{W m}^{-2}\hbox{K}^{-1}]$ is an effective heat transfer
coefficient and $\tilde{T}_{R}$ is an effective reference temperature. This
approximation assumes that the surface water temperature
$\tilde{T}_{\tilde{z}=\tilde{H}}$ and the atmospheric temperature
$\tilde{T}_{Air}$ are close
$(|\tilde{T}_{\tilde{z}=\tilde{H}}-\tilde{T}_{Air}|\ll\tilde{T}_{Air}$, in
Kelvin). This model can be similarly written to include the effects of
evaporation, which modify $\gamma$ and $\tilde{T}_{R}$, but not the functional
form of (3). We do not include them here as they increase the complexity of
the equations, but we refer the interested reader to Leppäranta [13]. For the
purposes of this paper, we will assume that $\gamma$ and $\tilde{T}_{R}$ are
given.

###  II.2 Equations of motion

The difference between the reference temperature $\tilde{T}_{R}$ and the
bottom temperature $\tilde{T}_{\tilde{z}=0}$ can be used to define a
characteristic velocity
$U_{0}^{2}={gH\alpha\left(\tilde{T}_{\tilde{z}=0}-\tilde{T}_{R}\right)}$,
where $\alpha$ is the assumed constant thermal expansion coefficient, $g$ is
the acceleration due to gravity, and $H$ is the water depth. We scale length
by $H$, velocity ($\mathbf{u}=(u,w)$) by $U_{0}$, pressure by
$\rho_{0}U_{0}^{2}$, and time by the advective timescale $\tau_{0}=H/U_{0}$.
Temperature is similarly scaled as

$\displaystyle\Theta=\frac{\tilde{T}-\tilde{T}_{R}}{\tilde{T}_{\tilde{z}=0}-\tilde{T}_{R}}.$
(5)

Under this non-dimensionalization, the equations of motion are

$\displaystyle\left(\frac{\partial}{\partial
t}+\mathbf{u}\cdot\nabla\right)\mathbf{u}=-\nabla
P+\Theta\mathbf{\hat{k}}+\sqrt{\frac{\Pr}{\hbox{Ra}_{0}}}\nabla^{2}\mathbf{u},$
(6) $\displaystyle\left(\frac{\partial}{\partial
t}+\mathbf{u}\cdot\nabla\right)\Theta=\frac{1}{\sqrt{\hbox{Ra}_{0}\
\Pr}}\nabla^{2}\Theta,$ (7) $\displaystyle\nabla\cdot\mathbf{u}=0.$ (8)

These equations contain two non-dimensional parameters: the Rayleigh number
($\hbox{Ra}_{0}$) and the Prandtl number (Pr).

$\displaystyle\hbox{Ra}_{0}=g\alpha\left(\tilde{T}_{\tilde{z}=0}-\tilde{T}_{R}\right)\frac{H^{3}}{\kappa\nu},\qquad\hbox{Pr}=\frac{\nu}{\kappa}.$
(9)

Here, $\nu$ is the kinematic viscosity and $\kappa$ is the thermal
diffusivity.

We model the top and bottom boundaries as impermeable that are free-slip (no
tangential stress) at the surface and no-slip at the bottom. The boundary
conditions are then written

$\displaystyle\frac{\partial u}{\partial z}=0,\quad$ $\displaystyle w=0,$
$\displaystyle\frac{\partial\Theta}{\partial z}$
$\displaystyle=-\beta\Theta,\qquad\hbox{ where }\beta=\frac{\gamma
H}{\rho_{0}c_{p}\kappa},$ $\displaystyle z=1,$ (10) $\displaystyle u=0,\quad$
$\displaystyle w=0,$ $\displaystyle\Theta$ $\displaystyle=1,$ $\displaystyle
z=0.$ (11)

For simplicity, the initial condition was set to $\Theta(z)=1$. Note that the
numerical solver instantaneously corrects the initial condition to match the
boundary condition. We perform a series of numerical simulations with
different values of $\hbox{Ra}_{0}$ and $\beta$. In all cases, we fix $\Pr=9$
as the approximate value for heat in water ($\Pr\approx 9.0$ at 11.5∘C and
atmospheric pressure at sea level).

###  II.3 Numerical Methods

We solve the system of equations (6)–(8) with Dedalus [2], using
pseudospectral spatial derivatives (Chebyshev polynomials in the vertical and
Fourier modes in the horizontal) and a second-order Runge-Kutta time-stepping
scheme. Numerical convergence was verified with grid resolution studies and by
ensuring that there were at least 8 grid points within the top boundary layer.

In total, 44 numerical simulations were run (see table 1). The simulations
were two-dimensional and run for a sufficiently long time that the system
reached a quasi-steady-state (final time was selected with $t/\tau\gtrsim 4$
as defined below). The domain width was 4 times the depth in all cases.

$\hbox{Ra}_{0}$ | $\beta$ | Nx | Nz | $\delta_{top}$ | $N_{B}$ | $t_{0}$ | $\epsilon$ | $t_{max}$ | $t_{max}/\tau$ |  Nu0 | $\Theta_{1}$ | $\sigma(\Theta_{1})$ | $\mathcal{K}_{0}$ | Re
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$10^{6}$ | $2^{0}$ | 256 | 128 | 0.034 | 15 | 15.7 | 0.058 | 1700 | 4.8 | 6.3 | 0.9 | 0.02 | 1.1$\times 10^{-3}$ | 15
$10^{6}$ | $2^{1}$ | 256 | 128 | 0.030 | 14 | 11.7 | 0.115 | 1300 | 4.4 | 7.4 | 0.8 | 0.03 | 1.7$\times 10^{-3}$ | 19
$10^{6}$ | $2^{2}$ | 256 | 128 | 0.030 | 14 | 9.3 | 0.230 | 1100 | 4.4 | 8.4 | 0.7 | 0.04 | 2.8$\times 10^{-3}$ | 24
$10^{6}$ | $2^{3}$ | 256 | 128 | 0.028 | 14 | 7.1 | 0.460 | 1300 | 6.2 | 8.9 | 0.5 | 0.05 | 4.3$\times 10^{-3}$ | 30
$10^{6}$ | $2^{4}$ | 256 | 128 | 0.025 | 13 | 6.5 | 0.920 | 1200 | 6.8 | 9.4 | 0.4 | 0.05 | 6.1$\times 10^{-3}$ | 36
$10^{6}$ | $2^{5}$ | 256 | 128 | 0.026 | 13 | 4.5 | 1.841 | 1100 | 7.4 | 10.3 | 0.2 | 0.05 | 7.8$\times 10^{-3}$ | 41
$10^{6}$ | $2^{6}$ | 256 | 128 | 0.026 | 13 | 3.7 | 3.681 | 1100 | 8.8 | 10.2 | 0.1 | 0.04 | 9.1$\times 10^{-3}$ | 44
$10^{6}$ | $2^{7}$ | 256 | 128 | 0.024 | 13 | 3.6 | 7.362 | 1000 | 9.7 | 10.3 | 0.1 | 0.03 | 9.7$\times 10^{-3}$ | 46
$10^{6}$ | $2^{8}$ | 256 | 128 | 0.026 | 13 | 3.4 | 14.724 | 1000 | 9.7 | 10.3 | 0.0 | 0.02 | 10.0$\times 10^{-3}$ | 47
$10^{6}$ | $2^{9}$ | 256 | 128 | 0.026 | 13 | 3.4 | 29.448 | 1000 | 9.7 | 10.3 | 0.0 | 0.01 | 1.0$\times 10^{-2}$ | 48
$10^{6}$ | Fixed | 256 | 128 | 0.025 | 13 | 2.7 | $\infty$ | 1000 | 9.7 | 10.5 | 0.0 | 0.00 | 1.1$\times 10^{-2}$ | 48
$10^{7}$ | $2^{0}$ | 256 | 128 | 0.025 | 13 | 18.4 | 0.028 | 3600 | 5.5 | 10.4 | 0.9 | 0.01 | 9.4$\times 10^{-4}$ | 45
$10^{7}$ | $2^{1}$ | 256 | 128 | 0.022 | 12 | 15.0 | 0.056 | 3000 | 5.5 | 11.9 | 0.9 | 0.02 | 1.7$\times 10^{-3}$ | 61
$10^{7}$ | $2^{2}$ | 256 | 128 | 0.019 | 11 | 10.8 | 0.113 | 2400 | 5.2 | 13.6 | 0.8 | 0.03 | 3.0$\times 10^{-3}$ | 81
$10^{7}$ | $2^{3}$ | 256 | 128 | 0.016 | 10 | 7.3 | 0.225 | 2000 | 5.2 | 15.3 | 0.7 | 0.04 | 4.8$\times 10^{-3}$ | 103
$10^{7}$ | $2^{4}$ | 256 | 128 | 0.016 | 10 | 6.2 | 0.451 | 1600 | 4.9 | 16.8 | 0.5 | 0.05 | 7.8$\times 10^{-3}$ | 131
$10^{7}$ | $2^{5}$ | 256 | 128 | 0.016 | 10 | 4.7 | 0.901 | 1600 | 5.8 | 17.8 | 0.4 | 0.06 | 1.1$\times 10^{-2}$ | 156
$10^{7}$ | $2^{6}$ | 512 | 128 | 0.014 | 10 | 3.7 | 1.803 | 1500 | 6.5 | 18.5 | 0.2 | 0.06 | 1.3$\times 10^{-2}$ | 173
$10^{7}$ | $2^{7}$ | 512 | 128 | 0.014 | 10 | 3.2 | 3.606 | 1450 | 7.5 | 19.0 | 0.1 | 0.04 | 1.6$\times 10^{-2}$ | 186
$10^{7}$ | $2^{8}$ | 512 | 128 | 0.013 | 9 | 2.9 | 7.212 | 1400 | 8.7 | 19.3 | 0.1 | 0.03 | 1.7$\times 10^{-2}$ | 193
$10^{7}$ | $2^{9}$ | 512 | 128 | 0.014 | 9 | 2.7 | 14.423 | 1450 | 9.0 | 19.2 | 0.0 | 0.02 | 1.7$\times 10^{-2}$ | 196
$10^{7}$ | Fixed | 256 | 128 | 0.013 | 9 | 2.4 | $\infty$ | 1450 | 9.0 | 19.2 | 0.0 | 0.00 | 1.8$\times 10^{-2}$ | 199
$10^{8}$ | $2^{0}$ | 256 | 128 | 0.016 | 10 | 22.9 | 0.014 | 7000 | 5.8 | 17.3 | 0.9 | 0.01 | 9.1$\times 10^{-4}$ | 142
$10^{8}$ | $2^{1}$ | 256 | 128 | 0.013 | 9 | 15.9 | 0.028 | 5400 | 5.3 | 19.9 | 0.9 | 0.01 | 1.7$\times 10^{-3}$ | 195
$10^{8}$ | $2^{2}$ | 512 | 128 | 0.009 | 8 | 11.9 | 0.055 | 4800 | 5.6 | 22.4 | 0.8 | 0.02 | 2.9$\times 10^{-3}$ | 254
$10^{8}$ | $2^{3}$ | 512 | 256 | 0.006 | 13 | 7.9 | 0.110 | 4100 | 5.7 | 28.6 | 0.8 | 0.03 | 4.9$\times 10^{-3}$ | 329
$10^{8}$ | $2^{4}$ | 512 | 256 | 0.006 | 13 | 6.9 | 0.221 | 3300 | 5.5 | 32.2 | 0.7 | 0.04 | 8.1$\times 10^{-3}$ | 423
$10^{8}$ | $2^{5}$ | 512 | 256 | 0.006 | 12 | 4.9 | 0.442 | 2700 | 5.3 | 35.7 | 0.5 | 0.06 | 1.2$\times 10^{-2}$ | 521
$10^{8}$ | $2^{6}$ | 512 | 256 | 0.006 | 12 | 3.9 | 0.883 | 2300 | 5.4 | 38.7 | 0.4 | 0.06 | 1.7$\times 10^{-2}$ | 610
$10^{8}$ | $2^{7}$ | 512 | 256 | 0.006 | 12 | 2.9 | 1.766 | 2300 | 6.4 | 40.1 | 0.2 | 0.06 | 2.1$\times 10^{-2}$ | 682
$10^{8}$ | $2^{8}$ | 512 | 256 | 0.005 | 12 | 2.9 | 3.532 | 2100 | 7.0 | 41.5 | 0.1 | 0.05 | 2.4$\times 10^{-2}$ | 736
$10^{8}$ | $2^{9}$ | 512 | 256 | 0.006 | 12 | 2.9 | 7.064 | 2000 | 8.1 | 42.1 | 0.1 | 0.03 | 2.6$\times 10^{-2}$ | 765
$10^{8}$ | Fixed | 512 | 256 | 0.006 | 13 | 1.9 | $\infty$ | 2000 | 8.1 | 42.4 | 0.0 | 0.00 | 2.8$\times 10^{-2}$ | 795
$10^{9}$ | $2^{0}$ | 512 | 256 | 0.006 | 13 | 23.9 | 0.007 | 14500 | 6.5 | 27.2 | 1.0 | 0.00 | 7.3$\times 10^{-4}$ | 403
$10^{9}$ | $2^{1}$ | 512 | 256 | 0.005 | 11 | 16.9 | 0.014 | 11900 | 6.3 | 36.5 | 0.9 | 0.01 | 1.4$\times 10^{-3}$ | 550
$10^{9}$ | $2^{2}$ | 512 | 256 | 0.004 | 11 | 12.9 | 0.027 | 9300 | 5.9 | 42.4 | 0.9 | 0.01 | 2.5$\times 10^{-3}$ | 744
$10^{9}$ | $2^{3}$ | 1024 | 256 | 0.004 | 11 | 9.9 | 0.054 | 7000 | 5.3 | 48.4 | 0.9 | 0.02 | 4.3$\times 10^{-3}$ | 980
$10^{9}$ | $2^{4}$ | 1024 | 256 | 0.004 | 10 | 6.9 | 0.108 | 6300 | 5.6 | 54.2 | 0.8 | 0.03 | 7.3$\times 10^{-3}$ | 1273
$10^{9}$ | $2^{5}$ | 1024 | 256 | 0.004 | 10 | 4.9 | 0.216 | 5000 | 5.3 | 60.0 | 0.7 | 0.04 | 1.1$\times 10^{-2}$ | 1579
$10^{9}$ | $2^{6}$ | 1024 | 256 | 0.004 | 10 | 3.9 | 0.432 | 3900 | 4.9 | 65.9 | 0.5 | 0.05 | 1.6$\times 10^{-2}$ | 1897
$10^{9}$ | $2^{7}$ | 1024 | 256 | 0.004 | 10 | 2.9 | 0.865 | 3000 | 4.5 | 71.0 | 0.4 | 0.05 | 2.2$\times 10^{-2}$ | 2189
$10^{9}$ | $2^{8}$ | 1024 | 256 | 0.004 | 10 | 2.9 | 1.730 | 2400 | 4.3 | 74.9 | 0.2 | 0.05 | 2.7$\times 10^{-2}$ | 2441
$10^{9}$ | $2^{9}$ | 1024 | 256 | 0.004 | 10 | 1.9 | 3.460 | 1900 | 4.0 | 77.6 | 0.1 | 0.04 | 3.1$\times 10^{-2}$ | 2629
$10^{9}$ | Fixed | 1024 | 256 | 0.004 | 10 | 1.3 | $\infty$ | 1900 | 4.9 | 76.4 | 0.0 | 0.00 | 3.8$\times 10^{-2}$ | 2915
Table 1: A table of the parameters associated with each experiment. The
horizontal and vertical grid resolution was $Nx,Nz$, respectively. The upper
thermal boundary layer thickness was $\delta_{top}$, with $N_{B}$ grid points
resolving the boundary layer. The initial time of instability is $t_{0}$ as
defined as the time of minimum kinetic energy $\mathcal{K}$. The end time of
each numerical run was $t_{max}$, which is at least $\geq 4\tau$, where
$\tau=\tau_{\epsilon}$ is used for $\epsilon<5$ and $\tau=\tau_{\infty}$
otherwise (both defined below). The estimated steady state values of
$\hbox{Nu}_{0}$, $\Theta_{1}$ and $\mathcal{K}_{0}$ are provided. The standard
deviation of the surface temperature $\sigma(\Theta_{1})$ is also provided.
Finally, the Reynolds number (Re) based upon the RMS velocity is given. The
‘Fixed’ cases have prescribed $\Theta_{z=1}=0$, i.e. the limit of
$\beta\to\infty$. The horizontal dimension Lx$=4$ in all cases.

##  III Theory

The Nusselt number (Nu) is the ratio between the measured surface vertical
heat flux ($J$) and the purely diffusive heat flux. In our non-
dimensionalization, it is computed as

$\displaystyle\hbox{Nu}=\frac{J}{\Delta},$ (12)

The temperature difference, $\Delta=\Theta_{z=0}-\Theta_{z=1}$, between the
upper and lower boundaries is typically prescribed for classical Rayleigh-
Bénard convection. For a fixed Prandtl number, we will show below that Nu
scales with the effective Rayleigh number $\hbox{Ra}_{\Delta}$,

$\displaystyle\hbox{Nu}-1\sim
C\hbox{Ra}_{\Delta}^{p},\qquad\hbox{Ra}_{\Delta}=\hbox{Ra}_{0}\Delta,$ (13)

where the effective Rayleigh number, $\hbox{Ra}_{\Delta}$, is defined by the
dynamic temperature values at the top and bottom boundary. While debate
persists concerning the ‘ultimate’ limit of the value of $p$ [16], it is often
shown to have a value of $p\approx 0.3$ at moderate Rayleigh numbers (e.g.
Niemela and Sreenivasan [15] estimate that $p\approx\frac{1}{3}$ for
$\hbox{Ra}\approx 10^{10}-10^{12}$, Plumley and Julien [16] estimate that
$p\approx 0.322$ for $\hbox{Ra}\approx 10^{5}-10^{15}$, Chillà et al. [5]
estimate $p\approx 0.3$ for $\hbox{Ra}\approx 10^{9}-10^{12}$ ).

In steady state, the horizontally averaged
($\overline{\left(\cdot\right)}=\frac{1}{A}\int_{\partial
V}\left(\cdot\right)\ \mathrm{d}A^{\prime}$) temperature equation (7)
simplifies to

$\displaystyle\frac{\partial}{\partial
z}{}\left(\frac{1}{\sqrt{\hbox{Ra}_{0}\Pr}}\frac{\partial}{\partial
z}{\overline{\Theta}}-\overline{w\Theta}\right)=0,$ (14)

where $w$ is the vertical velocity. Integrating the above equation, we show
that the Nusselt number can be equivalently determined by the surface heat
flux, or the volume averaged vertical heat flux as

$\displaystyle\hbox{Nu}=\frac{1}{\Delta}\frac{\partial}{\partial
z}{\overline{\Theta}}\bigg{|}_{z=1}=\frac{\sqrt{\hbox{Ra}_{0}\
\Pr}}{\Delta}\langle w\Theta\rangle+1,$ (15)

where $\langle\left(\cdot\right)\rangle=\frac{1}{V}\int_{V}\left(\cdot\right)\
\mathrm{d}V$). In the results presented below, we will numerically evaluate
the Nusselt number based upon the surface gradient. We will use these
relationships below.

###  III.1 Diffusive solution

The simplest solution to equations (6)-(8) with boundary conditions (10) is
the diffusive (non-convecting) solution. That is, if we make the ansatz that

$\displaystyle\Theta=az+1\implies a=-\beta\left(a+1\right),$ (16)
$\displaystyle\Theta=-\frac{\beta}{\beta+1}z+1.$ (17)

Therefore, for the purely diffusive problem,

$\displaystyle\Theta_{z=1}=-\frac{\beta}{\beta+1}+1\approx\frac{1}{\beta}+\dots,\qquad\beta\to\infty$
(18)

The surface temperature $\Theta_{z=1}$ decreases with increasing $\beta$. For
large $\hbox{Ra}_{0}$, convection increases the vertical heat flux, resulting
in a nonlinear temperature profile. Nonetheless, this diffusive solution
highlights the approximate functional dependence of the surface temperature on
$\beta$.

##  IV Results

Figure 2: (a) The time series of the mean top and mean domain temperatures.
The red dashed line is an exponential fit to the data. Note that these
temperatures asymptote to statistically stationary values. (b) Mean
temperature profiles at the times identified by vertical lines in (a). The
data in (a) and (b) are from $\hbox{Ra}_{0}=10^{7},\ \beta=2^{3}$. (c) The
Nusselt number as a function of $\hbox{Ra}_{\Delta}$. The lines denote the
time evolution of $\hbox{Nu}-1$ over the convective period of the simulation,
while the squares denote the equilibrium values. (d) The asymptotic values of
$\Theta_{1}$ as a function of $\hbox{Ra}_{\Delta}$ and $\beta$. The standard
deviation of the surface temperature at the end of the simulation are included
as vertical grey bars. The data agrees well with (19) (solid black line).

In all of the simulations presented here, $\hbox{Ra}_{0},\beta$ are
sufficiently large that the system will become convectively unstable. This
convection increases the vertical transport of heat. After an initial
transient, both the average top temperature ($\Theta_{z=1}$) and mean domain
temperature ($\langle\Theta\rangle$) decay nearly exponentially to their
equilibrium values (see figure 2a). For the remainder of this paper, we will
denote the stationary value of the surface temperature as $\Theta_{1}$, which
was computed from the exponential fit (see the horizontal (black) dashed line
in figure 2a).

After convection begins, the mean temperature profiles have a consistent
structure; the centre of the domain is well mixed with sharp temperature
gradients near the boundaries (see figure 2b). The feedback between the
convective heat flux and the boundary condition results in a gradual, not
instantaneous, decrease in $\overline{\Theta}_{z=1}\to\Theta_{1}$. The upper
and lower temperature gradients increase accordingly. Note that the upper
temperature boundary layer is typically thinner than the bottom boundary layer
due to the different boundary conditions, leading to a mean temperature that
is less than the mean of the boundary values $\langle\Theta\rangle<\Delta/2$.

###  IV.1  What is the heat transfer rate at the water surface?

As the surface temperature changes, so do the top and bottom temperature
gradients, and, correspondingly, Nu. At the onset of convection, Nu is higher
than its asymptotic value ($\hbox{Nu}_{0}=\lim_{t\to\infty}\hbox{Nu}$). That
is, the value of Nu decreases over time, with $\hbox{Nu}_{0}$ collapsing onto
the curve (13), with $C=0.138\pm 0.007,\ p=0.308\pm 0.003$ (see figure 2c).
For the remainder of this paper, we will select $C=0.14,p=0.31$ as the power-
law coefficients for Nu.

We note that the observed power-law is less than the optimal value for free-
slip boundaries [$p\approx\frac{1}{3}$, see 20], but greater than the
simulated values with no-slip boundaries and fixed aspect ratio [$p\approx
0.28$, see 19]. The results of Clarté et al. [6] suggest that, in the absence
of rotation, $p\approx 0.29$. The current simulations with a no-slip bottom
boundary and a free-slip surface are within the range of previously reported
values.

###  IV.2 What is the equilibrium surface water temperature?

From the Nu scaling (13) and the boundary condition (10), we can write down an
implicit equation for the top boundary temperature in the large
$\hbox{Ra}_{\Delta}$ limit,

$\displaystyle\frac{C\hbox{Ra}_{0}^{p}}{\beta}\sim\frac{\Theta_{1}}{\left(1-\Theta_{1}\right)^{1+p}},\qquad\hbox{Ra}_{\Delta}\gg
1.$ (19)

Noting that, here, the mean heat flux $J\sim C\hbox{Ra}_{0}^{p}\Delta^{1+p}$.
In the limit, $J\to 1$, we recover the diffusive solution (18).

We find good agreement between the measured values of $\Theta_{1}$ and (19)
with the estimated Nu coefficients $C,p$ (See figure 2d). That is, the ratio
$\beta/\hbox{Ra}_{0}^{p}$ is the fundamental parameter that determines the
equilibrium surface temperature. As discussed in the introduction, it is only
reasonable to assume a fixed surface temperature when
$\beta\gg\hbox{Ra}_{0}^{p}$.

For finite $\beta$, the surface temperature fluctuates. The standard deviation
of surface temperature at the end of the simulations (averaged in space and
time) are denoted by the grey errorbars in figure 2d and are provided in table
1. For both large and small values of $C\hbox{Ra}_{0}^{p}/\beta$, the surface
temperature is constrained and the temperature variations are minimal. Larger
variations of $\approx 0.06$ are consistently found when
$C\hbox{Ra}_{0}^{p}/\beta=O(1)$, and are primarily produced by the emergent
convective circulation cells (see figure 1). As the form of the surface
temperature is different in three-dimensions [3], we do not pursue this
further and will investigate the structure of the surface temperature
distribution in future work.

###  IV.3 How quickly does the system reach equilibrium?

The convection induced by surface cooling results in a rapid decrease in the
mean water temperature (figure 2a), which eventually reaches a steady state.
How long does it take to reach steady state? Figure 3(a) is a plot of the
temperature perturbation $\left(1-\langle\Theta\rangle\right)$ normalized by
its asymptotic value as a function of time for all cases. We find a clear
difference in the time to equilibrium between the different cases.

As mentioned previously, the internal water temperature is nearly uniform,
with strong gradients at the boundaries (figure 2b). Modelling this
temperature profile as piecewise linear, the mean heat equation (1) and
boundary condition (10) simplify to

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\langle\Theta\rangle=\frac{1}{\sqrt{\hbox{Ra}_{0}\
\Pr}}\left[\left(\frac{\Theta_{z=1}-\langle\Theta\rangle}{\delta_{top}}\right)-\left(\frac{\langle\Theta\rangle-\Theta_{z=0}}{\delta_{bottom}}\right)\right],\qquad\left(\frac{\Theta_{z=1}-\langle\Theta\rangle}{\delta_{top}}\right)=-\beta\Theta_{z=1}.$
(20)

Here, $\delta_{top},\delta_{bottom}$ are the top and bottom boundary layer
thicknesses, and $\Theta_{z=0}=~{}1$. We will attempt to simplify the equation
for $\langle\Theta\rangle$ in the large Rayleigh number limit.

We begin by assuming that each boundary layer separately satisfies a Nusselt-
Rayleigh number relationship. That is, the upper and lower boundary layers
satisfy,

$\displaystyle Nu_{u}$
$\displaystyle=\frac{1}{\langle\Theta\rangle-\Theta_{z=1}}\frac{\langle\Theta\rangle-\Theta_{z=1}}{\delta_{top}}\sim
C_{1}\left(Ra_{0}\left(\langle\Theta\rangle-\Theta_{z=1}\right)\right)^{p},$
(21) $\displaystyle Nu_{l}$
$\displaystyle=\frac{1}{\Theta_{z=0}-\langle\Theta\rangle}\frac{\Theta_{z=0}-\langle\Theta\rangle}{\delta_{bottom}}\sim
C_{2}\left(Ra_{0}\left(\Theta_{z=0}-\langle\Theta\rangle\right)\right)^{p},\qquad\hbox{Ra}_{0}\to\infty.$
(22)

As discussed in the introduction, we expect $C_{1}\neq C_{2}$ due to the
different upper and lower boundary conditions. Rearranging, we can simplify
equation (20) above to

$\displaystyle\phi\frac{\mathrm{d}}{\mathrm{d}t}\langle\Theta\rangle=\left(1-\langle\Theta\rangle\right)^{1+p}-\epsilon\Theta_{z=1},\qquad\left(\langle\Theta\rangle-\Theta_{z=1}\right)^{1+p}=\frac{\epsilon}{\zeta}\Theta_{z=1},$
(23)

where we have defined the parameters

$\displaystyle\phi=\frac{\sqrt{\hbox{Ra}_{0}\
\Pr}}{C_{2}\hbox{Ra}_{0}^{p}},\qquad\epsilon=\frac{\beta}{C_{2}\hbox{Ra}_{0}^{p}},\qquad\zeta=\frac{C_{1}}{C_{2}}.$
(24)

We consider two extreme limits for $\epsilon$. As discussed above, in the
$\epsilon=0$ limit, the surface boundary condition tends to an insulating
boundary condition. Conversely, the $\epsilon\to\infty$ limit fixes the
surface temperature, which is more similar to the classic Rayleigh-Bénard
problem. Investigating these two limits, we will estimate the timescale to
equilibrium that depends on $\phi$, $\epsilon$ and $\zeta$.

####  IV.3.1 Small $\epsilon\ll 1$

As we increase $\hbox{Ra}_{0}$ with fixed $\beta$, the parameter $\epsilon$
tends to 0. In fact, most of the simulations performed in this paper have
$\epsilon<1$. However, in the limit $\epsilon=0$, the insulating boundary
condition results in no temperature difference across the layer, and the
Nusselt number relationships (21)–(22) are invalid. In this analysis, we
assume that $\epsilon$ is small but sufficiently large that the system remains
convective.

We look for a perturbation series solution for $\langle\Theta\rangle$ and
$\Theta_{z=1}$. Scaling through, we determine a solution in small $\epsilon$
as

$\displaystyle\langle\Theta\rangle=1-\epsilon^{\frac{1}{1+p}}\left(x_{\epsilon,1}+\epsilon
x_{\epsilon,2}+\dots\right),\qquad\Theta_{z=1}=1-\epsilon^{\frac{1}{1+p}}\left(y_{\epsilon,1}+\epsilon
y_{\epsilon,2}+\dots\right),\qquad\epsilon\ll 1.$ (25)

For consistency, we further scale time as

$\displaystyle t=\phi\epsilon^{\frac{-p}{1+p}}\xi_{\epsilon}.$ (26)

At first order, (23) reduces for $x_{\epsilon,1}$ to,

$\displaystyle\frac{\mathrm{d}x_{\epsilon,1}}{\mathrm{d}\xi_{\epsilon}}=1-x_{\epsilon,1}^{1+p}.$
(27)

We plot the numerical solution to (27) in figure 3(b) along with a subset of
the numerical simulations (indicated by the squares in figure 3(c)). Here, we
select $C_{2}=0.24$ as an empirical parameter, without error bounds due to the
low number of data points. Note that the choice of scale for $\epsilon$
determines both the amplitude and time scale for $x_{\epsilon,1}$. For the
cases selected, the reduced model (27) provides a reasonable approximation.
Figure 3(c) plots the asymptotic value of $x_{\epsilon,1}$ for each of the
numerical simulations. While each different $\hbox{Ra}_{0}$ cases appear to
plateau for $\epsilon\ll 1$, the asymptotic values of $x_{\epsilon,1}$ are
significantly lower for $\hbox{Ra}_{0}=10^{6},10^{7}$, than for
$\hbox{Ra}_{0}=10^{8},10^{9}$. We suggest that this change in amplitude may
result from the corresponding Nu being closer to 1 in those cases and
therefore (21)–(22) are perturbed. Nevertheless, (26) provides a reasonable
estimate for the equilibrium timescale for $\epsilon\ll 1$ at large
$\hbox{Ra}_{0}$.

####  IV.3.2 Large $\epsilon\to\infty$

In the alternative limit, $\epsilon\to\infty$, the surface boundary condition
is isothermal at $\Theta_{z=1}=0$. We can then construct an asymptotic series
for $\langle\Theta\rangle$, $\Theta_{z=1}$ as

$\displaystyle\langle\Theta\rangle=1-\left(x_{\infty,0}+\frac{1}{\epsilon}x_{\infty,1}+\dots\right),\qquad\Theta_{z=1}=0+\frac{1}{\epsilon}y_{\infty,1}+\dots,\qquad\epsilon\to\infty.$
(28)

Again, we use the timescale

$\displaystyle t=\frac{\phi}{\zeta}\xi_{\infty}.$ (29)

Scaling time by $\zeta$ defines a timescale that is more consistent with (26)
as shown in figure 3 (i.e. more similar to an e-folding time). To leading
order, the heat budget (23) reduces to the simplified equation

$\displaystyle\frac{\mathrm{d}x_{\infty,0}}{\mathrm{d}\xi_{\infty}}=-\frac{1}{\zeta}\left(x_{\infty,0}\right)^{1+p}+\left(1-x_{\infty,0}\right)^{1+p}.$
(30)

In this large $\epsilon$ limit, the timescale and asymptotic value of
$x_{\infty,0}$ are determined by $\zeta$. We estimate $1/\zeta\approx 0.6$.
The solution to (30) is plotted in figure 3(d) along with a subset of the
numerical simulations (indicated by the squares in figure 3(e)). Figure 3(d)
also includes the cases with $\Theta_{z=1}$ fixed at 0. The equilibirum
perturbation temperature $(1-\langle\Theta\rangle_{0}$) decreases rapidly for
$\epsilon<5$ as shown in figure 3(e). Plotting the $\epsilon>5$ cases, the
solution to (30) agrees well with the numerical simulations in both amplitude
and timescale.

####  IV.3.3 Matching fits

In developing the asymptotic solutions, we have defined two timescales
($\tau$) that were given by

$\displaystyle\tau=\begin{cases}\tau_{\epsilon}=\phi\epsilon^{\frac{-p}{1+p}},&\epsilon\ll
1,\\\ \tau_{\infty}=\phi/\zeta,&\epsilon\to\infty,\\\ \end{cases}\qquad\hbox{
where }\epsilon=\frac{\beta}{C_{2}\hbox{Ra}_{0}^{p}}\quad\hbox{ and
}\phi=\frac{\sqrt{\hbox{Ra}_{0}\ \Pr}}{C_{2}\hbox{Ra}_{0}^{p}}.$ (31)

We suggest that we may be able to match these two solutions as
$\tau=\phi\left(\epsilon^{\frac{-p}{1+p}}+f(\epsilon)\frac{1}{\zeta}\right)$,
where $f(\epsilon)\to 0$ for $\epsilon\ll 0$, and $f(\epsilon)\to 1$ as
$\epsilon\to\infty$. Empirically, it appears that
$\tau=\phi\epsilon^{\frac{-p}{1+p}}$ provides a reasonable estimate for the
timescale in all simulated cases at finite $\epsilon$. Figure 3(f) is a plot
of the normalized temperature perturbation for all of the numerical
simulations with finite $\epsilon$. Further, we find that the asymptotic
solution where $\epsilon\ll 1$ ( i.e. (27)) provides a reasonable
approximation for the shape of the mean temperature timeseries when normalized
by their asymptotic values. Note that, contrary to the original derivation,
this includes the cases where $\epsilon>1$ ( and finite). Based upon this
surprising result, we suggest that $f(\epsilon)$ remains small even at
moderate values of $\epsilon$.

Figure 3: (a) Plot of the normalized temperature perturbation as a function of
time. (b) Plot of the scaled temperature perturbation as a function of scaled
time for $\epsilon\ll 1$ with $C_{2}=0.24$. (c) Scaled asymptotic value of the
temperature perturbations for all (non-fixed) cases as a function of
$\epsilon$. Only cases denoted by a square are plotted in (b). (d) Plot of the
temperature perturbation as a function of scaled time for $\epsilon\to\infty$
with $1/\zeta=0.6$. (e) Asymptotic value of the Temperature perturbation for
all (non-fixed) cases as a function of $\epsilon$. The cases plotted with a
square are plotted in (d) in addition to the cases with $\Theta_{z=1}$ fixed
at 0. (f) Plot of the normalized temperature perturbation for all cases at
finite $\epsilon$ with scaled time. The solution to (27) is included as a
dashed lines in (b) and (f), while the solution to (30) is plotted as a
dashed-dot line in (d). The timescales $\tau_{\epsilon}$ and $\tau_{\infty}$
are defined as in (31).

###  IV.4 How vigorous is the convection?

We now address our final question: How vigorous is the convection? We quantify
this by looking at the kinetic energy of the system
($\mathcal{K}=\frac{1}{2}\mathbf{u}\cdot\mathbf{u}$) in quasi-steady state. In
that state, the $\mathcal{K}$ equation reduces to a balance of the buoyant
production to dissipation. We have previously discussed the Nu scaling (13),
which we will use to scale $\langle w\Theta\rangle$ as in (15). Additionally,
we can use the typical scaling for the viscous dissipation rate in terms of a
turbulence length scale $l_{\epsilon}$. Writing these together, we have

$\displaystyle\langle w\Theta\rangle-\varepsilon=0,\qquad\langle
w\Theta\rangle\sim\frac{C\hbox{Ra}_{0}^{p}}{\sqrt{\hbox{Ra}_{0}\
\Pr}}\Delta^{p+1},\qquad\varepsilon\approx\frac{\mathcal{K}^{\frac{3}{2}}}{l_{\epsilon}},$
(32)

where
$\varepsilon=2\sqrt{\frac{\Pr}{\hbox{Ra}_{0}}}\langle\underline{\mathbf{e}}:\underline{\mathbf{e}}\rangle$
is the volume-averaged viscous dissipation rate based upon the rate of strain
tensor
$\underline{\mathbf{e}}=\frac{1}{2}\left(\nabla\mathbf{u}+\left(\nabla\mathbf{u}\right)^{T}\right)$.
We then solve for $\mathcal{K}$ as a function of $\hbox{Ra}_{0}$ and $\Delta$,

$\displaystyle\mathcal{K}=\left(l_{\epsilon}\frac{C\hbox{Ra}_{0}^{p}}{\sqrt{\hbox{Ra}_{0}\
\Pr}}\Delta^{p+1}\right)^{\frac{2}{3}}.$ (33)

To close this equation, we need to determine $l_{\epsilon}$.

Figure 4b is a plot of $l_{\epsilon}$, computed from the ratio of
$\mathcal{K}$ and $\epsilon$, as a function of $\hbox{Ra}_{\Delta}$. We find
that the data collapses well onto the curve,

$\displaystyle l_{\epsilon}=\left(5.3\pm 0.6\right)\times 10^{-4\
}\hbox{Ra}_{\Delta}^{0.475\pm 0.006}.$ (34)

While there is some variation in the $\mathcal{K}$ around the predicted mean
value, the scaling law (33) provides a reasonable estimate for the equilibrium
$\mathcal{K}$ in this system (see figure 4a). Much of the observed variance
about the horizontal line found in figure 4a could be removed by using a case
specific value of $l_{\epsilon}$, rather than using the general
parameterization (34). However, the case-specific value does not eliminate
oscillations about the mean.

Figure 4: (a) The scaled kinetic energy as a function of $\hbox{Ra}_{\Delta}$.
For this panel, $l_{\epsilon}$ is given by (34). The lines denote the time
evolution of $\mathcal{K}$ over the convective period of the simulation, while
the squares denote the mean steady state value. (b) Computed $l_{\epsilon}$ vs
$\hbox{Ra}_{\Delta}$.

The relative strength of the fluid inertia to viscosity can be quantified
through a Reynolds number Re. Using the predicted kinetic energy scaling, we
can define the Reynolds number as

$\displaystyle\hbox{Re}=\frac{\tilde{U}_{rms}H}{\nu}=\sqrt{\frac{2\hbox{Ra}_{0}}{\Pr}}\left(l_{\epsilon}\frac{C\hbox{Ra}_{0}^{p}}{\sqrt{\hbox{Ra}_{0}\
\Pr}}\Delta^{p+1}\right)^{\frac{1}{3}}.$ (35)

In the present paper, we did not vary the Prandtl number and we do not know
how $l_{\epsilon}$ changes with $\Pr$. Nevertheless, our kinetic energy
scaling suggests that
$\hbox{Re}\approx\Pr^{q}\hbox{Ra}_{0}^{0.60}\Delta^{0.44}$ for some constant
$q$ for large $\hbox{Ra}_{0}$. Ignoring the Prandlt number dependence, the
experimental work by Lam et al. [12] and Niemela and Sreenivasan [15] suggest
that $\hbox{Re}\sim\hbox{Ra}_{0}^{0.4}$. The results from nonlinear numerical
simulations by Clarté et al. [6] in a spherical shell with a similar boundary
condition appear to scale $\hbox{Re}\sim\hbox{Ra}_{0}^{0.5}$. However, the
present results are more consistent with the optimal numerical predictions by
Wen et al. [20] with $\hbox{Re}\sim\hbox{Ra}_{0}^{\frac{2}{3}}$. We are unsure
of the cause of this discrepancy, but one possibility is the two-dimensional
nature of the present simulations. Future work will investigate the cause of
this difference between studies.

##  V Conclusions

Aquatic systems are coupled to their atmospheric forcing. We have extended the
Rayleigh-Bénard problem to include this dynamic coupling. In particular, this
coupling results in an extra parameter $\beta$, a scaled effective thermal
conductivity that incorporates the effect of long-wave radiation, sensible
heat loss, and evaporative heat loss. This model reduces to the typical
Rayleigh-Bénard setup in the limit of $\beta\to\infty$. This extension to the
dominant theory provides a relatively simple model to translate the results of
Rayleigh-Bénard theory to environmental systems, such as lakes.

By defining an effective Rayleigh number,
$\hbox{Ra}_{\Delta}=\hbox{Ra}_{0}(1-\Theta_{z=1})$, we have answered our four
motivating questions:

  1. 1.

The surface heat transport is quantified through the Nusselt number Nu, which
we have shown to scale as

$\displaystyle\hbox{Nu}_{0}-1\sim C\hbox{Ra}_{\Delta}^{p},\qquad\hbox{ where
}C=0.138\pm 0.007,\ p=0.308\pm 0.003.$ (36)

  2. 2.

The equilibrium surface temperature ($\Theta_{1}$) is an implicit function of
$\hbox{Ra}_{0}$ and $\beta$, with

$\displaystyle\frac{C\hbox{Ra}_{0}^{p}}{\beta}\approx\frac{\Theta_{1}}{\left(1-\Theta_{1}\right)^{1+p}}.$
(37)

  3. 3.

The system rapidly cools to its equilibrium value on an approximate timescale
of

$\displaystyle\tau=\begin{cases}\phi\epsilon^{\frac{-p}{1+p}},&\epsilon\ll
1,\\\ \phi/\zeta,&\epsilon\to\infty,\\\ \end{cases}\qquad\hbox{ where
}\epsilon=\frac{\beta}{C_{2}\hbox{Ra}_{0}^{p}}\hbox{ and
}\phi=\frac{\sqrt{\hbox{Ra}_{0}\ \Pr}}{C_{2}\hbox{Ra}_{0}^{p}}.$ (38)

We estimate $C_{2}\approx 0.24$ and $1/\zeta\approx 0.6$.

  4. 4.

The kinetic energy of the induced convection scales as

$\displaystyle\mathcal{K}=\left(l_{\epsilon}\frac{C\hbox{Ra}_{0}^{p}}{\sqrt{\hbox{Ra}_{0}\
\Pr}}\Delta^{p+1}\right)^{\frac{2}{3}},\qquad\hbox{where}\quad
l_{\epsilon}=\left(5.3\pm 0.6\right)\times 10^{-4}\
\hbox{Ra}_{\Delta}^{0.475\pm 0.006}.$ (39)

While we have answered the questions put forth in this study, there remain
many open questions that have not been addressed. In particular, the present
simulations are two-dimensional. While outside the scope of this current
study, future work will investigate how three-dimensional effects modify the
present conclusions. Of particular interest would be the surface temperature
variations due to the convection. While the surface temperature variations
observed in this study were limited (see figure 2(d)), it would be interesting
to study these variations in three-dimensional simulations that will likely
have smaller scale structures at the water surface. In addition, the present
model assumes that the atmospheric conditions were not varying. For forcing
processes that change over a period much shorter than the equilibrium
timescale $\tau$, there may be some interplay between the various timescales
of interest.

We also note that there appears to be weak curvature in Nu as a function of
$\beta$. Busse and Riahi [3] highlighted that, near the critical Rayleigh
number, the structure of the flow is modified by the finite value of $\beta$.
Later, Ishiwatari et al. [11] showed that two-dimensional cells have different
preferred wavelengths in the two extreme cases of fixed temperature and fixed
flux. This has been further detailed by Clarté et al. [6] in nonlinear
simulations of spherical shells. A detailed investigation of how the change in
flow structure affects the curvature of Nu as a function of $\hbox{Ra}_{0}$
and $\beta$ is an interesting avenue for future research.

We conclude by providing an estimate for the value of $\beta$ and
$\hbox{Ra}_{0}$ in real lakes. Leppäranta [13] provides an estimate for
$\gamma\approx 20\ [\hbox{W m}^{-2}\hbox{K}^{-1}]$ for lakes in Southern
Finland. If we assume that the entire lake is well mixed to a depth of 10 m
with an atmospheric reference temperature of $\tilde{T}_{R}=283$ [K] and
bottom water temperature of $\tilde{T}_{\tilde{z}=0}=285$ [K] (e.g. in
September), we estimate that $\beta\approx 350\ (\approx 2^{8.5})$ and
$\hbox{Ra}_{0}=O\left(10^{13}\right)$. These ballpark estimates are higher,
but in a comparable range, to those considered here. While outside the scope
of this present paper, future research would verify the current scaling laws
in three-dimensions for these higher values of $\beta$ and $\hbox{Ra}_{0}$.

This manuscript discusses an extension to the usual Rayleigh-Bénard
configuration that accounts for common surface cooling processes. Our hope is
that this work will encourage a more direct comparison between Rayleigh-Bénard
theory and field measurements in lakes and oceans.

## Acknowledgement

This work was funded by the Natural Sciences and Engineering Research Council
of Canada (NSERC). I would like to thank Hugo Ulloa and Megan Davies Wykes for
their helpful feedback on this paper, along with three anonymous referees.

## References

  * Brown et al. [2005] E. Brown, A. Nikolaenko, D. Funfschilling, and G. Ahlers.  Heat transport in turbulent Rayleigh-Bénard convection: Effect of finite top- and bottom-plate conductivities.  _Physics of Fluids_ , 17(7):075108, July 2005\.  ISSN 1070-6631, 1089-7666.  doi: 10.1063/1.1964987.  URL http://aip.scitation.org/doi/10.1063/1.1964987.
  * Burns et al. [2020] K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown.  Dedalus: A flexible framework for numerical simulations with spectral methods.  _Physical Review Research_ , 2(2):023068, Apr. 2020.  ISSN 2643-1564.  doi: 10.1103/PhysRevResearch.2.023068.  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.023068.
  * Busse and Riahi [1980] F. H. Busse and N. Riahi.  Nonlinear convection in a layer with nearly insulating boundaries.  _Journal of Fluid Mechanics_ , 96(2):243–256, Jan. 1980.  ISSN 0022-1120, 1469-7645.  doi: 10.1017/S0022112080002091.  URL https://www.cambridge.org/core/product/identifier/S0022112080002091/type/journal_article.
  * Bénard [1901] H. Bénard.  Les tourbillons cellulaires dans une nappe liquide. - Méthodes optiques d’observation et d’enregistrement.  10(1):254–266, 1901.  doi: 10.1051/jphystap:0190100100025400.
  * Chillà et al. [2004] F. Chillà, M. Rastello, S. Chaumat, and B. Castaing.  Ultimate regime in Rayleigh–Bénard convection: The role of plates.  _Physics of Fluids_ , 16(7):2452–2456, July 2004\.  ISSN 1070-6631, 1089-7666.  doi: 10.1063/1.1751396.  URL http://aip.scitation.org/doi/10.1063/1.1751396.
  * Clarté et al. [2021] T. T. Clarté, N. Schaeffer, S. Labrosse, and J. Vidal.  The effects of a Robin boundary condition on thermal convection in a rotating spherical shell.  _Journal of Fluid Mechanics_ , 918:A36, July 2021.  ISSN 0022-1120, 1469-7645.  doi: 10.1017/jfm.2021.356.  URL https://www.cambridge.org/core/product/identifier/S0022112021003566/type/journal_article.
  * Doering [2020] C. R. Doering.  Turning up the heat in turbulent thermal convection.  117(18):3, 2020.  doi: 10.1073/pnas.2004239117.
  * Foster [1968] T. D. Foster.  Effect of Boundary Conditions on the Onset of Convection.  _Physics of Fluids_ , 11(6):1257, 1968.  ISSN 00319171.  doi: 10.1063/1.1692095.  URL https://aip.scitation.org/doi/10.1063/1.1692095.
  * Hitchen and Wells [2016] J. Hitchen and A. J. Wells.  The impact of imperfect heat transfer on the convective instability of a thermal boundary layer in a porous media.  _Journal of Fluid Mechanics_ , 794, May 2016.  ISSN 0022-1120, 1469-7645.  doi: 10.1017/jfm.2016.149.
  * Imboden and Wüest [1995] D. M. Imboden and A. Wüest.  Mixing Mechanisms in Lakes.  In A. Lerman, D. M. Imboden, and J. R. Gat, editors, _Physics and Chemistry of Lakes_ , pages 83–138. Springer Berlin Heidelberg, Berlin, Heidelberg, 1995.  ISBN 978-3-642-85134-6 978-3-642-85132-2.  doi: 10.1007/978-3-642-85132-2“˙4.
  * Ishiwatari et al. [1994] M. Ishiwatari, S.-I. Takehiro, and Y.-Y. Hayashi.  The effects of thermal conditions on the cell sizes of two-dimensional convection.  _Journal of Fluid Mechanics_ , 281:33–50, Dec. 1994.  ISSN 0022-1120, 1469-7645.  doi: 10.1017/S0022112094003022.  URL https://www.cambridge.org/core/product/identifier/S0022112094003022/type/journal_article.
  * Lam et al. [2002] S. Lam, X.-D. Shang, S.-Q. Zhou, and K.-Q. Xia.  Prandtl number dependence of the viscous boundary layer and the Reynolds numbers in Rayleigh-Bénard convection.  _Physical Review E_ , 65(6):066306, June 2002\.  ISSN 1063-651X, 1095-3787.  doi: 10.1103/PhysRevE.65.066306.  URL https://link.aps.org/doi/10.1103/PhysRevE.65.066306.
  * Leppäranta [2015] M. Leppäranta.  _Freezing of Lakes and the Evolution of their Ice Cover_.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2015.  ISBN 978-3-642-29080-0 978-3-642-29081-7.  doi: 10.1007/978-3-642-29081-7.
  * Lord Rayleigh [1916] O. Lord Rayleigh.  On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side.  _Philosophical Magazine_ , 32(192):17, 1916.
  * Niemela and Sreenivasan [2006] J. J. Niemela and K. R. Sreenivasan.  Turbulent convection at high Rayleigh numbers and aspect ratio 4.  _Journal of Fluid Mechanics_ , 557:411, June 2006.  ISSN 0022-1120, 1469-7645.  doi: 10.1017/S0022112006009669.  URL http://www.journals.cambridge.org/abstract_S0022112006009669.
  * Plumley and Julien [2019] M. Plumley and K. Julien.  Scaling Laws in Rayleigh‐Bénard Convection.  _Earth and Space Science_ , 6(9):1580–1592, Sept. 2019.  ISSN 2333-5084, 2333-5084.  doi: 10.1029/2019EA000583.  URL https://onlinelibrary.wiley.com/doi/abs/10.1029/2019EA000583.
  * Sparrow et al. [1964] E. M. Sparrow, R. J. Goldstein, and V. K. Jonsson.  Thermal instability in a horizontal fluid layer: effect of boundary conditions and non-linear temperature profile.  _Journal of Fluid Mechanics_ , 18(4):513–528, Apr. 1964.  ISSN 0022-1120, 1469-7645.  doi: 10.1017/S0022112064000386.  URL https://www.cambridge.org/core/product/identifier/S0022112064000386/type/journal_article.
  * Verzicco [2004] R. Verzicco.  Effects of nonperfect thermal sources in turbulent thermal convection.  _Physics of Fluids_ , 16(6):1965–1979, June 2004\.  ISSN 1070-6631, 1089-7666.  doi: 10.1063/1.1723463.  URL http://aip.scitation.org/doi/10.1063/1.1723463.
  * Waleffe et al. [2015] F. Waleffe, A. Boonkasame, and L. M. Smith.  Heat transport by coherent Rayleigh-Bénard convection.  _Physics of Fluids_ , 27(5):051702, May 2015.  ISSN 1070-6631, 1089-7666.  doi: 10.1063/1.4919930.  URL http://aip.scitation.org/doi/10.1063/1.4919930.
  * Wen et al. [2020] B. Wen, D. Goluskin, M. LeDuc, G. P. Chini, and C. R. Doering.  Steady Rayleigh–Bénard convection between stress-free boundaries.  _Journal of Fluid Mechanics_ , 905:R4, Dec. 2020.  ISSN 0022-1120, 1469-7645.  doi: 10.1017/jfm.2020.812.  URL https://www.cambridge.org/core/product/identifier/S0022112020008125/type/journal_article.
  * Wittenberg [2010] R. W. Wittenberg.  Bounds on Rayleigh–Bénard convection with imperfectly conducting plates.  _Journal of Fluid Mechanics_ , 665:158–198, Dec. 2010.  ISSN 0022-1120, 1469-7645.  doi: 10.1017/S0022112010003897.  URL https://www.cambridge.org/core/product/identifier/S0022112010003897/type/journal_article.