The instability of liquid films with temperature-dependent properties flowing down a heated incline

1.   Introduction

Thin liquid films with a free surface flowing down an inclined plane are subject to an interfacial instability brought about by inertial effects, which generates waves propagating along the surface. An illustration of this occurrence can be observed on rain water draining along sloped pavement. This phenomenon is actually encountered in many different settings where it has important consequences. In industrial processes interfacial instability is a crucial factor, for example, in coating applications and the effective operation of various equipment involving heat and mass transfer. In the environment, surface waves in lava and debris flows can coalesce into large bores with enough force to cause serious damage.

Due to the prevalence and significance of the instability of gravity-driven flows, researchers have shown considerable interest in predicting the conditions for the onset of instability, and determining the structure of the resulting waves. Groundbreaking experimental and theoretical investigations have been carried out by Kapitza ^[1] and Kapitza and Kapitza ^[2] in the late 1940s. They established a dimensionless group, now known as the ``Kapitza number", combining the surface tension and the viscosity of the fluid. They obtained a relationship between this quantity and the inertial effects, specified by the Reynolds number, describing the conditions for the onset of instability. A more systematic theoretical study, based on a linear stability analysis, was later undertaken by both Benjamin ^[3] and Yih ^[4] who found that the critical Reynolds number is actually proportional to the cotangent of the angle of inclination. Benney ^[5] exploited the relative thinness of the liquid layer and carried out an asymptotic analysis which yielded a single evolution equation for the thickness of the layer. Gjevik ^[6] improved the accuracy of the asymptotic expansions and obtained an equation for the evolution of the thickness, which is referred to as the Benney equation. This equation suffers from an artificial singularity for sufficiently large Reynolds numbers, however it proves to be useful for conducting weakly nonlinear analyses.

If the inclined fluid film is also heated then the flow is nonisothermal and an additional source of instability exists ^[7]. This is caused by the fact that temperature differences lead to nonuniformity in surface tension generating so-called Marangoni stresses, which can lead to fluid motion along the surface strong enough to destabilize the equilibrium flow. This phenomenon is referred to as thermocapillarity or the Marangoni effect. Kalliadasis et al. ^[8] employed an integral-boundary-layer approximation to study the stability of long-wavelength perturbations in film flow on a heated incline with negligible evaporation. This approximation allows for a nonlinear analysis, however it turns out to provide somewhat inaccurate predictions for the onset of instability when compared to experimental observations and the linear stability analysis of the full equations. Ruyer-Quil et al. ^[9] and Scheid et al. ^[10] rectified this shortcoming by employing a weighted residual method to obtain a reduction of the long-wave equations. Later, Trevelyan et al. ^[11] applied the weighted residual approach and studied both the case of a prescribed temperature at the bottom of the liquid film and the case of a prescribed heat flux. The effect of evaporation in heated inclined flow was taken into account by Lin ^[12].

In reality, in a nonisothermal flow the properties of the fluid vary due to changes in temperature. Kabova and Kuznelsov ^[13] have incorporated a temperature dependent viscosity into a steady flow model, but did not investigate the stability of the flow. Goussis and Kelly ^[14] and Hwang and Weng ^[15] have implemented an unsteady flow model for nonisothermal inclined film flow with variable viscosity and carried out stability analyses, however the temperature of the free surface is assumed to be prescribed. This is in fact a non-realistic condition, and furthermore it completely negates the Marangoni effect.

The reason why such little attention has been paid to temperature variation in the fluid properties in film flow is the assumption that the effect is negligible in comparison to that due to thermocapillarity. However, a systematic investigation of this issue has not been reported until the work of Ogden et al. ^[16]. They demonstrate that extremely weak variations in fluid properties must be assumed in order to prove that the effect is negligible. More precisely, the rate of variation must be of the same order of magnitude as the wavenumber of the very long perturbations that destabilize the flow. Pascal et al. ^[17] have developed a model with variation in all fluid properties: not just viscosity, but also mass density and thermal conductivity. The model captures thermocapillarity by implementing Newton's law of cooling at the surface of the liquid layer. Pascal et al. apply this model to investigate how the Marangoni effect enhances inertial instability. They perform a linear stability analysis, with asymptotic solutions being obtained under the assumption of perturbations with long wavelength (small wavenumber) and small rates of variation in the fluid properties, that are however asymptotically larger than the wavenumber of the perturbations. This investigation was subsequently revisited by D'Alessio et al. ^[18] who considered the more general case with no restriction on the rate of variation in the fluid properties. Just recently, Chhay et al. ^[19] have resorted to the long-wave approximation and developed a reduced dimesionality model which includes a temperature dependent viscosity. They obtained numerical simulations for the nonlinear evolution of the flow.

In this paper we use the full equations to study the stability of the flow of a liquid film with temperature dependent fluid properties down a heated inclined plane. We extend the work of D'Alessio et al. ^[18] in two important ways. First of all, our investigation is not restricted to perturbations of very long wavelength. Secondly, we investigate the coupling of the Marangoni and inertial effects like D'Alessio et al., however we also consider instability at low levels of inertia which is precipitated by the Marangoni effect alone. The paper is arranged as follows. In Section 2 we determine the equations governing the flow. In Section 3 we obtain an equilibrium solution and study its stability by means of a linear stability analysis. This involves solving an eigenvalue problem with Orr-Sommerfeld type equations, which we accomplished numerically. Then, in Section 4 we present some results and interpret our findings. Comparisons with nonlinear simulations are also made. Finally, in Section 5 we highlight conclusions drawn from our investigation.

2.   Governing equations

Here we consider the laminar gravity-driven flow of a thin liquid layer with Newtonian rheology down a heated surface inclined at an angle $\beta$ with the horizontal as depicted in Figure 1. The flow is assumed to be two dimensional and we set up a rectangular coordinate system with the $x$ axis along the incline and pointing in the downslope direction, and the $z$ axis pointing into the fluid. We denote by $u$ and $w$ the velocity components in the $x$ and $z$ directions respectively, and $z = h(x, t)$ is the position of the free surface. The temperature of the surrounding gas is $T_a$ and the inclined plane (substrate) is maintained at a higher temperature $T_b$.

Equations of motion for the flow are obtained from conservation of mass and momentum and in the general case of variable fluid properties and are given by ^[20]

where $\frac{D }{D t} = \frac{\partial }{\partial t} +u\frac{\partial}{\partial x} +w\frac{\partial}{\partial z}$ is the material derivative, $p$ is the pressure, $g$ is the acceleration due to gravity, $\rho$ is the mass density of the fluid and $\mu$ is the viscosity. Conservation of energy yields the equation

where $e$ is the internal energy, $T$ is the temperature of the fluid and $K$ is the thermal conductivity.

We assume the fluid properties to be temperature dependent, and specifically consider the following linear relations

where $\sigma$ is the surface tension and all the rates of variation: $\hat{\alpha}$, $\hat{\lambda}$, $\hat{\Lambda}$ and $\gamma$ are assumed to be nonnegative constants.

In order to obtain the governing equations for our investigation we apply the Boussinesq approximation by which we assume the mass density to be constant except when in gravitational terms. We thus obtain

where the kinematic viscosity is given by $\nu = \frac{\mu}{\rho_0} = \frac{\mu_0}{\rho_0}-\frac{\hat{\lambda }}{\rho_0} \left(T-T_a\right)$, and $c_p$ is the specific heat.

Next, we establish the conditions at the interfaces of the liquid layer: the fluid-gas interface at the free surface $z = h(x, t)$, and the fluid-solid interface at the bottom, $z = 0$. At $z = 0$ we apply the no-slip and no-penetration conditions, as well as the prescribed temperature

At the free surface we assume the ambient gas to be dynamically passive and the force exerted on the surface by the flow to be balanced by surface tension alone. Without loss of generality we take the constant ambient pressure to be zero, since only its gradient ends up being involved in the governing equations. The normal component of the force balance is then given by

while the tangential component is expressed as

where both equations are evaluated at $z = h(x, t)$.

Under an assumption of negligible evaporation or condensation, a kinematic condition for the free surface reads

The relation between the heat flux normal to the surface and the difference in the temperature of the liquid and the ambient gas can be expressed through Newton's law of cooling which is given by

where $\chi$ is the heat transfer coefficient.

The governing equations can be made dimensionless by applying an appropriate scaling. For the length scale we employ

which is the Nusselt thickness of an isothermal flow produced by a prescribed flow rate $Q$. We introduce the following transformation

where $\Delta T = T_b - T_a$ and $U = \frac{Q}{H}$.

Applying the above scaling in equations (2.1)–(2.9) and discarding the asterisks for notational convenience we get

and at the free surface, $z = h\left(x, t\right)$,

where $F = 1 + \left(\frac{\partial h}{\partial x}\right)^2 \; , \; G = 1 - \left(\frac{\partial h}{\partial x}\right)^2 \; ,$ while at $z = 0$ the conditions are

and

In the above equations the parameters are the Reynolds number, $Re = \frac{\rho_0 UH}{\mu_0}$, the Prandtl number, $Pr = \frac{ \mu_0c_p}{K_0}$, the Weber number, $We = \frac{\sigma_0}{\rho_0 U^2 H}$, the scaled heat transfer coefficient, $B = \frac{\chi H}{K_0}$ and the scaled rates of variation of the fluid properties: $M = \frac{\gamma \Delta T}{\rho_0 U^2 H}$, $\alpha = \frac{\hat{\alpha}\Delta T}{\rho_0}$, $\lambda = \frac{\hat{\lambda}\Delta T}{\mu_0}$, and $\Lambda = \frac{\hat{\Lambda}\Delta T}{K_0}$.

3.   Linear stability analysis

Our mathematical model given by equations (2.10)–(2.19) admits a solution that is independent of $x$ and $t$, which corresponds to a steady flow that is uniform in the streamwise direction. We study the linear stability of this solution, which we will refer to as the equilibrium flow.

In order to calculate the equilibrium solution we set the $x$ and $t$ derivatives to zero in the governing equations and set the thickness of the film to be $h\equiv 1$. Solving equation (2.10) together with condition (2.18) we find the equilibrium cross-stream velocity to be $w_s(z)\equiv 0$. We then obtain the following problem governing the equilibrium temperature distribution, $T_s(z)$, from equations (2.13), (2.17) and (2.19)

where $D\equiv \frac{d}{dz}$. Solving this problem we obtain

where

To get the equilibrium streamwise velocity, $u_s(z)$, we solve the following problem obtained from equations (2.11), (2.15) and (2.18)

and obtain

where

Now, from equations (2.12) and (2.14) we get

and

which can be solved to yield the equilibrium pressure distribution

In accordance with the procedure for linear stability analysis we consider the perturbed equilibrium solution expressed as

where $\tilde{u}, \tilde{w}, \tilde{p}, \tilde{T}$ and $\tilde{\eta}$ are the infinitesimal perturbations. Introducing this perturbed state into equations (2.10)–(2.13) and linearizing with respect to the perturbations we get

Transferring the boundary conditions at $z = 1 + \tilde{\eta}$ to $z = 1$ and linearizing yields the following conditions at $z = 1$

While at $z = 0$ the conditions are

Now we employ normal modes into the linearized perturbation equations expressed as

where $k$ is the wavenumber of the perturbation and $c$ is a complex quantity. The real part of $c$ is the phase velocity and the imaginary part multiplied by $k$ gives the temporal growth rate. In terms of the normal modes we get

with the conditions at $z = 1$ being

and the ones at $z = 0$ written as

We now eliminate $\hat{p}$ from the system and introduce the stream function, $\psi$, which is related to the velocity perturbations by

Introducing normal modes for the stream function expressed as

leads to the following Orr-Sommerfeld type equation

The interface conditions at $z = 1$ become

while at $z = 0$ we have

The problem given by equations (3.1)–(3.7) constitutes an eigenvalue problem with $c$ being the parameter that is to be assigned characteristic values. Solving for $c$ provides the growth rate of the perturbation with wavenumber $k$ for a given set of flow parameters. A positive value for the imaginary part of $c$ indicates that the perturbation amplitude grows in time, while if negative the perturbation is damped.

We will solve the eigenvalue problem numerically by using a Chebyshev spectral collocation method. We first transform the interval $z\in[0, 1]$ into $\xi\in[-1, 1]$ by letting $z = \frac{1}{2} \left(\xi +1\right)$. We then discretize the $\xi$ interval by the points

and expand the dependent variables in terms of the cardinal functions as

where

This leads to an algebraic system of the form $Ly = c\mathcal{M}y$ where

and $L$ and $\mathcal{M}$ are $(2N+3)\times(2N+3)$ matrices.

The eigenvalues of the system were calculated using the Matlab subroutine eig. For very small values of $k$ the solution is in excellent agreement with the asymptotic results as $k\rightarrow 0$ obtained by D'Alessio et al. ^[18]. Our numerical method allows us to determine the instability of perturbations of any wavenumber. Furthermore, it turns out that the parameters $We$, $M$ and $B$ are implicitly dependent on the Reynolds number. We thus introduce new parameters $Ka$, $Ma$ and $Bi$ which are independent of the Reynolds number and related to the previous ones by

where $Ka$, $Ma$ and $Bi$ are defined as

In terms of these control parameters we are able to capture not only the inertial instability of the flow, but also the pure thermocapillary instability which occurs at low inertia.

4.   Results

We now present and interpret results from our analysis. We first point out that since the variation with temperature of the density and viscosity is assumed to be governed by

and since the maximum of $T$ is $1$, the parameters $\alpha$ and $\lambda$ must be restricted to the interval $[0, 1)$.

We begin by considering the state of neutral stability, which refers to the set of parameter values for which a particular perturbation experiences zero growth in time. In Figure 2 we show neutral stability curves in the $(Re, k)$ plane. These curves delineate regions in the plane where perturbations are amplified or damped. As it can be seen, for sufficiently small Marangoni numbers the curves consist of two branches. The one on the left separates a region of perturbation amplification for the smaller Reynolds numbers from one where perturbations are damped. This kind of instability is referred to as the S mode and is generated by perturbations in the temperature of the free surface of the film. These cause nonuniformity in surface tension which generates Marangoni stresses which pull fluid along the surface toward regions where surface tension is stronger. This fluid motion, if amplified, destabilizes the equilibrium flow. The fluid motion induced by the Marangoni stresses is impeded by inertia, therefore, increasing the Reynolds number, which measures the level of inertia, plays a stabilizing role.

The right hand branch of the neutral stability curve corresponds to the amplification of perturbations in the position of the surface of the film and is referred to as the H mode of instability. In this case inertia is a destabilizing agent. It decelerates the flow ahead of the crests of surface waves and accelerates it ahead of the troughs. Fluid is thus pulled from beneath the troughs and piles up under the crests, and as a result wave amplitudes are amplified. Consequently, the region of instability for the H mode lies to the right of the neutral stability curve where the Reynolds number is larger. Now, the proximity of the troughs to the heated substrate weakens surface tension and as a result fluid is pulled by Marangoni stresses toward the crests, which also works to amplify the waves. Therefore, the Marangoni effect augments the H-mode instability.

Under flow conditions for which all perturbations are damped the equilibrium flow is stable. For sufficiently weak thermocapillarity, measured by the Marangoni number, there is an interval of Reynolds numbers for which the flow is stable. The endpoints of this interval are the critical Reynolds numbers for the onset of instability in the equilibrium flow and they correspond to the intercepts of the neutral stability curve with the $k = 0$ axis, i.e. the threshold of instability for the S and H modes of infinitely long perturbations. As the Marangoni number is increased both modes are destabilized and the two branches eventually merge. Beyond this point, the equilibrium flow is unstable for all Reynolds numbers, i.e. all levels of inertia. We can thus conclude that since the instability of infinitely long perturbations generate the instability of the flow, the asymptotic results as $k\rightarrow 0$ for the critical Reynolds numbers obtained in ^[18] accurately describe the instability of the flow due to the H mode, and the general conclusions drawn there are correct and consistent with the interpretation of the results made in our present investigation. However, here, in presenting the results we are also mindful of the critical Marangoni number for which the two modes merge. For larger values the $Re_{crit}$ results obtained in ^[18] do not actually represent a threshold for the instability of the flow.

We now focus on determining how various parameters affect the critical Reynolds number for the onset of instability in the equilibrium flow. To this end we construct, as exemplified in Figure 3, the graph of $Re_{crit}$ as a function of $Ma$ for different parameter values. The upper branch of these curves indicates the Reynolds number for the onset of the H mode, while the lower branch is related to the S mode. The region inside the curve corresponds to the interval of Reynolds numbers for which the flow is stable.

We begin by investigating how the variable density affects stability. The relevant parameter in this regard is $\alpha$, the scaled rate of change with temperature. Now, the temperature is a function of $z$, i.e. varies with depth in the fluid layer. So, increasing $\alpha$ lowers the density at every depth, but at the same time amplifies its gradient with $z$. We expect these two outcomes to cause competing instability mechanisms. A decrease in density results in less inertia in the flow which stabilizes the H mode and destabilizes the S mode. A greater variation with $z$ on the other hand magnifies the stratification of the fluid layer and thus causes greater density differences along the surface if its height is not uniform due to wavy perturbations. As a result, we get an enhanced momentum variation along the surface which generates impulses creating surges and destabilizing the H mode as well as the S mode. In addition, the stratification is top heavy and consequently the crests of waves have a higher flow rate which acts to increase the amplitude and contribute to instability. In Figure 3 we present $Re_{crit}$ vs. $Ma$ curves for different values of $\alpha$. Clearly, as $\alpha$ is increased the interval of $Re$ for which the S mode is unstable widens which confirms the expectation spelled out above that the S mode is destabilized by a variable density.

Focusing on the H mode, the results in Figure 3 demonstrate that for sufficiently small $Ma$ values the critical Reynolds number increases with $\alpha$ indicating that the stabilizing mechanism is dominant. For the larger Marangoni numbers however, the upper branches of some of the curves cross and interchange position specifying that $Re_{crit}$ changes from an increasing to a decreasing function of $\alpha$. This is more clearly shown in Figure 4 which contains the graph of the critical Reynolds number for the H mode as a function of $\alpha$ for some of the larger $Ma$ values from Figure 3. We see here that $Re_{crit}$ is in fact a decreasing function of $\alpha$ for certain ranges of this parameter. In these cases the destabilizing mechanism associated with varying density surpasses the stabilizing one. We thus infer that the Marangoni effect couples with the destabilizing action of variable density and if sufficient makes it dominant.

We next look at the effect of a variable viscosity, which is measured by $\lambda$, the scaled rate of variation with temperature. Increasing this parameter lowers the viscosity at every point in the flow. The results in Figure 5 indicate that doing this destabilizes the H mode. The explanation is that reducing viscosity leads to faster flow rates. This accentuates the acceleration/deceleration action of inertia associated with variation in the thickness of the fluid layer. The results also indicate that the effect on the S mode is negligible unless thermocapillarity is sufficiently strong. In this case the S mode is destabilized possibly due to the fact that a lower viscosity facilitates fluid motion along the surface.

We next turn to investigating how heat transfer across the free surface, quantified by the Biot number, affects stability. If $Bi = 0$ the free surface is thermally insulated and consequently at equilibrium the temperature of the fluid, including the surface, is constant ($T_s(z)\equiv 1$). Surface tension is then uniform and the Marangoni effect is neutralized. Similarly, as $Bi\rightarrow\infty$ the temperature of the surface approaches that of the ambient gas and again the surface tension is constant. As the Biot number is increased from zero, thermocapillarity becomes more significant and there is a critical point where the Marangoni effect is maximized and the flow conditions are most unstable. This is illustrated in Figure 6 which shows that as $Bi$ is increased the region of stability first shrinks, but eventually starts to grow.

It turns out that viscosity variation has an interesting qualitative effect on the evolution of $Re_{crit}$ with $Bi$ described above for the constant viscosity case. This is apparent by comparing the curves in Figure 6 to those in Figure 7 where $\lambda = 0.5$. It can be seen that with viscosity variation, for small Marangoni numbers $Re_{crit}$ for the H mode does not attain a minimum, but rather increases monotonically with $Bi$. In order to physically explain the correlation between heat transfer across the surface and variable viscosity, we point out that if $Bi = 0$, in which case $T_s(z)\equiv 1$, then the viscosity is reduced by the maximum amount. The flow rate is then maximized and as such with $Bi = 0$ we have the most unstable configuration from the point of view of inertial instability. In this regard then, increasing $Bi$ has a stabilizing effect. But doing this also results in greater surface temperature variation. However if the nature of the fluid is such that surface tension is only weakly dependent on temperature changes, signified by low $Ma$ values, the Marangoni instability is dominated by the stabilizing effect of reduced inertia.

Now we consider the effect of variation with temperature in thermal conductivity. Figure 8 contains $Re_{crit}$ vs. $Ma$ curves for different values of $\Lambda$ which measures the rate of change in thermal conductivity. The results indicate that for both the S and H modes, $Re_{crit}$ is a nonmonotonic function of this parameter, however the variation in $Re_{crit}$ is very small.

We next draw our attention to nonlinear effects and begin by discussing some approaches that can be used to account for nonlinear phenomena. One method involves solving the so-called Benney equation. Benney ^[5] proposed an asymptotic expansion procedure to derive a nonlinear evolution equation for the free surface of a fluid layer which has been used extensively by numerous researchers to investigate the stability of the flow. One of the first numerical studies using the Benney equation was conducted by Gjevik ^[6] to investigate finite amplitude effects. Another commonly used approach involves the integral-boundary-layer (IBL) formulation. The general idea with this technique is to first simplify the governing equations by discarding negligible terms based on a scaling argument. The simplified equations are then integrated across the fluid layer to eliminate the cross-flow variation. A parabolic velocity profile resulting from the uniform steady flow is inserted when necessary in order to perform the integration together with appropriate boundary conditions. The IBL model usually consists of a system of equations for the fluid thickness and flow rate. An improvement to the IBL approach was advanced by Ruyer-Quil and Manneville ^[21], ^[22] which involves multiplying the equations by carefully selected weight functions before integrating across the fluid layer. This yields a similar system of equations referred to as the weighted-residual model (WRM). As noted in the introduction, the advantage of the WRM is that it correctly predicts the critical Reynolds number for the onset of instability. Lastly, Chhay et al. ^[19] recently proposed a new asymptotic approach to model the heat transfer through a liquid film flowing down a vertical wall.

For our nonlinear analysis we have decided to implement a weighted-residual model. The details surrounding the numerical solution procedure are outlined in ^[23]. Here, our focus is on comparing the results from the nonlinear simulations with those obtained from the linear stability analysis. Simulations were conducted on a computational domain of length, $L$, subject to periodic boundary conditions. Although the domain length is arbitrary, it was observed that $L = 2$ was sufficiently large to trigger the long-wave instability. The evolution of the unsteady flow was computed by imposing small perturbations on the initial steady state solutions. By monitoring the growth of the disturbances as the perturbed solutions were marched in time we were able to estimate the onset of instability by incrementing the Reynolds number and noticing the first instance when the flow develops a permanent series of waves on the free surface. Contrasted in Figure 9 are $Re_{crit}$ values based on our linear stability analysis and our nonlinear simulations for the case having the following parameter values: $Pr = 7$, $B = \cot\!\beta = 1$, $\alpha = \lambda = \Lambda = 0.1$ and $We = 50$ over the interval $0 \leq M \leq 10$. From this we observe that linear theory is in close agreement with the fully nonlinear simulations.

5.   Conclusions

In this work we have considered the flow of a thin liquid film down a heated incline. We conducted a theoretical study of the stability of a steady equilibrium flow that is uniform in the streamwise direction. Fully nonlinear numerical simulations were also conducted and the agreement with the linear theory was good. Our mathematical model includes the Marangoni effect and captures the two related types of instability: The so-called H mode, occurring at larger Reynolds numbers, is the amplification caused directly by inertial factors and enhanced by the Marangoni effect. And the S mode which is observed at lower Reynolds numbers. Under these conditions, inertial effects are weak and surface waves related to flow perturbations are damped by gravity. If these effects are sufficiently weak however, significant fluid motion can be generated along the surface by the Marangoni effect, altering and destabilizing the equilibrium flow.

Our model also accounts for the dependence on temperature of the properties of the fluid. Besides changes in surface tension with temperature which leads to thermocapillarity, we have also included variation in mass density, viscosity and thermal conductivity. These quantities were taken to vary linearly with temperature, the rates of change being the control parameters.

The density was assumed to decrease with temperature, meaning that the fluid experiences thermal expansion, with the control parameter being the thermal expansion coefficient. Increasing this parameter lowers the density but magnifies the stratification of the fluid layer. Our results indicate that the S mode is destabilized due to the fact that the decrease in density facilitates the fluid motion along the surface driven by the Marangoni stresses. The effect of thermal expansion on the H mode is more complicated, as it triggers competing stability mechanisms. The enhanced stratification is a destabilizing factor, whereas the decrease in density diminishes the flow rate in the equilibrium flow which has a stabilizing effect. For cases where thermocapillarity is sufficiently strong, the destabilizing factors are the victors in this competition, and thermal expansion destabilizes the flow. If thermocapillarity is weak, the stabilizing factor is dominant.

In our investigation the viscosity of the fluid was also allowed to decrease with temperature. Increasing the rate of change allows the equilibrium temperature to induce a greater reduction in viscosity, which leads to faster flow rates and magnifies the inertial effects which destabilize the H mode. In addition, the varying viscosity was found to also cause a reversal in the destabilizing role played by the transfer of heat across the surface from the liquid to the ambient gas. When the heat transfer coefficient is small, and the rate of heat flow across the surface is low, an increase in this rate leads to a greater temperature variation on the surface and the flow is destabilized due to a stronger Marangoni effect. However, another consequence is a diminished reduction in viscosity in the bulk of the fluid layer which slows the equilibrium flow rate and has a stabilizing effect.

Acknowledgments

Financial support for this research was provided by the Natural Sciences and Engineering Research Council of Canada and the Faculty of Mathematics at the University of Waterloo.

Conflict of interest

The authors declare that there is no conflict of interest in this paper.