2. Problem formulation

We consider an incompressible Newtonian liquid layer of mean thickness $d^*$, constant dynamic viscosity $\mu ^*$ and thermal conductivity $k^*_{th}$ where superscript $(*)$ signifies a dimensional quantity. The layer is supported by a good conducting bottom wall and exposed to an inert gas at the free non-deformable surface, as shown schematically in figure 1. The density $\rho$ is assumed to be a linear function of temperature as follows:

(2.1)\begin{equation} \rho^*=\rho^*_0 [1-\alpha^* (T^*-T^*_0)], \end{equation}

where $\alpha ^*$ is the coefficient of thermal expansion and $T^*_0$ is the reference temperature. The entire system, consisting of the liquid layer, bottom wall and inert gas, is subjected to the VTG, $-\beta ^*$ and HTG, $-\eta ^*$, forming an OTG. Here, VTG and HTG are perpendicular and tangential to the bottom wall planes, respectively, as illustrated in figure 1. The imposed VTG is $\beta ^* = {q^* (T^*_w - T^*_\infty )}/({k^*_{th}+q^* d^*})$ where $q^*, T^*_w$ and $T^*_\infty$ are the convective heat transfer coefficient at the air–liquid interface, dimensional temperature of the bottom wall and the ambient gas, respectively. Note that in the absence of the HTG, the problem considered here reduces to the canonical Rayleigh–Bénard convection problem. Here, we assume a negative VTG implying $T^*_w > T^*_\infty$.

Figure 1. A schematic of the system considered here. The liquid layer is supported by a perfectly conducting bottom wall at $y=0$ and is subjected to an OTG. The dimensional VTG component is $-\beta ^*$, and the dimensional HTG component is $-\eta ^*$, each shown by an arrow. The HTG induces the Hadley circulation, while the imposed VTG leads to the density stratification.

We also assume the aspect ratio of the system to be small such that $d^*/L^*\ll 1$, where $L^*$ is the length of the container along the $x$ and $z$ directions. This allows us to analyse the stability of the core parallel flow existing far from the sidewalls and neglect the effect of sidewalls in agreement with the previous studies (Mercier & Normand Reference Mercier and Normand1996; Patne & Oron Reference Patne and Oron2022). The length, time, velocity, pressure and temperature are scaled by $d^*$, $d^{*2}/\kappa ^*$, $\kappa ^*/d^*$, $\mu ^* \kappa ^*/d^{*2}$ and $\beta ^* d^*$, respectively. The dimensionless VTG is ($-1$) and HTG is ($-\eta =-\eta ^*/\beta ^*$).

Let the scaled fluid velocity components be $\boldsymbol {v}=(v_x,v_y,v_z)$ with $v_j$ being the respective velocity component in the $j$ direction. On imposing the above assumptions in conjunction with the Boussinesq approximation, the dimensionless continuity and momentum conservation equations are given as

(2.2a)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}=0, \end{gather}

(2.2b)\begin{gather}\frac{1}{Pr} \left[ \partial_t \boldsymbol{v} + (\boldsymbol{v} \boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{v} \right] ={-} \boldsymbol{\nabla} p + RaT + \nabla^2 \boldsymbol{v}, \end{gather}

where $Pr={\mu ^*}/{\rho ^*_0 \kappa ^*}$ and $Ra={\rho ^*_0 \alpha ^* \beta ^* d^{*4} g^*}/{\mu ^*\kappa ^*}$ are the Prandtl and Rayleigh numbers, respectively. The gradient and Laplacian operators are $\boldsymbol {\nabla }=( \partial _ x,\partial _y,\partial _z )$ and $\nabla ^2 \equiv \partial ^2_ x+\partial ^2_ y+ \partial ^2_z$, respectively. Also, $\partial _j$ denotes the partial derivative with respect to $j=x,y,z,t$ and $p$ is the pressure. The scaled heat advection–diffusion equation is

(2.2c)\begin{equation} \partial_t T + \boldsymbol{v} \boldsymbol{\cdot} T = \nabla^{2} T. \end{equation}

Equations (2.2) are subjected to the following boundary conditions. At the bottom wall ($y=0$), we assume no-slip, an impermeable bottom wall and a linearly decreasing temperature in the $x$ direction. The thermal boundary condition arises on account of the consideration of a perfectly conducting bottom wall. Thus,

(2.3a)\begin{equation} v_x=0 ; \quad v_z=0 ; \quad v_y=0 ; \quad T=T_w - \eta x. \end{equation}

At the gas–liquid interface ($y=1$), we impose a kinematic boundary condition, continuity of heat flux and the vanishing tangential stresses, as follows:

(2.3b)\begin{gather} v_y=0, \end{gather}

(2.3c)\begin{gather}\partial_y v_x =0, \end{gather}

(2.3d)\begin{gather}\partial_y v_z =0, \end{gather}

(2.3e)\begin{gather}\partial_y T + Bi(T-T_\infty +\eta x) =0. \end{gather}

Here, $Bi={q^* d^*}/{k^*_{th}}$, $T_\infty$ and $q^*$ are the Biot number, the dimensionless temperature of the ambient gas and the heat transfer coefficient of the heat convection at the interface, respectively. It must be noted that the assumption of a non-deformable gas–liquid interface allows us to provide the boundary conditions in the above-simplified form.

Following the procedure of Mercier & Normand (Reference Mercier and Normand1996) and Patne & Oron (Reference Patne and Oron2022), the fully developed, steady-state core flow is

(2.4a)\begin{gather} \bar{v}_x= {-} \frac{\eta Ra}{48}(8y^3-15y^2+6y); \quad \bar{v}_y=0; \quad \bar{v}_z=0, \end{gather}

(2.4b)\begin{gather}\bar{p}=p_a + Ra \int_0^1 \bar{T}(x,y)\, {{\rm d}y}, \end{gather}

(2.4c)\begin{gather}\bar{T} = T_w -\eta x - \left[1 + \frac{Bi \, \eta^2 \, Ra}{320 (1 + Bi)} \right] y + \frac{\eta^2 \, Ra \, y^3 }{960} (20 -25y + 8 y^2). \end{gather}

A comparison of the above base state temperature with that of Patne & Oron (Reference Patne and Oron2022) shows that there is an extra term $({-Bi \, \eta ^2 \, Ra}/{320 (1 + Bi)}) y$ entering as a linear term in $y$. Additionally, similar to Patne & Oron (Reference Patne and Oron2022), there is a quintic polynomial term in $y$, namely $({\eta ^2 \, Ra \, y^3 }/{960}) (20 -25y + 8 y^2)$. These terms will be absent for $\eta =0$ and are functions of $y$ but not of $x$ (i.e. the imposed direction of the HTG), implying that these VTG terms are induced by the imposed HTG. Henceforth, the terms $({Bi \, \eta ^2 \, Ra}/{320 (1 + Bi)}) y$ and $({\eta ^2 \, Ra \, y^3 }/{960}) (20 -25y + 8 y^2)$ will be referred to as induced linear VTG and induced quintic VTG to differentiate from the imposed VTG and each other. Furthermore, the induced linear VTG term possesses the same sign as the imposed VTG, reinforcing the imposed VTG. In a similar setting, the previous studies of Patne, Agnon & Oron (Reference Patne, Agnon and Oron2020b, Reference Patne, Agnon and Oron2021a,Reference Patne, Agnon and Oronb) concerning the thermocapillary instabilities exhibited by a liquid layer predict the existence of an induced VTG which opposes the imposed VTG. Their analysis predicted a strong stabilising effect of the induced VTG on the instabilities caused by the imposed VTG. However, in the present case, the induced VTG reinforces the imposed VTG. This leads to anticipation that this extra VTG could lead to an interesting alteration of the liquid layer dynamics, which is indeed the case as discussed in § 4.

In the subsequent discussion, we will analyse the linear stability of the above base state. Thus, infinitesimally small perturbations are imposed on the base state (2.4). Squire's theorem does not apply to the system under consideration; thus, we consider three-dimensional perturbations. The governing equations (2.2) are linearised around the base state (2.4). In the linearised governing equations, we substitute normal modes of the form

(2.5)\begin{equation} f'(x,t)=\tilde{f}(y) \exp[{\rm i} (k x+m z)+s t]. \end{equation}

Here, $f'(x,t)$ is a perturbation to the dynamic quantity $f(x,t)$ and $\tilde {f}(y)$ is the corresponding eigenfunction in the Laplace–Fourier space. The parameters $k$ and $m$ are the wavenumbers in the $x$ and $z$ directions, respectively. The temporal stability of the base flow is determined by the complex frequency $s=s_r +i s_i$. The base flow is unstable if at least one eigenvalue possesses $s_r>0$. The linearised equations after substitution of the normal modes (2.5) lead to

(2.6a)\begin{gather} {\rm i}k\tilde{v}_x + D\tilde{v}_y + {\rm i}m \tilde{v}_z=0, \end{gather}

(2.6b)\begin{gather}\frac{1}{Pr} \left[ s \tilde{v}_x + {\rm i}k {\bar{v}}_x \tilde{v}_x + \tilde{v}_y D {\bar{v}}_x \right] ={-} {\rm i}k \tilde{p} + (D^2 - k^2 - m^2) \tilde{v}_x, \end{gather}

(2.6c)\begin{gather}\frac{1}{Pr} \left[ s \tilde{v}_y + {\rm i}k {\bar{v}}_x \tilde{v}_y \right] ={-} D \tilde{p} + (D^2 - k^2 - m^2) \tilde{v}_y + Ra\,\tilde{T}, \end{gather}

(2.6d)\begin{gather}\frac{1}{Pr} \left[ s \tilde{v}_z + {\rm i}k {\bar{v}}_x \tilde{v}_z \right] ={-} {\rm i}m \tilde{p} + (D^2 - k^2 - m^2) \tilde{v}_z, \end{gather}

(2.6e)\begin{gather}s \tilde{T} + {\rm i}k {\bar{v}}_x \tilde{T} + \partial_x \, {\bar{T}} \, \tilde{v}_x + \partial_y \, {\bar{T}} \, \tilde{v}_y= (D^2-k^2 - m^2) \tilde{T}, \end{gather}

where $D \equiv {{\rm d}}/{{\rm d}y}$.

The above (2.6) are to be solved using the following boundary conditions. At $y=0$, the boundary conditions become

(2.7a)\begin{equation} y=0: \quad \tilde{v}_x=0; \quad \tilde{v}_y=0; \quad \tilde{v}_z=0; \quad \tilde{T}=0. \end{equation}

At $y=1$, the assumption of a non-deformable interface leads to

(2.7b)\begin{gather} \tilde{v}_y=0, \end{gather}

(2.7c)\begin{gather}D\tilde{v}_x = 0, \end{gather}

(2.7d)\begin{gather}D\tilde{v}_z = 0, \end{gather}

(2.7e)\begin{gather}D\tilde{T} + Bi\tilde{T} =0. \end{gather}

Equations (2.6) along with boundary conditions (2.7) form an eigenvalue problem for $\omega$ for a specified value of the parameters $Ra, Bi, Pr, \eta$ and wavenumbers $k$ and $m$. To solve this eigenvalue problem, we employ the pseudospectral method, briefly explained below.

3. Numerical methodology

In the pseudospectral method, the eigenfunctions are expanded in the form of a series of Chebyshev polynomials

(3.1)\begin{equation} \tilde{f}(y)=\sum_{m=0}^{m=M} a_m T_m (y), \end{equation}

where $M$ is the highest degree of the polynomial in the series expansion and $T_m(y)$ are Chebyshev polynomials of degree $m$. The parameter $M$ is also termed as the number of collocation points. The series coefficients $a_m$ are the unknowns to be solved for.

To implement the pseudospectral method, we transform the fluid domain $0 \leqslant y \leqslant 1$ to $-1 \leqslant y \leqslant 1$ by using mapping $y \rightarrow 2y-1$. Upon discretisation using the above procedure, the eigenvalue problem becomes

(3.2)\begin{equation} \boldsymbol{Pe} + s \boldsymbol{Q}\boldsymbol{e}=0, \end{equation}

where $\boldsymbol {P}$ and $\boldsymbol {Q}$ are matrices derived following the discretisation procedure and $\boldsymbol {e}$ is the vector containing the coefficients of all series expansions.

The standard procedure to discretise the governing equations and boundary conditions using Chebyshev polynomials can be found in Trefethen (Reference Trefethen2000) and Schmid & Henningson (Reference Schmid and Henningson2001). The application of the pseudospectral method for similar problems can be found in Patne, Agnon & Oron (Reference Patne, Agnon and Oron2020a); Patne et al. (Reference Patne, Agnon and Oron2020b, Reference Patne, Agnon and Oron2021a,Reference Patne, Agnon and Oronb), Patne & Oron (Reference Patne and Oron2022), Patne & Chandarana (Reference Patne and Chandarana2023) and Patne (Reference Patne2024). To solve the eigenvalue problem (3.2), we use the eig MATLAB routine. The predicted eigenspectrum also contains numerical spurious eigenvalues. To filter out these modes, we execute the code for $N$ and $N+2$ collocation points, and the eigenvalues are compared with a specified tolerance, e.g. $10^{-4}$. To further ascertain the genuine eigenvalues, the number of collocation points is increased by $25$, and the variation of the eigenvalues is observed. If the eigenvalue does not change more than a specified precision, e.g. to the sixth significant digit, the same number of collocation points is used to determine the critical parameters of the system. In the present work, $N=60$ is found to be sufficient to achieve convergence and to determine the leading, most unstable eigenvalue within the investigated parameter range.

The numerical method employed here is validated as follows. The critical $Ra_c$ value for a liquid layer supported by a perfectly conducting bottom wall and heated from below (i.e. only VTG) and exposed to inert ambient gas is available in Drazin (Reference Drazin2002) (p. 99). Drazin (Reference Drazin2002) reports $Ra_c=1101$ and the critical wavenumber $k_c=2.682$. The corresponding problem in Drazin (Reference Drazin2002) differs due to the thermal boundary condition at the gas–liquid interface or free surface. He assumes a fixed temperature at the bottom wall and free surface as well. To validate the present numerical methodology, we change the boundary condition at the free surface to a fixed temperature instead of Newton's cooling law (2.7e) and substitute $\eta =0$. The resulting problem predicts $Ra_c=1101$ and $k_c=2.681$ in excellent agreement with Drazin (Reference Drazin2002), thereby validating the numerical methodology for $\eta ^* = 0$.

In the extreme limit of $\eta ^* \rightarrow \infty (\beta ^* =0)$, the present problem corresponds to the stability of Hadley circulation studied by Hart (Reference Hart1972), Gill (Reference Gill1974), Hart (Reference Hart1983), Laure & Roux (Reference Laure and Roux1989) with altered boundary conditions. Thus, to validate the present numerical methodology, we compare the results obtained using the pseudospectral method utilised here with those of Hart (Reference Hart1983) and Laure & Roux (Reference Laure and Roux1989) as follows.

Hart (Reference Hart1983) and Laure & Roux (Reference Laure and Roux1989) consider two cases based on the boundary conditions at the top and bottom of the liquid layer. Considering the present problem of having a free upper surface and a rigid bottom wall, we validate our results with their rigid-free case. Given the vast difference in the scaling and boundary conditions, we have directly utilised base state equations (5) and (6) and the linearised perturbation equations (10)–(12) of Hart (Reference Hart1983). Note that, for the parallel core flow far from the sidewalls and shallow cavity, the constant $a$ in his base state equations (5) and (6) is unity. Also, we have utilised equations (10)–(12) of Hart (Reference Hart1983) since we wish to reproduce select points from figure 3 of Laure & Roux (Reference Laure and Roux1989) for the spanwise (longitudinal) modes of Laure & Roux (Reference Laure and Roux1989). Regarding the boundary conditions, following Hart (Reference Hart1983) and Laure & Roux (Reference Laure and Roux1989), we assume an insulator at the bottom and a stress-free insulator at the upper free surface. Note that, in the present problem, we consider a perfectly conducting bottom wall and energy exchange at the free surface via convection.

Hart (Reference Hart1983) and Laure & Roux (Reference Laure and Roux1989) utilise the Grashof number, $Gr=Ra_H \, Pr$, to characterise the inception of instability where $Ra_H$ is the horizontal Rayleigh number based on the imposed dimensional HTG, $\eta ^*$. To validate our numerical methodology, we utilise three points digitally extracted from figure 3 of Laure & Roux (Reference Laure and Roux1989), viz., $(Pr, Gr_c)=(0.001,16450),(0.01,5310),(0.1,2054)$ where $Gr_c$ is the critical Grashof number. Upon using the pseudospectral method, we find for $Pr=0.001$, $Gr_c= 16455$, for $Pr=0.01$, $Gr_c = 5313$ and for $Pr=0.1$, $Gr_c \sim 2052$ in excellent agreement with figure 3 of Laure & Roux (Reference Laure and Roux1989) (curve for the spanwise or longitudinal modes) thereby validating the present numerical methodology for $\beta ^* = 0$.

Additionally, the present problem differs from the problem solved by Patne & Oron (Reference Patne and Oron2022) due to the boundary condition at the bottom wall, which also modifies the base state. Thus, if we switch the base state and boundary condition to that of Patne & Oron (Reference Patne and Oron2022), we should be able to reproduce their results, which is indeed the case. We obtained the variation of the eigenvalue $\omega ={\rm i}s$ with $\eta$ in the $\omega _r$–$\omega _i$ space at $Bi=10^{-5}, m=0, k=0.1, Ra=400$ and $Pr=10$ for the problem of Patne & Oron (Reference Patne and Oron2022) and compared with their figure 2. The comparison showed an excellent agreement thereby validating the present numerical procedure.

5. Perturbation energy budget analysis

In § 4, the stability analysis predicted new streamwise and spanwise modes for $Bi \neq 0$ and $Pr=7$. To understand the factors responsible for the existence of these modes, here we analyse the perturbation energy budget using the approach of Hu, He & Chen (Reference Hu, He and Chen2016), Hu et al. (Reference Hu, He, Chen and Liu2017), Patne et al. (Reference Patne, Agnon and Oron2021a) and Patne & Oron (Reference Patne and Oron2022). To accomplish this, the momentum conservation equations (2.2b) are linearised around an arbitrary base state to yield

(5.1)\begin{equation} \frac{1}{Pr} \frac{\partial \boldsymbol{v}^\prime}{\partial t} ={-}\boldsymbol{\nabla} p^\prime + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau}^\prime + Ra T' \boldsymbol{\nabla} y - \frac{1}{Pr} \left[ (\boldsymbol{v'} \boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{\bar{v}} + (\boldsymbol{\bar{v}} \boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{v'} \right], \end{equation}

where $\boldsymbol {\tau }$ is the stress tensor and a prime, $'$, implies a perturbation quantity. To obtain the perturbation energy equation, (5.1) is scalar multiplied by the perturbation velocity vector ${\boldsymbol {v}}^\prime$. The resulting equation is then integrated over the liquid layer to yield

(5.2)\begin{align} \frac{1}{Pr} \frac{\partial E}{\partial t} &={-} \frac{1}{2} \int \boldsymbol{\tau}^\prime : \boldsymbol{\dot{\gamma}^\prime} \, {\rm d}V + Ra \int T' \, v'_y \, {\rm d}V - \frac{1}{Pr} \int \boldsymbol{v'} \boldsymbol{\cdot} \boldsymbol{\bar{\dot{\gamma}} } \boldsymbol{\cdot} \boldsymbol{v'} \, {\rm d}V \nonumber\\ &\equiv{-} I_b + I_{RB} - I_R, \end{align}

(5.3)\begin{equation} E= \tfrac{1}{2}\int{\boldsymbol{v}}^\prime \boldsymbol{\cdot}{\boldsymbol{v}}^\prime\,{\rm d}V. \end{equation}

Here, the integrals $I_b, I_{RB}$ and $I_{R}$ are the bulk stress work, Rayleigh–Bénard integral or buoyancy work, and the Reynolds stress work (Drazin Reference Drazin2002), respectively, and ${\rm d}V$ is the volume element. The quantities $\boldsymbol {\dot {\gamma }} ^\prime = \boldsymbol {\nabla } \boldsymbol {v}^\prime + \boldsymbol {\nabla } {\boldsymbol {v}^\prime }^T$ and $\boldsymbol {\bar {\dot {\gamma }} }= \boldsymbol {\nabla } \boldsymbol {\bar {v}} + \boldsymbol {\nabla } \boldsymbol {\bar {v}}^T$ are the strain-rate tensors in the perturbed and base states, respectively. It must be noted that while obtaining the above equation, we have used the continuity equation, $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {v'}=0$. The pressure perturbation energy term vanishes at the bottom wall and free surface due to the vanishing normal velocity perturbations.

The bulk stress work $I_b = \frac {1}{2} \int \boldsymbol {\tau }^\prime : \boldsymbol {\dot {\gamma }^\prime } \, {\rm d}V$ represents viscous dissipation of the perturbation energy. For a Newtonian fluid, $\tau '_{ij} = \gamma '_{ij}$, thus, $I_b = \frac {1}{2} \boldsymbol {\tau }^\prime : \boldsymbol {\dot {\gamma }}^\prime = \frac {1}{2} \gamma '_{ij}\gamma '_{ji} \geqslant 0$. This leads to the immediate conclusion that the bulk stress work will always lead to a decrease in the perturbation energy for a Newtonian fluid. The same is not the case for viscoelastic fluids as predicted by Hu et al. (Reference Hu, He and Chen2016, Reference Hu, He, Chen and Liu2017) and Patne et al. (Reference Patne, Agnon and Oron2021a, Reference Patne, Agnon and Oron2020a) where $I_b$ could lead to new instabilities.

The Rayleigh–Bénard integral or buoyancy work $I_{RB}$ arises purely due to the buoyancy. Typically, for a liquid layer subjected to a purely negative VTG, causing unstable density stratification, $I_{RB}$ is positive, thereby leading to the VTG mode. To accomplish that, $I_{RB}$ must overcome the stabilising impact of the dissipation factors, such as viscous forces and thermal diffusion of the energy.

The third integral, viz., the Reynolds stress work ($I_{R}$), arises purely as a consequence of the imposed HTG. For a liquid layer subjected to a purely VTG, the base state is a vanishing velocity field. In that case $I_{R}=({1}/{Pr} )\int \boldsymbol {v'} \boldsymbol {\cdot } \boldsymbol {\bar {\dot {\gamma }} } \boldsymbol {\cdot } \boldsymbol {v'} \, {\rm d}V$ will vanish since $\boldsymbol {\bar {\dot {\gamma }} } = 0$. Thus, it does not contribute to the VTG mode. However, in the present case, due to the HTG component of the imposed OTG, the base state velocity in the $x$ direction is non-zero, thus $\boldsymbol {\bar {\dot {\gamma }} } \neq 0$. The integral $I_{R}$ represents a volume-averaged correlation between the perturbations in the horizontal and vertical components of the velocity field. This integral is also responsible for the energy exchange between the base state and the perturbation fields. The presence of ${1}/{Pr}$ in the integral $I_{R}$ indicates that for low $Pr$, this term will play a major role in determining the stability of the layer, which is indeed the case as discussed below. Furthermore, the base state strain-rate tensor $\boldsymbol {\bar {\dot {\gamma }} }$ directly depends on $\eta$. Thus, an increasing $\eta$ will lead to an increase in $\boldsymbol {\bar {\dot {\gamma }} }$ thereby enhancing the importance of the Reynolds stress work in influencing the stability of the layer.

Before proceeding with the numerical calculations of the integrals, the following conclusions can be drawn by using the definitions of the integrals. The absolute values of the integrals may not provide the exact roles of the different terms. Thus, we normalise the integrals by the absolute value of $I_{RB}$. From (5.2), we have

(5.4)\begin{equation} \frac{I_b}{I_{RB}} = \frac{1}{2Ra} \frac{ \displaystyle\int \boldsymbol{\tau}^\prime : \boldsymbol{\dot{\gamma}^\prime } \, {\rm d}V}{ \displaystyle\int T' \, v'_y \, {\rm d}V}. \end{equation}

As explained later in § 6, for the buoyancy convection to exist, the destabilising forces, i.e. Raleigh–Bénard work term, must overcome the viscous dissipation forces represented by $I_b$. From above expression, as $Ra \rightarrow \infty$, ${I_b}/{I_{RB}} \rightarrow 0$. Physically, this implies that as $Ra$ increases, the Rayleigh–Bénard work term overpowers the viscous dissipation term to result in an instability, which is indeed the case for $\eta =0$.

The other ratio is

(5.5)\begin{equation} \frac{I_R}{I_{RB}} = \frac{1}{Pr \, Ra} \frac{ \displaystyle\int \boldsymbol{v'} \boldsymbol{\cdot} \boldsymbol{\bar{\dot{\gamma}} } \boldsymbol{\cdot} \boldsymbol{v'} \, {\rm d}V}{ \displaystyle\int T' \, v'_y \, {\rm d}V}. \end{equation}

Only $\bar {\dot {\gamma }}_{xy}$ and $\bar {\dot {\gamma }}_{yx}$ components of the base state strain-rate are non-zero. Thus, the above equation using the base state velocity (2.4a) modifies to

(5.6)\begin{equation} \frac{I_R}{I_{RB}} ={-}\frac{\eta}{4 Pr} \frac{ \displaystyle\int (y-1) (4y-1) v'_x v'_y \, {\rm d}V}{ \displaystyle\int T' \, v'_y \, {\rm d}V}. \end{equation}

This implies that an increase in $Ra$ will not directly affect the value of this ratio. Instead, it is the strength of the imposed HTG $\eta$ that will directly affect this ratio. Additionally, the presence of $Pr$ in the denominator implies that cases with the $Pr<1$ will have a better chance of being destabilised by the Reynolds stress work $I_R$. This is indeed the case from figures 8 and 9 where $Pr<1$ leads to the destabilisation of the VTG mode unlike in the case of $Pr=7$. The volume-averaged correlation $\int (y-1) (4y-1) v'_x v'_y \, {\rm d}V$ must be positive to contribute to the growth of the perturbations via (5.2).

We focus on the new streamwise mode and try to understand the impact of various terms in (5.2) in introducing the new modes. The values of the integrals normalised by the Rayleigh–Bénard integral for a set of parameters corresponding to the new streamwise mode in an unstable regime are given in table 1. The bulk stress work ($I_b$) is positive and thus will lead to the dissipation of the perturbation energy. The Reynolds stress work ($I_R$), on the other hand, is always negative, implying that it will lead to the enhancement of the perturbation energy. For $\eta =1$, both normalised $I_b$ and $I_R$ are less than unity; thus, the Rayleigh–Bénard work term ($I_{RB}$) is larger than both the terms, indicating the major role of $I_{RB}$ in introducing the new mode of instability. It must be noted that the integral $I_{RB}$ could achieve such a high value as a direct result of the presence of the induced linear VTG term. This term, as noted in § 4, counters the stabilising influence of the induced quintic VTG term even for $\eta <1$ (see figure 6). The base-state velocity (2.4a) is proportional to $\eta$; thus, at low $\eta$, the integral $I_R$ would have a lower contribution. As $\eta$ increases, the value of $I_R$ also increases, and for $\eta =100$, it is two times larger than $I_{RB}$. Thus, while $I_{RB}$ plays a vital role in introducing the new mode, $I_R$ is the one which drives the new mode at high $\eta$.

Table 1. Values of the bulk stress work ($I_b$) and Reynolds stress work ($I_R$) normalised by the buoyancy work ($I_{RB}$) in the energy balance equation (5.2) for the streamwise new mode at $Bi=10, k=7, m=0$ and $Pr=7$. In each case, the streamwise new mode is unstable. As expected, the normalised bulk stress work is positive, implying its stabilising influence. The normalised Reynolds work term is negative, thereby having a destabilising influence. The Rayleigh–Bénard work term is found to be positive for the above combination of parameters.

It must be noted that the perturbation energy analysis carried out here can only provide qualitative information about perturbation energy sources/sinks that may lead to destabilisation (Drazin Reference Drazin2002). However, the correct critical parameter values may be obtained only using the general linear stability analysis discussed in § 4.

6. Physical mechanism

From figure 4, the VTG mode is completely stabilised by the imposed HTG for vanishing $Bi$ and $Pr>1$. While for $Pr<1$ and $Bi \sim 0$, the imposed HTG may have a stabilising effect at the beginning, but eventually, it leads to destabilisation of the VTG mode (see figures 8 and 9). For $Bi \neq 0$, new streamwise and spanwise modes appear as shown in figures 6 and 12. These modes do not exist when a liquid layer is subjected to a purely VTG but arise due to an interaction between the imposed VTG and HTG. This physical manifestation of the stabilising/destabilising effect of the imposed HTG remains to be understood, which is precisely the aim of this section.

A liquid layer subjected to a purely negative VTG (i.e. the bottom wall is at a higher temperature than the ambient gas) exhibits Rayleigh–Bénard convection due to the lighter liquid at the bottom and heavy liquid at the top owing to the temperature dependence of the liquid density. One can describe this convection phenomenon in terms of liquid movement as follows. Consider a liquid parcel located at $y=y_0$. Assuming the density of the liquid is a linearly decreasing function of the temperature, the parcel located at a higher level $y=y_0+\Delta y$ will have a higher density than that of the liquid layer present at $y=y_0$. Here $\varDelta \ll 1$ and a positive quantity. Due to the imposed gravity field, the liquid parcel at $y=y_0$ will try to move to $y=y_0+\Delta y$. The viscous forces and thermal diffusion of the energy from the parcel to the surrounding liquid will oppose this movement. If the imposed VTG strength is strong enough, the parcel will overcome the viscous and thermal energy diffusion resistance and move to $y=y_0+\Delta y$, resulting in the buoyancy convection. For the convenience of discussion, this section will be divided into two subsections: (i) $Bi \ll 1$ and (ii) $Bi \neq 0$.

6.1. For $Bi \ll 1$

6.1.1. For $Pr>1$

From figure 4, the VTG mode is stabilised for $\eta > 0.22$. Patne & Oron (Reference Patne and Oron2022) discussed the stabilisation of the longwave mode that originates due to the imposed VTG. Their analysis indicates the presence of an opposing VTG induced by the imposed HTG, which leads to the predicted stabilisation. From base state temperature (2.4c), at $Bi \rightarrow 0$, there is a linear term in $y$ with coefficient ($-1$) and a quintic polynomial term, viz., induced quintic VTG term. The average total VTG can then be obtained by differentiating base state temperature (2.4c) with respect to $y$ and integrating it over the fluid domain in the $y$ direction to obtain

(6.1)\begin{equation} \left(\frac{\partial {\bar{T}} }{\partial y} \right)_{avg} = \frac{\eta^2 Ra}{320}-1, \end{equation}

where the subscript $avg$ indicates the average temperature gradient. The positive term ${\eta ^2 Ra}/{320}$ opposes the imposed VTG, thereby weakening it.

The physical manifestation of this stabilisation can be understood as follows. As discussed above, for the buoyancy instability to exist, the imposed VTG must be strong enough to overcome the resistance offered by the viscous forces and thermal diffusion away from the fluid parcel. The quintic polynomial term induced by the HTG weakens the imposed VTG, thus adding more resistance to the development of buoyancy instability. The weakened VTG fails to cause instability, thus the predicted stabilisation. At a sufficiently high $\eta$, the average temperature gradient will become positive, thereby removing the possibility of the existence of buoyancy instability.

6.1.2. For $Pr<1$

As discussed in § 4, for $Pr=0.1$ and $Pr=0.01$ (see figures 8 and 18), the imposed HTG fails to stabilise the VTG mode even for $Bi \ll 1$. The failure to stabilise the VTG mode could be understood from the perturbation energy analysis discussed in § 5 as follows.

For $Bi \ll 1$, the induced linear VTG term does not affect the stability of the layer. Thus, the Rayleigh–Bénard work term $I_{RB}$ of (5.2), which plays a major role in introducing the new mode for $Bi \neq 0$, will not be of much relevance. Furthermore, the same $I_{RB}$ term fails to counter the stabilising effect of the imposed HTG through the induced quintic VTG term for $Pr=7$; thus, it is expected that such a term will not be the stability-determining term. In fact, at a sufficiently high $\eta$ from (6.1), owing to the induced quintic VTG term, the VTG will reverse the sign, resulting in a stable density stratification. In that case, the $I_{RB}$ term will have a stabilising effect. The bulk stress work term $I_b$ always leads to the viscous dissipation of the perturbation energy. The only remaining term that can lead to the growth of the perturbation energy is the Reynolds work term $I_R$. This term must overcome the stabilising effect of $I_b$ and $I_{RB}$ to result in instability. From figures 8 and 18, the $I_{R}$ term succeeds in this endeavour, thereby leading to the predicted results.

From perturbation energy equation (5.2), $I_R$ has a multiplier $1/Pr$. Thus, as $Pr$ decreases, there will be a corresponding increase in the $I_R$ value. The base state velocity (2.4a) also has a multiplier $\eta Ra$. Thus, an increasing $\eta$ and $Ra$ and a decreasing $Pr$ will exacerbate the destabilisation caused by $I_R$, which succeeds for $Pr=0.01$. For $Pr=0.1$, from figure 8, due to a moderate value of $Pr$, the $I_{R}$ term is weak at low $\eta$; thus, the VTG mode is stabilised. However, as $\eta$ and $Ra$ increase, the $I_{R}$ term gains enough strength to counter the stabilising influence of $I_b$ and $I_{RB}$ resulting in the predicted curve.

6.2. For $Bi \neq 0$

From figures 6 and 12, for $Bi \neq 0$, new streamwise and spanwise modes emerge. As noted in § 4, these modes originate due to the presence of the induced linear VTG term in the base-state temperature (2.4c), namely $({-Bi \, \eta ^2 \, Ra}/{320 (1 + Bi)} )y$. The physical mechanism by which the induced linear VTG term achieves this can be explained as follows.

The average total base state VTG for $Bi \neq 0$ is

(6.2)\begin{equation} \left(\frac{\partial {\bar{T}} }{\partial y} \right)_{avg} = \frac{\eta^2 Ra}{320} - \frac{\eta^2 \, Bi \, Ra}{320(1+Bi)}-1. \end{equation}

For low $Bi$, the term ${\eta ^2 \, Bi \, Ra}/{320(1+Bi)}$ will not be significant. However, as $Bi \rightarrow \infty$, the same term modifies to ${\eta ^2 Ra}/{320}$ which is exactly equal to the stabilising term. Thus, the induced linear VTG term can counter the effect of the induced quintic VTG term, provided that $Bi$ is sufficiently high. Once this induced quintic VTG term effect nullification occurs, the imposed VTG introduces the unstable density stratification, leading to the convection. At face value, we are returning to the liquid layer system subjected to a purely negative VTG. It turns out that the liquid layer is hiding a much more interesting dynamics, as described below.

To extract the additional information, we plot the variation of the base state temperature gradient obtained by differentiating (2.4c) with respect to $y$. The result of such an exercise is shown in figure 19 in the parametric regime where the new streamwise mode is unstable. Unlike in the case of a liquid layer subjected to a purely VTG, where ${\partial \bar {T}}/{\partial y} <0$ throughout the liquid layer, there are two regions with ${\partial \bar {T}}/{\partial y} <0$. Thus, there are two regions which could potentially lead to the existence of the new mode. This raises the question, which of these two regions is responsible for the existence of the predicted new streamwise mode?

Figure 19. The base state temperature gradient variation with $y$ at $Bi=10, \eta =10$ and $Ra=900$ where the streamwise new mode of instability is unstable. There are two regions with ${{\rm d}\bar {T}}/{{\rm d}y} <0$, i.e. with unstable density stratification.

The analysis of Patne & Oron (Reference Patne and Oron2022) for a similar problem with a poorly conducting bottom wall predicted a new mode of instability at low $Pr$. Their analysis shows that the reverse flow shown in the schematic 1, near the bottom wall, could be responsible for the predicted new mode. In the present case, there is also the presence of reverse flow, which plays a role in causing an unstable stratification near the bottom wall shown in figure 19. However, as pointed out previously, the term $({-Bi \, \eta ^2 \, Ra}/{320 (1 + Bi)}) y$ could be responsible for the existence of the new mode of instability. This term contains the Biot number $Bi$, which enters through the thermal boundary condition at the free surface at $y=1$, and the second region with unstable stratification exists near the free surface. These arguments indicate that the unstably stratified region near the free surface is responsible for the existence of the new mode.

The requirement of sufficiently high $Bi={q d}/{k_{th}}$ for the existence of the new mode also indicates that there is a need for efficient heat transfer from the liquid layer to the ambient gas, which then will result in the formation of a strong negative temperature gradient. As $Bi$ increases, the corresponding negative temperature gradient will also become severe, thereby exacerbating the instability and lowering $Ra_c$. This does not imply that $Bi \rightarrow \infty$ will lead to $Ra_c \rightarrow 0$ for sufficiently high $\eta$. For $Bi \gg 1$, it can be eliminated from the term ${-Bi \, \eta ^2 \, Ra}/{320 (1 + Bi)}$ by taking the limit $Bi \rightarrow \infty$ to obtain ${\eta ^2 Ra}/{320}$ thereby achieving a finite value. Thus, as $Bi$ increases beyond $10$, it will not lead to a much change in $Ra_c$ at high $\eta$. From perturbation energy budget analysis in § 5, an increasing $\eta$ also enhances the role of the Reynolds stress term, which dominates at high $\eta$ and the induced linear VTG term plays a supporting role.

7. Summary

In the present study, we carry out the linear stability analysis of a liquid layer supported by a perfectly conducting bottom wall and the other non-deformable free surface exposed to an inert ambient gas. The bottom wall, liquid layer and ambient gas are subjected to an OTG. The liquid density is assumed to be a linearly decreasing function of temperature. The imposed OTG consists of a negative VTG component, which leads to the canonical Rayleigh–Bénard convection and a HTG component, which leads to Hadley circulation. Using normal mode analysis, the linearised perturbation equations are reduced to a set of ordinary differential equations, becoming an eigenvalue problem in terms of growth rate $s$. We employ the pseudospectral method (using Chebyshev polynomials) to solve the resulting eigenvalue problem. The numerical predictions are further explained using the perturbation energy budget analysis and physical arguments.

In the absence of the imposed HTG (i.e. $\eta =0$), instability modes exist due to the imposed VTG. All the VTG modes of instability, viz., the streamwise ($k \neq 0$ and $m=0$), spanwise ($k=0$ and $m \neq 0$) and oblique modes ($k \neq 0$ and $m \neq 0$) become unstable at the same value of Rayleigh number $Ra =669$, thus the critical Rayleigh number $Ra_c =669$. Upon switching on the HTG, it leads to an induced linear VTG term ($({-Bi \, \eta ^2 \, Ra}/{320 (1 + Bi)}) y$) and an induced quintic VTG term ($({\eta ^2 \, Ra \, y^3 }/{960}) (20 -25y + 8 y^2)$).

For low $Bi$, the induced linear VTG term can be neglected. For $Pr>1$, the induced quintic VTG term has a strong stabilising effect on the VTG modes. This stabilisation leads to an island of instability in the $k\unicode{x2013}Ra$ and $\eta \unicode{x2013}Ra_c$ planes for $Pr>1$. The physical arguments reveal that the induced quintic VTG term weakens the imposed VTG, which fails to overcome the resistance offered by the viscous dissipation and thermal energy diffusion, thereby increasing $Ra_c$. The spanwise VTG mode exhibits a larger island of instability than the streamwise VTG mode. Thus, the spanwise mode is the dominant mode of instability for $Pr>1$ and low $Bi$.

For $Pr=0.1$ and $Bi \ll 1$, the HTG stabilises the VTG mode for $\eta <0.15$ and the opposite effect for $\eta > 0.15$. In this case, the stabilisation occurs through the induced quintic VTG term. However, as $\eta$ increases beyond $0.15$, the Reynolds stress term, which facilitates the energy transfer from the base flow to the perturbations, becomes relevant. This term then leads to the subsequent destabilisation. At sufficiently high $\eta$, the $Ra_c$ versus $\eta$ curve exhibits $Ra_c \sim 1/\eta$ scaling. For $Pr=0.01$ and $Bi \ll 1$, the HTG has a destabilising effect for an arbitrary value of $\eta$ on the streamwise VTG mode. For the spanwise VTG mode, the HTG exhibits a destabilising effect except for low $\eta$. For $Pr=0.01$ and $Bi \ll 1$, the streamwise mode is the dominant mode of instability.

For $Bi \neq 0$ and $Pr>1$, the induced linear VTG term becomes important. The perturbation energy analysis reveals that the same term reinforces the imposed VTG and helps cancel the stabilising effect of the induced quintic VTG term. This leads to the existence of new modes of instability. These new modes originate from the high $Ra$ part of the neutral stability curves for the VTG mode as $\eta$ increases. Thus, the liquid layer will exhibit instability at an arbitrary $\eta$ in contrast with the case for low $Bi$. The variation of the base state temperature gradient with $y$ reveals the existence of two zones with unstable density stratification, the first near the bottom wall and the second near the free surface. The presence of the Biot number ($Bi$) in the induced linear VTG term indicates that the unstable zone near the free surface could be the main culprit behind the existence of the new modes.

For $Pr=7$, the new spanwise mode becomes more unstable at a lower $Ra$ than the new streamwise mode; thus, the spanwise mode is the dominant mode of instability. At $Pr=0.01$, an increasing $Bi$ stabilises the streamwise mode at low $\eta$ while at high $\eta$, curves in the $\eta \unicode{x2013}Ra_c$ plane for various $Bi$ merge. However, the spanwise mode at $Pr=0.01$ is stabilised by the increasing $Bi$ for the whole range of $\eta$. Thus, the streamwise mode is the dominant mode of instability for $Pr<1$, irrespective of the value of $Bi$.

New streamwise and spanwise modes for $Pr>1$ and finite $Bi$, and VTG mode for $Pr<1$ exhibit the scaling $Ra_c \sim 1/\eta$ in the limit of large $\eta$. Thus, $Ra_c \rightarrow 0$ as $\eta \rightarrow \infty$ for the above-mentioned modes. This implies the existence of an unstable base flow even if the imposed VTG is negligible. The results predicted here could be readily recast in terms of the Rayleigh number based on the imposed HTG, $Ra_H$. By definition, the Rayleigh number based on the imposed VTG and $Ra_H$ are related by $Ra_H = \eta Ra$.

The present analysis thus shows the possibility of using an OTG to intensify the mixing by manipulating buoyancy instabilities. Industrial applications could involve non-Newtonian fluids; thus, the present study could be extended to non-Newtonian fluids. Also, the instabilities predicted in the present study could be experimentally confirmed.