Penetrative Convection at High Rayleigh Numbers

TL;DR

This study investigates penetrative convection at high Rayleigh numbers, introducing a new stability parameter, and demonstrates how it influences flow characteristics and heat transfer in a fluid between two plates.

Contribution

The paper presents a new non-dimensional parameter, $\Lambda$, and explores its impact on penetrative convection through simulations, extending understanding beyond classical Rayleigh-Bénard convection.

Findings

Nusselt number increases as $\Lambda$ decreases for fixed Rayleigh number.

Flow characteristics are highly sensitive to the stability parameter $\Lambda$.

In the limit $\Lambda ightarrow 0$, the system behaves like classical Rayleigh-Bénard convection.

Abstract

We study penetrative convection of a fluid confined between two horizontal plates, the temperatures of which are such that a temperature of maximum density lies between them. The range of Rayleigh numbers studied is $Ra = \left[10^6, 10^8 \right]$ and the Prandtl numbers are $Pr = 1$ and $11.6$. An evolution equation for the growth of the convecting region is obtained through an integral energy balance. We identify a new non-dimensional parameter, $\Lambda$, which is the ratio of temperature difference between the stable and unstable regions of the flow; larger values of $\Lambda$ denote increased stability of the upper stable layer. We study the effects of $\Lambda$ on the flow field using well-resolved lattice Boltzmann simulations, and show that the characteristics of the flow depend sensitively upon it. For the range $\Lambda = \left[0.01, 4\right]$, we find that for a fixed $Ra$ the Nusselt number, $Nu$, increases with decreasing $\Lambda$. We also investigate the effects of $\Lambda$ on the vertical variation of convective heat flux and the Brunt-V\"{a}is\"{a}l\"{a} frequency. Our results clearly indicate that in the limit $\Lambda \rightarrow 0$ the problem reduces to that of the classical Rayleigh-B\'enard convection.

Equations

$$\rho ={\rho }_0\left[1-\alpha {\left(T-T_M\right)}^2\right],$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}=-\frac{1}{{\rho }_0}\nabla p+g\alpha {\left(T-T_M\right)}^2\mathbf{k}+\nu {\nabla }^2\mathbf{u},$$

$$\frac{\partial T}{\partial t}+\mathbf{u}\cdot \nabla T=\kappa {\nabla }^{2}T.$$

$$\theta =\frac{T-T_M}{\Delta T}.$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}=Pr\left(-\nabla p+Ra{\theta }^2\mathbf{k}+{\nabla }^2\mathbf{u}\right),$$

$$\frac{\partial \theta }{\partial t}+\mathbf{u}\cdot \nabla \theta ={\nabla }^2\theta ,$$

$$Ra=\frac{g\alpha {\left(\Delta T\right)}^2{L}_z^3}{\nu \kappa }\text{and}Pr=\frac{\nu }{\kappa }$$

$$\frac{\partial \theta }{\partial t}+\frac{\partial }{\partial x}\left(u\theta \right)+\frac{\partial }{\partial z}\left(w\theta \right)=\frac{{\partial }^2\theta }{\partial x^2}+\frac{{\partial }^2\theta }{\partial z^2}.$$

$$\frac{\partial \overline{\theta }}{\partial t}+\frac{\partial }{\partial z}\left(\overline{w^{\prime }{\theta }^{\prime }}\right)=\frac{{\partial }^2\overline{\theta }}{\partial z^2},$$

$$\overline{\Psi }=\frac{1}{L_x}{\int }_0^{L_x}\Psi \mathrm{d}x$$

$${\int }_0^{h^{-}}\frac{\partial \overline{\theta }}{\partial t}\mathrm{d}z=-{\left(\overline{w^{\prime }{\theta }^{\prime }}-\frac{\partial \overline{\theta }}{\partial z}\right)|}_{z=h^{-}}+{\left(\overline{w^{\prime }{\theta }^{\prime }}-\frac{\partial \overline{\theta }}{\partial z}\right)|}_{z=0}.$$

$${\int }_0^{h^{-}}\frac{\partial \overline{\theta }}{\partial t}\mathrm{d}z=-{\left(\overline{w^{\prime }{\theta }^{\prime }}-\frac{\partial \overline{\theta }}{\partial z}\right)|}_{z=h^{-}}-{\frac{\partial \overline{\theta }}{\partial z}|}_{z=0}.$$

$${\left(\overline{w^{\prime }{\theta }^{\prime }}-\frac{\partial \overline{\theta }}{\partial z}\right)|}_{z=h^{-}}=-{\frac{\partial \overline{\theta }}{\partial z}|}_{z=h^{+}}.$$

$${\int }_0^{h^{-}}\frac{\partial \overline{\theta }}{\partial t}\mathrm{d}z={\frac{\partial \overline{\theta }}{\partial z}|}_{z=h^{+}}-{\frac{\partial \overline{\theta }}{\partial z}|}_{z=0}.$$

$$\overline{\theta }(z,t)=\{\begin{array}{ll}{\overline{\theta }}_1=\left({\theta }_{\text{mixed}}-{\theta }_H\right)\frac{z}{{\delta }_1}+{\theta }_H;& \text{if }0\le z\le {\delta }_1,\\ {\overline{\theta }}_2={\theta }_{\text{mixed}};& \text{if }{\delta }_1\le z\le h^{-}-{\delta }_2,\\ {\overline{\theta }}_3=\left(\frac{h^{-}-z}{{\delta }_2}\right){\theta }_{\text{mixed}};& \text{if }{h}^{-}-{\delta }_2\le z\le h^{-}.\end{array}$$

$${\int }_0^{h^{-}}\frac{\partial \overline{\theta }}{\partial t}\mathrm{d}z={\int }_0^{{\delta }_1}\frac{\partial {\overline{\theta }}_1}{\partial t}\mathrm{d}z+{\int }_{{\delta }_1}^{h^{-}-{\delta }_2}\frac{\partial {\overline{\theta }}_2}{\partial t}\mathrm{d}z+{\int }_{h^{-}-{\delta }_2}^{h^{-}}\frac{\partial {\overline{\theta }}_3}{\partial t}\mathrm{d}z.$$

$${\int }_0^{h^{-}}\frac{\partial \overline{\theta }}{\partial t}\mathrm{d}z={\theta }_{\text{mixed}}\frac{\mathrm{d}h}{\mathrm{d}t},$$

$${\theta }_{\text{mixed}}\frac{\mathrm{d}h}{\mathrm{d}t}={\frac{\partial \overline{\theta }}{\partial z}|}_{z=h^{+}}-{\frac{\partial \overline{\theta }}{\partial z}|}_{z=0}.$$

$${\frac{\partial \overline{\theta }}{\partial z}|}_{z=h^{+}}=\frac{{\theta }_C-{\theta }_M}{1-h}=-\frac{{\theta }_0}{1-h},$$

$$-{\frac{\partial \overline{\theta }}{\partial z}|}_{z=0}=-{\frac{\partial \overline{T}}{\partial z}|}_{z=0}{\left(\frac{\Delta T}{L_z}\right)}^{-1}=-{\frac{\partial \overline{T}}{\partial z}|}_{z=0}{\left(\frac{\Delta T_1}{h}\frac{h}{L_z}\frac{\Delta T}{\Delta T_1}\right)}^{-1}=\frac{Q}{h}{\theta }_H,$$

$$\frac{\mathrm{d}h}{\mathrm{d}t}=-\frac{\Lambda }{\gamma \left(1-h\right)}+\frac{1}{\gamma }\frac{Q}{h},$$

$$\Lambda =\frac{{\theta }_M-{\theta }_C}{{\theta }_H-{\theta }_M}=\frac{{\theta }_0}{{\theta }_H}.$$

$$\frac{\mathrm{d}h}{\mathrm{d}t}=-\frac{\Lambda }{\gamma \left(1-h\right)},$$

$$\frac{\mathrm{d}h}{\mathrm{d}t}=\frac{1}{\gamma }\frac{Q}{h},$$

$$h(t)=\sqrt{h_0^2+\left(\frac{2Q}{\gamma }\right)t},$$

$$\frac{\Lambda }{\gamma \left(1-h_s\right)}=\frac{1}{\gamma }\frac{Q}{h_s},$$

$$h_s=\frac{Q}{\Lambda +Q}.$$

$$Nu=-{{\left(\frac{\partial {\overline{\theta }}_T}{\partial z}\right)}_T|}_{z=0}.$$

Full text

I Introduction

Penetrative convection refers to situations where a gravitationally unstable layer of fluid advances into a stable layer of fluid (Veronis, 1963; Spiegel, 1972). The motion of the fluid in the unstable layer is typically driven by a source of heat. Penetrative convection is relevant in both astrophysical and geophysical settings (e.g., Veronis, 1963; Saunders, 1962; Zahn, 1977), with typical examples of the former being the interaction between convective and radiative zones in stars Brummell et al. (2002); Hanasoge et al. (2016) and of the latter being the destruction of the near-ground stable layer in the atmosphere due to radiative heating from the ground Deardorff et al. (1969); Tennekes (1973) and the deepening of the upper ocean mixed layer due to surface cooling or formation of sea ice Farmer (1975); Turner (1979).

For concreteness we study penetrative convection in water, which has a density maximum at $T_{M}=4$ *∘*C. If the upper surface of a column of water is maintained at a temperature below $T_{M}$ and the lower surface is maintained at a temperature above $T_{M}$, the layer of fluid with temperature below $T_{M}$ is stably stratified and the layer with temperature above $T_{M}$ is unstably stratified. As the value of $Ra$ for the unstable layer increases, convection will be initiated, which then leads to the entrainment of the fluid from the stable layer and hence growth of the convecting region.

The first stability analysis of penetrative convection was carried out by Veronis Veronis (1963), who considered a column of water the bottom of which is maintained at [math] *∘*C and the top of which is maintained at a temperature greater than $4$ *∘*C, along with stress-free conditions for velocity. From a linear stability analysis of the Boussinesq equations he found that as the temperature of the upper boundary increases, the critical Rayleigh number ($Ra_{c}$) for the unstable layer decreases from its value for the classical Rayleigh-Bénard problem, reaching a minimum before attaining an asymptotic value. Veronis Veronis (1963) argued that this behavior of $Ra_{c}$ is due to three competing factors: (1) The presence of a stable layer relaxes the upper boundary condition, thereby allowing the flow in the unstable region to reach an “optimum” state. As the thickness of the stable layer increases with the top-plate temperature, higher values of the temperature are preferred; (2) The number of cells in the vertical increases with increasing temperature, with the cell in the stable layer deriving its energy from the flow in the unstable layer. Hence, to minimize this energy loss, lower values of the top-plate temperature are preferred; (3) The available potential energy increases up to a top-plate temperature of $8$ *∘*C, and does not change with any further increase in the temperature, thereby clearly favoring a top-plate temperature of $8$ *∘*C. A combination of these three factors results in $Ra_{c}$ attaining a minimum at $6.7$ *∘*C. Veronis Veronis (1963) also discovered that convection could set in at subcritical values of Rayleigh number, because any finite-amplitude disturbance that mixes water layers above and below the level of maximum density leads to the creation of a deeper unstable layer, thereby favoring onset of convection.

The first experimental study of penetrative convection was by Townsend Townsend (1964), who examined turbulent natural convection over a layer of ice. The bottom surface of the tank was ice covered and the upper free surface was maintained at a temperature of $25$ *∘*C. We estimate that the $Ra$ in his experiments, based on the total depth of the cell, was about $4.36\times 10^{8}$, which is well into the turbulent regime. His key observations were:

The amplitude of temperature fluctuations was largest close to the base of the stable layer. 2. 2.

He released dye into the stable region, some of which was entrained into the convecting region to reveal the existence of elongated plume structures that extended from the base of the lower layer to the base of the stable layer.

Townsend Townsend (1964) attributed the large amplitude of the temperature fluctuations to the generation of internal gravity waves in the stable layer. These waves were generated at the interface between stable and unstable regions by the random impingement of plumes originating at the bottom surface. A systematic measurement of the heat flux could not be made due to heat loss from the sidewalls.

Deardorff et al. Deardorff et al. (1969) took a different approach to study the dynamics of penetrative convection. Using water as the working fluid, and a temperature range far from the temperature of maximum density, their initial condition was one of stable stratification. Convection ensued once the temperature of the bottom plate was increased. The motivation of this configuration was to understand the lifting of the inversion layer due to heating of the ground, and thus the central focus was to understand the evolution of the convecting layer. Their theoretical model predicted that the thickness of the convecting layer grows diffusively ($\propto\sqrt{t}$, where $t$ is time) when the heat flux from the bottom plate was assumed to be constant. However, when a constant temperature was imposed at the bottom plate, they derived a modified evolution equation whose results were in agreement with measurements. The best fit to their theoretical solution gave the growth of the layer as $\sim t^{0.41}$ (Figure 11 of Deardorff et al. (1969)). This indicates that the results for constant temperature and constant flux conditions are not substantially dissimilar, at least for the growth of the layer in this configuration. This is also supported by the fact that the heat transport in Rayleigh-Bénard convection is the same for constant temperature and constant flux conditions (Johnston and Doering, 2009). Similar theoretical models have been constructed by Tennekes Tennekes (1973) and Mahrt and Lenschow Mahrt and Lenschow (1976) to study the evolution of the convective layer. The model of Mahrt and Lenschow Mahrt and Lenschow (1976) is obtained by integrating the equations of motion in the convecting layer, and it reduces to that of Tennekes Tennekes (1973) when shear generation by turbulence is neglected.

Penetrative convection is also important in the study of stars. A typical star is comprised of three regions: an inner radiative zone, an outer convective zone, and the tachocline, which is a transition layer between the radiative and convective zones, and is stably stratified Spiegel and Zahn (1992). Cold plumes from the outer convective zone penetrate the upper layers of the tachocline generating internal gravity waves, which are thought to play an important role in the turbulent transport of momentum in the tachocline Brummell et al. (2002); Dintrans et al. (2005); Lecoanet et al. (2015). Hence, a detailed study of penetrative convection is necessary for the understanding of the coupling between these different zones and the effects of that coupling on the magnetic field of the star.

The situation studied here bears resemblance to penetrative convection in an internally heated fluid, where the fluid, which is bound by horizontal surfaces maintained at equal temperatures, is non-uniformly Whitehead and Chen (1970) or uniformly Goluskin and Spiegel (2012); Goluskin and van der Poel (2016) heated. The presence of the heat source leads to the generation of an unstable upper layer and a stable bottom layer. The relevant questions for this setting are Goluskin and van der Poel (2016): (1) how does the heat flux vary with the strength of the heat source? and (2) how does the mean temperature of the fluid vary with the strength of the heat source? An important distinction from our work is that, due to the asymmetry introduced by the heat source, the heat flux at the top and bottom surfaces are not equal in the stationary state. Also, the dependence of the heat flux on the heat source differs in two and three dimensions Goluskin and van der Poel (2016). Additionally, we note here that Chen and Whitehead Chen and Whitehead (1968) had previously used the idea of non-uniform heating of the fluid layer to study finite-amplitude motions in the classical Rayleigh-Bénard convection.

By modeling the fluid density with a piecewise linear phenomenological equation of state, with the same fixed linear increase in the unstable region and an arbitrary variable linear decrease (characterized by a parameter $\mathcal{S}=\left[2^{-8},2^{8}\right]$) in the stable region, Couston et al. (2017) numerically studied the flow in a similar geometry. Because of the quantitative difference between our equation of state and their parameter $\mathcal{S}$, we are not in a position to make a quantitative comparison to our work. However, we note that in this and a related study Couston et al. (2017); Lecoanet et al. (2015) they found (a) the convective region to be similar to that of the classical Rayleigh-Bénard convection for $\mathcal{S}\geq 100$ and (b) gravity waves at the density interface, heuristically as do we under certain conditions.

Both the astrophysical and geophysical settings in which penetrative convection is important offer a wide range of complications, such as rotation, not part of our study. However, in the spirit of the original paper of Veronis Veronis (1963), we have found further basic fluid mechanical processes free from the ravages of such complications and these are of interest to study in their own right.

In this paper we consider penetrative convection in a fluid that has a density maximum at a temperature between two horizontal plates. We derive an evolution equation for the thickness of the convecting layer by integrating the heat equation in the unstable layer and by using a form for the horizontally averaged temperature field based on our previous studies of turbulent Rayleigh-Bénard convection (Toppaladoddi et al., 2015a). We then compare the theory with the results from high-resolution numerical simulations for large Rayleigh numbers. Finally, we discuss the effects of boundary conditions on the flow field and on the heat transport.

II Equations of Motion

Figure 1 is a schematic of the system studied here. The width and depth of the domain are $L_{x}$ and $L_{z}$, respectively, the depth of the convecting layer is $h$, the bottom (top) plate is maintained at a temperature $T_{H}$ ($T_{C}$), and the fluid has a density maximum at a temperature $T_{M}$. The temperatures are such that $T_{C}<T_{M}<T_{H}$.

The fluid considered here is water, described well with the following equation of state Veronis (1963):

[TABLE]

showing that the fluid has a maximum density $\rho_{0}$ when $T=T_{M}$. Making the Boussinesq approximation, the equations of motion are

[TABLE]

[TABLE]

and

[TABLE]

Here, $\boldsymbol{u}(\boldsymbol{x},t)$ is the velocity field, $p(\boldsymbol{x},t)$ is the pressure field, $g$ is acceleration due to gravity, $\alpha$ is the coefficient of thermal expansion, $\boldsymbol{k}$ is the unit vector along the vertical, $\nu$ is kinematic viscosity, $T(\boldsymbol{x},t)$ is the temperature field, and $\kappa$ is the thermal diffusivity of the fluid.

To non-dimensionalize Eqs. (2) – (4), we choose $L_{z}$ as the length scale, $\Delta T=T_{H}-T_{C}$ as the temperature scale, $U_{0}=\kappa/L_{z}$ as the velocity scale, $\rho_{0}\,\nu\,\kappa/L_{z}$ as the pressure scale, and $t_{0}=L_{z}^{2}/\kappa$ as the time scale. We also introduce the non-dimensional temperature $\theta$ as

[TABLE]

Using these scales, but retaining the pre-scaled notation, save for the temperature field, we obtain

[TABLE]

[TABLE]

and

[TABLE]

where

[TABLE]

are the Rayleigh and Prandtl numbers, respectively. Hence, in non-dimensional units we have $\theta(z=0)=\theta_{H}$, $\theta_{M}=0$, and $\theta(z=1)=\theta_{C}=-\theta_{0}$, where $\theta_{0}>0$.

For velocity, the boundary conditions at the top and bottom surfaces are no-slip and no-penetration; and we assume periodicity in the horizontal direction.

III Numerical Scheme and Validation

We use the Lattice Boltzmann Method Benzi et al. (1992); Chen and Doolen (1998); Succi (2001); Shan (1997) to study penetrative convection for large Rayleigh number. The code developed has been extensively tested against results from spectral methods for shear and buoyancy driven flows Toppaladoddi et al. (2015b, a, 2017). The buoyancy force is introduced into the lattice Boltzmann equation using the scheme of Guo et al. Guo et al. (2002).

The code has also been validated against the results of Blake et al. (1984) for $\Gamma=L_{x}/L_{z}=2$, $T_{H}=8$ *∘*C, $T_{M}=3.98$ *∘*C, $T_{C}=0$ *∘*C, and $Pr=11.6$. Figure 2 shows the comparison of $Nu(Ra)$ with their simulations. Our values of $Nu$ are consistently lower than theirs, which we attribute to the low resolution of $22\times 42$ grid points used in their study; our resolution is an order of magnitude higher along both the horizontal and vertical directions.

We should note here that due to the presence of the stable layer, the time taken to reach a stationary state is much longer than in the classical Rayleigh-Bénard setting. The steady state thickness of the convecting layer is reached when the conductive heat flux in the stable layer is equal to the heat flux from the unstable layer Townsend (1964); Adrian (1975).

The results from numerical simulations presented in the following sections were obtained using $\Gamma\equiv L_{x}/L_{z}=2$ and $Pr=1$.

IV Results

IV.1 Analytical Results

IV.1.1 Evolution of the Convecting Layer

Here, using Eq. (8), we derive an evolution equation for the depth of the convecting layer, $h(t)$. The flow is assumed incompressible, and thus Eq. (8) can be written as

[TABLE]

Integrating along $x$ and assuming periodicity, we find

[TABLE]

where

[TABLE]

denotes the horizontal mean, and primes denote deviation from the horizontal means. Now, we integrate Eq. (11) along the vertical in the convecting region to find

[TABLE]

Owing to the no-penetration condition at $z=0$, Eq. (13) reduces to

[TABLE]

We assume that the dominant mode of heat transport in the stable layer is conduction, and by demanding the continuity of heat flux at the interface between the stable and unstable layers (e.g., Worster, 2004), we have

[TABLE]

Using condition (15) in Eq. (14), we find that

[TABLE]

To evaluate the integral on the left hand side of Eq. (16), we make the following assumptions about $\overline{\theta}(z,t)$:

The convecting layer consists of a well-mixed region that is bounded by boundary layers on its top and bottom surfaces. 2. 2.

The small boundary-layer thicknesses ($\delta_{1}$ and $\delta_{2}$) are assumed to be constants. The argument being that the boundary layers reach a stationary state much more rapidly than the well-mixed region.

Figure 3 shows the assumed profile for $\overline{\theta}(z,t)$.

Based on this, we write

[TABLE]

The integral in Eq. (16) can now be written as

[TABLE]

Assuming $\theta_{\text{mixed}}$ to be a constant, the integrals are easily evaluated to yield

[TABLE]

and hence, Eq. (16) becomes

[TABLE]

Moreover, we have

[TABLE]

and

[TABLE]

where $\Delta T_{1}=T_{H}-T_{M}$ and $Q$ ($>0$) is the non-dimensional heat flux delivered to the convecting region. Hence, we have the following evolution equation for the thickness of the convecting region,

[TABLE]

where $\theta_{\text{mixed}}=\gamma\,\theta_{H}$, with $0<\gamma<1$ a constant, and $\Lambda=\theta_{0}/\theta_{H}$. We note that in our approach the evolution equation has been obtained by assuming a profile for the mean temperature based on our quantitative understanding of the flow structure in classical turbulent Rayleigh-Bénard convection. Any need to parametrize the turbulent heat flux is circumvented by the requirement that the heat flux be continuous at the interface between the stable and unstable layers.

Our analysis also reveals that, in addition to $Ra$ and $Pr$, there is another governing parameter in the system, which is given by

[TABLE]

In general, the range of values $\Lambda$ can take is $[0,\infty)$. A large (small) value of $\Lambda$ indicates that the stable layer is strongly (weakly) stratified. The characteristics of the flow depend sensitively on the value of $\Lambda$, and hence this is a very important parameter in the description of penetrative convection.

In Eq. (22), there are different balances between the terms for different times, which is heuristically like the balances found in double-diffusive Worster (2004) and solidification problems Turner et al. (1986); Huppert and Worster (1992). Let $\mathcal{T}_{t}$ be the time at which the initial transients decay and $\mathcal{T}_{g}$ be the time beyond which the flow reaches a stationary state. The convective layer evolves in the following three stages:

Transient state: $0\leq t\leq\mathcal{T}_{i}$

The dominant balance during this time period is:

[TABLE]

which implies that the convective layer shrinks. This is expected on the grounds that the flow responds to the bottom heat flux on a time scale of $\mathcal{O}(\mathcal{T}_{i})$, during which the second term of Eq. (22) is smaller. The value of $\mathcal{T}_{i}$ would depend on $Ra$ and $\Lambda$, and in general can be expected to decrease with increasing $Ra$ and decreasing $\Lambda$. 2. 2.

Growth: $\mathcal{T}_{i}<t\leq\mathcal{T}_{g}$

During this stage we have

[TABLE]

which implies that the thickness of the layer increases with time. For $Q$ constant, the solution to Eq. (25) is

[TABLE]

where $h_{0}$ is the thickness at $t=0$. Thus, our analysis recovers the result discussed above that the convective layer grows diffusively for constant heat flux when $Ra\gg 1$ Deardorff et al. (1969); Tennekes (1973); Farmer (1975). 3. 3.

Steady state: $t>\mathcal{T}_{g}$

In the final stage, the flow reaches a steady state and Eq. (22) becomes

[TABLE]

where $h_{s}$ is final thickness of the layer, which is

[TABLE]

For a fixed $Q$, when the upper layer is unstratified ($\Lambda=0$), Eq. (28) gives $h_{s}=1$ and the convective layer occupies the whole domain. In the opposite limit of very strong stratification of the upper layer ($\Lambda\rightarrow\infty$), we have $h_{s}\rightarrow 0$. Both of these limits are found in our simulations.

We note that the expression for effective $Nu$ ($Q$ in our notation) in the work of Moore & Weiss (Eq. (15) of Moore and Weiss (1973)) reduces to Eq. (28) after some algebraic manipulation.

IV.2 Numerical Results

IV.2.1 Thickness of the Convecting Layer

We compute the thickness of the convecting layer, $h(t)$, which is defined as the height at which $\overline{\theta}=0$. Figures 4 and 5 show the evolution of the convecting layer for $Ra=10^{7}$ and $\Lambda=2$ and $0.25$, respectively.

A fit to the region where $h(t)$ increases in time for $\Lambda=2$ in Figure 4 gives $h(t)\propto t^{0.23}$; whereas, for $\Lambda=0.25$ one obtains $h(t)\propto t$. This shows that the growth for the convecting layer is much faster when $\Lambda$ is small, which arises from two effects. Firstly, the initial thickness of the convecting layer is larger for $\Lambda=0.25$ than for $\Lambda=2$ (see Figures 4 and 5 for thickness at $t=0$). The convective motions are more vigorous in the former case, leading to faster growth. Secondly, for lower values of $\Lambda$, the developing convecting layer experiences little resistance in entraining fluid from the stable layer, which again leads to a faster growth.

Once the flow has reached a stationary state, we compute the averaged thickness, $h_{s}$. Figure 6 shows $h_{s}$ as a function of $\Lambda$ for $Ra=10^{7}$. The agreement between the theory and simulations is very good for small $\Lambda$, but decreases for large $\Lambda$. The large $\Lambda$ behavior arises from the suppression of convective motions, and hence mixing, in the interior of the unstable region. This leads to both conduction and convection becoming important throughout the unstable layer, and hence the temperature profile assumed in the theoretical analysis is no longer valid for these large $\Lambda$.

IV.2.2 Temperature Field

Figures 7 and 8 show the time evolution of the temperature field for $Ra=10^{7}$ and $\Lambda=2$ and $0.25$. These values of $\Lambda$ were chosen to clearly reveal the effects of the stable layer stratification on the flow characteristics.

In Figure 7, the plumes that are generated from the hot bottom plate do not penetrate the stable layer because the strength of the stratification. The fluid from the stable layer is entrained slowly, and the flow takes a very long time to reach a stationary state. This is also clearly seen in Figure 4, where the growth of the convecting layer is subdiffusive. Additionally, we observe internal gravity waves generated at the interface between the stable and unstable layers, as well as in the interior of the stable layer in Figures 7(b) – 7(d). The structure of the flow here is in qualitative agreement with the experimental observations of Townsend Townsend (1964), with $\Lambda\approx 5$.

In contrast to this, the plumes penetrate the stable layer when $\Lambda=0.25$. In fact, the temperature fields closely resemble those in the Rayleigh-Bénard problem, where the fluid has a linear equation of state. Hence we intuitively expect that the Rayleigh-Bénard problem is realized in the limit of $\Lambda\rightarrow 0$. Indeed, the rapid growth of the convecting layer, as seen in Figure 5, is partly due to the weak stratification of the stable layer.

The effects of $\Lambda$ can also be discerned by studying the temporally and horizontally averaged temperature profiles. To that end, figures 9 and 10 show $\overline{\theta}_{T}(z)$ for $Ra=10^{7}$ and $\Lambda=2$ and $\Lambda=0.25$, respectively. The temperature profile for $\Lambda=2$ is more asymmetric than for $\Lambda=0.25$. The stable layer is much thicker for $\Lambda=2$, which is seen by the linear profile extending from $z=1$ to $z=0.8$. However, for $\Lambda=0.25$ the top/bottom symmetry of the temperature profile closely resembles that from turbulent Rayleigh-Bénard convection, consistent with the argument that the penetrative convective flow approaches that of the classical Rayleigh-Bénard problem as $\Lambda\rightarrow 0$.

Additionally, Figure 11 shows the averaged temperature profile for $Ra=10^{7}$ and $\Lambda=4$, which is in qualitative agreement with the experiments of Adrian Adrian (1975) who had $\Lambda\approx 5$ – $6$.

IV.2.3 Metastability of Plume Patterns

Another interesting consequence of the presence of the stable layer is its effect on the dynamics of plume generation in the convective layer. In Rayleigh-Bénard convection, the flow, for a given $Ra$ and $\Gamma$, settles into a stationary state with a fixed number of convection rolls that transport heat from the bottom wall to the top Johnston and Doering (2009). However, in penetrative convection, we find that for large $\Lambda$ and certain $Ra$, the flow structures enter a metastable state.

Figures 12 and 13 show the evolution of temperature field for $Ra=5\times 10^{6}$ and $\Lambda=4$. Focussing on the number of upwelling plumes, we see that there are four plumes in Figure 12(b) and as $h$ increases this configuration becomes unstable and two of the four plumes merge in Figure 12(c) forming now a total of three plumes as seen in Figure 12(d).

As $h$ increases further, the new configuration becomes unstable and two of the three plumes merge, giving rise to a total of two plumes [Figure 13(a)]. With increasing time, two smaller plumes are generated which then merge with one of the two larger plumes [Figures 13(a) - 13(d)]. This cycle of generation and merger of plumes continues and the flow does not settle into a stationary state with respect to flow structure. It is interesting to note that qualitatively similar observations of plume merger and generation were made by Whitehead and Chen (1970) in their study of penetrative convection in internally heated fluid.

These merger events give rise to more energetic plumes that then impinge upon the stable layer. This is seen in Figure 13(d). However, because of the stability of the upper layer, the plumes only generate low-frequency oscillations, as seen in Figure 14.

IV.2.4 Heat Transport

The non-dimensional heat flux, $Nu$, from the lower surface to the upper surface can be obtained using

[TABLE]

We note here that only the choice of $L_{z}$ as the characteristic length scale and $\Delta T=T_{H}-T_{C}$ as the characteristic temperature scale gives $Nu=1$ when $Ra=0$. The simulations were run for sufficiently long times to obtain converged statistics to compute $Nu$.

Figure 15 shows the least-squares fits for $Nu$($Ra$, $\Lambda$) data for $Ra=\left[10^{6},10^{8}\right]$ and $\Lambda=\left[0.01,4\right]$. For each $\Lambda$, the relation between $Nu$ and $Ra$ is sought in terms of a power law: $Nu=A\times Ra^{\beta}$. Clearly, for a fixed value of $Ra$, $Nu$ increases with decreasing $\Lambda$. This is due to the fact the stability of the upper layer decreases as $\Lambda$ decreases, which in turn leads to more vigorous convective motions in the unstable layer and larger heat transport. For $\Lambda=4$, there is no appreciable convective motion even when $Ra=10^{6}$ and the heat transport in the entire domain is dominated by conduction.

Another quantity that is of interest is the convective heat flux, $Q_{c}=\left(\overline{w^{\prime}\,\theta^{\prime}}\right)_{T}$, and its variation with height. Deardorff et al. Deardorff et al. (1969) found that $Q_{c}$ remains positive in the convective region, but becomes negative near the interface due to entrainment of the fluid from the stable layer. Similar observations have also been made by Adrian Adrian (1975). For $Ra=10^{7}$, Figures 16 and 17 show how $Q_{c}$ changes as $\Lambda$ changes from 4 to 0.01, respectively.

It is clear that, for $\Lambda=0.01$, except for the boundary layers, $Q_{c}$ is constant in the unstable region. Hence, in this case, convective motions transport nearly all the heat. On the other hand, when $\Lambda=4$ convection is not the dominant mode of transport, even in the unstable region. This is reflected by the fact that $Q_{c}$ attains a maximum value, equal to $Nu$ for this case, in only a small region of the flow. Moreover, $Q_{c}$ changes sign again in the stable layer, which is due to the combined effects of entrainment of the fluid from the stable layer and the excitation of internal gravity waves. This is quantified by studying the height dependence of the Brunt-Väisälä frequency Turner (1979), which is defined in dimensional units as

[TABLE]

Figure 18 shows the height dependence of $\mathcal{N}^{2}$, scaled by the convective time scale $t_{c}=\sqrt{H/g\,\alpha\,\Delta T^{2}}$, for $Ra=10^{7}$ and $\Lambda=4$.

By definition in the stable region $\mathcal{N}^{2}>0$, and in the region where $0\leq\mathcal{N}^{2}\leq 1$ both entrainment and the internal gravity waves drive vertical motions of the fluid; internal gravity waves become dominant only for $z>0.73$, where $\mathcal{N}^{2}>1$.

In contrast, for $\Lambda=0.01$, the internal gravity waves play no appreciable role in generating vertical motions. This can be seen from Figure 19, where $\mathcal{N}^{2}<0$ in the entire domain, showing that convective motion of the fluid dominates.

V Conclusions

We have systematically studied penetrative convection of a fluid with a density maximum using both analytical and numerical tools. We derived an evolution equation for the growth of the convecting layer by integrating the heat equation in the convecting layer and by constructing the mean temperature profile based on our knowledge of the flow in turbulent Rayleigh-Bénard convection. In so doing, we have identified a new governing parameter, $\Lambda$, that measures the strength of the stratification of the stable upper layer and thereby exerts a controlling influence on the evolution of the underlying convecting layer. For a constant heat flux, $Q$, we recover the result from previous studies Deardorff et al. (1969); Tennekes (1973); Farmer (1975) that the convecting layer grows diffusively. The final steady thickness is shown to depend solely on the values of $Q$ and $\Lambda$.

In order to obtain an analytic equation for the evolution of $h(t)$, Eq. 22, we assumed that the heat transport in the stable layer is controlled by conduction. The veracity of this assumption is justified by the results of the analysis in the extreme limits of $\Lambda$, which provides the framework for the utility of such a simple approach.

High-resolution numerical simulations using the lattice Boltzmann method reveal that the growth of the convecting layer at a same $Ra$ depends sensitively on the value of $\Lambda$ – the smaller the value of $\Lambda$, the faster the convecting layer grows. The flow field was also found to depend sensitively on $\Lambda$. For larger values of $\Lambda$, the penetrative entrainment of the plumes by the stable upper layer is suppressed. However, for smaller $\Lambda$ entrainment into the stable layer is efficient and the flow rapidly reaches a stationary state. The temporally and horizontally averaged temperature profile for $\Lambda=4$ and $Ra=10^{7}$ was found to be in qualitative agreement with the temperature profile from the experiments of Adrian Adrian (1975).

We computed $Nu$ for $Ra=\left[10^{6},10^{8}\right]$ and $\Lambda=\left[0.01,4\right]$ and found that for a fixed $Ra$, as $\Lambda$ decreases, $Nu$ increases. This is consistent with the limit of $\Lambda\rightarrow 0$ in penetrative convection reducing to that of the classical Rayleigh-Bénard convection. For $\Lambda$ = $\left[0.01,4\right]$, power-laws were obtained for the data using a linear least-squares fit, giving the exponent $\beta$ in $Nu=A\times Ra^{\beta}$. Both $A$ and $\beta$ vary non-monotonically with $\Lambda$, but a consistent physical interpretation is only possible by studying the changes in $Nu$ and not $A$ or $\beta$ individually.

We conclude by noting that whilst the complexities of many of the astrophysical and geophysical settings in which penetrative convection is operative are not at play in our study, nonetheless key qualitative phenomena will not differ. Of principle relevance is the influence of rotation, which has the general effect of suppressing convection, as does stratification. Indeed, there is a direct mathematical analogy between rotating and stratified fluids, and under some conditions the analogy is exact Veronis (1970). Thus, because penetrative convection, as we have studied it here, couples a convective region with a strongly stratified region, we suggest that the analogy between rotation and stratification is of some use in considering the qualitative influence of rotation on our results.

Acknowledgements.

The authors acknowledge the support of the University of Oxford and Yale University, and the facilities and staff of the Yale University Faculty of Arts and Sciences High Performance Computing Center. S.T. acknowledges a NASA Graduate Research Fellowship. J.S.W. acknowledges NASA Grant NNH13ZDA001N-CRYO, Swedish Research Council grant no. 638-2013-9243, and a Royal Society Wolfson Research Merit Award for support.