Penetrative convection: heat transport with marginal stability assumption

Abstract

This paper investigates heat transport in penetrative convection with a marginally stable temporal-horizontal-averaged field or background field. Assuming that the background field is steady and is stabilised by the nonlinear perturbation terms, we obtain an eigenvalue problem with an unknown background temperature $\tau$ by truncating the nonlinear terms. Using a piecewise profile for $\tau$, we derived an analytical scaling law for heat transport in penetrative convection as $Ra\rightarrow \infty$: $Nu=(1/8)(1-T_M)^{5/3}Ra^{1/3}$ ($Nu$ is the Nusselt number; $Ra$ is the Rayleigh number and $T_M$ corresponds to the temperature at which the density is maximal). A conditional lower bound on $Nu$, under the marginal stability assumption, is then derived from a variational problem. All the solutions to the full system should deliver a higher heat flux than the lower bound if they satisfy the marginal stability assumption. However, data from the present direct numerical simulations and previous optimal steady solutions by Ding & Wu (J. Fluid Mech., vol. 920, 2021, A48) exhibit smaller $Nu$ than the lower bound at large $Ra$, indicating that these averaged fields are over-stabilised by the nonlinear terms. To incorporate a more physically plausible constraint to bound heat transport, an alternative approach, i.e. the quasilinear approach is invoked which delivers the highest heat transport and agrees well with Veronis's assumption, i.e. $Nu\sim Ra^{1/3}$ (Astrophys. J., vol. 137, 1963, p. 641). Interestingly, the background temperature $\tau$ yielded by the quasilinear approach can be non-unique when instability is subcritical.

1. Introduction

Penetrative convection refers to the phenomenon of an unstably stratified fluid layer advancing into a stably stratified fluid layer, which is ubiquitous in nature and industrial applications. Well-known examples include the stellar convection (Dintrans et al. 2005), water convection near $4\,^\circ {\rm C}$ (Townsend 1964; Couston & Siegert 2021; Olsthoorn, Tedford & Lawrence 2021) and ice growth (Wang et al. 2021a; Wang, Calzavarini & Sun 2021b). The penetrative convection can be either caused by a parabolic base temperature profile due to internal heating (Goluskin & Spiegel 2012; Goluskin & van der Poel 2016) or a nonlinear constitutive relationship between the liquid density and temperature (Veronis 1963; Musman 1968; Moore & Weiss 1973; Payne & Straughan 1987; Lecoanet et al. 2015; Toppaladoddi & Wettlaufer 2018; Wang et al. 2019, 2021a,b). For instance, the density of water exhibits a parabolic dependence on temperature near $4\,^\circ {\rm C}$ (higher than the icing point)

(1.1)\begin{equation} \rho=\rho_0\left(1-\alpha(T^*-T_M^*)^2\right), \end{equation}

where $T_M^*$ corresponds to the temperature at which the density is maximal. The thermal expansion coefficient $\alpha >0$ is in units of ${\rm K}^{-2}$, which defines a different Rayleigh number from that in the classical Rayleigh–Bénard convection. Unlike the classical Rayleigh–Bénard convection in which the flow is unstably stratified (density decreases with temperature linearly), the quadratic constitutive relation can result in a more complicated flow configuration: a stably stratified flow coexisting with an unstably stratified flow. However, the constitutive relation (1.1) deteriorates when the temperature $T$ becomes high and water's density shows a linear dependence on the temperature. Nevertheless, this paper will be constrained to the valid regime of the quadratic constitutive relationship and aims to provide understanding of the heat transport below ice.

An intriguing problem in turbulent penetrative convection is to estimate the global heat transfer between two parallel plates, which remains challenging. Over the past decades, however, there has been a relatively larger amount of research on the scaling law of heat transport in Rayleigh–Bénard convection (Ahlers, Grossmann & Lohse 2009). In the Rayleigh–Bénard convection, a classical scaling $Nu\sim Ra^{1/3}$ ($Nu$ being the Nusselt number and $Ra$ being the Rayleigh number) proposed by Malkus (1954) and an ultimate scaling $Nu\sim Ra^{1/2}$ (possibly with logarithmic corrections) proposed by Kraichnan (1962) are competing. It seems that experiments and direct numerical simulations (DNS) of turbulent convection support the classical scaling (Urban, Musilová & Skrbek 2011; Bouillaut et al. 2019; Iyer et al. 2020), but there are some competing claims of the ultimate scaling (He et al. 2012; Zhu et al. 2018). Indeed, when the internal heating effect is limited within the boundary layers, the classical scaling is dominant (Bouillaut et al. 2019). Either internal heating in the bulk or breakdown of the boundary layers can cause the transition to the ultimate scaling (Zhu et al. 2018; Bouillaut et al. 2019). Recent computations of steady solutions of the Boussinesq equations (no internal heating and the bounded walls are smooth) exhibit higher heat flux than all the existing numerical simulations or experimental observations, which are more consistent with the classical scaling (Sondak, Smith & Waleffe 2015; Waleffe, Boonkasame & Smith 2015; Wen et al. 2020a). In penetrative convection, Ding & Wu (2021) found that heat transfer by exact steady solutions is higher than in turbulent flows and it roughly follows the classical scaling $Ra^{1/3}$. They also derived an upper bound on heat transport in penetrative convection: $Nu\leqslant cRa^{1/2}$ (the best value of $c$ remains unknown). Yet, it is unknown if the $1/2$ scaling law can be achieved as $Ra\rightarrow \infty$ both in Rayleigh–Bénard convection and penetrative convection. It is worth mentioning that only quadratic constraints are incorporated in the upper bound problem, which yields bounds unrelated to that of the real physical system (Olson et al. 2021). Although Wen et al. (2022) found that the $1/2$ bound on heat transport in Rayleigh–Bénard convection can be lowered by adding the marginal stability assumption as a further constraint, their solution is not physically meaningful either. Very recently, Wen, Goluskin & Doering (2020b) computed exact steady solutions in Rayleigh–Bénard convection up to $Ra=10^{14}$ and found that the scaling approaches $1/3$ asymptotically. Thus, they raised the intriguing possibility that optimal steady rolls may transport more heat than turbulent convection as $Ra\rightarrow \infty$.

The $1/3$ scaling of heat transfer in Rayleigh–Bénard convection can be derived by an assumption of marginally stable boundary layers (Howard 1963). Therefore, it was thought that the horizontal-averaged field of the exact steady solutions is close to a marginally stable state (Sondak et al. 2015; Waleffe et al. 2015; Wen et al. 2020a,b; Ding & Wu 2021). In Waleffe et al. (2015) and Sondak et al. (2015), heat transport is actually slightly lower than $Ra^{1/3}$, which was thought to be attributed to small departures from the marginally stable configuration. For instance, unstably stratified layers were observed in the optimal steady solution. In fact, those exact steady solutions are dynamically unstable. For turbulent flows, an assumption of a ‘statistically’ marginally stable boundary layer might be more appropriate for the derivation of the $1/3$ scaling law. Malkus (1956) conceived the idea of ‘statistical’ stability of turbulent flows and he suggested that the turbulent mean field should be marginally stable and nonlinear momentum transport only play a stabilising role. Following Malkus's idea, a recent work by O'Connor, Lecoanet & Anders (2021) showed that a quasilinear approximation of the Boussinesq equation can exhibit a clear $1/3$ scaling, in which the nonlinear terms in the momentum equation were dropped and the mean temperature field was assumed to be driven by fluctuations. A temporal-horizontal averaged field was assumed to be marginally stable. A more interesting observation is that the quasilinear problem transports more heat than in turbulent flows. This implies that the truncated nonlinear terms are ‘stabilising’ turbulent flows and are ‘suppressing’ heat transport. However, there is a ‘dip’ in the mean temperature profile (the temperature gradient reverses) which contradicts the hypothesis of Malkus (1956).

In penetrative convection, the $1/3$ scaling can also be derived by assuming the boundary layers are marginally stable (Veronis 1963). However, current existing numerical studies of turbulent penetrative convection (Lecoanet et al. 2015; Toppaladoddi & Wettlaufer 2018; Wang et al. 2019) showed that the heat transport is below the $1/3$ scaling. Although Ding & Wu (2021) claimed that their exact steady solutions can achieve the $1/3$ scaling, an examination of the averaged field showed that it is not marginally stable. In fact, the shape of the marginally stable field remains an unknown question, e.g. is there also a ‘dip’ structure in the mean temperature field? In addition, how the scaling depends on the thickness of the stably stratified layer is not clear. In his paper, Malkus (1956) conjectured that the turbulent flow can be highly degenerate if the mean field is marginally stable, i.e. different turbulent flows. A recent study of two-dimensional plane Poiseuille flow confirmed Malkus's conjecture (Markeviciute & Kerswell 2021), which demonstrates that symmetries play important roles in turbulence: turbulence with a different symmetry can be different (the mean field is different consequently) (Xie, Ding & Xia 2018). Inspired by Markeviciute & Kerswell (2021), it is therefore natural to ask if the marginally stable mean field is unique in penetrative convection, i.e. for a given parameter, are there different marginally stable mean fields? A relevant study is Wen et al. (2022), in Rayleigh–Bénard convection, which demonstrated that the marginally stable mean field can be non-unique: the mean field derived from a variational problem is different from the quasilinear problem. However, it is questionable as to whether non-unique fields can be derived from the same approach or not, e.g. the variational problem or the quasilinear approach. Another motivation of this study is to move beyond quadratic functionals and incorporate more physically plausible constraint to bound heat transport in penetrative convection, e.g. deriving a conditional lower bound under the marginal stability constraint. Moreover, it will also be interesting to examine if the nonlinear terms are ‘stabilising’ turbulent penetrative convection and ‘suppressing’ heat transport in penetrative convection as in the Rayleigh–Bénard convection (O'Connor et al. 2021). We hope to provide insights into the above questions in this paper.

The rest of the paper is organised as follows. Section 2 formulates the problem mathematically and the linearised problem of a two-dimensional flow is obtained. A piecewise background field is used to test the stability problem and derive an analytical scaling law for penetrative convection under the marginal stability assumption. Furthermore, a variational problem which delivers a conditional lower bound on the heat flux is constructed in § 3. Section 4 presents the study of a quasilinear approach to penetrative convection, which demonstrates that the marginally stable background field can be non-unique. Direct numerical simulations are performed to assess the effect of truncated nonlinear terms in § 5. A conclusion is given in § 6.

2. Mathematical formulation

2.1. Governing equations

By introducing the liquid depth $d$ as the length scale, $\Delta T=T_{bottom}^*-T_{top}^*$ ($T_{top}^*$ and $T_{bottom}^*$ are temperatures of the top wall and bottom walls) as the temperature scale, the free-fall velocity $\sqrt {\alpha g(\Delta T)^2 d}$ as the velocity scale, $d/\sqrt {\alpha g(\Delta T)^2 d}$ as the time scale and $\rho _0\alpha g(\Delta T)^2 d$ as the pressure scale, we start with the following dimensionless equations for penetrative convection:

(2.1)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}

(2.2)\begin{gather}\partial_t\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}= - \boldsymbol{\nabla}p+ \sqrt{\frac{Pr}{Ra}}\nabla^2\boldsymbol{u}+(T-T_M)^2\boldsymbol{e}_z, \end{gather}

(2.3)\begin{gather}\partial_tT+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}T=\frac{1}{\sqrt{PrRa}}\nabla^2T, \end{gather}

where $\boldsymbol{u}$ is the velocity of the fluid, $\boldsymbol{e}_z$ is unit vector in z-direction, $T=(T^*-T_{top})/\Delta T$ is the dimensionless temperature ($T^*$ is the dimensional temperature), $Ra:=\beta g(\Delta T)^2d^3/\nu \kappa$ ($\beta$ being the thermal expansion coefficient and g being gravitational acceleration) is the Rayleigh number and $Pr:=\nu /\kappa$ is the Prandtl number ($\nu$ being the kinematic viscosity and $\kappa$ being the thermal diffusivity); $T_M=(T_M^*-T_{top}^*)/\Delta T$ describes the density maximal effect and $T_M=0,\,1$ indicates that the maximal density locates at the top/bottom boundary respectively. When $T_M=1$, the flow is stable for all $Ra$.

We assume that there is no slip on the walls. Therefore, boundary conditions for $\boldsymbol {u}$ and $T$ are

(2.4a–c)\begin{equation} \boldsymbol{u}|_{z=0,1}=\boldsymbol{0},\quad T|_{z=0}=1,\quad T|_{z=1}=0. \end{equation}

We decompose the flow into a ‘steady’ background state and a perturbation state

(2.5a,b)\begin{equation} {\boldsymbol{u}}=\bar{\boldsymbol{u}}+\boldsymbol{u}',\quad T=\tau+\theta'. \end{equation}

One may also assume that $(\bar {\bullet })$ is the ensemble average of $(\bullet )$. In the present work, we shall consider that the background state is defined using the following temporal-horizontal average:

(2.6a,b)\begin{equation} \bar{\boldsymbol{u}}=\lim_{t^*\rightarrow\infty}\frac{1}{2Lt^*}\int^{t^*}_0\int^L_{ - L}\boldsymbol{u}\,{\rm d}\kern0.06em x\,{\rm d}t,\quad \tau=\lim_{t^*\rightarrow\infty}\frac{1}{2Lt^*}\int^{t^*}_0\int^L_{ - L}T\,{\rm d}\kern0.06em x\,{\rm d}t. \end{equation}

If turbulence is ergodic, then the ensemble average should be identical to the temporal-spatial average. However, recent study suggests that there exist multiple states of turbulence (Xie et al. 2018; Wang et al. 2020). For a marginally stable field, the present study also suggests that the background field is steady by using a bifurcation analysis and branch continuation. But it is not guaranteed that the marginally stable field is always steady and one may relax the steadiness assumption and consider a time-periodic background flow. Since we will assess the stability of the background field, it is sufficient to consider a two-dimensional problem and the governing equations (2.1)–(2.3) are restated using Reynolds’ decomposition

(2.7)\begin{gather} \sqrt{\frac{Pr}{Ra}}\frac{{\rm d}^2\bar{u}}{{\rm d}z^2}=\partial_z\overline{w'{u'}}, \end{gather}

(2.8)\begin{gather}\frac{{\rm d}\bar{p}}{{\rm d}z}=\overline{{\theta'}^2}-\partial_z\overline{{w'}^2}, \end{gather}

(2.9)\begin{gather}\frac{1}{\sqrt{PrRa}}\frac{{\rm d}^2\tau}{{\rm d}z^2}=\partial_z\overline{w'{\theta'}}, \end{gather}

(2.10)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}'=0 \end{gather}

(2.11)\begin{gather}\partial_t\boldsymbol{u}'+\bar{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}'+\boldsymbol{u}' \boldsymbol{\cdot}\boldsymbol{\nabla}\bar{\boldsymbol{u}}+\boldsymbol{\nabla}p'- \sqrt{\frac{Pr}{Ra}}\nabla^2\boldsymbol{u}'-2(\tau-T_M)\theta'\boldsymbol{e}_z=\boldsymbol{\mathcal{T}}, \end{gather}

(2.12)\begin{gather}\partial_t\theta'+\bar{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}\theta'+ \frac{{\rm d}\tau}{{\rm d}z}w'-\frac{1}{\sqrt{PrRa}}\nabla^2\theta'=\mathcal {S}, \end{gather}

where the nonlinear terms are $\boldsymbol {\mathcal {T}}=\partial _z\overline {w'{u'}}\boldsymbol {e}_x-\boldsymbol {u}'\boldsymbol {{\cdot }} \boldsymbol {\nabla }{\boldsymbol {u}}'+{\theta '}^2\boldsymbol {e}_z-(\overline {{\theta '}^2}-\partial _z\overline {{w'}^2})\boldsymbol {e}_z$ and $\mathcal {S}=\partial _z\overline {w'\theta '}-\boldsymbol {u}'\boldsymbol {{\cdot }}\boldsymbol {\nabla }{\theta }'$. The set of governing equations shares many similar features to the resolvent analysis when using a forcing term to replace the nonlinear terms, which assumes the mean flow profile to predict the structure of the accompanying fluctuation fields (McKeon 2017). Following the original idea of Malkus (1956), we assume that the background field is marginally stable and is stabilised by the nonlinear terms, whence the nonlinear terms $\boldsymbol {\mathcal {T}}$ and $\mathcal {S}$ are dropped from the governing equations. And for the minimal choice of the background velocity field, we set $\bar {\boldsymbol {u}}=\boldsymbol {0}$, which yields the linear stability problem

(2.13)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}'=0, \end{gather}

(2.14)\begin{gather}\partial_t\boldsymbol{u}'+\boldsymbol{\nabla}p'-\sqrt{\frac{Pr}{{Ra}}}\nabla^2\boldsymbol{u}'-2( \tau-T_M)\theta'\boldsymbol{e}_z=\boldsymbol{0}, \end{gather}

(2.15)\begin{gather}\partial_t\theta'-\frac{1}{\sqrt{PrRa}}\nabla^2\theta'+\frac{{\rm d}\tau}{{\rm d}z}w'=0. \end{gather}

Actually, only $\bar {\boldsymbol {u}}=\boldsymbol {0}$ is found in the present work because branch continuation is used. However, it does not indicate that the background flow should be zero. If an alternative method is used, e.g. time-stepping method (O'Connor et al. 2021; Wen et al. 2022), or there is another branch disconnected from the present solution of $\bar {\boldsymbol {u}}=\boldsymbol {0}$, one may find a non-zero background field $\bar {\boldsymbol {u}}$, e.g. a mean ‘wind’ flow $\bar {\boldsymbol {u}}=u(z)\boldsymbol {e}_x\neq \boldsymbol {0}$. However, we will leave this possibility for future exploration. In the Rayleigh–Bénard convection, the time derivative terms can be safely dropped because exchange of stability always holds for arbitrary $\tau$ (Wen et al. 2022). In the penetrative convection, exchange of stability can be established for the trivial base state $(\bar {\boldsymbol {u}},\tau )=(\boldsymbol {0},1-z)$. For an arbitrary $\tau$, the eigenvalue can be complex, i.e. no exchange of stability (we tested this numerically but failed to prove it rigorously). Nevertheless, when using the bifurcation analysis with branch continuation, we observe that the marginally stable state is always neutral (i.e. the eigenvalue is always zero). Hence, in what follows, we drop the time derivative terms in the linear stability problem. For neutral stability, we rescale the pressure by $p'/\sqrt {Pr}\rightarrow p'$ and temperature $\theta '$ by $\theta '/\sqrt {Pr}\rightarrow \theta '$. Then, the Prandtl number will not appear in the governing equations of the neutral stability problem.

To characterise the heat flux, a Nusselt number defined as follows:

(2.16)\begin{equation} Nu:=\langle|\boldsymbol{\nabla} T|^2\rangle=\int^1_0\left(\frac{{\rm d}\tau}{{\rm d}z}\right)^2\,{\rm d}z+\langle|\boldsymbol{\nabla} \theta'|^2\rangle, \end{equation}

where $\langle (\bullet )\rangle =\int ^1_0\overline {(\bullet )}\,{\rm d}z$ is the volume-temporal average.

2.2. A piecewise background temperature

To examine the stability of (2.13)–(2.15), a simple and straightforward approach is to assume a piecewise $\tau$ such that we only need to adjust the thickness of boundary layer ($\epsilon$) to make it marginally stable. But the piecewise $\tau$ does not satisfy the equations of motion nor is it driven by the nonlinear fluctuation terms in (2.9). The choice of $\tau$ was inspired by previous numerical simulations that indicated that the averaged temperature field presents a similar structure as the piecewise $\tau$. In this section, we consider the case of $Ra\rightarrow \infty$ (namely, the fluid depth $d\rightarrow \infty$). Hence, the boundary layer thickness $\epsilon$ is used as the length scale. For simplicity, we remove the top plane since $d\rightarrow \infty$. Therefore, we assume the following background temperature profile:

(2.17)\begin{equation} \tau=\begin{cases} 1-(1-\bar{T})z, & 0\leqslant z\leqslant1 \\ \bar{T}={\rm const}. & 1< z<\infty. \end{cases} \end{equation}

Clearly, unlike the Rayleigh–Bénard convection in which the mixing region temperature $\bar {T}$ was given (see Currie 1967; Kerr 2016), $\bar {T}$ in the penetrative convection is an unknown parameter and is to be determined.

Eliminating the pressure and using the normal mode analysis $(\boldsymbol {u}',\theta ')=(\hat {\boldsymbol {u}},\hat {\theta })\exp ({\rm i}kx)$, we obtain the linearised equations

(2.18)\begin{gather} \frac{1}{\sqrt{Ra_\epsilon}}(D^2-k^2)^2\hat{w}-2k^2(\tau-T_M)\hat\theta=0, \end{gather}

(2.19)\begin{gather}\frac{1}{\sqrt{Ra_\epsilon}}(D^2-k^2)\hat{\theta}-\frac{{\rm d}\tau}{{\rm d}z}\hat{w}=0, \end{gather}

where k represents wavenumber and $Ra_\epsilon =\beta g(\Delta T)^2\epsilon ^3/\nu \kappa$.

There is an analytical solution in $1< z<\infty$

(2.20a,b)\begin{equation} \hat{w}=(c_0+c_1z+c_2z^2){\rm e}^{ - kz},\quad \hat{\theta}=\frac{4c_2}{(\bar{T}-T_M)\sqrt{Ra_\epsilon}}{\rm e}^{ - kz}. \end{equation}

The solution in $0\leqslant z\leqslant 1$ is solved numerically by requiring $\hat {w}$ and $\hat \theta$ to be continuous across the interface $z=1$

(2.21) \begin{equation} \hat{w}|^{1^+}_{1^-}=\frac{{\rm d}^i \hat{w}}{{\rm d} z^i}|^{1^+}_{1^-}=\hat{\theta}|^{1^+}_{1^-}=\frac{{\rm d}\hat{\theta}}{{\rm d} z}|^{1^+}_{1^-}=0,\quad i=1,2,3. \end{equation}

Numerical result shows that the critical Rayleigh number $Ra_\epsilon$ is achieved around $k\rightarrow 0$ as demonstrated in figure 1. And $Ra_\epsilon$ exhibits a strong dependence on the temperature in the mixing region $\bar {T}$, e.g. $Ra_\epsilon$ is minimal at $\bar {T}=0.5$ for $T_M=0$. To explore how $Ra_\epsilon$ varies with $\bar {T}$, we set the wavenumber $k$ at a small value $k=10^{-4}$ (for $k=0$, the eigenvalue problem is numerically singular) and the results are illustrated in figure 2. It is interesting that the minimal $Ra_\epsilon$ is achieved around $\bar {T}=(1+T_M)/2$ and the minimal $Ra_\epsilon$ varies like $64/(1-T_M)^2$. The penetrative convection is very similar to the classical Rayleigh–Bénard convection when $T_M=0$ (Wang et al. 2019; Ding & Wu 2021). Optimal steady solutions indicate the temperature in the mixing region is around $0.5(T^*_{bottom}+T^*_{top})$, which coincides with the stability analysis in this section. This may imply that the temperature in the bulk region tends to adjust its value such that the system is least stable and more heat can be transferred.

Figure 1. The contour plot of $Ra_\epsilon$ in the $(\bar {T},k)$ plane: (a) $T_M=0$ and (b) $T_M=0.7$.

Figure 2. The critical Rayleigh number $Ra_\epsilon$ computed numerically at $k=10^{-4}$ vs the temperature $\bar {T}$ of the mixing region. Solid lines are obtained by fixing $T_M$ as indicated in the plot. The dashed line illustrates that the minimal Rayleigh number $Ra_\epsilon =64/(1-T_M)^2$ can be achieved at $\bar {T}=(T_M+1)/2$, which agrees well with the numerical solution.

To account for the singular limit of $k\rightarrow 0$, an asymptotic solution is sought. We expand the solution in the region of $z\in [0,1]$ as

(2.22a–c)\begin{equation} \hat{w}=\hat{w}_0(z)+k\hat{w}_1(z)+\dots, \quad \hat{\theta}=\hat{\theta}_0(z)+k\hat{\theta}_1(z)\dots,\quad Ra_\epsilon=Ra_0+kRa_1+\dots, \end{equation}

and the leading-order solution is obtained by using the boundary condition at $z=0$

(2.23a,b)\begin{equation} \hat{w}_0=C_0z^2+C_1z^3,\quad \hat{\theta}_0= - \sqrt{Ra_0}(1-\bar{T})\left(\frac{C_0z^4}{12}+\frac{C_1z^5}{20}+C_2z\right). \end{equation}

Using the conditions in (2.21), we obtain

(2.24a–d)\begin{equation} c_2=C_0,\quad c_0+c_1=0,\quad C_1=0,\quad C_2= - C_0/3\neq0, \end{equation}

and

(2.25)\begin{equation} Ra_0=\frac{16}{(1-\bar{T})(\bar{T}-T_M)}. \end{equation}

One can also solve the problem up to the first-order approximation and obtain the solution for $Ra_1$ (Kerr 2016). However, we will not pursue this since the numerical simulation indicates that $Ra_1>0$. Therefore, the critical Rayleigh number is achieved at $k=0$, in contrast to Currie (1967), which showed that $k=O(1)$ in the Rayleigh–Bénard convection when the top wall is retained. However, if the top wall is moved away, Kerr (2016) showed that the critical Rayleigh number is also achieved at $k=0$ in the Rayleigh–Bénard convection, which is a common feature between penetrative convection and Rayleigh–Bénard convection. This common feature indicates that large-scale convection can transfer more heat than small-scale motions.

An optimal Rayleigh number $Ra_{opt}$ occurring at $\bar {T}=(1+{T}_M)/2$ can be found from (2.25)

(2.26)\begin{equation} Ra_\epsilon\geqslant Ra_{opt}=\frac{64}{(1-T_M)^2}. \end{equation}

The flow is linearly stable for $Ra_\epsilon < Ra_{opt}$. Note that the numerical solutions in figure 2 are in excellent agreement with the analytical solutions in (2.25) and (2.26) with relative errors smaller than $0.1\,\%$. However, most numerical simulations and the recent optimal steady solutions show that $\bar {T}$ is dependent on the Prandtl number, which also deviates from $(1+T_M)/2$ when $T_M$ is large (e.g. $T_M=0.5,0.7$) (Wang et al. 2019; Ding & Wu 2021). Perhaps, in these previous works, $Ra$ is not high enough or the truncated nonlinear terms are responsible for this deviation.

Under the marginal stability assumption, the piecewise temperature gives the following scaling law of heat transport in penetrative convection as $Ra\rightarrow \infty$:

(2.27)\begin{equation} Ra_\epsilon=2^{ - 3}RaNu^{ - 3}(1-T_M)^3=\frac{64}{(1-T_M)^2}\Rightarrow Nu=\frac{1}{8}(1-T_M)^{5/3}Ra^{1/3}, \end{equation}

where $Nu=(1-T_M)d/(2\epsilon )$. In Ding & Wu (2021), they assumed that the temperature in the mixing region is linearly dependent on $T_M$ and the bottom wall such that they derived (2.27) (the coefficient was not given). The present work justified their assumption and determined the coefficient $1/8$ using the piecewise background temperature analysis. Currie (1967) derived $Ra_\epsilon =32$ in the classical Rayleigh–Bénard convection using the piecewise background temperature, which is half of the present work ($Ra_\epsilon =64$ when $T_M=0$). It implies that the penetrative convection will show a different heat transport scaling law from the classical Rayleigh–Bénard convection ($Nu\sim 0.315Ra^{1/3}$) at very large $Ra$, although their heat transport is very close for $Ra\leqslant 10^9$ (see Wang et al. 2019; Ding & Wu 2021). It is worth pointing out that the prefactor $1/8$ depends on the choice of the background temperature $\tau$. For instance, in Kerr (2016), an error function yields $Ra_\epsilon =\sqrt {{\rm \pi} }$ and $Nu\sim 0.330Ra^{1/3}$ can be derived. Therefore, if one chooses a different background temperature profile for penetrative convection, a different prefactor would be expected but the $Ra^{1/3}$ should remain invariant.

3. Conditional lower bound: variational problem

In the previous section, we assumed a given shape for the background profile $\tau$, in this section we solve for the background profile $\tau$ as the numerical solution to a variational problem. In Wen et al. (2022), it is shown that a conditional lower bound on the heat transfer in Rayleigh–Bénard convection can be derived based on the marginal stability assumption. It is interesting that the conditional lower bound also exhibits the $1/3$ scaling. In this section, we follow Wen et al. (2022) and aim to derive a conditional lower bound on heat transfer in penetrative convection. The neutral stability problem of an unknown $\tau$ is stated as follows:

(3.1)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}'=0, \end{gather}

(3.2)\begin{gather}-\boldsymbol{\nabla}p'+\frac{1}{\sqrt{Ra}}\nabla^2\boldsymbol{u}'+2(\tau-T_M)\theta'\boldsymbol{e}_z=\boldsymbol{0}, \end{gather}

(3.3)\begin{gather}\frac{1}{\sqrt{Ra}}\nabla^2\theta'-\frac{{\rm d}\tau}{{\rm d}z}w'=0. \end{gather}

Here, we can also conclude that the scaling law is independent of $Pr$ under the marginal stability assumption.

Using (2.16), we can conclude that

(3.4a,b)\begin{equation} Nu\geqslant Nu_l,\quad Nu_l=\int^1_0\left(\frac{{\rm d}\tau}{{\rm d}z}\right)^2\,{\rm d}z. \end{equation}

It will be interesting to extremise $Nu_l$ such that it delivers a lower bound of heat transport. Therefore, we construct the following Lagrangian:

(3.5)\begin{align} \mathscr{L}&:=\int^1_0\left(\frac{{\rm d}\tau}{{\rm d}z}\right)^2\,{\rm d}z- \left\langle\boldsymbol{u}^{\dagger}\boldsymbol{\cdot}(-\boldsymbol{\nabla}p'+\frac{1}{\sqrt{Ra}}\nabla^2\boldsymbol{u}'+2( \tau-T_M)\theta'\boldsymbol{e}_z)\right\rangle\nonumber\\ &\quad -\left\langle\theta^{\dagger}\left(\frac{1}{\sqrt{Ra}}\nabla^2\theta'-\frac{{\rm d}\tau}{{\rm d}z}w'\right)\right\rangle-\langle p^{\dagger}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}'\rangle, \end{align}

wherein the marginal stability is imposed as a constraint; $\boldsymbol {u}^{\dagger}$, $\theta ^{\dagger}$ and $p^{\dagger}$ are adjoint variables of $\boldsymbol {u}'$, $\theta '$ and $p'$. It should be stressed that the conditional lower bound is not rigorously derived from the equations of motion. Therefore, the solutions of the full system can either deliver a higher heat flux or a lower heat flux than the conditional lower bound.

Variation of $\mathscr {L}$ gives the following Euler–Lagrange equations:

(3.6)\begin{gather} \delta\mathscr{L}/\delta p^{\dagger}:=\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}'=0, \end{gather}

(3.7)\begin{gather}\delta\mathscr{L}/\delta \boldsymbol{u}^{\dagger}:= - \boldsymbol{\nabla}p'+\frac{1}{\sqrt{Ra}}\nabla^2\boldsymbol{u}'+2(\tau-T_M)\theta'\boldsymbol{e}_z=0, \end{gather}

(3.8)\begin{gather}\delta\mathscr{L}/\delta \theta^{\dagger}:=\frac{1}{\sqrt{Ra}}\nabla^2\theta'-\frac{{\rm d}\tau}{{\rm d}z}w'=0, \end{gather}

(3.9)\begin{gather}\delta\mathscr{L}/\delta p:=\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^{\dagger} = 0, \end{gather}

(3.10)\begin{gather}\delta\mathscr{L}/\delta \boldsymbol{u}:= - \boldsymbol{\nabla}p^{\dagger} + \frac{1}{\sqrt{Ra}}\nabla^2\boldsymbol{u}^{\dagger} - \frac{{\rm d}\tau}{{\rm d}z} \theta^{\dagger}\boldsymbol{e}_z=0, \end{gather}

(3.11)\begin{gather}\delta\mathscr{L}/\delta \theta:=\frac{1}{\sqrt{Ra}}\nabla^2\theta^{\dagger} + 2(\tau-T_M)w^{\dagger} = 0, \end{gather}

(3.12)\begin{gather}\delta\mathscr{L}/\delta \tau:=2\frac{{\rm d}\tau^2}{{\rm d}z^2}+2\overline{\theta' w^{\dagger}}+\frac{{\rm d}\overline{w'\theta^{\dagger}}}{{\rm d}z}=0. \end{gather}

It is noted that there is a close relationship between the adjoint variables ($\boldsymbol {u}^{\dagger},\theta ^{\dagger},p^{\dagger}$) and the primitive variables ($\boldsymbol {u}',\theta ',p'$) in the Rayleigh–Bénard convection (Wen et al. 2022), which significantly simplifies the Euler–Lagrange system. However, in the penetrative convection, such a relationship is lost due to the nonlinear constitutive relation (1.1), which complicates the calculation of the Euler–Lagrange equations, i.e. we have to solve a coupled system between the primitive variables and adjoint variables. The Euler–Lagrange equations may admit non-unique solutions when the instability of the base state $(\bar {\boldsymbol {u}},\tau )=(\boldsymbol {0},1-z)$ is subcritical, i.e. solutions can exist when $Ra< Ra_c$ ($Ra_c$ is the critical Rayleigh number for the base state). However, to find a ‘lower’ bound, we only seek solutions which produce the lowest $Nu_l$ for $Ra>Ra_c$. Then, a parametric continuation is used to track the solutions of the Euler–Lagrange equations by continuously increasing $Ra$.

To solve the Euler–Lagrange equations, we expand the variables according to

(3.13a,b)\begin{equation} \begin{pmatrix} u' \\ w' \\ p' \\ \theta' \\ \end{pmatrix} = \begin{pmatrix} \sum_{1}^K u_m(z)\sin(k_mx) \\[6pt] \sum_{1}^K w_m(z)\cos(k_mx) \\[6pt] \sum_{1}^K p_m(z)\cos(k_mx) \\[6pt] \sum_{1}^K \theta_m(z)\cos(k_mx) \\ \end{pmatrix},\quad \begin{pmatrix} u^{\dagger} \\ w^{\dagger} \\ p^{\dagger} \\ \theta^{\dagger} \\ \end{pmatrix} = \begin{pmatrix} \sum_{1}^K u^{\dagger}_m(z)\sin(k_mx) \\[6pt] \sum_{1}^K w^{\dagger}_m(z)\cos(k_mx) \\[6pt] \sum_{1}^K p^{\dagger}_m(z)\cos(k_mx) \\[6pt] \sum_{1}^K \theta^{\dagger}_m(z)\cos(k_mx) \\ \end{pmatrix}. \end{equation}

Here, we assumed that there is a mirror symmetry in the solution (Ding & Wu 2021). Under mirror symmetry, it is easy to prove that $\overline {u'w'}=0$ and $\bar {\boldsymbol {u}}=0$. Although the mirror symmetry may be broken at high $Ra$ (see DNS below), we do not pursue this possibility in the present work because our numerical approach does not find an asymmetric solution of the Euler–Lagrange equations. Hence, the amplitude of the triangle mode satisfies the following equations:

(3.14)\begin{gather} k_mu_m+Dw_m=0, \end{gather}

(3.15)\begin{gather}k_mp_m+\frac{1}{\sqrt{Ra}}(D^2-k^2_m)u_m=0, \end{gather}

(3.16)\begin{gather}-Dp_m+\frac{1}{\sqrt{Ra}}(D^2-k^2_m)w_m+2(\tau-T_M)\theta_m=0, \end{gather}

(3.17)\begin{gather}\frac{1}{\sqrt{Ra}}(D^2-k^2_m)\theta_m-\frac{{\rm d}\tau}{{\rm d}z}w_m=0, \end{gather}

(3.18)\begin{gather}k_mu_m^{\dagger} + Dw_m^{\dagger} = 0, \end{gather}

(3.19)\begin{gather}k_mp_m^{\dagger} + \frac{1}{\sqrt{Ra}}(D^2-k^2_m)u_m^{\dagger} = 0, \end{gather}

(3.20)\begin{gather}-Dp_m^{\dagger} + \frac{1}{\sqrt{Ra}}(D^2-k^2_m)w_m^{\dagger} - \frac{{\rm d}\tau}{{\rm d}z}\theta^{\dagger}_m=0, \end{gather}

(3.21)\begin{gather}\frac{1}{\sqrt{Ra}}(D^2-k^2_m)\theta_m^{\dagger} + 2(\tau-T_M)w^{\dagger}_m=0, \end{gather}

(3.22)\begin{gather}2\frac{{\rm d}\tau^2}{{\rm d}z^2}+\sum^K_1{\theta_m w^{\dagger}_m}+ \frac{1}{2}\sum_1^K\frac{{\rm d}{w_m\theta^{\dagger}_m}}{{\rm d}z}=0. \end{gather}

Here, the notation $D:={\rm d}/{\rm d}z$ is the differential operator with respect to $z$.

The optimal wavenumber $k_m$ for the $m{{\rm th}}$ marginal mode is sought by adding the following equation:

(3.23)\begin{equation} \delta\mathscr{L}/\delta k_m:=\int^1_0u^{\dagger}_m p_m-\frac{2k_m}{\sqrt{Ra}}(u^{\dagger}_m u_m+w^{\dagger}_m w_m+\theta^{\dagger}_m\theta_m)+p^{\dagger}_m u_m\,{\rm d}z=0. \end{equation}

It is crucial to ensure all the modes are marginally stable such that a true lower bound can be yielded. A well-honed technique, i.e. Newton algorithm, has been developed in the study of upper bound problems (Plasting & Kerswell 2003; Ding & Kerswell 2020). A complementary method, which uses a time-stepping algorithm, also succeeded in obtaining the upper bound of heat transfer in Rayleigh–Bénard convection (Wen et al. 2020a). The advantage of the time-stepping algorithm is that it does not require a bifurcation analysis nor a branch continuation technique. Recently, the time-stepping algorithm was utilised to solve the quasilinear problem by O'Connor et al. (2021) and then the lower bound on heat transport in Rayleigh–Bénard convection (Wen et al. 2022). We should point out that, unlike the upper bound problem, the time-stepping algorithm does not guarantee the global optimal solution in either the quasilinear problem or the conditional lower bound problem. Also, it requires quite a long time to reach to a steady state and the time-stepping algorithm cannot be used to find an unstable ‘steady’ state. Indeed, we have known that the stability of the trivial state $(\bar {\boldsymbol {u}},\tau )=(\boldsymbol {0},1-z)$ is subcritical when $T_M\geqslant 0.4$ and $Pr=1$ (Ding & Wu 2021). Therefore, the present work chooses the Newton algorithm and branch continuation to track the solutions for the conditional lower bound problem and the quasilinear problem in § 4. To check if there is a new unstable mode emerging as $Ra$ increases, we have to solve (2.13)–(2.15) and examine if there is a positive eigenvalue $\lambda$, which is introduced into the problem by the normal mode analysis: $(\boldsymbol {u}',\theta ')=(\hat {\boldsymbol {u}},\hat {\theta })\exp ({\rm i}kx+\lambda t)$.

Figure 3 shows the dispersive relation of a marginally stable state at $Ra=10^9$. For $T_M=0$, we have detected $7$ marginally stable modes, while there is only $1$ marginally stable mode for $T_M=0.7$. Another interesting feature is that the eigenvalue $\lambda$ is real for all wavenumbers $k$ when $T_M=0$, while it is complex for $40< k<90$ when $T_M=0.7$ (in Rayleigh–Bénard convection, all the leading eigenvalues are real O'Connor et al. 2021; Wen et al. 2022). Although the mode with a complex eigenvalue is linearly stable, it might be important when the fully nonlinear system is considered. The bifurcation diagram of the critical modes is illustrated in figure 4. For small $T_M$ ($T_M\leqslant 0.1$), the penetrative convection is very similar to the Rayleigh–Bénard convection, and the bifurcation diagram exhibits a hierarchical structure (new unstable modes emerge as $Ra$ increases). When $T_M>0.1$, only a single critical mode is detected using our algorithm, e.g. $T_M=0.3$, $0.7$. This demonstrates that the thickness of the stable layer plays an important role in selecting the marginal modes.

Figure 3. The dispersive relation of a marginally stable state obtained using the variational problem at $Ra=10^9$.

Figure 4. The bifurcation diagram of the wavenumber $k_m$, (a–d) $T_M=0,0.1,0.3,0.7$. Each line stands for the wavenumber $k_m$ of the $m{{\rm th}}$ mode vs the Rayleigh number $Ra$. In (a,b), large wavenumber modes appear in pairs (adjacent solid and dashed lines). In (c,d), only a single critical mode is detected.

To illustrate the critical modes at different values of $T_M$, two typical cases of $T_M=0$ and $T_M=0.3$ are plotted in figures 5 and 6. For $T_M=0$, there are seven marginal modes: a central mode with six wall modes (the wall mode either concentrates near the top wall or the bottom wall, and they usually appear in pairs). The wavenumber of the central mode is of order $O(1)$, which seems to persist as $Ra\rightarrow \infty$, and the wavenumber of other modes scales as $Ra^{1/3}$ (see figure 4). Equation (2.27) shows that the thickness of a marginally stable boundary layer $\epsilon$ scales as $Ra^{-1/3}$. Therefore, if the convection cell is confined within the boundary layer, then its size is proportional to the boundary layer thickness and we can immediately obtain $k=2{\rm \pi} /\epsilon \sim Ra^{1/3}$. Like the long-wave analytical solution in § 2.2, the central mode, which occupies the whole domain, represents the large-scale motion in the system. However, for $T_M=0.3$, there is only a single marginal mode which always concentrates near the bottom wall (the corresponding wavenumber also scales as $Ra^{1/3}$). Such an interesting feature indicates that convection is confined within the boundary layer near the bottom wall and the outside fluid is stationary. Figure 7 demonstrates that the background temperature profile does have two boundary layers when $T_M\lesssim 0.1$ and it only has a single boundary layer when $T_M>0.1$. Nevertheless, we observe that the background temperature profile is not physical: the heat flux is not conserved – the heat flux at the bottom wall is not equal to that at the top wall. This non-physical feature is caused by the non-physical driving force in (3.22). Therefore, we can conclude that the variational problem delivers a non-physical bound for the penetrative convection, which, however, delivers a physically plausible solution in the classical Rayleigh–Bénard convection (Wen et al. 2022). If the solution of the full system satisfies the marginal stability assumption, it should deliver a higher heat flux than the present conditional lower bound. If heat flux in the full system is lower than the present conditional lower bound, the corresponding mean field should violate the marginal stability assumption, i.e. the averaged field is not marginally stable.

Figure 5. The eigenvectors of the marginally stable state obtained from the variational problem: $T_M=0$ and $Ra=10^9$. In (c,d), zoom-in insets are shown.

Figure 6. The eigenvectors of the marginally stable state obtained from the variational problem at different Rayleigh numbers and $T_M=0.3$.

Figure 7. The marginally stable temperature profile obtained from the variational problem at $Ra=10^9$.

4. Quasilinear problem

In § 3, the background temperature derived from the variational problem is non-physical because of the non-physical driving force. Now, we examine the marginal stability assumption of $\tau$ that is yielded from the quasilinear problem, i.e. the background temperature field is driven by the perturbation field (see (2.9))

(4.1)\begin{equation} \frac{{\rm d}^2\tau}{{\rm d} z^2}=\frac{\sqrt{Ra}}{2}\sum_1^K\frac{{\rm d} (w_m\theta_m)}{{\rm d}z}. \end{equation}

Note that the Prandtl number $Pr$ can also be scaled out from (4.1) by introducing $\theta '/\sqrt {Pr}\rightarrow \theta '$. The global heat transfer is conserved when the background temperature field is driven by the physical perturbations. Using (4.1), we can expect that if there is a boundary layer near the bottom wall, there should be a complementary boundary layer near the top wall. Using the Newton algorithm, it is more tedious to find the solution of (3.14)–(3.21) and (4.1) than the variational problem in § 3, because there is no equation for the wavenumber $k_m$ in the quasilinear problem. In O'Connor et al. (2021), both harmonic and sub-harmonic waves were included in their quasilinear problem, i.e. $k_m=mk_c/2$ ($m\in \mathbb {N}^+$ and $k_c$ is the critical wavenumber of primary instability). It indicates that we should consider an infinitely long domain to allow those harmonic and sub-harmonic waves. In Wen et al. (2022), a finite domain of size $L=2{\rm \pi} /k_c$ was considered and $k_m=mk_c$ ($m\in \mathbb {N}^+$). Comparison between O'Connor et al. (2021) and Wen et al. (2022) indicates that the domain size's effect is negligible. Hence, we follow Wen et al. (2022) and fix the domain size in the present work, which also allows us to perform DNS with initial conditions seeded by the quasilinear problem.

The bifurcation diagram of $k_m$ is illustrated in figure 8, indicating that the quasilinear problem undergoes a hierarchical bifurcation as $Ra$ increases for all $T_M$, which is an obviously different feature from the lower bound problem. The central mode $m=1$ is always critical for all $T_M$ in the quasilinear problem. The largest two wavenumbers exhibit a $Ra^{1/3}$ scaling, as demonstrated in figure 8, which implies that there are two boundary layers in the system. Thicknesses of the two boundary layers are proportional to $Ra^{-1/3}$ and convection cells should be confined within them. Figure 9 demonstrates that the two modes $m=23$, $25$ are concentrated in the near-wall regions, which are due to the instability of the boundary layers. Figure 9 also indicates that larger-scale convection contains more kinetic energy.

Figure 8. The bifurcation diagram of the wavenumber $k_m$ of the quasilinear problem, (a–d) $T_M=0,0.1,0.3,0.7$.

Figure 9. The eigenvectors of the marginally stable state obtained from the quasilinear problem: $T_M=0$ and $Ra=10^9$.

Figure 10 depicts the profiles of the background temperature of the quasilinear problem at $Ra=10^9$. It demonstrates that the boundary layers become thicker and the asymmetry of $\tau$ becomes more conspicuous as $T_M$ increases (the top boundary layer is thicker than the bottom boundary layer). The mean temperature in the bulk region increases but the mean temperature is below $(T_M+1)/2$. It may be caused by the constraint of the top wall or the relatively low Rayleigh numbers ($Ra\ll \infty$). A similar feature to the classical Rayleigh–Bénard convection is that there are also two dips in the profile of background temperature. However, as $T_M$ increases, the protruding part of the lower ‘dip’ becomes smaller, while it remains significant in the top ‘dip’. The ‘dip’ structure, however, is not significant in the optimal horizontal-averaged temperature profiles. In addition, figure 10 indicates that the quasilinear problem predicts a thinner boundary layer than the optimal steady solution by Ding & Wu (2021), and it seems that the optimal steady solution yields a mean temperature profile that is more consistent with the Malkus (1954) hypothesis. Unlike the penetrative convection, the ‘dip’ structures are significant in both the steady optimal solutions (Sondak et al. 2015) and the quasilinear problem (O'Connor et al. 2021) of the classical Rayleigh–Bénard convection. The predicted mean averaged temperature in the bulk region is a bit lower than the optimal horizontal-averaged temperature due to the consequence of thinner boundary layers. It also indicates that the quasilinear problem predicts a higher heat flux than the optimal steady solution.

Figure 10. The marginally stable temperature profiles (lines) obtained from the quasilinear problem at $Ra=10^9$. Symbols are the optimal horizontal-averaged temperature obtained from Ding & Wu (2021) at $Ra=10^9$, $Pr=1$.

Furthermore, we plot the Nusselt number predicted by the quasilinear approach and the conditional lower bound in figure 11, and the Nusselt number is rescaled by $(1-T_M)^{-5/3}$. It is interesting that the quasilinear approach shows that $Nu\sim 0.17(1-T_M)^{5/3}Ra^{1/3}$ when $T_M\geqslant 0.3$ (lines are almost overlapping). But for small $T_M$, the quasilinear approach shows that heat transport is slightly lower than $0.17(1-T_M)^{5/3}Ra^{1/3}$, which, perhaps, is because the Rayleigh number is not large enough (however, they are very close to $(1/8)(1-T_M)^{5/3}Ra^{1/3}$). An interesting feature is that the line can be multi-valued when the instability of the base state $(\bar {\boldsymbol {u}},\tau )=(\boldsymbol {0},1-z)$ is subcritical (see the solid lines for $T_M=0.5$, $0.7$ in figure 11). This indicates that the background temperature field is non-unique, although the lower branch is linearly unstable (Ding & Wu 2021). Note that the subcritical instability nature of penetrative convection indicates that the retained nonlinear term $\partial _z\overline {w'\theta '}$ destabilises the base state $(\bar {\boldsymbol {u}},\tau )=(\boldsymbol {0},1-z)$, thus contradicting Malkus's assumption. At high $Ra$, our algorithm did not find a double-branch solution but this does not mean the solution of the quasilinear problem is unique. The lower bound roughly follows $(1-T_M)^{5/3}Ra^{1/3}$ when $T_M\leqslant 0.3$, and it deviates from $(1-T_M)^{5/3}Ra^{1/3}$ when $T_M$ is large, which perhaps is also because of the relatively low $Ra$. Both the quasilinear problem and the lower bound indicate that $Nu\sim Ra^{1/3}$.

Figure 11. The Nusselt number $Nu$ vs the Rayleigh number $Ra$; ‘QL’ stands for quasilinear problem (solid lines) and ‘VS’ is for the variational lower bound problem (dashed lines).

5. Direct numerical simulations

In previous sections, different choices of the background field $\tau$ were examined and the ensuing turbulent heat transport under the marginal stability assumption was examined. Now we investigate the heat transport and flow fields using DNS. To check the influence of the truncated nonlinear terms in the quasilinear problem, DNS are performed by initialising simulations with the quasilinear states. Aiming at understanding how mirror symmetry influences the turbulent state, two different sets of simulations were performed: (i) with a mirror symmetry (ii) without a mirror symmetry. Here, we only consider a two-dimensional study. If there are two different statistical states emerging from the two sets of numerical simulations, this indicates that the turbulent state is degenerate in penetrative convection. The spatial discretisation was accomplished using a Fourier-spectral method. The axial direction was discretised with a Fourier transform (Frigo & Johnson 2005) and the 2/3 dealiasing rule was implemented. A finite difference method with non-uniform grid clustering near the walls (Chebyshev points) was used in the wall-normal direction. The time discretisation was accomplished with a second-order Adams–Bashforth–Crank–Nicolson method. Details of the numerical method can be found in Willis (2017) and Ding & Wu (2021).

To illustrate the nonlinear evolution of the quasilinear state, we consider two cases of $T_M=0$, $0.7$. The Nusselt number is monitored during the numerical simulation, as shown in figures 12 and 14, which demonstrates that the quasilinear steady state is nonlinearly unstable. It also indicates that the Nusselt numbers of the symmetric and asymmetric cases are identical and the flow becomes steady when the Rayleigh number is low (e.g. $Ra=10^7$). However, as $Ra$ increases to $10^8$, we observe that the time series of $Nu(t)$ start to deviate from each other as $t$ increases and both cases are chaotic. This suggests that the system can allow different chaotic solutions, i.e. the ‘turbulent’ state can be degenerate in penetrative convection. It also implies that the background temperature field can be non-unique at high $Ra$ since we can obtain different averaged fields using the simulation data. Figures 12 and 14 also demonstrate that spontaneous symmetry breaking occurs at high Rayleigh number in the system. To illustrate the spontaneous symmetry breaking, we plot the flow field and thermal field in figures 13 and 15. For $T_M=0$, only two large convection cells are observed in both the symmetric simulation and asymmetric simulation. However, in the asymmetric simulation, the hot plume near the bottom wall and the cold plume near the top wall start to oscillate when symmetry breaks, causing a mean ‘wind’ flow in the system. The mean ‘wind’ flow is time periodic as indicated by the vortex incline angle. As demonstrated in figure 12, the oscillating asymmetric flow generally transfers less heat than the symmetric flow. However, for $T_M=0.7$, we observe that the asymmetric flow can transfer more heat than the symmetric flow, as demonstrated in figure 14, in contrast to the case of $T_M=0$. Unlike the case of $T_M=0$, we observe that the initial field develops into a two-layer structure which undergoes the spontaneous symmetry breaking at $t\approx 160$. The two-layer structure is unstable, which evolves into a symmetric four-convection-cell structure or an asymmetric two-convection-cell structure as demonstrated in figure 15. Note that the width of the symmetric four-convection-cell structure is approximately half of the asymmetric two-convection-cell, indicating that large-scale motion can transfer more heat.

Figure 12. The evolution of Nusselt number with time $t$ at $Pr=1$, $T_M=0$. The solid (dashed) lines represent simulation results with (without) mirror symmetry.

Figure 13. The instantaneous flow and thermal fields computed at $Pr=1$, $T_M=0$ and $Ra=10^9$.

Figure 14. The evolution of Nusselt number with time $t$ at $Pr=1$, $T_M=0.7$. The dashed (solid) lines represent simulation results with (without) mirror symmetry.

Figure 15. The instantaneous flow and thermal fields computed at $Pr=1$, $T_M=0.7$ and $Ra=10^{10}$.

The long-time average of Nusselt number $Nu$ is rescaled by factor $(1-T_M)^{-5/3}Ra^{-1/3}$ ($Nu$ is averaged within $t\in [1000,2000]$) and plotted against $Ra$ as shown in figure 16. For $T_M=0.5$, $0.7$, we observe that heat transport via optimal steady solution (Ding & Wu 2021) follows $(1-T_M)^{5/3}Ra^{1/3}$ asymptotically. And the optimal steady solution transfers more heat than the two-dimensional ‘turbulent’ states (both symmetric and asymmetric simulations deliver lower $Nu$ than the optimal steady solution). The heat transport by two-dimensional ‘turbulent’ states and the optimal steady state is bounded by the quasilinear solution and the lower bound. However, for $T_M=0$, 0.1, 0.3, heat transport by the two-dimensional ‘turbulent’ states and the optimal steady state is lower than the lower bound when the Rayleigh number is large. This indicates that the averaged thermal field of the two-dimensional ‘turbulent’ states is not marginally stable. In this sense, it contradicts the hypothesis of Malkus (1956). However, figure 16 implies that the background field is over-stabilised by the nonlinear terms at high $Ra$ and the global heat transfer is consequently reduced, which agrees well with Malkus's assumption. It should be also stressed that the heat transfer computed using the full system is less than the quasilinear problem, similar to the findings by Wen et al. (2022) in Rayleigh–Bénard convection. Although many previous works claimed that penetrative convection is almost identical to the Rayleigh–Bénard convection when $T_M=0$, it is questionable whether the $Nu\sim Ra^{1/2}$ scaling can be observed in penetrative convection as $Ra\rightarrow \infty$.

Figure 16. The Nusselt number $Nu$ vs the Rayleigh number $Ra$. ‘OP’ stands for optimal solution in Ding & Wu (2021). ‘DNSS’ (‘DNSA’) is for the direct numerical simulation with (without) mirror symmetry. Solid lines are quasilinear approximation and dashed lines are the variational lower bound for various values of $T_M$: blue lines for $T_M=0$; red lines for $T_M=0.1$; green lines for $T_M=0.3$; magenta lines for $T_M=0.5$; cyan lines for $T_M=0.7$.

6. Conclusion

The present study investigates heat transport in penetrative convection by assuming a marginally stable temporal-horizontal-averaged field or background field. A linear stability problem is built around the background field by truncating the nonlinear perturbation terms. A minimal choice of the background field is then considered: $(\bar {\boldsymbol {u}},\bar {T})=(\boldsymbol {0},\tau (z))$, leading to an eigenvalue problem of an unknown background temperature $\tau$. Three different approaches were used to determine the profile of the marginally stable temperature $\tau$: a piecewise profile, a variational problem and a quasilinear approach.

Using a piecewise temperature profile, we derived that large convection cells transfer more heat and the optimal temperature in the mixing region is the average value of the bottom wall temperature $T_{bottom}$ and the liquid temperature corresponding to the maximal density $T_M$. An analytical scaling law for heat transfer using the piecewise background temperature is obtained: $Nu=(1/8)(1-T_M)^{5/3}Ra^{1/3},\quad Ra\rightarrow \infty$. For $T_M=0$, although previous work claimed that penetrative convection is almost identical to the Rayleigh–Bénard convection (Wang et al. 2019), under the marginal stability assumption, the present work indicates that they are different as $Ra\rightarrow \infty$ and their difference is more conspicuous when $T_M>0$.

A variational problem is then constructed to derive a conditional lower bound for heat transfer under the marginal stability constraint. When $T_M\leqslant 0.1$, the variational problem shows that the marginally stable modes undergo a hierarchical bifurcation as $Ra$ increases, which yields a background temperature profile of two boundary layer structures: one located near the bottom wall and the other near the top wall. For $T_M>0.1$, only a single mode is detected for all $Ra$ and the background temperature only has one boundary layer structure nested to the bottom wall. Unlike the Rayleigh–Bénard problem, in which the variational problem can yield a physically meaningful background temperature profile (Wen et al. 2022), the present study shows that the conditional lower bound is not physically meaningful because the global heat transfer is not conserved. One potential approach to derive physical meaningful bound on heat transport in penetrative convection (either upper bound or lower bound) is the auxiliary function method (Tobasco, Goluskin & Doering 2018), which deserves future exploration. In the present study, to incorporate more physical information, we invoked the quasilinear approach and solved the ensuing problem using a Newton method with branch continuation. For all $T_M$, we observed that the marginal modes undergo a hierarchical bifurcation as $Ra$ increases and the background temperature has two boundary layer structures. It is interesting that, when the base state $(\boldsymbol {0},1-z)$ is subcritically unstable, the background temperature profile can be non-unique as $Ra< Ra_{c}$ ($Ra_c$ being the critical Rayleigh number of the base state), and the sophisticated time-stepping algorithm as developed in O'Connor et al. (2021) and Wen et al. (2022) cannot be used therein. For $T_M\geqslant 0.3$, quasilinear approach agrees well with the analytical scaling but with a slightly lower coefficient: $Nu\approx 0.17(1-T_M)^{5/3}Ra^{1/3}$. For $T_M<0.3$, we observe that the scaling yielded from the quasilinear approach is not in good agreement with the analytical scaling, which perhaps is because $Ra$ is not large enough. It is clear that the marginally stable background field $\tau$ is non-unique as it can either be constructed using the variational problem or the quasilinear approach. But it should be stressed that, although $\tau$ is different in shape, the $Ra^{1/3}$ scaling remains the same.

Last, to test the effect of truncated nonlinear terms, DNS were performed by utilising the quasilinear problem as initial condition. The DNS show that the nonlinear terms are over-stabilising the assumed marginally stable field at high $Ra$, which suppresses heat transfer in penetrative convection. Thus, our present study is in excellent agreement with the hypothesis of Malkus (1956) in the turbulent flow regime, whence it is questionable if the $Nu\sim Ra^{1/2}$ scaling can be observed in penetrative convection as $Ra\rightarrow \infty$. A possible way is to break the thermal boundary layer using rough walls. However, it should be pointed out that the fluid is always stably stratified near the top wall in penetrative flow, which should impede the transition into the ultimate regime. The mirror symmetry was either imposed or removed in the DNS, aiming at examining the existence of different turbulent states. We demonstrated that spontaneous symmetry breaking can occur at high $Ra$ (e.g. $T_M=0.7$, $Ra=10^{10}$), leading to different statistics of the nonlinear system. Such a phenomenon indicates that the marginally stable background field can be non-unique at high $Ra$ and the turbulent states in penetrative convection are degenerate. For instance, a horizontal ‘wind’ flow was observed in direct numerical simulation without mirror symmetry, which indicates that a time-periodic background flow might also satisfy the marginal stability assumption. More interestingly, the Prandtl number's influence on heat transfer can be examined using the unsteady background field.