Convective Instabilities in Two Liquid Layers

We perform linear stability calculations for horizontal fluid bilayers, taking into account both buoyancy effects and thermocapillary effects in the presence of a vertical temperature gradient. To help understand the mechanisms driving the instability, we have performed both long-wavelength and short-wavelength analyses. The mechanism for the large wavelength instability is complicated, and the detailed form of the expansion is found to depend on the Crispation and Bond numbers. The system also allows a conventional Rayleigh-Taylor instability if heavier fluid overlies lighter fluid, and the long-wavelength analysis describes this case as well. In addition to the asymptotic analyses for large and small wavelengths, we have performed numerical calculations using materials parameters for a benzene-water system.

effects and thermocapillary effects in the presence of a vertical temperature gradient. The problem depends on a number of dimensionless parameters, including the ratio of material properties in each layer, the Rayleigh number (dimensionless temperature gradient), the Marangoni number (dimensionless force due to surface tension gradients), the Bond number (dimensionless buoyancy force at the interface), and the Crispation number (a measure of interface deformability). We find that in some situations the onset of instability can be due to a overstable mode that oscillates in time, with a finite critical wavenumber. In other circumstances the instability is stationary in time with a finite critical wavenumber. The stationary instabilities can extend to small wavenumbers, though the critical instability usually does not occur at zero wavenumber. To help understand the mechanisms driving the instability we have performed both long-wavelength and short-wavelength analyses. The mechanism for the large wavelength instability is complicated, and the detailed form of the expansion is found to depend on the Crispation and Bond numbers. The system allows a conventional Rayleigh-Taylor instability if heavier fluid overlies lighter fluid, and this case is included in the results as well. In addition to the asymptotic analyses for large and small wavenumbers, we have performed numerical calculations using materials parameters for a benzenewater system.
The paper is organized as follows. Governing equations are given in the next section. The numerical procedure is described briefly in Sec. 3. Linear stability results for various cases are presented next, including a comparison of the numerical results with large and small wavenumber expansions. A discussion is presented in Sec. 5, followed by conclusions in Sec. 6. An appendix contains a summary of the expansion results.

Equations
We consider a semi-infinite horizontal bilayer system, with vertical heating across the layers. The unperturbed upper layer (denoted by α) extends over the interval 0 < z < H α , and the unperturbed lower layer (denoted by β) extends over the interval -H β < z < 0. Without loss of generality we consider linear stability results for a two-dimensional system. The horizontal coordinate extends from -∞ < x < ∞, and the velocity u has components in the x and z directions denoted by u and w, respectively.

Governing Equations in the Bulk
In each layer, we consider the Boussinesq equations Here p is the pressure, µ is the dynamic viscosity, T is the temperature, κ is the thermal diffusivity, g is the gravitational acceleration, t is the time, and ẑ is the unit vector in the z-direction (anti-parallel to gravity). We assume µ and κ are uniform in each layer, and also assume the density ρ is uniform in all terms except the gravitational term, where the density is given by in each layer. Here ρ − is the density in each layer at the reference temperature T R , and the thermal expansion coefficient η is assumed to be uniform in each layer.

Boundary Conditions
The upper boundary at z = H α and the lower boundary at z = -H β are assumed to be isothermal with noslip boundary conditions. The temperature is continuous across the interface, where T = T α -T β denotes the temperature jump across the interface. The tangential velocity is assumed to satisfy the no-slip condition where t is any unit vector tangent to the interface. The stress boundary condition is where n is a unit normal vector to the interface, T j k = -pδ j k + µ(∂u j /∂x k + ∂u k /∂x j ) is the stress tensor, γ is the surface tension, K is the curvature, υ n is the normal velocity of the interface, and ∇s is the surface gradient.
Here our sign convention is that the curvature K is defined to be positive for a spherical inclusion of β phase. For example, in two dimensions with an inter- ] 1/2 , and the surface gradient of the temperature-dependent surface energy γ = γ(T) is given by (8) where γ T = dγ/dT and t is the unit tangent vector to the interface in the direction of increasing x. Here h xx indicates the second derivative of h, etc.
The interface is a material surface, so that we have Under these conditions the left hand side of Eq. (7) vanishes. The continuity of heat flux gives (11) where k is the thermal conductivity and ∂T/∂n = erature field in each layer.

Base State
We linearize about a quiescent base state (also indicated by bars). The thermal field is (12) in the α layer, and (13) in the β layer, where T E is the unperturbed interface temperature. The temperature gradients in the base state satisfy (14) The pressure field in the base state is hydrostatic, with

Dimensionless Parameters and Governing Equations
Following the treatment by Zeren and Reynolds [9], we make the equations dimensionless based on a length scale given by the total depth d = H α + H β , a time scale based on the thermal time d 2 /κ β , a velocity scale κ β /d, a temperature scale G β d, and a pressure scale ν β κ β ρ − β /d 2 . These scales introduce the dimensionless parameters (16) (17) (18) and the geometrical parameter = -H β /H α . Here ν α = µ α / ρ − α and ν β = µ β / ρ − β are the kinematic viscosities in each layer, Pr is the Prandtl number, Ra is the Rayleigh number, Cr is the Crispation number, Bo is the Bond number, and Ma is the Marangoni number. The minus sign in the definition of Ma is conventional, since for most materials, as for the benzene-water system, we have γ T < 0.
We assume a horizontal wavenumber a and a temporal growth rate σ = σ r + iσ i . The perturbed quantities (indicated by tildes) then satisfy )ˆ, Critical conditions are often determined experimentally by varying the temperature gradient across the system. The temperature gradient G β appears in the dimensionless parameters Ma and Ra. Zeren and Reynolds [9] introduce the parameter (34) which is independent of the temperature gradient, and perform calculations by varying Ma for a fixed value of Γ.

Numerical Implementation
We solve the eigenvalue problem that governs the linear stability of the system by using two complementary procedures. In the first approach, the equations are discretized using pseudo-spectral Chebyshev collocation, and the resulting generalized matrix eigenvalue problem is solved using the package RGG from the EISPACK software library [10]. For a discretization with N degrees of freedom, this routine produces approximations to the first N eigenvalues of the system. The second approach is to use the two-point boundary value solver BVSUP [11], coupled with the root finder SNSQ [12], both from the SLATEC library [13], to implement a method described by Keller [14] to solve the eigenvalue problem. The BVSUP procedure provides a very accurate solution for a given eigenmode provided a good enough initial estimate is available for the root-finding procedure. The pseudospectral method is efficient for small values of N, and is well-suited for searching parameter space to detect real and complex eigenvalues. Rather than performing fine grid calculations with the pseudospectral procedure, however, the coarse grid results from the pseudospectral method are often used as initial guesses for the BVSUP code. Continuation from previous solutions also is used once an eigenmode has been identified.
The BVSUP software works in a single domain, so we have mapped the two layers to a common domain by setting (35) so that 0 < z -< H α in each layer. We then have To simplify the treatment of the problem, we also introduce an auxiliary ordinary differential equation in zfor the interface h , by setting (37) which allows us to avoid eliminating h as an unknown from the interface boundary conditions.

Results
In this section we present numerical results for the linear stability of the bilayer system. We provide a comparison with previous results of Zeren and Reynolds [9], who consider benzene overlying water. We use the thermophysical parameters given in their Table II for a temperature of 16°C, and also consider their case d * b = H β /(H α + H β ) = 0.4. Zeren and Reynolds consider both the case of heating from above (positive Ma), where the main driving force is the Marangoni effect and buoyancy effects are expected to be stabilizing, and the case of heating from below (negative Ma), where both buoyancy effects and Marangoni effects can produce instabilities.
If the system is heated from above, Zeren and Reynolds compute a positive critical Marangoni number of 1486, with a critical wavenumber of a = 2.6; they consider only stationary modes with σ i = 0. For a = 2.6 we find Ma = 1468.4635, in fair agreement with their results. We find a critical wavenumber to three digits of a = 2.66, with Ma = 1466.8951. Ferm and Wollkind [15] also obtained similar agreement in a comparison with Zeren and Reynolds, although they chose different thermophysical properties for the benzene-water system.
Our computed neutral stability curves for this case are shown in Fig. 1. As discussed by Zeren and Reynolds, the critical mode is mainly driven by Marangoni convection. We show the two stationary modes that are given in the second figure of Zeren and Reynolds; the convection pattern in these modes has a single cell in each layer, with negligible surface deformation. The convection is more concentrated near the interface for the higher wavenumber mode. We also find an oscillatory mode (σ i ≠ 0) at intermediate wavenumbers and Marangoni numbers. The oscillatory mode has a minimum near a = 5.5, with Ma = 22379.593 and σ i = ± 165.1. The surface deflection for this mode is negligible. The critical branch of the stationary mode does not extend to very small wavenumbers, but instead shows a large-Ma asymptote as the wavenumber tends to the value a = 0.12.
If the system is heated from below, Zeren and Reynolds compute a negative critical Marangoni number of Ma = -6068, with a critical wavenumber of a = 9. For a = 9.0 we find Ma = -6014.4082. We show our computed marginal stability curves for this case in Fig. 2. We find that the critical mode is actually oscillatory, with Ma = -3146.7277 for a = 4.6 and σ i = ± 64.93. The surface deflection for this mode is negligible. This mode exhibits a single convective cell in the upper layer, and two vertically-stacked cells in the lower layer of comparable strength. This mode merges with another stationary branch that has a longwavelength asymptote; σ i tends to zero on the oscilla-tory branch as the modes merge. On the long-wavelength branch, for a = 0.001 we find Ma ≈ -19318.282. Analytic results for long-wavelength modes are described in the appendix; for this case the asymptotic result is Ma = -19318.06, in excellent agreement with the numerical results.
There is a single convection cell in each layer for the long wavelength mode, with a significant interface deflection. Examination of the eigenmode shows that there is down-flow in the layer beneath the elevated portion of the perturbed interface. As discussed by Scriven and Sternling [16] for the single layer case, Marangoni modes and buoyant modes can be distinguished by the direction of the vertical flow in the layer beneath an interface elevation: there is downflow beneath elevations for Marangoni flow, and upflow beneath elevations for buoyancy-driven flow. This is consistent with the Zeren and Reynold's interpretation of the long-wavelength mode as driven by Marangoni effects. The other two stationary modes in Fig. 2 with minima near a = 10 have negligible surface deflection. For both modes the lower layer is nearly isothermal with a weak unicellular flow, and the upper layer exhibits two vertically-stacked convective cells for the mode with Ma ≈ -6000, and three vertically-stacked convective cells for the mode with Ma ≈ -50,000. There is in fact an entire family of additional stationary buoyant modes, not shown in Fig. 2, with larger values of | Ma |, each containing multiple vertically-stacked cells in the two layers.

Discussion
For the case of heating from below, as shown in Fig. 2 there is a small wavenumber instability that asymptotically approaches a finite Marangoni number as the wavenumber a tends to zero. To examine this behavior in more detail, we have performed additional calculations with Γ = 0, which eliminates the effects of buoyancy. Results are shown in Fig. 3 for both heating from below (Ma < 0) and for heating from above (Ma > 0). Suppressing buoyancy eliminates the oscillatory modes that prevail at intermediate wavenumbers for Γ = 0.142. For heating from below with our parameter values the stationary mode has a small-wavenumber limiting behavior that is given approximately by the expression (38) (see Eq. (51) in the appendix). Note that in this approximation there is a pole at a = 0.335. This is in approximate agreement with the vertical asymptote that is obtained numerically for a = 0.268 for heating from either above or below. In fact, plotting instead 1/Ma versus a produces a single smooth curve crossing the y-axis with 1/Ma = 0 at a wavenumber a = 0.268.
To help further understand the low wavenumber instability that is observed numerically, we have also performed numerical computations and asymptotic expansions that illustrate the effects of the Bond number Bo and the Crispation number Cr for this mode. In Fig. 4 the upper curve corresponds to the data given in Table I. This curve asymptotically approaches a small wavenumber limit for Ma = -2.27 × 10 4 . The lower curve corresponds to setting Bo = 0. The symbols on the curves correspond to numerical results, and the curves themselves correspond to analytical results from the small wavenumber approximation given in the appendix. The small wavenumber results depend strongly on both Bo and Cr. For the parameter values given in Table I Here the coefficients c 1 , c 2 , c 4 , and c 6 in these expressions generally depend on the layer geometry and the remaining material constants and are given in the appendix; here we have evaluated these expressions for the values given in Table I.

Large Wavenumbers
The temperature and flow field for the Marangoni instability can be studied in the large wavenumber limit that is summarized in the appendix. We ignore buoyancy by setting Γ = 0, and also ignore interface deformation by setting Cr = 0, which the numerical results indicate are good approximations in this case. In the large wavenumber limit the flow is confined to the vicinity of the interface and the effect of the upper and lower boundaries at z = H α and z = -H β is negligible. To help visualize the flow we introduce a two-dimensional streamfunction of the temperature and streamfunction near the inter-face are shown in Fig. 6. Here we have exaggerated the size of the perturbation to emphasize the distortion of the isotherms near the interface. To make the plots easier to interpret we have assumed equal thermal conductivities (k * = 1), which equalizes the unperturbed temperature gradients in the two layers and facilitates comparison of the perturbed temperature fields in each layer. There is a temperature gradient along the interface (at z = 0), and the streamlines of the flow are along the interface. Over the single period of the flow shown in Fig. 6, there are four convective cells near the interface that alternately compress and expand the isotherms near the interface. For our parameters with κ * < 1, the Marangoni instability occurs for heating from above (Ma > 0) in the large wavenumber limit. If κ * > 1 the instabilities occur for heating from below (Ma < 0) instead.

Rayleigh-Taylor Instability
The two-layer system can exhibit the classical Rayleigh-Taylor instability if heavier fluid overlies lighter fluid. In the simplest case the Rayleigh-Taylor instability can be understood by a simple potential energy argument that balances the increased surface energy of a deformed interface y = h(x) against the change in the gravitational potential energy of the displaced fluid,  Table I. Dimensionless parameters for the benzene (α layer) and water (β layer) system at T = 16°C, from Ref. [9] ratio of densities ρ * = ρα /ρ - which can be seen as a factor in the normal stress balance boundary condition (29).
In the situation we have studied above, we have lighter benzene overlying heavier water, so the Rayleigh-Taylor instability does not occur. To study this instability for our system with a minimal change in notation, we temporarily choose to change the direction of gravity while keeping benzene and water in the original orientation, so that the water and benzene are unstably stratified with respect to gravity. We take G β < 0, so that buoyancy has a stabilizing effect on the system; the resulting sign conventions produce Ma < 0, Ra > 0, Γ < 0, and Bo < 0. In Fig. 7 we show the corresponding numerical results for the Rayleigh-Taylor instability. The dotted curve in Fig. 7 shows the curve -Bo = a 2 /(1ρ*) that holds for Ra = Ma = 0. The solid curve shows numerical results for Ra = 1.0 × 10 5 and Ma = 0. The stabilizing effect of buoyancy is evident at small wavenumbers, where the system is then stable if | Bo | is sufficiently small. The marginal stability curve asymptotes at small wavenumbers to the value (44) in our case, which follows from Eq. (45) in the appendix. The dashed curve shows numerical results for Ra = 0 and Ma = -1.0 × 10 5 . The destabilizing Marangoni effect in this configuration is evident at small wavenumbers. and there is a cut-off wavenumber at a = 0.333 where Bo = 0; at lower wavenumbers this mode merges into the Marangoni instability with positive Bond numbers. This value of the cut-off wavenumber is too large to be quantitatively approximated by the small wavenumber expansion given by Eq. (51) in the appendix. For the benzene-water system with Ra/Ma = -0.142, numerical results are qualitatively similar to those for the case with Ra = 0.
Our small wavenumber expansion of the instability, as summarized in the appendix, corrects and extends a similar analysis given by Zeren and Reynolds [9], who note that at leading order in a only the interface position and temperature field are perturbed, with no leading order velocity profiles. However, the asymptotic values for the critical Marangoni number do depend on the viscosity of the fluid, so flow does have an effect on the stability of the system, even at small wavenumbers. The small wavenumber analysis proceeds by an expansion in powers of a 2 , and the critical Marangoni number is determined at a later point in the expansion where the O(a 2 ) flow perturbation is significant. We also note that the expansion procedure is based on a limit of small but non-zero wavenumbers. In the special one-dimensional case a = 0 the perturbed interface is flat and mass conservation restricts its location to z = 0. In a sense this represents a discontinuity in the problem in the limit a → 0, since in the expansion procedure a shift in the interface location is obtained at leading order. This is somewhat analogous to the dependence of a perturbed interface of the form z = H 0 cos a x, which conserves mass for a ≠ 0 (i.e., the average interface position is zero), but represents a shifted flat interface at z = H 0 for a = 0. In general, for real experiments the smallest allowable value of a will be limited by the finite lateral extent of the system, but whether or not a perturbation with a = 0 is physically relevant depends on the specific problem under consideration.

Conclusions
We have performed linear stability calculations for horizontal fluid bilayers, taking into account both Marginal stability curves for a water-benzene system that is unstably stratified with respect to gravity. The dotted curve (Ma = Ra = 0) corresponds to the classical Rayleigh-Taylor instability in the absence of buoyancy, given by -Bo = a buoyancy effects and thermocapillary effects. We consider the case of a lighter layer overlying a heavier layer, so that the base state is stably stratified in this sense. We find that the system can be linearly unstable to either heating from above (Ma > 0) or below (Ma < 0).
The mechanism for the large wavelength (small a) instability is complicated (see Figs. 4 and 5). The detailed form of the expansion, viz, the exponent n in the leading order expansion Ma ∼ a n depends on the Crispation and Bond numbers.

Appendix
Here we consider the limits of large and small wavenumbers. We introduce the dimensionless layer widths

Small Wavenumbers
By performing long-wavelength asymptotics for Bo ≠ 0 and Cr ≠ 0 we find the leading order result This expression is in good agreement with the numerical results for small wavenumbers as shown in Fig. 2. For Ra = 0, this expression agrees with that given by Zeren and Reynolds [9] up to a difference in sign; their expression also agrees with their numerical results for Ra = 0 if the sign of their expression is changed. For Ra ≠ 0 our expression gives results which differ in both magnitude and sign from that given by Zeren and Reynolds. We note that the leading order result in Eq. (45) is proportional to Bo and inversely proportional to Cr, so that the special cases Bo = 0 or Cr = 0 require additional attention; we consider these cases in more detail below for Ra = 0.
For the values given in Table I,  If we set Ra = 0 the linear eigenfunctions for the general problem can be written down explicitly, and the dispersion relation for σ = 0 can be evaluated symbolically in closed form, though the resulting expression is lengthy and difficult to interpret. However, the result can be expanded for small wavenumbers, and produces expressions that are not too unwieldy.
The general expression for σ = 0 and Ra = 0 can be found by expanding in the wavenumber a to give

Large Wavenumbers
We next consider the limit of large wavenumbers for a stationary mode. For our system the numerical results suggest that buoyancy effects and interface deformation are unimportant in this limit, so we also consider the formal limit of small Crispation number, Cr → 0, along with Ra = 0. For Ra = 0, the governing equations for the velocity field are decoupled from the thermal field, which simplifies the analysis. For Cr → 0 the dimensionless form of the normal momentum balance, then reduces to (58) For Bo (ρ* -1) -a 2 ≠ 0, we conclude that the interface deformation vanishes, h ∼ = 0.
In the limit of large wavenumber the disturbances are concentrated near the interface and the effects of the outer boundaries are insignificant. The appropriate solution can then be computed by applying decay conditions in an unbounded domain as z → ± ∞. The vertical components of the velocity field are given by The interface boundary conditions determine the remaining six constants A α , B α , C α , A β , B β , and C β .
The velocity boundary conditions require A α = A β = 0, and B α = B β , and the solution to the thermal problem then gives (67) The remaining boundary condition gives the dispersion relation, For the values in Table I