3.1 Steady solution

The mass conservation equations admit steady solutions that satisfy

(3.0a,b ) $$\begin{eqnarray}\boldsymbol{q}_{\boldsymbol{l}}={\mathcal{Q}}\boldsymbol{e}_{\boldsymbol{x}},\quad \boldsymbol{q}=\boldsymbol{e}_{\boldsymbol{x}},\end{eqnarray}$$

along with (3.1). These equations can be solved analytically to yield the two-dimensional steady solutions

(3.1) $$\begin{eqnarray}(H_{s},h_{s})=(A,a)(1-x)^{1/4}.\end{eqnarray}$$

The constants $A$ and $a$ depend on the dimensionless parameters through the algebraic conditions given in (A 1) and (A 2) in appendix A.

3.2 Linear perturbation equations

We investigate the linear stability of the steady basic state of § 3.1 by introducing small disturbances $h_{p}$ , $H_{p}\ll 1$ . These give rise to the following linearised perturbation fluxes:

(3.2) $$\begin{eqnarray}\displaystyle \boldsymbol{q}_{\boldsymbol{l}\boldsymbol{p}} & = & \displaystyle -\left[\frac{{\mathcal{M}}}{3}h_{s}^{3}({\mathcal{D}}\unicode[STIX]{x1D735}h_{p}+\unicode[STIX]{x1D735}H_{p})+\frac{{\mathcal{M}}}{2}h_{s}^{2}(H_{s}-h_{s})\unicode[STIX]{x1D735}H_{p}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.{\mathcal{M}}h_{s}^{2}h_{p}({\mathcal{D}}\unicode[STIX]{x1D735}h_{s}+\unicode[STIX]{x1D735}H_{s})+\frac{{\mathcal{M}}}{2}(h_{s}^{2}(H_{p}-h_{p})+2h_{s}h_{p}(H_{s}-h_{s}))\unicode[STIX]{x1D735}H_{s}\right],\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle \boldsymbol{q}_{\boldsymbol{p}} & = & \displaystyle -\left[\frac{1}{3}(H_{s}-h_{s})^{3}\unicode[STIX]{x1D735}H_{p}+\frac{{\mathcal{M}}}{2}h_{s}^{2}(H_{s}-h_{s})({\mathcal{D}}\unicode[STIX]{x1D735}h_{p}+\unicode[STIX]{x1D735}H_{p})\right.\nonumber\\ \displaystyle & & \displaystyle +\,{\mathcal{M}}h_{s}(H_{s}-h_{s})^{2}\unicode[STIX]{x1D735}H_{p}++(H_{s}-h_{s})^{2}(H_{p}-h_{p})\unicode[STIX]{x1D735}H_{s}\nonumber\\ \displaystyle & & \displaystyle \frac{{\mathcal{M}}}{2}(2h_{s}h_{p}(H_{s}-h_{s})+h_{s}^{2}(H_{p}-h_{p}))({\mathcal{D}}\unicode[STIX]{x1D735}h_{s}+\unicode[STIX]{x1D735}H_{s})\nonumber\\ \displaystyle & & \displaystyle +\left.{\mathcal{M}}(h_{p}(H_{s}-h_{s})^{2}+2h_{s}(H_{s}-h_{s})(H_{p}-h_{p}))\unicode[STIX]{x1D735}H_{s}\vphantom{\frac{{\mathcal{M}}}{2}}\right],\end{eqnarray}$$

and mass conservation equations

(3.4a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}h_{p}}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}_{\boldsymbol{l}\boldsymbol{p}},\quad \frac{\unicode[STIX]{x2202}(H_{p}-h_{p})}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}_{\boldsymbol{p}}\end{eqnarray}$$

for the perturbed lower and upper layers, respectively. We impose the boundary conditions

(3.5a,b ) $$\begin{eqnarray}\boldsymbol{q}_{\boldsymbol{l}\boldsymbol{p}}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{x}}=0,\quad \boldsymbol{q}_{\boldsymbol{p}}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{x}}=0\quad (x=0),\end{eqnarray}$$

(3.6a,b ) $$\begin{eqnarray}h_{p}=0,\quad H_{p}=0\quad (x=1).\end{eqnarray}$$

That is, it is required that the source fluxes of both fluids remain fixed and that both fluids drain at the nose.

We search for normal-mode solutions proportional to $\text{e}^{\unicode[STIX]{x1D70E}t+\text{i}ky}$ with amplitudes $\tilde{h}_{p}(x)$ , $\tilde{H}_{p}(x)$ . Here, $k$ is the transverse wavenumber of the disturbances and $\unicode[STIX]{x1D70E}$ is their (complex) growth rate. After lengthy algebra, introducing the notation $\unicode[STIX]{x1D709}=1-x$ and dropping tildes, the mass conservation equations (3.4a,b ) can be written compactly in matrix form:

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D70E}\unicode[STIX]{x1D63E}\boldsymbol{v}=\left(\unicode[STIX]{x1D709}^{3/4}\unicode[STIX]{x1D640}\boldsymbol{v}^{\prime }-{\textstyle \frac{1}{4}}\unicode[STIX]{x1D709}^{-1/4}\unicode[STIX]{x1D641}\boldsymbol{v}\right)^{\prime }-k^{2}\unicode[STIX]{x1D709}^{3/4}\unicode[STIX]{x1D640}\boldsymbol{v},\end{eqnarray}$$

and the boundary conditions reduce to

(3.8a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D709}^{3/4}\unicode[STIX]{x1D640}\boldsymbol{v}^{\prime }-{\textstyle \frac{1}{4}}\unicode[STIX]{x1D709}^{-1/4}\unicode[STIX]{x1D641}\boldsymbol{v}=0\quad (x=0),\quad \boldsymbol{v}=0\quad (x=1),\end{eqnarray}$$

where the prime denotes differentiation with respect to $x$ . Here, $\boldsymbol{v}=[h_{p},H_{p}]^{\text{T}}$ , and the entries of the matrices $\unicode[STIX]{x1D63E}$ , $\unicode[STIX]{x1D640}$ and $\unicode[STIX]{x1D641}$ are given in (A 3)–(A 5) in appendix A. It will be convenient to define $\unicode[STIX]{x1D645}=\unicode[STIX]{x1D641}+\unicode[STIX]{x1D640}$ and to consider the transformation

(3.9a-c ) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D63E}}=\unicode[STIX]{x1D63E}\boldsymbol{T},\quad \tilde{\unicode[STIX]{x1D640}}=\unicode[STIX]{x1D640}\boldsymbol{T},\quad \tilde{\unicode[STIX]{x1D645}}=\unicode[STIX]{x1D645}\boldsymbol{T},\end{eqnarray}$$

where

(3.10) $$\begin{eqnarray}\boldsymbol{T}=\left(\begin{array}{@{}cc@{}}1 & 0\\ 1 & 1\end{array}\right).\end{eqnarray}$$

This transformation is equivalent to reformulating the problem in terms of layer thicknesses $[h,H-h]^{\text{T}}$ rather than layer heights $[h,H]^{\text{T}}$ . Some key properties of these matrices will prove useful in the later sections, including that the trace and determinant of $\tilde{\unicode[STIX]{x1D640}}$ and $\tilde{\unicode[STIX]{x1D645}}$ are strictly positive and that $\tilde{\unicode[STIX]{x1D63E}}$ is the identity matrix.

We note that (3.7) has a singular point at the front $x=1$ ( $\unicode[STIX]{x1D709}=0$ ), which poses a hindrance in our stability calculations. In order to make progress, we analyse the behaviour of the perturbations near this singular point asymptotically in § 3.3.

3.3 Asymptotic solution near the singular front

We explore the asymptotic behaviour of the perturbations in an inner region of size $\unicode[STIX]{x1D716}\ll 1$ about $\unicode[STIX]{x1D709}=0$ and define an inner variable $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}/\unicode[STIX]{x1D716}$ . In terms of this inner variable, equation (3.7) becomes

(3.11) $$\begin{eqnarray}\unicode[STIX]{x1D70E}\unicode[STIX]{x1D716}^{5/4}\unicode[STIX]{x1D63E}\boldsymbol{v}+k^{2}\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D701}^{3/4}\unicode[STIX]{x1D640}\boldsymbol{v}=\frac{\text{d}}{\text{d}\unicode[STIX]{x1D701}}\left(\unicode[STIX]{x1D701}^{3/4}\unicode[STIX]{x1D640}\frac{\text{d}\boldsymbol{v}}{\text{d}\unicode[STIX]{x1D701}}+\frac{1}{4}\unicode[STIX]{x1D701}^{-1/4}\unicode[STIX]{x1D641}\boldsymbol{v}\right).\end{eqnarray}$$

Expanding the solution as a series $\boldsymbol{v}=\boldsymbol{v}_{0}+\unicode[STIX]{x1D716}\boldsymbol{v}_{1}+\cdots \,$ gives, at leading order in $\unicode[STIX]{x1D716}$ ,

(3.12) $$\begin{eqnarray}\frac{\text{d}}{\text{d}\unicode[STIX]{x1D701}}\left(\unicode[STIX]{x1D701}^{3/4}\unicode[STIX]{x1D640}\frac{\text{d}\boldsymbol{v}_{0}}{\text{d}\unicode[STIX]{x1D701}}+\frac{1}{4}\unicode[STIX]{x1D701}^{-1/4}\unicode[STIX]{x1D641}\boldsymbol{v}_{\boldsymbol{ 0}}\right)=0,\end{eqnarray}$$

which has general solution

(3.13) $$\begin{eqnarray}\boldsymbol{v}_{\mathbf{0}}=\unicode[STIX]{x1D701}^{1/4}\boldsymbol{A}_{\boldsymbol{ 0}}+\unicode[STIX]{x1D701}^{-\boldsymbol{M}}\boldsymbol{B}_{\boldsymbol{ 0}},\end{eqnarray}$$

obtained after direct integration and the use of an integrating factor, where $\unicode[STIX]{x1D648}=(1/4)\unicode[STIX]{x1D640}^{-1}\unicode[STIX]{x1D641}$ and $\boldsymbol{A}_{\mathbf{0}}$ , $\boldsymbol{B}_{\mathbf{0}}\in \mathbb{R}^{2}$ are constant vectors. The matrix $\unicode[STIX]{x1D701}^{-\boldsymbol{M}}=\exp (-\log (\unicode[STIX]{x1D701})\unicode[STIX]{x1D648})$ is defined in terms of its Taylor series. It is possible to show (see (B 3)–(B 4) in appendix B) that $\unicode[STIX]{x1D648}$ is positive definite for all parameter values. This implies that, for non-zero choices of $\boldsymbol{B}_{\mathbf{0}}$ , the second term on the right-hand side of (3.13) tends to infinity as $\unicode[STIX]{x1D701}\rightarrow 0$ , contradicting the zero-thickness boundary condition (3.8b ). This can be seen by changing basis to canonical form. Therefore, $\boldsymbol{B}_{\mathbf{0}}=\mathbf{0}$ , and so the leading-order solution is simply proportional to $\unicode[STIX]{x1D709}^{1/4}$ .

3.4 Stability results

In this section, we show that the system is stable to all perturbations, whatever the choice of the three dimensionless parameters ${\mathcal{M}}$ , ${\mathcal{D}}$ and ${\mathcal{Q}}$ .

After defining $\boldsymbol{w}=\unicode[STIX]{x1D709}^{-1/4}\boldsymbol{v}$ , equation (3.7) becomes

(3.14) $$\begin{eqnarray}\unicode[STIX]{x1D70E}\unicode[STIX]{x1D709}^{1/4}\unicode[STIX]{x1D63E}\boldsymbol{w}+k^{2}\unicode[STIX]{x1D709}\unicode[STIX]{x1D640}\boldsymbol{w}=\left(\unicode[STIX]{x1D709}\unicode[STIX]{x1D640}\boldsymbol{w}^{\prime }-{\textstyle \frac{1}{4}}\unicode[STIX]{x1D645}\boldsymbol{w}\right)^{\prime }.\end{eqnarray}$$

Multiplying on the left by $\boldsymbol{w}^{\dagger }\unicode[STIX]{x1D645}^{\dagger }$ and integrating gives

(3.15) $$\begin{eqnarray}\unicode[STIX]{x1D70E}\int _{0}^{1}\unicode[STIX]{x1D709}^{1/4}\boldsymbol{w}^{\dagger }\boldsymbol{J}^{\dagger }\unicode[STIX]{x1D63E}\boldsymbol{w}\,\text{d}x+k^{2}\int _{0}^{1}\unicode[STIX]{x1D709}\boldsymbol{w}^{\dagger }\unicode[STIX]{x1D645}^{\dagger }\unicode[STIX]{x1D640}\boldsymbol{w}\,\text{d}x=\int _{0}^{1}\boldsymbol{w}^{\dagger }\unicode[STIX]{x1D645}^{\dagger }\left(\unicode[STIX]{x1D709}\unicode[STIX]{x1D640}\boldsymbol{w}^{\prime }-{\displaystyle \frac{1}{4}}\unicode[STIX]{x1D645}\boldsymbol{w}\right)^{\prime }\,\text{d}x.\end{eqnarray}$$

After integrating by parts and noting that $\unicode[STIX]{x1D645}^{\dagger }\boldsymbol{J}$ is symmetric, the right-hand side becomes

(3.16) $$\begin{eqnarray}\boldsymbol{w}^{\dagger }\unicode[STIX]{x1D645}^{\dagger }\left.\left(\unicode[STIX]{x1D709}\unicode[STIX]{x1D640}\boldsymbol{w}-\frac{1}{4}\unicode[STIX]{x1D645}\boldsymbol{w}\right)\right|_{0}^{1}-\int _{0}^{1}\boldsymbol{w}^{\prime \dagger }\unicode[STIX]{x1D645}^{\dagger }\unicode[STIX]{x1D640}\boldsymbol{w}^{\prime }\unicode[STIX]{x1D709}\,\text{d}x+\left.\frac{1}{8}\boldsymbol{w}^{\dagger }\unicode[STIX]{x1D645}^{\dagger }\unicode[STIX]{x1D645}\boldsymbol{w}\right|_{0}^{1}.\end{eqnarray}$$

Noting that $\unicode[STIX]{x1D709}\unicode[STIX]{x1D640}\boldsymbol{w}-(1/4)\unicode[STIX]{x1D645}\boldsymbol{w}$ is proportional to $(\boldsymbol{q}_{\boldsymbol{l}\boldsymbol{p}}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{x}},\boldsymbol{q}_{\boldsymbol{p}}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{x}})^{\text{T}}$ , using the boundary conditions and noting from § 3.3 that $\boldsymbol{w}$ is regular at the nose $x=1$ , we find that

(3.17) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D70E}\int _{0}^{1}\unicode[STIX]{x1D709}^{1/4}\boldsymbol{w}^{\dagger }\boldsymbol{J}^{\dagger }\unicode[STIX]{x1D63E}\boldsymbol{w}\,\text{d}x+k^{2}\int _{0}^{1}\unicode[STIX]{x1D709}\unicode[STIX]{x1D66C}^{\dagger }\unicode[STIX]{x1D645}^{\dagger }\unicode[STIX]{x1D640}\boldsymbol{w}\,\text{d}x+\int _{0}^{1}\boldsymbol{w}^{\prime \dagger }\unicode[STIX]{x1D645}^{\dagger }\unicode[STIX]{x1D640}\boldsymbol{w}^{\prime }\unicode[STIX]{x1D709}\,\text{d}x\nonumber\\ \displaystyle & & \displaystyle \quad =\left.-\frac{1}{8}\boldsymbol{w}^{\dagger }\unicode[STIX]{x1D645}^{\dagger }\unicode[STIX]{x1D645}\boldsymbol{w}\right|_{x=1}-\left.\frac{1}{8}\boldsymbol{w}^{\dagger }\unicode[STIX]{x1D645}^{\dagger }\unicode[STIX]{x1D645}\boldsymbol{w}\right|_{x=0}\leqslant 0.\end{eqnarray}$$

We wish to show that $\unicode[STIX]{x1D645}^{\dagger }\unicode[STIX]{x1D63E}$ and $\unicode[STIX]{x1D645}^{\dagger }\unicode[STIX]{x1D640}$ are both positive definite as this implies that each of the integrals on the left-hand side of (3.17) is strictly positive for non-zero perturbations $\boldsymbol{w}$ . This implies that $\unicode[STIX]{x1D70E}\leqslant 0$ , else, if $\unicode[STIX]{x1D70E}>0$ then the left-hand side of (3.17) is strictly positive while the right hand side is negative – a contradiction. It is sufficient to show that $\tilde{\unicode[STIX]{x1D645}}^{\dagger }\tilde{\unicode[STIX]{x1D63E}}$ and $\tilde{\unicode[STIX]{x1D645}}^{\dagger }\tilde{\unicode[STIX]{x1D640}}$ are positive definite. After some algebra and using the fact that $A>0$ and $A>a$ , it is possible to show that $\tilde{\unicode[STIX]{x1D640}}$ and $\tilde{\unicode[STIX]{x1D645}}$ have all strictly positive entries, hence the trace of both matrices is strictly positive, and that their determinants are also strictly positive. Therefore, $\tilde{\unicode[STIX]{x1D640}}$ and $\tilde{\unicode[STIX]{x1D645}}$ are positive definite. Since $\tilde{\unicode[STIX]{x1D640}}$ and $\tilde{\unicode[STIX]{x1D645}}$ have strictly positive entries, so does the product $\tilde{\unicode[STIX]{x1D645}}^{\dagger }\tilde{\unicode[STIX]{x1D640}}$ , and in particular, the trace of this product is strictly positive. Using this and the fact that $\det (\tilde{\unicode[STIX]{x1D645}}^{\dagger }\tilde{\unicode[STIX]{x1D640}})=\det (\tilde{\unicode[STIX]{x1D645}})\det (\tilde{\unicode[STIX]{x1D640}})>0$ , we conclude that $\tilde{\unicode[STIX]{x1D645}}^{\dagger }\tilde{\unicode[STIX]{x1D640}}$ , and hence $\unicode[STIX]{x1D645}^{\dagger }\unicode[STIX]{x1D640}$ , is positive definite. See appendix C for details, in particular (C 5) and (C 10) and the surrounding text. Note also that $\tilde{\unicode[STIX]{x1D63E}}\equiv \unicode[STIX]{x1D644}$ is the identity matrix, therefore $\tilde{\unicode[STIX]{x1D645}}^{\dagger }\tilde{\unicode[STIX]{x1D63E}}=\tilde{\unicode[STIX]{x1D645}}^{\dagger }$ , and hence $\unicode[STIX]{x1D645}^{\dagger }\unicode[STIX]{x1D63E}$ is positive definite as well. This concludes the argument that $\unicode[STIX]{x1D70E}\leqslant 0$ . This conclusion holds for all values of the dimensionless parameters ${\mathcal{M}}$ , ${\mathcal{D}}$ and ${\mathcal{Q}}$ . That is, it is not possible to find a choice of parameters for which the system is linearly unstable to small perturbations.

3.5 One-layer gravity currents with drainage

Although the framework above was developed for lubricated, two-layer currents, it may also be applied to a single-fluid current by taking ${\mathcal{Q}}=0$ . The shape of the steady-state current remains qualitatively the same, except that $a=0$ and, by (A 2), $A=12^{1/4}$ . There are no free dimensionless parameters in this case. Although the matrices used in the previous sections are no longer well defined, we can instead straightforwardly show that the perturbation $H_{p}$ obeys the perturbation equation (3.7) with $\unicode[STIX]{x1D63E}$ replaced by the scalar $3/A^{3}$ , $\unicode[STIX]{x1D640}$ replaced by $1$ and $\unicode[STIX]{x1D641}$ replaced by $3$ . As these are all strictly positive scalars, by following the analysis of § 3.4 it is immediate that the one-layer gravity current with drainage is linearly stable as well, as expected.

4.1 Mechanism of instability

To illustrate the instability mechanism, consider the case ${\mathcal{M}}>1$ in which a less viscous fluid intrudes into a more viscous fluid in the unlubricated region, depicted in figure 3(b). In this case, surface slopes, and so pressure gradients, are larger in the unlubricated region than in the lubricated region. If the front is perturbed and the less viscous fluid finds itself in the unlubricated region, it will become subject to higher hydrostatic pressure gradients and so will advance forward, causing an instability. The scenario is reversed in the opposite case ${\mathcal{M}}<1$ , depicted in figure 3(a). In this case, the jump in slope is reversed so that surface slopes, and so pressure gradients, are smaller in the unlubricated region than in the lubricated region, and if the front is perturbed so that the more viscous, underlying fluid gets displaced past the position of the basic-state front, then it will become subject to lower hydrostatic pressure gradients and become suppressed. That is, a viscosity contrast brings with it a change in pressure gradient, which triggers the instability.

There is a fundamental similarity between this instability and the Saffman–Taylor instability in that both are caused by discontinuities in a driving pressure gradient. However, in this case, the underlying pressure gradients originate from changes in free-surface slope and are related to hydrostatic pressure rather than dynamic pressure.

We proceed with a formal justification of the physical mechanism in the following sections.

4.2 Perturbations resisted by vertical shear

We use the non-dimensionalisation (2.2), for which we use the length scale

(4.1) $$\begin{eqnarray}{\mathcal{L}}=\left(\frac{gq_{0}^{3}T_{0}^{4}}{\unicode[STIX]{x1D708}}\right)^{1/5},\end{eqnarray}$$

where $T_{0}$ is the time delay between the initiation of flow of the two layers, as in the experiments of Kowal & Worster (Reference Kowal and Worster2015). The shear-dominated theory of § 2 is supplemented by the following mass conservation equations in the no-slip region (Huppert Reference Huppert1982):

(4.2a,b ) $$\begin{eqnarray}\boldsymbol{q}=-\frac{1}{3}H^{3}\unicode[STIX]{x1D735}H,\quad \frac{\unicode[STIX]{x2202}H}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}\quad (x\geqslant x_{L}),\end{eqnarray}$$

the matching condition

(4.3a,b ) $$\begin{eqnarray}[H]_{-}^{+}=0\quad (x=x_{L})\end{eqnarray}$$

at the lubrication front, and the flux condition and kinematic condition

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle [\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{n}_{\boldsymbol{L}}+\boldsymbol{q}_{\boldsymbol{l}}\boldsymbol{\cdot }\boldsymbol{n}_{\boldsymbol{L}}]^{-}=[\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{n}_{\boldsymbol{L}}]^{+}\quad (x=x_{L}), & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{n}_{\boldsymbol{L}}\boldsymbol{\cdot }\dot{\boldsymbol{x}}_{\boldsymbol{L}}=\lim _{x\rightarrow x_{L}}\boldsymbol{n}_{\boldsymbol{L}}\boldsymbol{\cdot }\boldsymbol{q}_{\boldsymbol{l}}/h. & \displaystyle\end{eqnarray}$$

Here, we take ${\mathcal{D}}=0$ .

4.3 Linearisation

We investigate the linear stability of the system by introducing small disturbances of size $\unicode[STIX]{x1D716}\ll 1$ (distinct from $\unicode[STIX]{x1D716}$ used in the previous section) and write $h=h_{0}+\unicode[STIX]{x1D716}h_{1}$ , $H=H_{0}+\unicode[STIX]{x1D716}H_{1}$ , $x_{L}=x_{L0}+\unicode[STIX]{x1D716}x_{L1}$ and so on. These give rise to the following linearised perturbation fluxes:

(4.6) $$\begin{eqnarray}\boldsymbol{q}_{\boldsymbol{l}\mathbf{1}}=-\left[\left(-\frac{{\mathcal{M}}}{2}h_{0}^{2}h_{1}+{\mathcal{M}}h_{0}h_{1}H_{0}+\frac{{\mathcal{M}}}{2}h_{0}^{2}H_{1}\right)\unicode[STIX]{x1D735}H_{0}+\left(-\frac{{\mathcal{M}}}{6}h_{0}^{3}+\frac{{\mathcal{M}}}{2}h_{0}^{2}H_{0}\right)\unicode[STIX]{x1D735}H_{1}\right],\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle \boldsymbol{q}_{\mathbf{1}}+\boldsymbol{q}_{\boldsymbol{l}\mathbf{1}} & = & \displaystyle -\left[\vphantom{\left(\frac{1-{\mathcal{M}}}{3}(H_{0}-h_{0})^{3}+\frac{{\mathcal{M}}}{3}H_{0}^{3}\right)}((1-{\mathcal{M}})(H_{0}-h_{0})^{2}(H_{1}-h_{1})+{\mathcal{M}}H_{0}^{2}H_{1})\unicode[STIX]{x1D735}H_{0}\right.\nonumber\\ \displaystyle & & \displaystyle +\left.\left(\frac{1-{\mathcal{M}}}{3}(H_{0}-h_{0})^{3}+\frac{{\mathcal{M}}}{3}H_{0}^{3}\right)\unicode[STIX]{x1D735}H_{1}\right]\end{eqnarray}$$

in the lubricated region and

(4.8) $$\begin{eqnarray}\boldsymbol{q}_{\mathbf{1}}=-\left[H_{0}^{2}H_{1}\unicode[STIX]{x1D735}H_{0}+{\textstyle \frac{1}{3}}H_{0}^{3}\unicode[STIX]{x1D735}H_{1}\right]\end{eqnarray}$$

in the unlubricated region. First-order mass conservation equations in the lubricated region are

(4.9a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}h_{1}}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}_{\boldsymbol{l}\mathbf{1}},\quad \frac{\unicode[STIX]{x2202}(H_{1}-h_{1})}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}_{\mathbf{1}},\end{eqnarray}$$

and in the unlubricated region

(4.10) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}H_{1}}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}_{\mathbf{1}}.\end{eqnarray}$$

The condition of flux continuity to first order in $\unicode[STIX]{x1D716}$ can be expressed as

(4.11) $$\begin{eqnarray}\left[x_{L1}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(\boldsymbol{q}_{\boldsymbol{l}\mathbf{0}}+\boldsymbol{q}_{\mathbf{0}})+(\boldsymbol{q}_{\boldsymbol{l}\mathbf{1}}+\boldsymbol{q}_{\mathbf{1}})\right]^{-}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{ x}}=\left[x_{L1}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\boldsymbol{q}_{\mathbf{0}}+\boldsymbol{q}_{\mathbf{1}}\right]^{+}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{ x}}\quad (x=x_{L0}).\end{eqnarray}$$

Height continuity to first order gives

(4.12) $$\begin{eqnarray}\left[x_{L1}\frac{\unicode[STIX]{x2202}H_{0}}{\unicode[STIX]{x2202}x}+H_{1}\right]^{-}=\left[x_{L1}\frac{\unicode[STIX]{x2202}H_{0}}{\unicode[STIX]{x2202}x}+H_{1}\right]^{+}\quad (x=x_{L0}),\end{eqnarray}$$

and the kinematic condition to first order gives

(4.13) $$\begin{eqnarray}{\dot{x}}_{L1}=\frac{1}{h_{0}}\left(x_{L1}\frac{\unicode[STIX]{x2202}\boldsymbol{q}_{\boldsymbol{l}\mathbf{0}}}{\unicode[STIX]{x2202}x}+\boldsymbol{q}_{\boldsymbol{l}\mathbf{1}}\right)\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{x}}-\frac{\boldsymbol{q}_{\boldsymbol{l}\mathbf{0}}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{x}}}{h_{0}^{2}}\left(x_{L1}\frac{\unicode[STIX]{x2202}h_{0}}{\unicode[STIX]{x2202}x}+h_{1}\right).\end{eqnarray}$$

We wish to perform a local analysis in an inner region of size $\unicode[STIX]{x1D6FF}\ll x_{L}$ near the lubrication front in which the surface heights of the two layers in the lubricated region as well as that of the unlubricated region are linear to first order in $\unicode[STIX]{x1D6FF}$ . This can be attained by rescaling

(4.14a-c ) $$\begin{eqnarray}x-x_{L0}=\unicode[STIX]{x1D6FF}(\tilde{x}-x_{L0}),\quad y=\unicode[STIX]{x1D6FF}{\tilde{y}},\quad t-1=\unicode[STIX]{x1D6FF}^{2}\tilde{t},\end{eqnarray}$$

assuming $\tilde{x}$ , ${\tilde{y}}$ , $\tilde{t}=\mathit{O}(1)$ as $\unicode[STIX]{x1D6FF}\rightarrow 0$ , and expanding to first order in $\unicode[STIX]{x1D6FF}$ . Such an approach can be justified by identifying $\unicode[STIX]{x1D6FF}$ with $k^{-1}$ . Such an approach also ensures that the linearised quantities are independent of time, which, effectively, imposes a frozen-time approximation. We additionally impose that all perturbations are local to the lubrication front and decay away from it. We present a more detailed analysis involving the full region of the flow and perturbation equations involving the full, space-dependent basic state in the companion paper (Kowal & Worster Reference Kowal and Worster2019).

In what follows, we split the domain of the flow into two asymptotic regions: the above-mentioned inner region near the lubrication front, composed of both the lubricated and no-slip domains, and two outer regions away from the lubrication front, which also consist of both the lubricated and no-slip domains.

Evaluating basic-state quantities at $t=1$ defines the dimensional time scale ${\mathcal{T}}$ . We denote the basic-state surface heights in the inner region to first order in $\unicode[STIX]{x1D6FF}$ by

(4.15a,b ) $$\begin{eqnarray}h_{0}=a+b(x-x_{L}),\quad H_{0}=A+B(x-x_{L})\quad (x\leqslant x_{L0})\end{eqnarray}$$

for the lubricated region and

(4.16) $$\begin{eqnarray}H_{0}=\unicode[STIX]{x1D6FC}+\unicode[STIX]{x1D6FD}(x-x_{L})\quad (x\geqslant x_{L0})\end{eqnarray}$$

for the unlubricated region, as depicted in the schematic diagram of figure 4.

Figure 4.

Schematic of the basic state in the inner region, in which the surface heights are locally linear.

By searching for normal-mode solutions in the frame of the lubrication front and writing

(4.17) $$\begin{eqnarray}(h_{1}(x,y,t),H_{1}(x,y,t),x_{L1}(y,t))=({\hat{h}}_{1}(x),{\hat{H}}_{1}(x),\hat{x}_{L1})\exp (\unicode[STIX]{x1D70E}t+\text{i}ky),\end{eqnarray}$$

where $k$ is the transverse wavenumber of the disturbances and $\unicode[STIX]{x1D70E}$ is their (complex) growth rate, we find that the amplitudes satisfy the following mass conservation equations (after dropping hats):

(4.18) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D70E}h_{1}-Vh_{1}^{\prime }=(a_{l1}H_{0}^{\prime }+a_{l2}H_{1}^{\prime })^{\prime }-a_{l3}k^{2}H_{1}\quad (x\leqslant x_{L0}), & \displaystyle\end{eqnarray}$$

(4.19) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D70E}H_{1}-VH_{1}^{\prime }=(a_{1}H_{0}^{\prime }+a_{2}H_{1}^{\prime })^{\prime }-a_{3}k^{2}H_{1}\quad (x\leqslant x_{L0}) & \displaystyle\end{eqnarray}$$

for the lubricated region and

(4.20) $$\begin{eqnarray}\unicode[STIX]{x1D70E}H_{1}-VH_{1}^{\prime }=\left[H_{0}^{2}H_{1}H_{0}^{\prime }+\frac{1}{3}H_{0}^{3}H_{1}^{\prime }\right]^{\prime }-\frac{k^{2}}{3}H_{0}^{3}H_{1}\quad (x\geqslant x_{L0})\end{eqnarray}$$

for the unlubricated region, where $a_{li},a_{i}$ , for $i=1,\ldots ,3$ , are cubic functions of the basic-state layer thicknesses and the viscosity ratio ${\mathcal{M}}$ , which are given in (S.1)–(S.9) of the supplementary material available at https://doi.org/10.1017/jfm.2019.321. Here, $V={\dot{x}}_{L0}$ is the velocity of the lubrication front, which is related to $A$ , $B$ , $a$ and $b$ through its dependence on the three dimensionless parameters. Of interest to us are the equations (4.18)–(4.20) up to first order in $\unicode[STIX]{x1D6FF}$ (see (D 1)–(D 3) of appendix D for details).

The first-order flux continuity boundary condition reduces to

(4.21) $$\begin{eqnarray}\left\{\bar{c}_{1}x_{L1}+\bar{c}_{2}H_{1}+\bar{c}_{3}h_{1}+\bar{c}_{4}\frac{\unicode[STIX]{x2202}H_{1}}{\unicode[STIX]{x2202}x}\right\}^{-}=\left\{x_{L1}\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D6FD}^{2}+\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D6FD}H_{1}+\frac{1}{3}\unicode[STIX]{x1D6FC}^{3}\frac{\unicode[STIX]{x2202}H_{1}}{\unicode[STIX]{x2202}x}\right\}^{+}\end{eqnarray}$$

for $x=x_{L0}$ , where $\bar{c}_{1},\ldots ,\bar{c}_{4}$ are functions of the viscosity ratio ${\mathcal{M}}$ , the basic-state thicknesses $a,A$ and slopes $b,B$ , and are given in (S.43)–(S.46) of the supplementary material. Similarly, the condition for continuity of height becomes

(4.22) $$\begin{eqnarray}[x_{L1}B+H_{1}]^{-}=[x_{L1}\unicode[STIX]{x1D6FD}+H_{1}]^{+}\end{eqnarray}$$

for $x=x_{L0}$ . The kinematic condition becomes

(4.23) $$\begin{eqnarray}\unicode[STIX]{x1D70E}x_{L1}={\mathcal{M}}\left(c_{8}x_{L1}+c_{9}h_{1}+c_{10}H_{1}+c_{11}\frac{\unicode[STIX]{x2202}H_{1}}{\unicode[STIX]{x2202}x}\right)\end{eqnarray}$$

for $x=x_{L0}^{-}$ , where $c_{8},\ldots ,c_{11}$ are functions of the layer thicknesses and their derivatives, and are given in (S.47)–(S.50) of the supplementary material.

The equations (4.18)–(4.20) up to first order in $\unicode[STIX]{x1D6FF}$ along with the boundary conditions (4.21)–(4.23) form a system of linear differential equations with constant coefficients, with solutions of the form

(4.24a ) $$\begin{eqnarray}\displaystyle & (h_{1},H_{1})=(\tilde{h}_{1}^{-},\tilde{H}_{1}^{-})\exp ((x-x_{L0})\unicode[STIX]{x1D706}^{-})\quad (x\leqslant x_{L0}), & \displaystyle\end{eqnarray}$$

(4.24b ) $$\begin{eqnarray}\displaystyle & H_{1}=\tilde{H}_{1}^{+}\exp ((x-x_{L0})\unicode[STIX]{x1D706}^{+})\quad (x\geqslant x_{L0}), & \displaystyle\end{eqnarray}$$

$\tilde{h}_{1}^{-}$

$\tilde{H}_{1}^{-}$

$\tilde{H}_{1}^{+}$

$\unicode[STIX]{x1D70E}$

$x_{L1}$

$\tilde{h}_{1}^{-}$

$\tilde{H}_{1}^{+}$

$\unicode[STIX]{x1D706}^{-}$

$\unicode[STIX]{x1D706}^{+}$

$\tilde{H}_{1}^{-}\equiv 1$

$\unicode[STIX]{x1D706}^{-}>0$

$\unicode[STIX]{x1D706}^{+}<0$

Expressing the mass conservation equations in the lubricated region in matrix form gives

(4.25) $$\begin{eqnarray}(\unicode[STIX]{x1D64B}-(\unicode[STIX]{x1D70E}-V\unicode[STIX]{x1D706}^{-})\unicode[STIX]{x1D644})[\tilde{h}_{1}^{-},\tilde{H}_{1}^{-}]^{\text{T}}=\mathbf{0},\end{eqnarray}$$

where $\unicode[STIX]{x1D64B}$ is a $2\times 2$ matrix with entries that are a function of $\unicode[STIX]{x1D706}^{-}$ , $k$ , $a,A,b,B$ and ${\mathcal{M}}$ . The entries of $\unicode[STIX]{x1D64B}$ are given explicitly in (S.51)–(S.54) of the supplementary material. Requiring non-trivial solutions gives that $\unicode[STIX]{x1D706}^{-}$ satisfies

(4.26) $$\begin{eqnarray}\det (\unicode[STIX]{x1D64B}-(\unicode[STIX]{x1D70E}-V\unicode[STIX]{x1D706}^{-})\unicode[STIX]{x1D644})=0,\end{eqnarray}$$

which is a cubic equation for $\unicode[STIX]{x1D706}^{-}$ . Recalling that $\tilde{H}_{1}^{-}=1$ , without loss of generality, it follows that $\tilde{h}_{1}^{-}$ is given by the corresponding solution to (4.26), that is,

(4.27) $$\begin{eqnarray}\tilde{h}_{1}^{-}=\frac{(a(b(6B-3a\unicode[STIX]{x1D706}^{-}+6A\unicode[STIX]{x1D706}^{-})+a(6B\unicode[STIX]{x1D706}^{-}+(a-3A)(k-\unicode[STIX]{x1D706}^{-})(k+\unicode[STIX]{x1D706}^{-}))){\mathcal{M}})}{(3B(-2Ab+a^{2}\unicode[STIX]{x1D706}^{-}+2a(b-B-A\unicode[STIX]{x1D706}^{-})){\mathcal{M}}+6\unicode[STIX]{x1D70E}-6V\unicode[STIX]{x1D706}^{-})}.\end{eqnarray}$$

Solving the mass conservation equation within the unlubricated region subject to decay conditions gives $\unicode[STIX]{x1D706}^{+}$ explicitly as

(4.28) $$\begin{eqnarray}\unicode[STIX]{x1D706}^{+}=-\frac{1}{2\unicode[STIX]{x1D6FC}^{3}}(6\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D6FD}+3V+\sqrt{4\unicode[STIX]{x1D6FC}^{3}(3\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}^{2}+\unicode[STIX]{x1D6FC}^{3}k^{2}+3\unicode[STIX]{x1D70E})+36\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D6FD}V+9V^{2}}).\end{eqnarray}$$

The flux continuity condition specifies the dispersion relation for $\unicode[STIX]{x1D70E}$ in the form

(4.29) $$\begin{eqnarray}\unicode[STIX]{x1D701}_{1}x_{L1}+\unicode[STIX]{x1D701}_{2}\tilde{H}_{1}^{-}+\unicode[STIX]{x1D701}_{3}\tilde{h}_{1}^{-}=-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D6FC}^{2}(3\unicode[STIX]{x1D6FD}\tilde{H}_{1}^{+}+\unicode[STIX]{x1D6FC}\tilde{H}_{1}^{+}\unicode[STIX]{x1D706}^{+}+3\unicode[STIX]{x1D6FD}^{2}x_{L1}),\end{eqnarray}$$

where $\unicode[STIX]{x1D701}_{1},\ldots ,\unicode[STIX]{x1D701}_{3}$ are functions of $a,A,b,B,{\mathcal{M}}$ and $\unicode[STIX]{x1D706}^{-}$ and are given in (S.55)–(S.57) of the supplementary material. This is an implicit equation for $\unicode[STIX]{x1D70E}$ , which appears via $\unicode[STIX]{x1D706}^{-}$ and $\unicode[STIX]{x1D706}^{+}$ . The $\unicode[STIX]{x1D706}^{-}$ and $\unicode[STIX]{x1D706}^{+}$ appear linearly in the dispersion relation and are solutions of cubic and quadratic equations, respectively. More than one solution for $\unicode[STIX]{x1D706}^{-}$ and $\unicode[STIX]{x1D706}^{+}$ , and hence for $\unicode[STIX]{x1D70E}$ , are possible and we are interested in the solution for which the real part of $\unicode[STIX]{x1D70E}$ is largest. Parameters $\tilde{H}_{1}^{+}$ and $x_{L1}$ are found by using the height continuity and kinematic conditions to be

(4.30) $$\begin{eqnarray}\tilde{H}_{1}^{+}=x_{L1}(B-\unicode[STIX]{x1D6FD})+\tilde{H}_{1}^{-}\end{eqnarray}$$

and

(4.31) $$\begin{eqnarray}x_{L1}=\frac{1}{\unicode[STIX]{x1D701}_{6}}(\unicode[STIX]{x1D701}_{4}\tilde{h}_{1}^{-}+\unicode[STIX]{x1D701}_{5}\tilde{H}_{1}^{-}),\end{eqnarray}$$

where $\unicode[STIX]{x1D701}_{4},\ldots ,\unicode[STIX]{x1D701}_{6}$ are functions of $a,A,b,B,{\mathcal{M}},\unicode[STIX]{x1D70E}$ and $\unicode[STIX]{x1D706}^{-}$ and are given in (S.58)–(S.60) of the supplementary material.

Figure 5.

Asymptotic (dashed) and numerical (solid) solutions of the equations of § 4 for the growth rates $\unicode[STIX]{x1D70E}$ against wavenumber $k$ for perturbations resisted dominantly by vertical shear. Parameter values used are ${\mathcal{M}}=10$ and ${\mathcal{Q}}=0.1$ .

In this section, we have determined the dispersion relation for $\unicode[STIX]{x1D70E}$ in terms of quantities such as the perturbations to the surface heights and leading edge position. Equations in this section are solved numerically and the solutions for $\unicode[STIX]{x1D70E}$ are shown in figure 5. In order to formulate a stability criterion, we simplify these using large-wavenumber asymptotics in the following section.

4.4 Large- $k$ asymptotics

In what follows, we expand in series of $1/k$ for $k$ large and write $X=X_{0}k+X_{1}+\cdots \,$ , where $X=(\unicode[STIX]{x1D70E},\unicode[STIX]{x1D706}^{\pm },\tilde{H}_{1}^{\pm },\tilde{h}_{1}^{-},x_{L1})$ . Using (4.26) gives the leading-order governing equation in the lubricated region

(4.32) $$\begin{eqnarray}{\textstyle \frac{1}{6}}k^{3}[(\unicode[STIX]{x1D706}_{0}^{-})^{2}-1][\unicode[STIX]{x1D706}_{0}^{-}-\tilde{\unicode[STIX]{x1D706}}_{0}^{-}]=O(k^{2})\end{eqnarray}$$

giving $\unicode[STIX]{x1D706}_{0}^{-}=\pm 1$ or $\unicode[STIX]{x1D706}_{0}^{-}=\tilde{\unicode[STIX]{x1D706}}_{0}^{-}$ , where

(4.33) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D706}}_{0}^{-}=-\frac{2\unicode[STIX]{x1D70E}_{0}\tilde{\unicode[STIX]{x1D706}}_{n0}}{\tilde{\unicode[STIX]{x1D706}}_{d0}}\end{eqnarray}$$

and $\tilde{\unicode[STIX]{x1D706}}_{n0},\tilde{\unicode[STIX]{x1D706}}_{d0}$ are functions of the layer thicknesses, slopes, viscosity ratio ${\mathcal{M}}$ and $V$ . Explicit expressions are given in (S.61)–(S.62) of the supplementary material. It can be deduced from the latter expression that $\tilde{\unicode[STIX]{x1D706}}_{0}^{-}\rightarrow 0$ as ${\mathcal{M}}\rightarrow 0$ . Taking $\unicode[STIX]{x1D706}_{0}^{-}=1$ gives the next-order equation outlined in the supplementary material, which gives the next-order contribution

(4.34) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{1}^{-}=\frac{\unicode[STIX]{x1D706}_{n1}^{-}}{\unicode[STIX]{x1D706}_{d1}^{-}},\end{eqnarray}$$

where $\unicode[STIX]{x1D706}_{n1}^{-},\unicode[STIX]{x1D706}_{d1}^{-}$ are given in (S.64)–(S.65) of the supplementary material. The relation (4.27) between $\tilde{h}_{1}^{-}$ and $\tilde{H}_{1}^{-}$ , given without loss of generality that $\tilde{H}_{0}^{-}=1$ , yields

(4.35) $$\begin{eqnarray}\tilde{h}_{10}^{-}k+\tilde{h}_{11}^{-}+O(1/k)=\frac{\unicode[STIX]{x1D701}_{7}}{\unicode[STIX]{x1D701}_{8}},\end{eqnarray}$$

where $\unicode[STIX]{x1D701}_{7},\unicode[STIX]{x1D701}_{8}$ are given in (S.66)–(S.67) of the supplementary material. Here $\unicode[STIX]{x1D701}_{7},\unicode[STIX]{x1D701}_{8}$ are functions of $k$ , among other quantities. Expanding in series of $1/k$ gives that $\tilde{h}_{10}^{-}=0$ and

(4.36) $$\begin{eqnarray}\tilde{h}_{11}^{-}=\left[\frac{a{\mathcal{M}}(6Ab-2a^{2}\unicode[STIX]{x1D706}_{1}^{-}+a(-3b+6(B+A\unicode[STIX]{x1D706}_{1}^{-})))}{3a(a-2A)B{\mathcal{M}}+6(\unicode[STIX]{x1D70E}_{0}-V)}\right].\end{eqnarray}$$

The governing equation (4.28) for the unlubricated region gives that

(4.37) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}_{0}^{+}k+\unicode[STIX]{x1D706}_{1}^{+}O(1/k) & = & \displaystyle \frac{1}{2}\left[-\frac{6\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6FC}}-\frac{3V}{\unicode[STIX]{x1D6FC}^{3}}-\sqrt{4\unicode[STIX]{x1D6FC}^{6}k^{2}+12\unicode[STIX]{x1D6FC}^{3}\unicode[STIX]{x1D70E}_{0}k+O(1)}\right]\nonumber\\ \displaystyle & = & \displaystyle -k-\frac{3}{2}\left[2\frac{\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6FC}}+\frac{\unicode[STIX]{x1D70E}_{0}+V\unicode[STIX]{x1D706}^{+}}{\unicode[STIX]{x1D6FC}^{3}}\right]+O(1/k),\end{eqnarray}$$

and so $\unicode[STIX]{x1D706}_{0}^{+}=-1$ and $\unicode[STIX]{x1D706}_{1}^{+}=-(3/2)[2(\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D6FC})+((\unicode[STIX]{x1D70E}_{0}+V\unicode[STIX]{x1D706}^{+})/\unicode[STIX]{x1D6FC}^{3})]$ . The kinematic condition gives that

(4.38) $$\begin{eqnarray}\displaystyle & \displaystyle x_{L10}=0, & \displaystyle\end{eqnarray}$$

(4.39) $$\begin{eqnarray}\displaystyle & \displaystyle x_{L11}=-\frac{2(A-a)^{3}+3a(2A-a)(A-a)\{}{6a\unicode[STIX]{x1D70E}_{0}}. & \displaystyle\end{eqnarray}$$

The next-order contribution, $x_{L12}$ , which will be needed in determining the first-order correction to the dispersion relation (note $x_{L11}$ was needed in determining the dispersion relation at leading order), is given in (S.68) of the supplementary material. The height continuity condition gives

(4.40) $$\begin{eqnarray}\displaystyle \tilde{H}^{+} & = & \displaystyle x_{L1}(B-\unicode[STIX]{x1D6FD})+\tilde{H}^{-}\nonumber\\ \displaystyle & = & \displaystyle x_{L11}(B-\unicode[STIX]{x1D6FD})+\tilde{H}_{11}^{-}+x_{L12}(B-\unicode[STIX]{x1D6FD})\frac{1}{k}+O(1/k^{2})\end{eqnarray}$$

and so

(4.41a-c ) $$\begin{eqnarray}\tilde{H}_{10}^{+}=0,\quad \tilde{H}_{11}^{+}=x_{L11}(B-\unicode[STIX]{x1D6FD})+\tilde{H}_{11}^{-},\quad \tilde{H}_{12}^{+}=x_{L12}(B-\unicode[STIX]{x1D6FD}).\end{eqnarray}$$

The leading-order flux condition reduces to

(4.42) $$\begin{eqnarray}\tilde{H}_{11}^{+}=[(1-{\mathcal{M}})(1-a/A)^{3}+{\mathcal{M}}]\frac{\tilde{H}_{11}^{-}}{\unicode[STIX]{x1D706}_{0}^{+}},\end{eqnarray}$$

which, by comparison to (4.41b ), gives

(4.43) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{0}=\unicode[STIX]{x1D6FE}(a,A,{\mathcal{M}})(B-\unicode[STIX]{x1D6FD}),\end{eqnarray}$$

where

(4.44) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}(a,A,{\mathcal{M}})=\frac{A-a}{6a}\cdot \frac{2(A-a)^{2}+3a(2A-a){\mathcal{M}}}{({\mathcal{M}}+1)+(1-a/A)^{3}(1-{\mathcal{M}})}.\end{eqnarray}$$

4.5 Instability threshold

The formula (4.43) shows that to leading order, the growth rate behaves as

(4.45) $$\begin{eqnarray}\unicode[STIX]{x1D70E}\sim -k\unicode[STIX]{x1D6FE}(a,A,{\mathcal{M}})\left[\frac{\unicode[STIX]{x2202}H_{0}}{\unicode[STIX]{x2202}x}\right]_{-}^{+},\end{eqnarray}$$

at $x=x_{L0}$ , as $k\rightarrow \infty$ . That is, the growth rate varies linearly with $k$ , as $k\rightarrow \infty$ , which signifies that there is no stabilising mechanism at large wavenumbers under the physics considered so far, similarly to the original result of Saffman & Taylor (Reference Saffman and Taylor1958) for the onset of viscous fingering in a porous medium without the stabilising influence of surface tension. We explore possible stabilising mechanisms in §§ 5 and 6 of this paper and in the companion paper. Since $\unicode[STIX]{x1D6FE}(a,A,{\mathcal{M}})$ in this expression is a positive constant, one can deduce that there are positive growth rates at large wavenumbers precisely when

(4.46) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}H_{0}}{\unicode[STIX]{x2202}x}^{-}>\frac{\unicode[STIX]{x2202}H_{0}}{\unicode[STIX]{x2202}x}^{+}\quad (x=x_{L0}),\end{eqnarray}$$

that is, when the upper surface slope is gentler in the lubricated region than in the unlubricated region: $|\unicode[STIX]{x2202}H_{0}^{-}/\unicode[STIX]{x2202}x|<|\unicode[STIX]{x2202}H_{0}^{+}/\unicode[STIX]{x2202}x|$ , since the surface slopes are negative.

To understand the conditions under which precisely this break in slope occurs, the condition of continuity of flux for the basic state $(\boldsymbol{q}_{\mathbf{0}}+\boldsymbol{q}_{\boldsymbol{l}\mathbf{0}})^{-}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{x}}=\boldsymbol{q}_{\mathbf{0}}^{+}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{x}}$ , written explicitly as

(4.47) $$\begin{eqnarray}\left\{-\frac{1}{3}[(1-{\mathcal{M}})(H_{0}-h_{0})^{3}+{\mathcal{M}}H_{0}^{3}]\frac{\unicode[STIX]{x2202}H_{0}}{\unicode[STIX]{x2202}x}\right\}^{-}=\left\{-\frac{1}{3}H_{0}^{3}\frac{\unicode[STIX]{x2202}H_{0}}{\unicode[STIX]{x2202}x}\right\}^{+},\end{eqnarray}$$

may be used to relate the upstream and downstream upper surface slopes as

(4.48) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}=[(1-{\mathcal{M}})(1-a/A)^{3}+{\mathcal{M}}]B,\end{eqnarray}$$

showing that $B$ and $\unicode[STIX]{x1D6FD}$ have the same sign (negative), since the prefactor on the right-hand side is positive, and that the condition $B>\unicode[STIX]{x1D6FD}$ is equivalent to

(4.49) $$\begin{eqnarray}[(1-{\mathcal{M}})(1-a/A)^{3}+{\mathcal{M}}]>1,\end{eqnarray}$$

that is,

(4.50) $$\begin{eqnarray}({\mathcal{M}}-1)[1-(1-a/A)^{3}]>0,\end{eqnarray}$$

or, equivalently,

(4.51) $$\begin{eqnarray}{\mathcal{M}}>1.\end{eqnarray}$$

Therefore, the jump in slope is completely determined by the viscosity ratio between the two fluids, which further determines the stability of the front at large wavenumbers.

4.6 Higher-order contributions

In order to determine an asymptotic solution, such that the difference between it and the full numerical solution decays to zero as $k\rightarrow \infty$ , we examine the next order. This provides visual agreement between asymptotic and numerical results for moderate wavenumbers. To the next order, the flux condition gives

(4.52) $$\begin{eqnarray}\tilde{H}_{12}^{+}=\unicode[STIX]{x1D719}_{1}\tilde{h}_{11}^{-}+\unicode[STIX]{x1D719}_{2}\tilde{H}_{11}^{-}+\unicode[STIX]{x1D719}_{3}\tilde{H}_{11}^{+}+\unicode[STIX]{x1D719}_{4}x_{L11},\end{eqnarray}$$

where $\unicode[STIX]{x1D719}_{1},\ldots ,\unicode[STIX]{x1D719}_{4}$ are given in (S.69)–(S.72) of the supplementary material. Relating this to (4.41c ) gives the next-order correction for the growth rate

(4.53) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{1}=\unicode[STIX]{x1D719}_{5}\tilde{h}_{11}^{-}+\unicode[STIX]{x1D719}_{6}\tilde{H}_{11}^{-}+\unicode[STIX]{x1D719}_{7}\tilde{H}_{11}^{+}+\unicode[STIX]{x1D719}_{8}x_{L11},\end{eqnarray}$$

where the coefficients $\unicode[STIX]{x1D719}_{5},\ldots ,\unicode[STIX]{x1D719}_{8}$ are given in (S.73)–(S.76) of the supplementary material.

The asymptotic solution

(4.54) $$\begin{eqnarray}\unicode[STIX]{x1D70E}\sim \unicode[STIX]{x1D70E}_{0}k+\unicode[STIX]{x1D70E}_{1}\quad (k\rightarrow \infty ),\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{0}$ and $\unicode[STIX]{x1D70E}_{1}$ are derived above, is compared against the numerical solution for $\unicode[STIX]{x1D70E}$ against $k$ in figure 5. The asymptotic solution is shown to converge towards the numerical solution (solid curve in figure 5) as $k$ increases, with good agreement obtained for $k\gtrsim 0.5$ .

It is seen from the asymptotic solution that the growth rates increase unboundedly as the wavenumber increases. Therefore, the limit in which vertical shear provides the dominant resistance to the flow of the perturbations provides no stabilising mechanism for large wavenumbers, indicating that there is another physical mechanism in our experiments that stabilises the flow for large enough wavenumbers. To explore one such possible mechanism, we consider the effect of horizontal shear as the dominant resistance to the flow of the perturbations in the next section. Horizontal shear stresses become important (or dominate) when the wavelength of the perturbations is comparable to (or much smaller than) the thickness of both layers of fluid. The next section concerns the latter limit.

5.1 Model equations

5.1.1 Unlubricated region

To investigate the stability of the flow, we introduce small disturbances of size $\unicode[STIX]{x1D716}\ll 1$ and search for normal-mode solutions with transverse wavenumber $k$ and growth rate $\unicode[STIX]{x1D70E}$ of the form $X=X_{0}+\unicode[STIX]{x1D716}X_{1}+\cdots \,$ , where $X_{1}=\hat{X}_{1}\exp (\text{i}ky+\unicode[STIX]{x1D70E}t)$ and $X=(H(x,y,t),\boldsymbol{u}(x,y,z,t),\boldsymbol{q}(x,y,t))$ . With the assumption that the perturbations to the flow are resisted dominantly by horizontal shear stresses, the gradients in $z$ are negligible and so the perturbation to the velocity obeys

(5.1) $$\begin{eqnarray}\boldsymbol{u}_{\mathbf{1}}=-\frac{1}{k^{2}}\unicode[STIX]{x1D735}H_{1}.\end{eqnarray}$$

In contrast with the previous analysis, here the velocity involves a $1/k^{2}$ contribution, reflecting the resistence to the flow by horizontal shear stresses. We will use the convention that the $y$ -component of the gradient operator means either $\text{d}/\text{d}y$ or $\text{i}k$ , interchangeably.

The volume flux, per unit width, can be expressed as

(5.2) $$\begin{eqnarray}\displaystyle \boldsymbol{q} & = & \displaystyle \int _{0}^{H}\boldsymbol{u}\,\text{d}z=\int _{0}^{H_{0}+\unicode[STIX]{x1D716}H_{1}}(\boldsymbol{u}_{\mathbf{0}}+\unicode[STIX]{x1D716}\boldsymbol{u}_{\mathbf{1}})\,\text{d}z\end{eqnarray}$$

(5.3) $$\begin{eqnarray}\displaystyle & = & \displaystyle \int _{0}^{H_{0}}\boldsymbol{u}_{\mathbf{0}}\,\text{d}z+\unicode[STIX]{x1D716}\left[\left.\int _{0}^{H_{0}}\boldsymbol{u}_{\mathbf{1}}\,\text{d}z+H_{1}\boldsymbol{u}_{\mathbf{0}}\right|_{z=H_{0}}\right]+O(\unicode[STIX]{x1D716}^{2})\end{eqnarray}$$

so that the perturbation to the volume flux, per unit width, is

(5.4) $$\begin{eqnarray}\boldsymbol{q}_{\mathbf{1}}=\left.\int _{0}^{H_{0}}\boldsymbol{u}_{\mathbf{1}}\,\text{d}z+H_{1}\boldsymbol{u}_{\mathbf{0}}\right|_{z=H_{0}}=\int _{0}^{H_{0}}\boldsymbol{u}_{\mathbf{1}}\,\text{d}z-\frac{1}{2}H_{1}H_{0}^{2}\unicode[STIX]{x1D735}H_{0}\end{eqnarray}$$

since the basic-state velocity for a classical current is primarily along the $x$ -direction with $u_{0}=-(1/2)z(2H_{0}-z)\unicode[STIX]{x2202}H_{0}/\unicode[STIX]{x2202}x$ . Substituting the expression (5.1) for the perturbation to the velocity gives

(5.5) $$\begin{eqnarray}\boldsymbol{q}_{\mathbf{1}}=-\frac{1}{k^{2}}H_{0}\unicode[STIX]{x1D735}H_{1}-\frac{1}{2}H_{1}H_{0}^{2}\unicode[STIX]{x1D735}H_{0}.\end{eqnarray}$$

This differs from (4.8) in the $1/k^{2}$ term to reflect that the flow of the perturbations is resisted by horizontal shear. The time evolution of the perturbations, in the frame of the lubrication front, is governed by the mass conservation equation for the next order in $\unicode[STIX]{x1D716}$ :

(5.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}H_{1}-V\frac{\unicode[STIX]{x2202}H_{1}}{\unicode[STIX]{x2202}x} & = & \displaystyle -\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}_{\pmb{{\it 1}}}\nonumber\\ \displaystyle & = & \displaystyle -\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(\boldsymbol{q}_{\mathbf{1}}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{x}})-\text{i}k\boldsymbol{q}_{\mathbf{1}}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{y}}.\end{eqnarray}$$

5.1.2 Lubricated region

We introduce small perturbations of size $\unicode[STIX]{x1D716}\ll 1$ and investigate the stability of the flow by searching for normal-mode solutions with transverse wavenumber $k$ and growth rate $\unicode[STIX]{x1D70E}$ . We expand in series of $\unicode[STIX]{x1D716}$ , $Y=Y_{0}+\unicode[STIX]{x1D716}Y_{1}$ where $Y_{1}={\hat{Y}}_{1}\exp (\text{i}ky+\unicode[STIX]{x1D70E}t),$ and

(5.7) $$\begin{eqnarray}Y=[H(x,y,t),h(x,y,t),\boldsymbol{u}(x,y,z,t),\boldsymbol{u}_{\boldsymbol{l}}(x,y,z,t),\boldsymbol{q}(x,y,t),\boldsymbol{q}_{\boldsymbol{l}}(x,y,t)].\end{eqnarray}$$

As the perturbations to the flow are resisted dominantly by horizontal shear stresses, the gradients in $z$ are negligible and so

(5.8) $$\begin{eqnarray}\boldsymbol{u}_{\mathbf{1}}=-\frac{1}{k^{2}}\unicode[STIX]{x1D735}H_{1}\end{eqnarray}$$

for the upper layer and

(5.9) $$\begin{eqnarray}\boldsymbol{u}_{\boldsymbol{l}\mathbf{1}}=-\frac{{\mathcal{M}}}{k^{2}}\unicode[STIX]{x1D735}H_{1}\end{eqnarray}$$

for the lower layer. There is a small boundary layer at the interface between the two fluids and at the base, in which vertical shear stresses are important and resolve the no-slip boundary condition as well as the continuity of velocity and shear stress at the interface. However, the effect of this boundary layer is negligible in a depth-averaged description of the flow when the wavenumber of the perturbations is large, specifically when $k\gg 1/h$ . As an alternative to scaling arguments, this can also be deduced from considering the full calculation with $k=O(1)$ so that both vertical and transverse shear contribute to resisting the perturbations to the flow, and considering the limit $kh\rightarrow \infty$ .

The volume flux, per unit width, of the upper layer fluid is

(5.10) $$\begin{eqnarray}\displaystyle \boldsymbol{q} & = & \displaystyle \int _{h}^{H}\boldsymbol{u}\,\text{d}z=\int _{h_{0}+\unicode[STIX]{x1D716}h_{1}}^{H_{0}+\unicode[STIX]{x1D716}H_{1}}(\boldsymbol{u}_{\mathbf{0}}+\unicode[STIX]{x1D716}\boldsymbol{u}_{\mathbf{1}})\,\text{d}z\nonumber\\ \displaystyle & = & \displaystyle \int _{h_{0}}^{H_{0}}\boldsymbol{u}_{\mathbf{0}}\,\text{d}z+\unicode[STIX]{x1D716}\left[\left.\int _{h_{0}}^{H_{0}}\boldsymbol{u}_{\boldsymbol{ 1}}\,\text{d}z+H_{1}\boldsymbol{u}_{\mathbf{0}}\right|_{z=H_{0}}\left.-h_{1}\boldsymbol{u}_{\mathbf{0}}\right|_{z=h_{0}}\right]+O(\unicode[STIX]{x1D716}^{2})\end{eqnarray}$$

and in the lower layer is

(5.11) $$\begin{eqnarray}\displaystyle \boldsymbol{q}_{\boldsymbol{l}} & = & \displaystyle \int _{0}^{h}\boldsymbol{u}_{\boldsymbol{l}}\,\text{d}z=\int _{0}^{h_{0}+\unicode[STIX]{x1D716}h_{1}}(\boldsymbol{u}_{\boldsymbol{l}\mathbf{0}}+\unicode[STIX]{x1D716}\boldsymbol{u}_{\boldsymbol{l}\mathbf{1}})\,\text{d}z\nonumber\\ \displaystyle & = & \displaystyle \int _{0}^{h_{0}}\boldsymbol{u}_{\boldsymbol{l}\mathbf{0}}\,\text{d}z+\unicode[STIX]{x1D716}\left[\left.\int _{0}^{h_{0}}\boldsymbol{u}_{\boldsymbol{l}\mathbf{1}}\,\text{d}z+h_{1}\boldsymbol{u}_{\boldsymbol{l}\mathbf{0}}\right|_{z=h_{0}}\right]+O(\unicode[STIX]{x1D716}^{2}).\end{eqnarray}$$

By using expressions for the velocity fields for the two layers at zeroth order in $\unicode[STIX]{x1D716}$ given Kowal & Worster (Reference Kowal and Worster2015), the perturbation to the volume flux in the upper layer can be found to be

(5.12) $$\begin{eqnarray}\displaystyle \boldsymbol{q}_{\mathbf{1}} & = & \displaystyle -\frac{1}{k^{2}}(H_{0}-h_{0})\unicode[STIX]{x1D735}H_{1}-\frac{1}{2}(H_{1}-h_{1})[\!{\mathcal{M}}h_{0}^{2}\nonumber\\ \displaystyle & & \displaystyle +\,2{\mathcal{M}}h_{0}(H_{0}-h_{0})\!]\unicode[STIX]{x1D735}H_{0}-\frac{1}{2}H_{1}(H_{0}-h_{0})^{2}\unicode[STIX]{x1D735}H_{0}\end{eqnarray}$$

and in the lower layer to be

(5.13) $$\begin{eqnarray}\boldsymbol{q}_{\boldsymbol{l}\mathbf{1}}=-\frac{1}{k^{2}}{\mathcal{M}}H_{0}\unicode[STIX]{x1D735}H_{1}-\frac{1}{2}h_{1}[{\mathcal{M}}h_{0}^{2}+2{\mathcal{M}}h_{0}(H_{0}-h_{0})]\unicode[STIX]{x1D735}H_{0}.\end{eqnarray}$$

These differ from (4.6) and (4.7) in the $1/k^{2}$ term, as before. Conservation of mass at first order in $\unicode[STIX]{x1D716}$ in the two layers is described by

(5.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}(H_{1}-h_{1})-V(H_{1}^{\prime }-h_{1}^{\prime }) & = & \displaystyle -\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}_{\mathbf{1}}\nonumber\\ \displaystyle & = & \displaystyle -\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(\boldsymbol{q}_{\mathbf{1}}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{x}})-\text{i}k\boldsymbol{q}_{\mathbf{1}}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{y}}\end{eqnarray}$$

for the upper layer and

(5.15) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}h_{1}-Vh_{1}^{\prime } & = & \displaystyle -\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}_{\boldsymbol{l}\mathbf{1}}\nonumber\\ \displaystyle & = & \displaystyle -\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(\boldsymbol{q}_{\boldsymbol{l}\mathbf{1}}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{x}})-\text{i}k\boldsymbol{q}_{\boldsymbol{l}\mathbf{1}}\boldsymbol{\cdot }\boldsymbol{e}_{\boldsymbol{y}}\end{eqnarray}$$

for the lower layer. Here, we take the basic state to be quasi-steady with respect to the lubrication front.

5.2 Local stability analysis: results

Figure 6.

Asymptotic (dashed) and numerical (solid) solutions for the growth rates $\unicode[STIX]{x1D70E}$ against wavenumber $k$ for perturbations resisted dominantly by horizontal shear. Parameter values used are ${\mathcal{M}}=10$ and ${\mathcal{Q}}=0.1$ .

We locally approximate the basic-state surface heights linearly near the lubrication front and summarise the resulting set of equations in the lubricated and unlubricated regions along with jump conditions in (E 1)–(E 9) of appendix E. The method used to derive these local equations, valid in the vicinity of the lubrication front, follows that of the previous section on perturbations resisted dominantly by vertical shear. We solve these equations explicitly to obtain the dispersion relation in § E.0.2 and also perform an asymptotic analysis for large wavenumbers $k$ . The asymptotics are outlined in § E.0.3 of appendix E. The main message to take away can be captured concisely in the limit of ${\mathcal{M}}\rightarrow \infty$ , in which the large- $k$ behaviour of $\unicode[STIX]{x1D70E}$ is given by the simple relation

(5.16) $$\begin{eqnarray}\unicode[STIX]{x1D70E}=\frac{\unicode[STIX]{x1D711}_{8}{\mathcal{M}}^{2}+\unicode[STIX]{x1D711}_{9}{\mathcal{M}}}{\unicode[STIX]{x1D711}_{10}}+O(1/k,1/{\mathcal{M}}),\end{eqnarray}$$

where $\unicode[STIX]{x1D711}_{8},\ldots ,\unicode[STIX]{x1D711}_{10}$ are functions of the basic state thicknesses, slopes and lubrication front velocity. Explicit expressions can be found in (S.113)–(S.115) of the supplementary material. The important message is that $\unicode[STIX]{x1D70E}$ tends to a constant at large wavenumbers. This, in fact, is the case not only in the limit ${\mathcal{M}}\rightarrow \infty$ , but holds for any value of ${\mathcal{M}}$ and indicates that horizontal shear significantly stabilises the flow at large wavenumber.

The large- $k$ asymptotic solution for the growth rate $\unicode[STIX]{x1D70E}$ is shown in figure 6, where it is additionally seen that the full numerical solution of the equations (E 1)–(E 9) described in appendix E for the growth rate converges towards the asymptotic result.

This analysis reveals that horizontal shear provides a desirable stabilising mechanism for large wavenumbers. In the following section, we combine the physical ideas of § 4 and the current section to determine the behaviour at intermediate wavenumbers.

6.1 Model equations

6.1.1 Unlubricated region

With the lubrication approximation, the Stokes equations for the fluid flow in the unlubricated region are given by

(6.1) $$\begin{eqnarray}0=-\unicode[STIX]{x1D735}p-\boldsymbol{e}_{\boldsymbol{z}}+\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}z^{2}}\right)\boldsymbol{u}\end{eqnarray}$$

in dimensionless variables. A vertical balance of forces gives that the pressure is given hydrostatically by

(6.2) $$\begin{eqnarray}p=H-z,\end{eqnarray}$$

whereas a horizontal balance gives

(6.3) $$\begin{eqnarray}\unicode[STIX]{x1D735}H=\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}z^{2}}\right)\boldsymbol{u}.\end{eqnarray}$$

To investigate the stability of the flow, we introduce small disturbances of size $\unicode[STIX]{x1D716}\ll 1$ and search for normal-mode solutions with transverse wavenumber $k$ and growth rate $\unicode[STIX]{x1D70E}$ of the form $X=X_{0}+\unicode[STIX]{x1D716}X_{1}+\cdots \,$ , where $X_{1}=\hat{X}_{1}\exp (\text{i}ky+\unicode[STIX]{x1D70E}t)$ and $X=(H(x,y,t),\boldsymbol{u}(x,y,z,t),\boldsymbol{q}(x,y,t))$ . The horizontal part of Stokes equations becomes

(6.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}H_{0}}{\unicode[STIX]{x2202}x}\boldsymbol{e}_{\boldsymbol{x}}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D735}H_{1}=\frac{\unicode[STIX]{x2202}^{2}u_{0}}{\unicode[STIX]{x2202}z^{2}}\boldsymbol{e}_{\boldsymbol{x}}+\unicode[STIX]{x1D716}\left(-k^{2}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}z^{2}}\right)\boldsymbol{u}_{\mathbf{1}},\end{eqnarray}$$

which has the following contribution at $O(\unicode[STIX]{x1D716})$ :

(6.5) $$\begin{eqnarray}\unicode[STIX]{x1D735}H_{1}=\left(-k^{2}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}z^{2}}\right)\boldsymbol{u}_{\mathbf{1}}.\end{eqnarray}$$

For intermediate wavenumbers $k$ , the effect of both horizontal and vertical shear stresses can become important and we retain both in the governing equations. The perturbed velocity is therefore exponential in $z$ and resolves the boundary layers that were previously neglected in the limit in which horizontal shear was dominant. The full details are given in (F 1)–(F 2) of appendix F, where the perturbed velocities and fluxes are given explicitly along with the mass conservation equations. We further apply local approximations valid in the vicinity of the lubrication front and match to the outer regions through decay conditions. This enables us to search for normal-mode solutions. The approach follows that of the previous two sections.

6.2 Lubricated region

Linearising the Stokes equations and looking for normal-mode solutions, while keeping the terms arising from both vertical and horizontal shear stresses, gives the following horizontal part of Stokes equations:

(6.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}H_{0}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D735}H_{1}=\frac{\unicode[STIX]{x2202}^{2}u_{0}}{\unicode[STIX]{x2202}z^{2}}\boldsymbol{e}_{\boldsymbol{x}}+\unicode[STIX]{x1D716}\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}z^{2}}\right)\boldsymbol{u}_{\mathbf{1}}, & \displaystyle\end{eqnarray}$$

(6.7) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{M}}(\unicode[STIX]{x1D735}H_{0}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D735}H_{1})=\frac{\unicode[STIX]{x2202}^{2}u_{l0}}{\unicode[STIX]{x2202}z^{2}}\boldsymbol{e}_{\boldsymbol{x}}+\unicode[STIX]{x1D716}\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}z^{2}}\right)\boldsymbol{u}_{\boldsymbol{l}\mathbf{1}} & \displaystyle\end{eqnarray}$$

for the two layers. This yields the following $O(\unicode[STIX]{x1D716})$ contributions:

(6.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}H_{1}=\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}z^{2}}\right)\boldsymbol{u}_{\mathbf{1}}, & \displaystyle\end{eqnarray}$$

(6.9) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{M}}\unicode[STIX]{x1D735}H_{1}=\left(\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}z^{2}}\right)\boldsymbol{u}_{\boldsymbol{l}\mathbf{1}} & \displaystyle\end{eqnarray}$$

for the two layers. As the perturbations to the flow are resisted by both vertical and horizontal shear stresses, the gradients in $z$ are comparable to the gradients in $y$ . We therefore retain both and by looking for normal-mode solutions obtain the following equations for the horizontal velocity:

(6.10) $$\begin{eqnarray}\unicode[STIX]{x1D735}H_{1}=\left(-k^{2}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}z^{2}}\right)\boldsymbol{u}_{\mathbf{1}}\end{eqnarray}$$

for the upper layer and

(6.11) $$\begin{eqnarray}\unicode[STIX]{x1D735}H_{1}=\left(-k^{2}+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}z^{2}}\right)\boldsymbol{u}_{\boldsymbol{l}\mathbf{1}}\end{eqnarray}$$

for the lower layer. This approach resolves the small boundary layer at the interface between the two fluids and at the base, in which vertical shear stresses are important. The limit $kh\rightarrow \infty$ simplifies the problem to that of the previous section, in which horizontal shear stresses provide the dominant resistance to the flow of the perturbations and in which the effect of this boundary layer is negligible in a depth-averaged description of the flow.

6.3 Local stability analysis: results

Due to the terms proportional to $-k^{2}$ , the horizontal velocities are given in terms of exponentials in $z$ (in a finite domain) and present difficulties in obtaining a concise lubrication model. The reader is referred to appendix F, specifically (F 1)–(F 2) and (F 7)–(F 10), where the details of the derivation of the model as well as the application of the local approximation and matching to the outer regions can be found.

Figure 7 shows the growth rates $\unicode[STIX]{x1D70E}$ as a function of the transverse wavenumber $k$ for representative parameter values, for the intermediate regime, considered in this section, in which the perturbations are resisted by both vertical and horizontal shear. Figure 7 also displays the growth rates for the small-wavenumber limit in which vertical shear provides the dominant resistance to the flow of the perturbations and also the large-wavenumber limit in which horizontal shear provides the dominant resistance to the flow of the perturbations. The intermediate regime, involving both vertical and horizontal shear stresses, is seen to converge to these small- and large-wavenumber limits in figure 7.

Figure 7.

Growth rates $\unicode[STIX]{x1D70E}$ against wavenumber $k$ for perturbations resisted by vertical shear only (dotted), horizontal shear only (dashed) and both vertical and horizontal shear (solid). This figure uses results presented in figure 6 (dashed) for $k\gg 1$ . Parameter values used are ${\mathcal{M}}=10$ and ${\mathcal{Q}}=0.1$ .

Figure 8.

Growth rates $\unicode[STIX]{x1D70E}$ against wavenumber $k$ for perturbations resisted by both vertical and horizontal shear. Parameter values used are ${\mathcal{M}}=10$ , 15, 20, 25, 30 and ${\mathcal{Q}}=0.1$ .

Apart from the large-wavenumber stabilisation analysed in the previous section, horizontal shear also stabilises the flow at small wavenumbers. Figure 8 shows representative growth rates as a function of the wavenumber, for the intermediate regime, involving both vertical and horizontal shear stresses, for varying viscosity ratios ${\mathcal{M}}$ . It is seen that small wavenumbers are the most unstable, and that the band of most unstable wavenumbers shifts towards smaller wavenumbers as the viscosity ratio increases. This does not necessarily imply that $k=0$ is the most unstable wavenumber for these viscosity ratios, as a non-local analysis, valid for small wavenumbers, would be necessary to examine this.