2.1 Governing equations

We will assume that the liquid layer has depth $d_{l}$ and lateral dimensions much larger than $d_{l}$ . Since in most experiments the system is covered by a horizontal plate to limit evaporation, we will assume that the gas layer also has a finite depth $d_{g}$ . To describe convection in this two-layer system, we will use a variation of the two-sided model originally introduced by Qin et al. (Reference Qin, Tuković and Grigoriev2014) for near-atmospheric conditions, and later extended by Qin & Grigoriev (Reference Qin and Grigoriev2015) to the limit when the gas phase is dominated by the vapour, rather than air. A version of the model interpolating between the two limits is summarized below. The momentum transport in the bulk is described, for both the liquid and the gas phase, by the Navier–Stokes equation

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D70C}(\unicode[STIX]{x2202}_{t}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u})=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\unicode[STIX]{x1D70C}(T,c_{a})\boldsymbol{g},\end{eqnarray}$$

where $u$ is the fluid velocity, $T$ is the fluid temperature, $p$ is the fluid pressure, $\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D707}$ are the fluid’s density and viscosity, respectively, $c_{a}$ is the concentration of air, and $\boldsymbol{g}$ is the gravitational acceleration. (The air is non-condensable, so $c_{a}=0$ in the liquid phase.) Following standard practice, we use the Boussinesq approximation, where the density is considered constant everywhere except in the last term on the right-hand side representing the buoyancy force. In particular, the mass conservation equation in both phases reduces to

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0.\end{eqnarray}$$

To account for buoyancy, the density of the liquid phase is assumed to depend linearly on the temperature

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{l}=\unicode[STIX]{x1D70C}_{l}^{0}[1-\unicode[STIX]{x1D6FD}_{l}(T-T_{0})],\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}_{l}^{0}$ is the density at the temperature $T_{0}$ describing global thermodynamic equilibrium, and $\unicode[STIX]{x1D6FD}_{l}=-\unicode[STIX]{x1D70C}_{l}^{-1}\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{l}/\unicode[STIX]{x2202}T$ is the coefficient of thermal expansion. Here and below, subscripts $l$ , $g$ , $v$ , $a$ and $i$ denote properties of the liquid and gas phase, vapour and air component, and the liquid–vapour interface, respectively. The index $0$ will be used throughout the paper to distinguish the equilibrium values from the non-equilibrium ones, where the choice is not obvious. In the gas phase

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{g}=\unicode[STIX]{x1D70C}_{a}+\unicode[STIX]{x1D70C}_{v},\end{eqnarray}$$

where both vapour and air are considered to be ideal gases

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{a,v}=\frac{p_{a,v}}{\bar{R}_{a,v}T},\end{eqnarray}$$

$\bar{R}_{a,v}=R/M_{a,v}$ , $R$ is the universal gas constant, and $M_{a,v}$ is the molar mass of air/vapour. Correspondingly, the total pressure in the gas is

(2.6) $$\begin{eqnarray}p_{g}=p_{a}+p_{v}.\end{eqnarray}$$

Mass transport in the gas phase is described by the standard conservation equation(s)

(2.7) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}n_{a,v}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{j}_{a,v}=0,\end{eqnarray}$$

for the number density of the two components

(2.8) $$\begin{eqnarray}n_{a,v}=\frac{\unicode[STIX]{x1D70C}_{a,v}}{m_{a,v}},\end{eqnarray}$$

where $m_{a,v}=M_{a,v}/N_{A}$ is the mass of one air/vapour molecule. The number density flux is given by

(2.9) $$\begin{eqnarray}\boldsymbol{j}_{a,v}=n_{a,v}\boldsymbol{u}-n_{g}D\unicode[STIX]{x1D735}c_{a,v}=n_{g}(\boldsymbol{u}c_{a,v}-D\unicode[STIX]{x1D735}c_{a,v}),\end{eqnarray}$$

where the first and the second term on the right-hand side represent the contributions due to advection and diffusion, respectively, $D$ is the binary diffusion coefficient, and

(2.10) $$\begin{eqnarray}c_{a,v}=\frac{n_{a,v}}{n_{g}}=\frac{p_{a,v}}{p_{g}}\end{eqnarray}$$

are the concentrations (or, more precisely, the molar fractions) of the two components. Therefore, the conservation equation(s) (2.7) can be rewritten as

(2.11) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}(n_{g}c_{a,v})+n_{g}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c_{a,v}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(n_{g}D\unicode[STIX]{x1D735}c_{a,v}).\end{eqnarray}$$

From the ideal gas law, the total number density in the gas phase

(2.12) $$\begin{eqnarray}n_{g}=\frac{p_{g}}{k_{B}T},\end{eqnarray}$$

where $k_{B}=R/N_{A}$ is the Boltzmann constant. As we showed previously (Qin et al. Reference Qin, Tuković and Grigoriev2015), spatial variation of $p_{g}$ can in practice be neglected. Furthermore, the largest temperature variation in relevant experimental and numerical studies is around $5\,\%$ (and typically much smaller than that), such that $n_{g}$ can be assumed constant. Consequently (2.11) reduces to an advection–diffusion equation for the two concentrations

(2.13) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}c_{a,v}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c_{a,v}=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(D\unicode[STIX]{x1D735}c_{a,v}).\end{eqnarray}$$

These two equations are equivalent, so either one can be used, since $c_{a}+c_{v}=1$ .

Finally, the transport of heat is also described using an advection–diffusion equation

(2.14) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}T+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T=\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FB}^{2}T,\end{eqnarray}$$

where $\unicode[STIX]{x1D6FC}=k/(\unicode[STIX]{x1D70C}C_{p})$ is the thermal diffusivity, $k$ is the thermal conductivity, and $C_{p}$ is the heat capacity, of the fluid.

2.2 Boundary conditions at the interface

The system of coupled evolution equations for the velocity, pressure, temperature and concentration fields should be solved in a self-consistent manner, subject to the boundary conditions describing the balance of momentum, heat and mass fluxes between the two phases. The phase change at the interface can be described using kinetic theory (Schrage Reference Schrage1953). As we have shown previously (Qin et al. Reference Qin, Tuković and Grigoriev2015), the choice of the phase change model has negligible effect on the results. The mass flux across the interface is given by (Klentzman & Ajaev Reference Klentzman and Ajaev2009)

(2.15) $$\begin{eqnarray}J=\frac{2\unicode[STIX]{x1D712}}{2-\unicode[STIX]{x1D712}}\unicode[STIX]{x1D70C}_{v}\sqrt{\frac{\bar{R}_{v}T_{i}}{2\unicode[STIX]{x03C0}}}\left[\frac{p_{l}-p_{g}}{\unicode[STIX]{x1D70C}_{l}\bar{R}_{v}T_{i}}+\frac{{\mathcal{L}}}{\bar{R}_{v}T_{i}}\frac{T_{i}-T_{s}}{T_{s}}\right],\end{eqnarray}$$

where $\unicode[STIX]{x1D712}$ is the accommodation coefficient, ${\mathcal{L}}$ is the latent heat of vaporization, and subscript $s$ denotes saturation values for the vapour. The dependence of the local saturation temperature $T_{s}$ on the partial pressure of vapour $p_{v}$ is described using the Clausius–Clapeyron equation for phase equilibrium

(2.16) $$\begin{eqnarray}\ln \frac{p_{v}}{p_{v}^{0}}=\frac{{\mathcal{L}}}{\bar{R}_{v}}\left[\frac{1}{T_{0}}-\frac{1}{T_{s}}\right],\end{eqnarray}$$

where $p_{v}^{0}$ is the saturation pressure of the vapour at the equilibrium temperature $T_{0}$ . The first term in (2.15) is proportional to the Young–Laplace pressure and can be ignored in this problem, since the interface is considered flat.

The mass flux balance at the interface can be expressed with the help of (2.7). In the reference frame of the interface, the mass flux of the vapour is given by

(2.17) $$\begin{eqnarray}\frac{J}{m_{v}}=\hat{\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{j}_{v}=n_{g}\hat{\boldsymbol{n}}\boldsymbol{\cdot }([\boldsymbol{u}_{g}-\boldsymbol{u}_{i}]c_{v}-D\unicode[STIX]{x1D735}c_{v}),\end{eqnarray}$$

where $\boldsymbol{u}_{i}$ is the velocity of the liquid at the interface and $\hat{\boldsymbol{n}}$ is the unit vector normal to the interface. Since air is non-condensable, its mass flux across the interface is zero:

(2.18) $$\begin{eqnarray}0=\hat{\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{j}_{a}=n_{g}\hat{\boldsymbol{n}}\boldsymbol{\cdot }([\boldsymbol{u}_{g}-\boldsymbol{u}_{i}]c_{a}-D\unicode[STIX]{x1D735}c_{a}).\end{eqnarray}$$

Since $c_{a}+c_{v}=1$ , these two relations can be solved yielding two of the boundary conditions for (2.1) and (2.7) in the gas phase

(2.19) $$\begin{eqnarray}\hat{\boldsymbol{n}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c_{v}=-\frac{c_{a}J}{m_{v}n_{g}D}\end{eqnarray}$$

and

(2.20) $$\begin{eqnarray}\hat{\boldsymbol{n}}\boldsymbol{\cdot }(\boldsymbol{u}_{g}-\boldsymbol{u}_{i})=\frac{J}{m_{v}n_{g}}.\end{eqnarray}$$

The heat flux balance is given by

(2.21) $$\begin{eqnarray}{\mathcal{L}}J=\hat{\boldsymbol{n}}\boldsymbol{\cdot }k_{g}\unicode[STIX]{x1D735}T_{g}-\hat{\boldsymbol{n}}\boldsymbol{\cdot }k_{l}\unicode[STIX]{x1D735}T_{l},\end{eqnarray}$$

where the advective contribution to the heat flux is negligible on both sides of the interface. Indeed, in the gas phase, conduction is the dominant contribution (Qin & Grigoriev Reference Qin and Grigoriev2015), while on the liquid side

(2.22) $$\begin{eqnarray}\hat{\boldsymbol{n}}\boldsymbol{\cdot }(\boldsymbol{u}_{l}-\boldsymbol{u}_{i})=\frac{J}{\unicode[STIX]{x1D70C}_{l}}.\end{eqnarray}$$

Since the liquid density is much greater than that of the gas, the left-hand side of (2.22) is very small compared with $\hat{\boldsymbol{n}}\boldsymbol{\cdot }(\boldsymbol{u}_{g}-\boldsymbol{u}_{i})$ and can be ignored.

The remaining boundary conditions for $\boldsymbol{u}$ and $T$ at the liquid–vapour interface are standard: the temperature is continuous

(2.23) $$\begin{eqnarray}T_{l}=T_{i}=T_{v}\end{eqnarray}$$

as are the tangential velocity components

(2.24) $$\begin{eqnarray}(\mathbb{I}-\hat{\boldsymbol{n}}\hat{\boldsymbol{n}})\boldsymbol{\cdot }(\boldsymbol{u}_{l}-\boldsymbol{u}_{g})=0.\end{eqnarray}$$

The stress balance

(2.25) $$\begin{eqnarray}(\unicode[STIX]{x1D72E}_{l}-\unicode[STIX]{x1D72E}_{g})\boldsymbol{\cdot }\hat{\boldsymbol{n}}=\hat{\boldsymbol{n}}\unicode[STIX]{x1D705}\unicode[STIX]{x1D70E}-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D735}_{s}T_{i}\end{eqnarray}$$

incorporates both the viscous drag between the two phases and the thermocapillary effect. Here $\unicode[STIX]{x1D72E}=\unicode[STIX]{x1D707}[\unicode[STIX]{x1D735}\boldsymbol{u}+(\unicode[STIX]{x1D735}\boldsymbol{u})^{\text{T}}]-p$ is the stress tensor, $\unicode[STIX]{x1D705}(=0)$ is the interfacial curvature, $\unicode[STIX]{x1D735}_{s}=(\mathbb{I}-\hat{\boldsymbol{n}}\hat{\boldsymbol{n}})\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ is the surface gradient, and $\unicode[STIX]{x1D6FE}=-\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}/\unicode[STIX]{x2202}T$ is the temperature coefficient of surface tension.

We will further assume that side/top/bottom walls are adiabatic, $\hat{\boldsymbol{n}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T=0$ , which is a good approximation for most experimental set-ups. The heat transport through the end walls at $x=0$ and $x=L_{x}$ is not treated explicitly in this study, but in principle a variety of boundary conditions, from constant temperature to constant flux to mixed boundary conditions can be accommodated. Standard no-slip boundary conditions $\boldsymbol{u}=0$ for the velocity and no-flux boundary conditions $\hat{\boldsymbol{n}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c_{a,v}=0$ for concentration apply on all the walls.

3.1 Boundary conditions

At the bottom of the liquid layer ( $\tilde{z}=-1$ ) and the top of the gas layer ( $\tilde{z}=A$ ), no-slip and adiabatic boundary conditions apply

(3.12) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\tilde{\boldsymbol{u}}=0,\\ \unicode[STIX]{x2202}_{\tilde{z}}\tilde{T}=0.\end{array}\right\}\end{eqnarray}$$

At the interface ( $\tilde{z}=0$ ), the temperature and velocity fields are continuous

(3.13) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\tilde{T}_{l}=\tilde{T}_{g}=\tilde{T}_{i},\\ \tilde{\boldsymbol{u}}_{l}=\tilde{\boldsymbol{u}}_{g}=\tilde{\boldsymbol{u}}_{i},\end{array}\right\}\end{eqnarray}$$

where for our choice of the origin of coordinate system $\tilde{T}_{i}=\tilde{x}$ . Since $\unicode[STIX]{x1D707}_{l}\gg \unicode[STIX]{x1D707}_{g}$ and typically $d_{g}>d_{l}$ , the viscous stress in the gas layer can be ignored, yielding a simplified expression for the shear stress balance at the interface

(3.14) $$\begin{eqnarray}\unicode[STIX]{x2202}_{\tilde{z}}\tilde{u} _{l,x}=-Re\unicode[STIX]{x2202}_{\tilde{x}}\tilde{T}_{l},\end{eqnarray}$$

where

(3.15) $$\begin{eqnarray}Re=\frac{Ma}{Pr}=\frac{\unicode[STIX]{x1D6FE}d_{l}^{2}}{\unicode[STIX]{x1D707}_{l}\unicode[STIX]{x1D708}_{l}}\unicode[STIX]{x1D70F}\end{eqnarray}$$

is the Reynolds number.

The heat flux balance (2.21) at the interface reduces to

(3.16) $$\begin{eqnarray}\unicode[STIX]{x2202}_{\tilde{z}}\tilde{T}_{l}=\frac{k_{g}}{k_{l}}\unicode[STIX]{x2202}_{\tilde{z}}\tilde{T}_{g}-\frac{V}{Ma}\tilde{J},\end{eqnarray}$$

where $k_{g}=c_{v}^{0}k_{v}+c_{a}^{0}k_{a}$ , and $\tilde{J}=Jd_{l}/(Dm_{v}n_{g})$ is the dimensionless mass flux. The dimensionless combination

(3.17) $$\begin{eqnarray}V=\frac{{\mathcal{L}}\unicode[STIX]{x1D6FE}d_{l}}{\unicode[STIX]{x1D6FC}_{l}\unicode[STIX]{x1D707}_{l}k_{l}}\frac{Dp_{g}^{0}}{\bar{R}_{v}T_{0}},\end{eqnarray}$$

or more precisely the ratio $V/Ma$ , describes the relative magnitude of the latent heat released (absorbed) at the interface due to condensation (evaporation) compared with the vertical heat flux in the liquid layer due to conduction. It should be noted that, although the product $Dp_{g}^{0}$ , and consequently $V$ , is a function of $T_{0}$ , it is independent of the gas pressure, and therefore $c_{a}^{0}$ .

Non-dimensionalizing (2.19) and (2.20), we obtain

(3.18) $$\begin{eqnarray}\unicode[STIX]{x2202}_{\tilde{z}}\tilde{c}_{v}=-c_{a}^{0}c_{v}^{0}\tilde{J},\end{eqnarray}$$

and

(3.19) $$\begin{eqnarray}\hat{\boldsymbol{n}}\boldsymbol{\cdot }(\tilde{\boldsymbol{u}}_{g}-\tilde{\boldsymbol{u}}_{i})=c_{v}^{0}K_{D}\tilde{J}.\end{eqnarray}$$

The mass flux vanishes at the top of the gas layer

(3.20) $$\begin{eqnarray}\unicode[STIX]{x2202}_{\tilde{z}}\tilde{c}_{v}=0.\end{eqnarray}$$

The base uniform flow corresponds to a vanishing mass flux at the interface, $\tilde{J}=0$ , which is also an assumption made by all one-sided models. The base flow is spatially uniform in the lateral direction(s) which, coupled with incompressibility, requires $\hat{\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{u}_{i}=0$ . This leads to a number of simplifications. In particular, (3.18) and (3.19) require that $\tilde{u} _{l,z}=\tilde{u} _{g,z}=0$ and $\unicode[STIX]{x2202}_{\tilde{z}}\tilde{c}_{v}=0$ at the interface. Furthermore, the heat flux at the interface should also vanish, so that (2.21) gives $\unicode[STIX]{x2202}_{\tilde{z}}\tilde{T}_{l}=\unicode[STIX]{x2202}_{\tilde{z}}\tilde{T}_{g}=0$ .

The boundary conditions for the stream function can be easily obtained from those for the velocities. At the bottom of the liquid layer and the top of the gas layer

(3.21) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\tilde{\unicode[STIX]{x1D713}}=0,\\ \unicode[STIX]{x2202}_{\tilde{x}}\tilde{\unicode[STIX]{x1D713}}=\unicode[STIX]{x2202}_{\tilde{z}}\tilde{\unicode[STIX]{x1D713}}=0.\end{array}\right\}\end{eqnarray}$$

For a uniform flow, the net flux through any vertical plane vanishes. Combined with conditions (2.24) and (2.25) this requires that at the interface

(3.22) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\tilde{\unicode[STIX]{x1D713}}_{l}=\tilde{\unicode[STIX]{x1D713}}_{g}=0,\\ \unicode[STIX]{x2202}_{\tilde{z}}\tilde{\unicode[STIX]{x1D713}}_{l}-\unicode[STIX]{x2202}_{\tilde{z}}\tilde{\unicode[STIX]{x1D713}}_{g}=0,\\ \unicode[STIX]{x2202}_{\tilde{z}}^{2}\tilde{\unicode[STIX]{x1D713}}_{l}+Re\unicode[STIX]{x2202}_{\tilde{x}}\tilde{T}_{l}=0.\end{array}\right\}\end{eqnarray}$$

3.2 Fluid flow and temperature in the liquid layer

For a uniform horizontal flow, where both $\tilde{\unicode[STIX]{x1D713}}_{l}$ and $\tilde{\unicode[STIX]{x1D713}}_{g}$ are functions of $\tilde{z}$ alone, we can look for solutions to (3.10) in the form

(3.23) $$\begin{eqnarray}\tilde{T}_{l}=\tilde{x}+\tilde{\unicode[STIX]{x1D703}}_{l}(\tilde{z}),\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D703}}_{l}(0)=0$ . With this choice, the system (3.10) reduces to coupled ordinary differential equations (ODEs)

(3.24) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}-\tilde{\unicode[STIX]{x1D713}}_{l}^{\prime \prime \prime \prime }+Gr=0,\\ Pr\,\tilde{\unicode[STIX]{x1D713}}_{l}^{\prime }-\tilde{\unicode[STIX]{x1D703}}_{l}^{\prime \prime }=0,\end{array}\right\}\end{eqnarray}$$

where prime stands for the derivatives with respect to the $\tilde{z}$ coordinate.

Solving the system (3.24) subject to the boundary conditions at the bottom and the free surface of the liquid layer, we find the steady-state solutions for the stream function

(3.25) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D713}}_{l}=Re\left[-\frac{\tilde{z}(\tilde{z}+1)^{2}}{4}+Bo_{D}\frac{\tilde{z}(\tilde{z}+1)^{2}(2\tilde{z}-1)}{48}\right],\end{eqnarray}$$

velocity

(3.26) $$\begin{eqnarray}\tilde{\boldsymbol{u}}_{l}=Re\left[-\frac{(\tilde{z}+1)(3\tilde{z}+1)}{4}+Bo_{D}\frac{(\tilde{z}+1)(8\tilde{z}^{2}+\tilde{z}-1)}{48}\right]\hat{\boldsymbol{x}},\end{eqnarray}$$

and temperature field

(3.27) $$\begin{eqnarray}\tilde{T}_{l}=\tilde{x}+Ma\left[-\frac{\tilde{z}^{2}(3\tilde{z}^{2}+8\tilde{z}+6)}{48}+Bo_{D}\frac{\tilde{z}^{2}(8\tilde{z}^{3}+15\tilde{z}^{2}-10)}{960}\right],\end{eqnarray}$$

describing the base flow. They agree with the analytical solutions originally obtained by Birikh (Reference Birikh1966) and later rederived by Kirdyashkin (Reference Kirdyashkin1984) and Villers & Platten (Reference Villers and Platten1987) using a one-sided model that ignores the effects of the gas phase.

The assumption that the interfacial temperature varies linearly has been widely used in previous studies, without much justification, for deriving the solutions (3.26) and (3.27) for the return flow underlying the stability analyses (Parmentier et al. Reference Parmentier, Regnier and Lebon1993; Mercier & Normand Reference Mercier and Normand1996; Priede & Gerbeth Reference Priede and Gerbeth1997) as well as in models of the adiabatic section of heat pipes (Ha & Peterson Reference Ha and Peterson1994; Suman & Kumar Reference Suman and Kumar2005; Markos, Ajaev & Homsy Reference Markos, Ajaev and Homsy2006), which assume unidirectional flow in the liquid phase. However, the validity of this assumption cannot be established by a one-sided model which ignores heat and mass transport in the gas phase. In fact, when $c_{a}^{0}$ becomes sufficiently low, the interfacial temperature profile becomes nonlinear (Li et al. Reference Li, Grigoriev and Yoda2014; Qin & Grigoriev Reference Qin and Grigoriev2015). Proper justification of the linearity assumption requires showing that it is consistent with a steady-state solution of the transport equations in the gas phase, which satisfies all of the boundary conditions at the free surface. We turn to this next.

3.3 Fluid flow, temperature and composition in the gas layer

The solutions for the velocity, temperature and composition of the gas phase can be found in the same way the solutions (3.26) and (3.27) were obtained for the liquid phase. We will start by finding the solution for the vapour concentration at the interface before deriving the solution in the bulk. Since there is no phase change, $T_{s}=T_{i}$ , according to (2.15), and $\unicode[STIX]{x2202}_{x}T_{s}=\unicode[STIX]{x1D70F}$ . Using the Clausius–Clapeyron relation (2.16) and neglecting the deviation of $p_{g}$ from the equilibrium value (3.8) we therefore find

(3.28) $$\begin{eqnarray}\unicode[STIX]{x2202}_{x}c_{v}=\frac{{\mathcal{L}}\unicode[STIX]{x1D70F}}{\bar{R}_{v}T_{i}^{2}}c_{v},\end{eqnarray}$$

where for $\unicode[STIX]{x1D70F}x\ll T_{0}$ we can replace $T_{i}$ with the equilibrium temperature $T_{0}$ . In general, the solution to this equation yields an exponential concentration profile for $c_{v}$ and hence a nonlinear interfacial temperature profile (Qin & Grigoriev Reference Qin and Grigoriev2015). However, when the deviation of $c_{v}$ from $c_{v}^{0}$ is small, the concentration profile at the interface becomes approximately linear

(3.29) $$\begin{eqnarray}\tilde{c}_{v}\approx \unicode[STIX]{x1D6FA}\tilde{x},\end{eqnarray}$$

where

(3.30) $$\begin{eqnarray}\unicode[STIX]{x1D6FA}=c_{v}^{0}\frac{{\mathcal{L}}d_{l}\unicode[STIX]{x1D70F}}{\bar{R}_{v}T_{0}^{2}}=c_{v}^{0}\frac{H}{V}Ma,\end{eqnarray}$$

and

(3.31) $$\begin{eqnarray}H=\frac{{\mathcal{L}}^{2}Dp_{g}^{0}}{\bar{R}_{v}^{2}T_{0}^{3}k_{l}}\end{eqnarray}$$

is another non-dimensional parameter, the meaning of which will become clear later. Again, since the product $Dp_{g}^{0}$ is independent of $c_{a}^{0}$ , so is $H$ .

Given the boundary conditions for $\tilde{T}_{g}$ and $\tilde{c}_{v}$ , we can look for solutions to these two fields in the gas layer of the form

(3.32) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\tilde{T}_{g}=\tilde{x}+\tilde{\unicode[STIX]{x1D703}}_{g}(\tilde{z}),\\ \tilde{c}_{v}=\unicode[STIX]{x1D6FA}[\tilde{x}+\tilde{\unicode[STIX]{x1D70D}}_{v}(\tilde{z})],\end{array}\right\}\end{eqnarray}$$

where

(3.33) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\tilde{\unicode[STIX]{x1D703}}_{g}(0)=0,\\ \tilde{\unicode[STIX]{x1D70D}}_{v}(0)=0.\end{array}\right\}\end{eqnarray}$$

For a uniform flow, $\tilde{\unicode[STIX]{x1D713}}_{g}$ should only depends on $\tilde{z}$ , so the system (3.11) reduces to

(3.34) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}-\tilde{\unicode[STIX]{x1D713}}_{g}^{\prime \prime \prime \prime }+\unicode[STIX]{x1D6EF}_{T}+\unicode[STIX]{x1D6EF}_{c}\unicode[STIX]{x1D6FA}=0,\\ \tilde{\unicode[STIX]{x1D713}}_{g}^{\prime }-K_{\unicode[STIX]{x1D6FC}}\tilde{\unicode[STIX]{x1D703}}_{g}^{\prime \prime }=0,\\ \tilde{\unicode[STIX]{x1D713}}_{g}^{\prime }-K_{D}\tilde{\unicode[STIX]{x1D70D}}_{v}^{\prime \prime }=0.\end{array}\right\}\end{eqnarray}$$

Solving these equations subject to the boundary conditions (3.20), (3.21), (3.22) and (3.33) at the interface and the top of the gas layer yields the steady-state solutions for the stream function

(3.35) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D713}}_{g}=-{\mathcal{R}}\left[\frac{\tilde{z}(\tilde{z}-A)^{2}}{4A^{2}}+{\mathcal{B}}\frac{\tilde{z}(\tilde{z}-A)^{2}(2\tilde{z}+A)}{48A^{3}}\right],\end{eqnarray}$$

velocity

(3.36) $$\begin{eqnarray}\tilde{\boldsymbol{u}}_{g}=-{\mathcal{R}}\left[\frac{(\tilde{z}-A)(3\tilde{z}-A)}{4A^{2}}+{\mathcal{B}}\frac{(\tilde{z}-A)(8\tilde{z}^{2}-A\tilde{z}-A^{2})}{48A^{3}}\right]\hat{\boldsymbol{x}},\end{eqnarray}$$

temperature

(3.37) $$\begin{eqnarray}\tilde{T}_{g}=\tilde{x}-\frac{{\mathcal{R}}}{K_{\unicode[STIX]{x1D6FC}}}\left[\frac{\tilde{z}^{2}(3\tilde{z}^{2}-8A\tilde{z}+6A^{2})}{48A^{2}}+{\mathcal{B}}\frac{\tilde{z}^{2}(8\tilde{z}^{3}-15A\tilde{z}^{2}+10A^{3})}{960A^{3}}\right],\end{eqnarray}$$

and vapour concentration in the gas phase

(3.38) $$\begin{eqnarray}\tilde{c}_{v}=\unicode[STIX]{x1D6FA}\tilde{x}-\unicode[STIX]{x1D6FA}\frac{{\mathcal{R}}}{K_{D}}\left[\frac{\tilde{z}^{2}(3\tilde{z}^{2}-8A\tilde{z}+6A^{2})}{48A^{2}}+{\mathcal{B}}\frac{\tilde{z}^{2}(8\tilde{z}^{3}-15A\tilde{z}^{2}+10A^{3})}{960A^{3}}\right].\end{eqnarray}$$

The parameters ${\mathcal{R}}$ and ${\mathcal{B}}$ are analogous to the Reynolds number and the dynamic Bond number, but incorporate the properties of both fluid layers:

(3.39) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\mathcal{R}}=Re\left(1+\frac{Bo_{D}}{12}\right)+\frac{A^{3}}{12K_{\unicode[STIX]{x1D708}}}(\unicode[STIX]{x1D6EF}_{T}+\unicode[STIX]{x1D6EF}_{\unicode[STIX]{x1D70D}}),\\ \displaystyle {\mathcal{B}}=-\frac{A^{3}}{{\mathcal{R}}K_{\unicode[STIX]{x1D708}}}(\unicode[STIX]{x1D6EF}_{T}+\unicode[STIX]{x1D6EF}_{\unicode[STIX]{x1D70D}}),\end{array}\right\}\end{eqnarray}$$

where we defined

(3.40) $$\begin{eqnarray}\unicode[STIX]{x1D6EF}_{\unicode[STIX]{x1D70D}}=\unicode[STIX]{x1D6EF}_{c}\unicode[STIX]{x1D6FA}=-\frac{(1-c_{a}^{0})(\bar{R}_{a}-\bar{R}_{v})}{\bar{R}_{a}-(\bar{R}_{a}-\bar{R}_{v})c_{a}^{0}}\frac{{\mathcal{L}}Bo_{D}}{\unicode[STIX]{x1D6FD}_{l}\bar{R}_{v}T^{2}}\frac{Ma}{Pr}.\end{eqnarray}$$

Note that the form of the analytical solutions (3.26)–(3.27) and (3.36)–(3.38) is different from that of the solutions derived by Qin et al. (Reference Qin, Tuković and Grigoriev2014) because (i) the buoyancy force caused by the variation in the composition of the gas is explicitly taken into account in the present analysis, (ii) the mass transport in the gas phase is described using the concentration rather than the mass density of the vapour, and (iii) a different choice of the origin of the coordinate system is made. Figure 2 compares the analytical solutions with the numerical ones (cf. figure 1) computed by Qin & Grigoriev (Reference Qin and Grigoriev2015). We find that the analytical solutions accurately describe the core region of the flow away from the end walls for sufficiently low $\unicode[STIX]{x1D70F}$ at which the base flow is stable.

Figure 2.

Comparison between the numerical and analytical solutions for $Ma=190$ and $c_{a}^{0}=0.96$ (atmospheric conditions) in the middle of the cell, $x=L_{x}/2$ . Shown are the normalized vertical profiles of (a) the horizontal velocity $u_{x}=u_{x}$ , (b) the variation $\unicode[STIX]{x1D6FF}T=T-T_{i}$ of the temperature, and (c) the variation $\unicode[STIX]{x1D6FF}c_{v}=c_{v}-c_{v,i}$ of vapour concentration. Open circles correspond to numerical results and solid lines correspond to the analytical solutions.

4.1 Diffusion-dominated case

While the non-dimensional equations (3.4) for the liquid layer can be used without modification, in order to obtain a simplified problem, which can generate useful physical insight, we will start by making several approximations in the treatment of heat and mass transport on the gas side. The relative contribution of advection and diffusion to the mass and heat transport in the gas layer are described by the Péclet numbers $Pe_{m}$ and $Pe_{t}$ , respectively. Since the relevant length scale here is the wavelength of the convective pattern, which is comparable to the depth of the liquid layer, and the velocity scale is determined by the interfacial velocity $u_{i}=\unicode[STIX]{x1D708}_{l}\tilde{u} _{l}(0)/d_{l}$ (cf. (3.26)), we have

(4.1) $$\begin{eqnarray}Pe_{m}=\frac{|u_{i}|d_{l}}{D}=\frac{12+Bo_{D}}{48}\frac{Re}{K_{D}},\end{eqnarray}$$

and

(4.2) $$\begin{eqnarray}Pe_{t}=\frac{|u_{i}|d_{l}}{\unicode[STIX]{x1D6FC}_{g}}=\frac{12+Bo_{D}}{48}\frac{Re}{K_{\unicode[STIX]{x1D6FC}}}.\end{eqnarray}$$

In typical experiments (Villers & Platten Reference Villers and Platten1992; Riley & Neitzel Reference Riley and Neitzel1998; Li et al. Reference Li, Grigoriev and Yoda2014) both Péclet numbers are at most $O(1)$ , so we can drop the advective terms in the transport equations (2.13) and (2.14) in the gas phase. Furthermore, if we are only interested in the transition thresholds for stationary instability, we can also ignore the time derivatives, such that the system (3.11) simplifies to

(4.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\tilde{\unicode[STIX]{x1D713}}_{g}=0,\\ (\unicode[STIX]{x2202}_{\tilde{x}}^{2}+\unicode[STIX]{x2202}_{\tilde{z}}^{2})\tilde{T}_{g}=0,\\ (\unicode[STIX]{x2202}_{\tilde{x}}^{2}+\unicode[STIX]{x2202}_{\tilde{z}}^{2})\tilde{c}_{v}=0.\end{array}\right\}\end{eqnarray}$$

Note that neglecting the transient dynamics of the gas phase changes neither the critical Marangoni number nor the critical wavelength of a (stationary) pattern.

The base solution $\tilde{\unicode[STIX]{x1D713}}_{l0}$ , $\tilde{T}_{l0}$ describing the uniform return flow in the liquid layer is given by (3.25)–(3.27). The base solution to (4.3) describing the gas layer in the diffusion-dominated case can be obtained by setting ${\mathcal{R}}=0$ in (3.37)–(3.38), which is equivalent to setting $\tilde{\boldsymbol{u}}_{g}=0$ ( $\tilde{\unicode[STIX]{x1D713}}_{g}=0$ ) and gives:

(4.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\tilde{T}_{g0}=\tilde{x},\\ \tilde{c}_{v0}=\unicode[STIX]{x1D6FA}\tilde{x}.\end{array}\right\}\end{eqnarray}$$

Although the base state satisfies the adiabatic boundary condition at the interface, perturbations in the temperature will give rise to heat and mass flux across the interface and, consequently, through the gas phase. Perturbed solutions can be written in the form of Fourier integrals

(4.5) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{\unicode[STIX]{x1D713}}_{m}=\tilde{\unicode[STIX]{x1D713}}_{m0}+\int _{-\infty }^{\infty }\tilde{\unicode[STIX]{x1D713}}_{mq}(\tilde{z})\text{e}^{\text{i}qx+\unicode[STIX]{x1D70E}_{q}t}\,\text{d}q,\\ \displaystyle \tilde{T}_{m}=\tilde{T}_{m0}+\int _{-\infty }^{\infty }\tilde{\unicode[STIX]{x1D703}}_{mq}(\tilde{z})\text{e}^{\text{i}qx+\unicode[STIX]{x1D70E}_{q}t}\,\text{d}q,\\ \displaystyle \tilde{c}_{v}=\tilde{c}_{v0}+\unicode[STIX]{x1D6FA}\int _{-\infty }^{\infty }\tilde{\unicode[STIX]{x1D70D}}_{vq}(\tilde{z})\text{e}^{\text{i}qx+\unicode[STIX]{x1D70E}_{q}t}\,\text{d}q,\end{array}\right\}\end{eqnarray}$$

where $m=\{l,g\}$ denotes the phase, $q$ is the wavenumber describing the variation in the horizontal direction, and $\unicode[STIX]{x1D70E}_{q}$ is the temporal growth rate. The functions $\tilde{\unicode[STIX]{x1D713}}_{mq}(\tilde{z})$ , $\tilde{\unicode[STIX]{x1D703}}_{mq}(\tilde{z})$ , and $\tilde{\unicode[STIX]{x1D70D}}_{vq}(\tilde{z})$ define the vertical profile for the perturbations in, respectively, the stream function $\tilde{\unicode[STIX]{x1D713}}$ , temperature $\tilde{T}$ , and vapour concentration $\tilde{c}_{v}$ with wavenumber $q$ .

Rewriting (4.3) in terms of the perturbations, we obtain

(4.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{\unicode[STIX]{x1D713}}_{gq}=0,\\ \displaystyle \tilde{\unicode[STIX]{x1D703}}_{gq}^{\prime \prime }=q^{2}\tilde{\unicode[STIX]{x1D703}}_{gq},\\ \displaystyle \tilde{\unicode[STIX]{x1D70D}}_{vq}^{\prime \prime }=q^{2}\tilde{\unicode[STIX]{x1D70D}}_{vq}.\end{array}\right\}\end{eqnarray}$$

Temperature continuity at the interface requires

(4.7) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D703}}_{gq}(0)=\tilde{\unicode[STIX]{x1D703}}_{lq}(0).\end{eqnarray}$$

At the top of the gas layer we have

(4.8) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{\unicode[STIX]{x1D703}}_{gq}^{\prime }(A)=0,\\ \displaystyle \tilde{\unicode[STIX]{x1D70D}}_{vq}^{\prime }(A)=0.\end{array}\right\}\end{eqnarray}$$

The solution of (4.6) satisfying these boundary conditions is

(4.9) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D703}}_{gq}(\tilde{z})=\frac{\cosh [q(\tilde{z}-A)]}{\cosh (qA)}\tilde{\unicode[STIX]{x1D703}}_{gq}(0).\end{eqnarray}$$

The solution for the perturbation in the vapour concentration is analogous,

(4.10) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D70D}}_{vq}(\tilde{z})=\frac{\cosh [q(\tilde{z}-A)]}{\cosh (qA)}\tilde{\unicode[STIX]{x1D70D}}_{vq}(0),\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D70D}}_{vq}(0)$ can be related to the perturbation in the interfacial temperature via the Clausius–Clapeyron relation (2.16)

(4.11) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D70D}}_{vq}(0)=\tilde{\unicode[STIX]{x1D703}}_{gq}(0).\end{eqnarray}$$

Fourier transforming (3.18) yields

(4.12) $$\begin{eqnarray}\tilde{J}_{q}=-\frac{\unicode[STIX]{x1D6FA}}{c_{a}^{0}}\tilde{\unicode[STIX]{x1D70D}}_{vq}^{\prime }(0),\end{eqnarray}$$

where $\tilde{J}_{q}$ is the Fourier coefficient of $\tilde{J}(\tilde{x})$ . Subtracting the base solutions from the heat balance (3.16), Fourier transforming the result, and using (4.12) gives the following relation

(4.13) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D703}}_{lq}^{\prime }(0)=\frac{k_{g}}{k_{l}}\tilde{\unicode[STIX]{x1D703}}_{gq}^{\prime }(0)+\frac{\unicode[STIX]{x1D6FA}}{c_{a}^{0}}\frac{V}{Ma}\tilde{\unicode[STIX]{x1D70D}}_{vq}^{\prime }(0).\end{eqnarray}$$

With the help of (4.9), (4.10) and (4.11) this can rewritten in the form of Newton’s cooling law

(4.14) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D703}}_{lq}^{\prime }(0)=-Bi_{q}\tilde{\unicode[STIX]{x1D703}}_{lq}(0),\end{eqnarray}$$

where we introduced a wavenumber-dependent analogue of the Biot number

(4.15) $$\begin{eqnarray}Bi_{q}=q\tanh (qA)\left[\frac{k_{g}}{k_{l}}+\frac{1-c_{a}^{0}}{c_{a}^{0}}H\right],\end{eqnarray}$$

which will be referred to as the Biot coefficient below to highlight the fact that it is a function of $q$ and $c_{a}^{0}$ , not a constant. The prefactor $\tanh (qA)$ describes the effect of finite thickness of the gas layer. For $d_{g}\gg d_{l}$ ( $A\gg 1$ ), it approaches unity and (4.15) reduces to the expression derived by Chauvet, Dehaeck & Colinet (Reference Chauvet, Dehaeck and Colinet2012) in the context of stability of a volatile liquid layer in the presence of a vertical temperature gradient. The first term describes the effect of thermal conduction through the gas layer, while the second term describes the effect of latent heat released/absorbed at the interface as a result of condensation/evaporation.

By linearizing the governing equations (3.10) around the base state (3.25)–(3.27), we obtain the evolution equations for the perturbations $\tilde{\unicode[STIX]{x1D713}}_{lq}$ and $\tilde{\unicode[STIX]{x1D703}}_{lq}$ in the liquid layer

(4.16) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{\unicode[STIX]{x1D6FB}}_{q}^{2}\tilde{\unicode[STIX]{x1D6FB}}_{q}^{2}\tilde{\unicode[STIX]{x1D713}}_{lq}+\text{i}qC_{1}(\tilde{z})Re\tilde{\unicode[STIX]{x1D713}}_{lq}^{\prime \prime }-\text{i}qC_{2}(\tilde{z})Re\tilde{\unicode[STIX]{x1D713}}_{lq}-\text{i}qGr\tilde{\unicode[STIX]{x1D703}}_{lq}=\unicode[STIX]{x1D70E}_{q}\tilde{\unicode[STIX]{x1D6FB}}_{q}^{2}\tilde{\unicode[STIX]{x1D713}}_{lq},\\ Pr^{-1}\tilde{\unicode[STIX]{x1D6FB}}_{q}^{2}\tilde{\unicode[STIX]{x1D703}}_{lq}+\text{i}qC_{1}(\tilde{z})Re\tilde{\unicode[STIX]{x1D703}}_{lq}-\text{i}qC_{3}(\tilde{z})Ma\tilde{\unicode[STIX]{x1D713}}_{lq}-\tilde{\unicode[STIX]{x1D713}}_{lq}^{\prime }=\unicode[STIX]{x1D70E}_{q}\tilde{\unicode[STIX]{x1D703}}_{lq},\end{array}\right\}\end{eqnarray}$$

where we defined $\tilde{\unicode[STIX]{x1D6FB}}_{q}^{2}=\unicode[STIX]{x2202}_{\tilde{z}}^{2}-q^{2}$ and

(4.17) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle C_{1}(\tilde{z})=\frac{(\tilde{z}+1)(3\tilde{z}+1)}{4}-Bo_{D}\frac{(\tilde{z}+1)(8\tilde{z}^{2}+\tilde{z}-1)}{48},\\ \displaystyle C_{2}(\tilde{z})=q^{2}C_{1}(\tilde{z})+\frac{3}{2}-Bo_{D}\frac{8\tilde{z}+3}{8},\\ \displaystyle C_{3}(\tilde{z})=\frac{\tilde{z}(\tilde{z}+1)^{2}}{4}-Bo_{D}\frac{\tilde{z}(\tilde{z}+1)^{2}(2\tilde{z}-1)}{48}.\end{array}\right\}\end{eqnarray}$$

This is a system of ODEs which is fourth order with respect to $\tilde{\unicode[STIX]{x1D713}}_{lq}$ and second order with respect to $\tilde{\unicode[STIX]{x1D703}}_{lq}$ and hence needs a total of six boundary conditions. These boundary conditions are:

(4.18) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{\unicode[STIX]{x1D713}}_{lq}(-1)=0,\\ \tilde{\unicode[STIX]{x1D713}}_{lq}^{\prime }(-1)=0,\\ \tilde{\unicode[STIX]{x1D703}}_{lq}^{\prime }(-1)=0,\end{array}\right\}\end{eqnarray}$$

at the bottom of the liquid layer and

(4.19) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{\unicode[STIX]{x1D713}}_{lq}(0)=0,\\ \tilde{\unicode[STIX]{x1D713}}_{lq}^{\prime \prime }(0)=-\text{i}qRe\tilde{\unicode[STIX]{x1D703}}_{lq}(0),\end{array}\right\}\end{eqnarray}$$

and

(4.20) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D703}}_{lq}^{\prime }(0)=-Bi_{q}\tilde{\unicode[STIX]{x1D703}}_{lq}(0),\end{eqnarray}$$

at the free surface. In the subsequent discussion we will refer to the system (4.16) with the boundary conditions (4.18)–(4.20) as the enhanced one-sided model, since it incorporates the effect of the gas phase entirely through boundary conditions at the interface.

4.2 Transient dynamics in the gas layer

A more accurate description of the instability, stationary or oscillatory, can be obtained by restoring the time dependence of the temperature and composition of the gas phase. This corresponds to replacing the Laplace equations (4.3) with

(4.21) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x2202}_{\tilde{t}}\tilde{T}_{g}=K_{\unicode[STIX]{x1D6FC}}\tilde{\unicode[STIX]{x1D6FB}}^{2}\tilde{T}_{g},\\ \displaystyle \unicode[STIX]{x2202}_{\tilde{t}}\tilde{c}_{v}=K_{D}\tilde{\unicode[STIX]{x1D6FB}}^{2}\tilde{c}_{v}.\end{array}\right\}\end{eqnarray}$$

The corresponding equations for the perturbations are

(4.22) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{\unicode[STIX]{x1D713}}_{gq}=0,\\ \displaystyle \tilde{\unicode[STIX]{x1D703}}_{gq}^{\prime \prime }-q^{2}\tilde{\unicode[STIX]{x1D703}}_{gq}=\unicode[STIX]{x1D70E}_{q}K_{\unicode[STIX]{x1D6FC}}^{-1}\tilde{\unicode[STIX]{x1D703}}_{gq},\\ \displaystyle \tilde{\unicode[STIX]{x1D70D}}_{vq}^{\prime \prime }-q^{2}\tilde{\unicode[STIX]{x1D70D}}_{vq}=\unicode[STIX]{x1D70E}_{q}K_{D}^{-1}\tilde{\unicode[STIX]{x1D70D}}_{vq}.\end{array}\right\}\end{eqnarray}$$

These equations should be solved subject to the boundary conditions (4.7), (4.8), (4.11), (4.18), (4.19) and

(4.23) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D703}}_{lq}^{\prime }(0)=-\frac{k_{g}}{k_{l}}\bar{\unicode[STIX]{x1D703}}_{gq}^{\prime }(0)-\frac{1-c_{a}^{0}}{c_{a}^{0}}H\bar{\unicode[STIX]{x1D70D}}_{vq}^{\prime }(0),\end{eqnarray}$$

which follows from the heat flux balance (3.16) and the mass flux balance (4.12) and replaces (4.14).

Since it involves an incomplete description of the gas layer (the effects of advection are not included), this model based on the governing equations (4.16) and (4.22) can be called ‘one-and-a-half-sided.’ Note that, since both $K_{\unicode[STIX]{x1D6FC}}$ and $K_{D}$ are typically large compared with unity, the terms on the right-hand side of (4.22) are small. Dropping them would lead to the enhanced one-sided model derived in § 4.1.

4.3 The effect of advection in the gas layer

The effect of advection can also be incorporated for the range of $c_{a}^{0}$ where the base solution (3.35)–(3.38) is valid. By setting $\tilde{\unicode[STIX]{x1D713}}_{gq}\neq 0$ we ensure that all of the transport mechanisms in the gas phase are accounted for and both layers are treated using an equally comprehensive model. The linearized evolution equations for the perturbations $\tilde{\unicode[STIX]{x1D713}}_{gq}$ , $\tilde{\unicode[STIX]{x1D70D}}_{vq}$ and $\tilde{\unicode[STIX]{x1D703}}_{gq}$ in the gas phase are

(4.24) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle K_{\unicode[STIX]{x1D708}}\tilde{\unicode[STIX]{x1D6FB}}_{q}^{2}\tilde{\unicode[STIX]{x1D6FB}}_{q}^{2}\tilde{\unicode[STIX]{x1D713}}_{gq}+\text{i}q\tilde{C}_{1}(\tilde{z}){\mathcal{R}}\tilde{\unicode[STIX]{x1D713}}_{gq}^{\prime \prime }-\text{i}q\tilde{C}_{2}(\tilde{z}){\mathcal{R}}\tilde{\unicode[STIX]{x1D713}}_{gq}-\text{i}q\unicode[STIX]{x1D6EF}_{T}\tilde{\unicode[STIX]{x1D703}}_{gq}-\text{i}q\unicode[STIX]{x1D6EF}_{\unicode[STIX]{x1D70D}}\tilde{\unicode[STIX]{x1D70D}}_{vq}=\unicode[STIX]{x1D70E}_{q}\tilde{\unicode[STIX]{x1D6FB}}_{q}^{2}\tilde{\unicode[STIX]{x1D713}}_{gq},\\ \displaystyle K_{\unicode[STIX]{x1D6FC}}\tilde{\unicode[STIX]{x1D6FB}}_{q}^{2}\tilde{\unicode[STIX]{x1D703}}_{gq}+\text{i}q\tilde{C}_{1}(\tilde{z}){\mathcal{R}}\tilde{\unicode[STIX]{x1D703}}_{gq}-\text{i}q\tilde{C}_{3}(\tilde{z}){\mathcal{R}}K_{\unicode[STIX]{x1D6FC}}^{-1}\tilde{\unicode[STIX]{x1D713}}_{gq}-\tilde{\unicode[STIX]{x1D713}}_{gq}^{\prime }=\unicode[STIX]{x1D70E}_{q}\tilde{\unicode[STIX]{x1D703}}_{gq},\\ \displaystyle K_{D}\tilde{\unicode[STIX]{x1D6FB}}_{q}^{2}\tilde{\unicode[STIX]{x1D70D}}_{vq}+\text{i}q\tilde{C}_{1}(\tilde{z}){\mathcal{R}}\tilde{\unicode[STIX]{x1D70D}}_{vq}-\text{i}q\tilde{C}_{3}(\tilde{z}){\mathcal{R}}K_{D}^{-1}\tilde{\unicode[STIX]{x1D713}}_{gq}-\tilde{\unicode[STIX]{x1D713}}_{gq}^{\prime }=\unicode[STIX]{x1D70E}_{q}\tilde{\unicode[STIX]{x1D70D}}_{vq},\end{array}\right\}\qquad\end{eqnarray}$$

where

(4.25) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{C}_{1}(\tilde{z})=\frac{(\tilde{z}-A)(3\tilde{z}-A)}{4A^{2}}+{\mathcal{B}}\frac{(\tilde{z}-A)(8\tilde{z}^{2}-A\tilde{z}-A^{2})}{48A^{3}},\\ \displaystyle \tilde{C}_{2}(\tilde{z})=q^{2}\tilde{C}_{1}(\tilde{z})+\frac{3}{2A^{2}}+{\mathcal{B}}\frac{8\tilde{z}-3A}{8A^{3}},\\ \displaystyle \tilde{C}_{3}(\tilde{z})=\frac{\tilde{z}(\tilde{z}-A)^{2}}{4A^{2}}+{\mathcal{B}}\frac{\tilde{z}(\tilde{z}-A)^{2}(2\tilde{z}+A)}{48A^{3}}.\end{array}\right\}\end{eqnarray}$$

The boundary conditions are given by (4.7), (4.8), (4.11), (4.18), (4.19) and (4.23). In addition, the boundary conditions for $\bar{\unicode[STIX]{x1D713}}_{gq}(\bar{z})$ follow from (3.21) and (3.22):

(4.26) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \tilde{\unicode[STIX]{x1D713}}_{gq}(A)=0,\\ \tilde{\unicode[STIX]{x1D713}}_{gq}^{\prime }(A)=0,\end{array}\right\}\end{eqnarray}$$

at the top of the gas layer and

(4.27) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\tilde{\unicode[STIX]{x1D713}}_{gq}(0)=0,\\ \tilde{\unicode[STIX]{x1D713}}_{gq}^{\prime }(0)=\tilde{\unicode[STIX]{x1D713}}_{lq}^{\prime }(0),\end{array}\right\}\end{eqnarray}$$

at the free surface. In the subsequent discussion we will refer to the system (4.16) and (4.24) with the appropriate boundary conditions as the two-sided model, since it treats both the liquid and the gas phase with the same level of detail.

4.4 Different modes of instability

The boundary value problem, whether it is described by the system (4.16) by itself or combined with (4.22) or (4.24), has a spectrum of eigenvalues $\unicode[STIX]{x1D70E}_{q}^{n}$ and the corresponding eigenfunctions $\{\unicode[STIX]{x1D713}_{lq}^{n},\unicode[STIX]{x1D703}_{lq}^{n},\ldots \}$ , $n=0,1,\ldots \,$ , for a given complex wavenumber $q=k+\text{i}s$ and Marangoni number $Ma$ . The stability of the base flow is determined by the eigenvalue $\unicode[STIX]{x1D70E}_{q}^{0}=\unicode[STIX]{x1D705}_{q}+\text{i}\unicode[STIX]{x1D714}_{q}$ with the largest real part $\unicode[STIX]{x1D705}_{q}$ and the character of the instability (oscillatory or stationary) is determined by the imaginary part $\unicode[STIX]{x1D714}_{q}$ .

So far we have only focused on the boundary conditions in the vertical (confined) direction. In both experiment and simulations the fluid layers are also bounded in the horizontal (extended) directions $0<x<L_{x}$ and $0<y<L_{y}$ . Comparison of various experimental results suggests that the spanwise aspect ratio $\unicode[STIX]{x1D6E4}_{y}=L_{y}/d_{l}$ is not particularly important: aside from the orientation of hydrothermal waves at low $Bo_{D}$ , the pattern that emerges above the threshold of primary instability is basically the same for both $\unicode[STIX]{x1D6E4}_{y}=O(1)$ (Villers & Platten Reference Villers and Platten1992; De Saedeleer et al. Reference De Saedeleer, Garcimartín, Chavepeyer, Platten and Lebon1996; Garcimartín et al. Reference Garcimartín, Mukolobwiez and Daviaud1997; Li et al. Reference Li, Grigoriev and Yoda2014) and $\unicode[STIX]{x1D6E4}_{y}\gg 1$ (Riley & Neitzel Reference Riley and Neitzel1998; Burguete et al. Reference Burguete, Mukolobwiez, Daviaud, Garnier and Chiffaudel2001). This is not very surprising since the boundaries in the spanwise direction play a passive role – they tend to suppress, rather than enhance convection.

The boundaries in the streamwise direction are less trivial and play a crucial role in selecting the convection pattern for any $\unicode[STIX]{x1D6E4}_{x}=L_{x}/d_{l}$ . The corresponding boundary conditions are of the enhancing type (Cross & Greenside Reference Cross and Greenside2009): buoyancy, which plays a significant role for $Bo_{D}=O(1)$ , drives an upward flow near the hot end wall, generating a single localized convection roll. Similarly, buoyancy drives a downward flow near the cold end wall, generating another, much weaker, convection roll. These rolls appear for $Ma$ well below the critical value $Ma_{c}$ at which convection pattern emerges away from the end walls. The rolls are steady over a wide range of $Ma$ and have a horizontal extent of $O(d_{l})$ . For $\unicode[STIX]{x1D6E4}_{x}\gg 1$ and sufficiently low $Ma$ , the flow in the core region between these two convection rolls, i.e. for $O(d_{l})<x<L_{x}-O(d_{l})$ , is described well by the base flow solution (3.26) and (3.36), as mentioned previously (cf. figure 1). Riley & Neitzel (Reference Riley and Neitzel1998) refer to this state as a steady unicellular flow (SUF). Since the perturbation strength is always finite at the end walls, defining the threshold of the primary instability at finite $\unicode[STIX]{x1D6E4}_{x}$ requires some care.

4.4.1 Convective instability

The dominant mode of primary instability in laterally unbounded systems (infinite $\unicode[STIX]{x1D6E4}_{x}$ ) tends to be convective. The threshold of convective instability corresponds to $\unicode[STIX]{x1D705}_{q}$ being negative over all $q$ -real, except for $q=k_{c}$ , where $k_{c}$ is the (real) critical wavenumber which defines the wavelength $\unicode[STIX]{x1D706}/d_{l}=2\unicode[STIX]{x03C0}/k_{c}$ of the pattern that develops above threshold. Therefore, the critical wavenumber $k_{c}$ and the critical Marangoni number $Ma_{c}$ are given by the solution of the two equations

(4.28) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D705}_{q}|_{q_{c}}=0,\\ \displaystyle \left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D705}_{q}}{\unicode[STIX]{x2202}q}\right|_{q_{c}}=0,\end{array}\right\}\end{eqnarray}$$

over $q$ -real. Generally, the critical frequency $\unicode[STIX]{x1D714}_{c}=\unicode[STIX]{x1D714}_{q}|_{q_{c}}\neq 0$ , so the resulting pattern will be oscillatory, with disturbances amplified in the frame of reference moving with the group velocity (Huerre Reference Huerre, Batchelor, Moffatt and Worster2000)

(4.29) $$\begin{eqnarray}v_{g}=\left.-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}_{q}}{\unicode[STIX]{x2202}q}\right|_{q_{c}}.\end{eqnarray}$$

In laterally bounded systems, such as the fluid volumes with finite aspect ratios $\unicode[STIX]{x1D6E4}_{x}$ discussed in this paper, convective instabilities should not be observed in the absence of a localized oscillatory source of the right temporal frequency ( $\unicode[STIX]{x1D714}\approx \unicode[STIX]{x1D714}_{c}$ ). An example of such a source could be an oscillating convection roll near the left (right) end wall for $v_{g}>0\,(v_{g}<0)$ driven by a completely different mechanism (e.g. buoyancy due to cooling or heating at the end wall). While such oscillating convection rolls have indeed been observed in experiment (De Saedeleer et al. Reference De Saedeleer, Garcimartín, Chavepeyer, Platten and Lebon1996; Garcimartín et al. Reference Garcimartín, Mukolobwiez and Daviaud1997; Li et al. Reference Li, Grigoriev and Yoda2014) and simulations (Qin et al. Reference Qin, Tuković and Grigoriev2014), they arise after convection pattern in the core region of the flow has already been established. Hence, this is a secondary, rather than a primary, instability. We should only expect to see a primary convective instability in systems with very large streamwise aspect ratio $\unicode[STIX]{x1D6E4}_{x}$ where infinitesimal disturbances have sufficient time $\unicode[STIX]{x1D6E4}_{x}/v_{g}$ to develop.

4.4.2 Absolute instability

Absolute instabilities arise when infinitesimal disturbances grow in the laboratory frame. Since they do not need a finite-amplitude source, they can arise both in laterally unbounded systems as well as in bounded systems with moderate $\unicode[STIX]{x1D6E4}_{x}$ . Absolute instability requires that the group velocity $v_{g}=0$ at onset, yielding a system of three real equations for three real unknown $k_{c}=\text{Re}(q_{c})$ , $s_{c}=\text{Im}(q_{c})$ , and $Ma_{c}$ , which can be rewritten as a system of two equations (Huerre Reference Huerre, Batchelor, Moffatt and Worster2000)

(4.30) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D705}_{q}|_{q_{c}}=0,\\ \displaystyle \left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{q}^{0}}{\unicode[STIX]{x2202}q}\right|_{q_{c}}=0,\end{array}\right\}\end{eqnarray}$$

the second of which is complex, for $q_{c}$ -complex and $Ma_{c}$ .

The resulting pattern will generally be oscillatory, $\unicode[STIX]{x1D714}_{c}\neq 0$ , just like in the case of convective instability. Although its group velocity vanishes, the phase velocity $v_{p}=-\unicode[STIX]{x1D714}_{c}/k_{c}$ of the pattern will be non-zero. A good example of such an absolute instability is hydrothermal waves (Smith & Davis Reference Smith and Davis1983a ,Reference Smith and Davis b ) that are observed in the limit of low $Bo_{D}$ (Daviaud & Vince Reference Daviaud and Vince1993; Riley & Neitzel Reference Riley and Neitzel1998; Burguete et al. Reference Burguete, Mukolobwiez, Daviaud, Garnier and Chiffaudel2001). In the following, where necessary to avoid confusion, we will use the superscript $a$ to denote the critical values corresponding to absolute instability.

4.4.3 Stationary instability

The threshold of a stationary instability corresponds to both the real and imaginary part of $\unicode[STIX]{x1D70E}_{q}^{0}$ vanishing

(4.31) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D705}_{q}|_{q_{c}}=0,\\ \displaystyle \unicode[STIX]{x1D714}_{q}|_{q_{c}}=0.\end{array}\right\}\end{eqnarray}$$

These two equations should be solved for the critical wavenumber $k_{c}$ and the critical Marangoni number $Ma_{c}$ , which both depend on the spatial attenuation rate $s$ . In the following, where necessary to avoid confusion, we will use the superscript $s$ to denote the critical values corresponding to stationary instability. In particular, we will write $Ma_{c}^{s}(s)$ to explicitly denote the dependence on $s$ .

In laterally infinite systems the transition threshold is well defined and corresponds to a spatially uniform stationary pattern ( $s=0$ ). Such a pattern will be referred to as the steady multicellular (SMC) state below. For $Ma<Ma_{c}^{s}(0)$ the convection pattern arises due to forcing at one of the boundaries (convection rolls driven by buoyancy, cf. figure 1(c,d)), rather than spontaneously, and generally does not extend across the entire length of the system. In this case the absolute value of $s=\text{Im}(q_{c})$ defines the spatial attenuation of the perturbation (or the number of convection rolls that can be detected in the liquid layer), while $k_{c}=\text{Re}(q_{c})$ determines the wavelength $\unicode[STIX]{x1D706}$ of the pattern (or the distance between convection rolls). For $s\neq 0$ the pattern appears below $Ma_{c}^{s}(0)$ and is localized near the right (hot) end wall for $s<0$ or near the left (cold) end wall for $s>0$ . As we will see below, $\text{d}Ma_{c}^{s}/\text{d}s>0$ , which means that the pattern always appears near the hot end wall (cf. figure 1 d) and we can only consider $s\leqslant 0$ .

When $\text{e}^{-|s|\unicode[STIX]{x1D6E4}_{x}}\ll 1$ , the pattern extends over a region of length $O(2d_{l}/|s|)<L_{x}$ . We will refer to this as partial multicellular (PMC) state following Qin & Grigoriev (Reference Qin and Grigoriev2015) and Li et al. (Reference Li, Grigoriev and Yoda2014). In this analysis, we set $s_{PMC}=-1$ to define the transition from SUF to PMC at $Ma_{c}^{s}(-1)<Ma_{c}^{s}(0)$ , which roughly corresponds to one extra roll nucleating near the (hot) end wall. As $Ma$ is increased beyond this value, the convection pattern extends further away from the end wall, until it eventually covers the entire fluid layer. In practice this roughly corresponds to $\text{e}^{-|s|\unicode[STIX]{x1D6E4}_{x}}\approx 0.1$ or $s_{SMC}=\pm 2/\unicode[STIX]{x1D6E4}_{x}$ . This is the definition that has to be used for $Pr\lesssim 1$ (i.e. for gases or metallic liquids featuring sensitive dependence of $Ma_{c}^{s}$ on $s$ for small $|s|$ , as shown below). For typical (non-metallic) liquids with $Pr=O(10)$ , we can set $s_{SMC}=0$ , as in the case of laterally unbounded layers, as long as $\unicode[STIX]{x1D6E4}_{x}\gtrsim 20$ . This is the case in both the experimental (Riley & Neitzel Reference Riley and Neitzel1998; Li et al. Reference Li, Grigoriev and Yoda2014) and numerical studies (Qin et al. Reference Qin, Tuković and Grigoriev2014; Qin & Grigoriev Reference Qin and Grigoriev2015) used for comparison below.

Note that, for $\unicode[STIX]{x1D70E}_{q}^{0}=0$ , equations (4.22) reduce to (4.6) even if $K_{\unicode[STIX]{x1D6FC}}$ and $K_{D}$ are not large. Hence the two formulations always give identical predictions for the threshold of stationary instability (even though they predict different growth rates for modes with wavenumbers $q$ different from the critical one). Therefore, in what follows we will only focus on the limits when either both $Pe_{m}$ and $Pe_{t}$ are vanishingly small and stability is described by (4.16) alone, or when they are both $O(1)$ , so (4.16) and (4.24) have to be solved simultaneously.