2. Theoretical development

Consider a thin film of viscous fluid (fluid 1) of dynamic viscosity $\mu ^{-}$, density $\rho$ and thickness $h(x,y,t)$ spreading under gravity over a surface $z=0$ and inclined at an angle $\theta$ to the horizontal in the semi-infinite region $-\infty < x\leq x_L(y,t)$. The point $x=x_L(y,t)$ is an intrusion front, where the viscous fluid makes contact with another thin film of viscous fluid (fluid 2) of different dynamic viscosity $\mu ^{+}$, density $\rho$ and thickness $H(x,y,t)$ occupying the region $x_L(y,t)\leq x \leq x_N(y,t)$. As depicted in the schematic of figure 1, the thickness of fluid 2 vanishes at the leading edge and fluid 1 is supplied at constant flux in the far field (as $x\to -\infty$), which leads to a uniform thickness in the far field.

Figure 1. Schematic diagram of the profile (a) and plan view (b) of a thin film of viscous fluid of dynamic viscosity $\mu ^{-}$ intruding into another fluid of (different) dynamic viscosity $\mu ^{+}$ down an inclined plane.

We assume that the length scale of each thin film of viscous fluid is much larger than its thickness, and that the thickness gradients may be of the order of the surface slope. We assume that the effects of inertia and surface tension are negligible and formulate a mathematical model of the flow by applying lubrication theory.

Within the limits of lubrication theory, it is the depth-averaged flux that dictates the behaviour of the flow to leading order. Integrating a parabolic velocity profile for such a flow, which satisfies no-slip at the base, yields a non-zero ($z$-independent) velocity at all points and a non-zero flux. This may lead to a vertical front, such as the intrusion front here. Similar vertical fronts occur, for example, for lubricated viscous gravity currents, in which a thin film of viscous fluid intrudes underneath another thin film of viscous fluid of the same density (Kowal & Worster Reference Kowal and Worster2015, Reference Kowal and Worster2019b) or for thin films of a single viscous fluid down an inclined plane before being smoothed off by surface tension (Huppert Reference Huppert1982a). The former becomes unstable to a similar mechanism to that investigated here. The latter becomes unstable to a different kind of frontal instability of a wavelength chosen by a balance with surface tension.

The study most relevant to this work involves lubricated viscous gravity currents (Kowal & Worster Reference Kowal and Worster2015), which are examined both with and without a density difference. The kinematic propagation is occurring at the same time as the intrusion front becomes perpendicular to the underlying plane, while the intrusion front, as well as each vertical slice either side of the intrusion front, propagates kinematically at a rate determined locally by hydrostatic forces, which are directly associated with the upper-surface slope. That rate is the same on either side of the intrusion front, preventing the intruding fluid from overtaking the ‘pre-existing’ fluid.

A vertical force balance within the viscous gravity current upstream and downstream of the intrusion front yields a hydrostatic pressure within each viscous fluid, given by

(2.1)\begin{gather} p^{-} = \rho g \cos\theta(h-z) \quad \text{for}\ x< x_L, \end{gather}

(2.2)\begin{gather}p^{+} = \rho g \cos\theta (H-z) \quad \text{for}\ x_L< x< x_N, \end{gather}

respectively, from which it is possible to formulate and integrate (twice) the horizontal force balance, leading to the the depth-integrated flux, per unit width, given by

(2.3)\begin{gather} \boldsymbol{q^{-}}=\frac{\rho g \cos\theta}{3\mu^{-}}h^{3}(-\boldsymbol{\nabla} h + \tan\theta \boldsymbol{e}_x) \quad \text{for}\ x< x_L, \end{gather}

(2.4)\begin{gather}\boldsymbol{q^{+}}=\frac{\rho g \cos\theta}{3\mu^{+}}H^{3}(-\boldsymbol{\nabla} H + \tan\theta \boldsymbol{e}_x) \quad \text{for}\ x_L< x< x_N, \end{gather}

where $\boldsymbol {e}_x=(1,0,0)$ and $\boldsymbol {e}_y=(0,1,0)$ are the unit basis vectors in the $x$- and $y$-directions, respectively, and the gradient operator is given by $\boldsymbol {\nabla }=({\partial }/{\partial x},\ \partial /\partial y,\ \partial /\partial z)$. The evolution of the thicknesses of both layers of fluid is determined via the mass conservation equations

(2.5)\begin{equation} \frac{\partial h}{\partial t} ={-} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{q^{-}} = \frac{\rho g\cos\theta}{3\mu^{-}}\left[ \boldsymbol{\nabla}\boldsymbol{\cdot} \left(h^{3}\boldsymbol{\nabla} h\right) - \tan\theta \frac{\partial}{\partial x}\left(h^{3}\right)\right], \end{equation}

for $x< x_L$, and

(2.6)\begin{equation} \frac{\partial H}{\partial t} ={-} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{q^{+}} = \frac{\rho g \cos\theta}{3\mu^{+}}\left[ \boldsymbol{\nabla}\boldsymbol{\cdot} \left(H^{3}\boldsymbol{\nabla} H\right) - \tan\theta \frac{\partial}{\partial x}\left(H^{3}\right)\right], \end{equation}

for $x_L< x< x_N$. The intruding fluid is supplied at constant flux in the far field, so that

(2.7)\begin{equation} \boldsymbol{q^{-}\cdot e_x} \to q_\infty \quad \text{as}\ x\to-\infty. \end{equation}

The normal component of the flux and the thickness of the current are both continuous across the intrusion front, so that

(2.8)\begin{equation} \boldsymbol{q^{-}\cdot n_L} = \boldsymbol{q^{+}\cdot n_L} \quad \text{at}\ x=x_L \end{equation}

and

(2.9)\begin{equation} h=H \quad \text{at}\ x=x_L, \end{equation}

where $\boldsymbol n_L$ is the outward unit vector normal to the intrusion front $x=x_L$. The depth-integrated flux vanishes at the leading edge, giving

(2.10)\begin{equation} \boldsymbol{q^{+}\cdot n_N}=0\quad \text{at}\ x=x_N, \end{equation}

where $\boldsymbol n_N$ is the outward unit vector normal to the leading front $x=x_N$. Equivalently, this yields the vanishing frontal thickness condition

(2.11)\begin{equation} H=0\quad \text{at}\ x=x_N. \end{equation}

A fixed volume $V$ of fluid is released ahead of the intrusion front, so that

(2.12)\begin{equation} \int_{x_L}^{x_N} H \,\mathrm{d}\kern0.7pt x = V \end{equation}

and the leading edge $x=x_N$ evolves kinematically, so that

(2.13)\begin{equation} \dot x_N = \lim_{x\to x_N} \frac{\boldsymbol{q\cdot n_N}}{H}. \end{equation}

Similarly, it may be deduced that

(2.14)\begin{equation} \dot x_L = \lim_{x\to x_L} \frac{\boldsymbol{q\cdot n_L}}{H}. \end{equation}

Although the last equation is redundant, we list it here for completeness.

2.1. Asymptotics near the nose

Upon inspection of the governing equations, it follows that there is a stress singularity at the nose $x=x_N$, as a result of the diverging gradient of the thickness of the current near the nose. Such a singularity occurs also for single-layer flows (see Huppert Reference Huppert1982b), irrespective of any viscosity contrast. Near the nose, we have $|\boldsymbol {\nabla } H|\gg 1$, and so the dominant balance of the governing equation (2.6) becomes ${\partial H}/{\partial t} \sim ({\rho g \cos \theta }/{3\mu ^{+}})[ \boldsymbol {\nabla }\boldsymbol {\cdot } (H^{3}\boldsymbol {\nabla } H)],$ akin to rescaled flow over a horizontal substrate. Looking for solutions of the form $H\sim A(x_N-x)^{p}$ yields

(2.15)\begin{equation} H \sim \left[ \frac{9\mu^{+} \dot x_N}{\rho g\cos\theta} \left(1+\left(\frac{\partial x_N}{\partial y}\right)^{2}\right)^{{-}1} (x_N-x)\right]^{1/3}. \end{equation}

Similar asymptotic solutions have been obtained in the two-dimensional analyses of Huppert (Reference Huppert1982b) and Kowal & Worster (Reference Kowal and Worster2015) and the analysis of Kowal & Worster (Reference Kowal and Worster2019b), for perturbations about an axisymmetric base flow. This balance is independent of the extent of the viscous fluid ahead of the intrusion front, as well as any quantities upstream of the intrusion front.

As discussed by Mathunjwa & Hogg (Reference Mathunjwa and Hogg2006), in the context of the stability of the Barenblatt–Pattle solution, by Kowal & Worster (Reference Kowal and Worster2019a,Reference Kowal and Worsterb) in the context of the stability of lubricated viscous gravity currents, and by others, a singularity of this type poses problems when examining the linear stability of the underlying systems. As can be seen by considering small perturbations to the frontal position and linearising, the perturbed frontal thickness itself is singular. Possible remedies include transforming the dependent variable to ensure that the gradient of the new variable is non-singular (see Mathunjwa & Hogg Reference Mathunjwa and Hogg2006), adopting the method of strained coordinates (see Grundy & McLaughlin Reference Grundy and McLaughlin1982) or pinning down the front (see Kowal & Worster Reference Kowal and Worster2019a,Reference Kowal and Worsterb). We remedy this by adopting the latter approach.

2.3. Dimensionless equations

We note that the far-field condition (2.7) implies that the current reaches towards a steady thickness, $h\to h_\infty$ as $x\to -\infty$, where

(2.16)\begin{equation} h_\infty = \left(\frac{3\mu^{-} q_\infty}{\rho g \sin\theta} \right)^{1/3}. \end{equation}

There is, therefore, a one-to-one mapping between the far-field flux and the far-field thickness. We scale against this thickness, and change variables as follows:

(2.17a–c)\begin{gather} (H,h,z)= h_\infty (H^{*},h^{*},z^{*}),\quad (x,y) = (h_\infty \cot\theta) (\eta+x_N^{*},\zeta),\quad t=\left(\frac{3\mu^{-}\cos\theta}{\rho g h_\infty \sin^{2}\theta}\right)\tau, \end{gather}

(2.18)\begin{gather}(x_N,x_L) = (h_\infty \cot\theta) (x_N^{*},x_L^{*}). \end{gather}

The asterisks $^{*}$ denote dimensionless variables. In this way, the leading front, $x=x_N$, corresponds to $\eta =0$ and the intrusion front, $x=x_L$, corresponds to $\eta =\eta _L\equiv x_L-x_N$.

After dropping asterisks, for convenience, the transformed, dimensionless fluxes reduce to

(2.19)\begin{gather} \boldsymbol{q^{-}}= h^{3}\left(-\frac{\partial h}{\partial \eta} + 1\right) \boldsymbol{e_x} + h^{3}\left( -\frac{\partial h}{\partial \zeta} + \frac{\partial x_N}{\partial \zeta}\frac{\partial h}{\partial \eta} \right) \boldsymbol{e_y} \quad \text{for}\ \eta<\eta_L, \end{gather}

(2.20)\begin{gather}\boldsymbol{q^{+}}= \mathcal{M}^{{-}1}H^{3}\left(-\frac{\partial H}{\partial \eta} + 1\right)\boldsymbol{e_x} + \mathcal{M}^{{-}1}H^{3}\left(-\frac{\partial H}{\partial \zeta} + \frac{\partial x_N}{\partial \zeta}\frac{\partial H}{\partial \eta}\right) \boldsymbol{e_y} \quad \text{for}\ \eta>\eta_L. \end{gather}

The governing mass conservation equations become

(2.21)\begin{gather} \frac{\partial h}{\partial \tau} - \dot x_N \frac{\partial h}{\partial \eta} ={-} \frac{\partial}{\partial \eta} \left(q^{-}_x\right) - \left( \frac{\partial}{\partial \zeta} - \frac{\partial x_N}{\partial \zeta}\frac{\partial}{\partial \eta} \right) \left( q^{-}_y\right), \quad \text{for}\ \eta<\eta_L, \end{gather}

(2.22)\begin{gather}\frac{\partial H}{\partial \tau} - \dot x_N \frac{\partial H}{\partial \eta} ={-} \frac{\partial}{\partial \eta} \left(q^{+}_x\right) - \left( \frac{\partial}{\partial \zeta} - \frac{\partial x_N}{\partial \zeta}\frac{\partial}{\partial \eta} \right) \left( q^{+}_y\right) \quad \text{for}\ \eta>\eta_L, \end{gather}

where the subscripts $_x$ and $_y$ denote the $x$ and $y$ vector components, respectively, and the overdot denotes partial differentiation with respect to $\tau$. Similarly, the boundary conditions become

(2.23)\begin{gather} q^{-}_x \to 1 \quad \text{as}\ \eta\to-\infty, \end{gather}

(2.24)\begin{gather}q^{-}_x - \left( \frac{\partial x_L}{\partial \zeta} \right) q^{-}_y= q^{+}_x - \left( \frac{\partial x_L}{\partial \zeta} \right) q^{+}_y \quad \text{at}\ \eta=\eta_L, \end{gather}

(2.25)\begin{gather}h=H \quad \text{at}\ \eta=\eta_L, \end{gather}

(2.26)\begin{gather}H=0 \quad \text{at}\ \eta=0. \end{gather}

The global mass conservation condition reduces to

(2.27)\begin{equation} \int_{\eta_L}^{0} H \,\mathrm{d}\eta = \mathcal{V} \end{equation}

and the kinematic condition reduces to

(2.28)\begin{equation} \dot x_N = \lim_{\eta\to 0} H^{{-}1}\left(q^{+}_x - \frac{\partial x_N}{\partial \zeta}q^{+}_y\right) \left( 1 + \left( \frac{\partial x_N}{\partial \zeta} \right)^{2} \right)^{{-}1/2}. \end{equation}

It is not necessary to use a kinematic condition for the intrusion front $x_L$, as the global condition (2.27) is sufficient to specify $x_L$.

The asymptotic solution near the leading front is given, in dimensionless form, by

(2.29)\begin{equation} H \sim \left[- 3\mathcal{M} \dot x_N \eta \left(1+\left(\frac{\partial x_N}{\partial \zeta}\right)^{2}\right)^{{-}1}\right]^{1/3} \quad \text{as}\ \eta\to 0^{-}. \end{equation}

The solutions are determined by two dimensionless parameters: the viscosity ratio $\mathcal {M}$ and the fixed volume $\mathcal {V}$ of viscous fluid ahead of the intrusion front. Explicitly, these dimensionless parameters are given by

(2.30a,b)\begin{equation} \mathcal{M}=\frac{\mu^{+}}{\mu^{-}}, \quad \text{and} \quad \mathcal{V}=\frac{V\tan\theta}{h^{2}_\infty}. \end{equation}

2.4. Travelling-wave solution

The basic state, which we wish to perturb to examine its stability, takes the form of a $y$-independent travelling-wave solution $\boldsymbol {\varPsi }_0=(h_0(\eta ),\ H_0(\eta ),\ \eta _{L0},\ x_{N0})$ of unit speed in the frame of the advancing front. Here, the subscript $_0$ refers to quantities corresponding to the basic state. Without loss of generality, the origin is set to coincide with the leading front of the current, so that $\eta =\eta _{N0}=0$ corresponds to $x=x_{N0}$. As the travelling-wave solution is of unit speed (equal to the far-field speed), both fronts evolve as $\dot x_{N0}=\dot x_{L0}=1$, so that $x_{N0}=\tau$ and $x_{L0}=\tau +\eta _{L0}$. That is, the velocity of the travelling-wave solution is determined fully by the upstream, far-field flux.

The depth-integrated fluxes per unit width reduce to

(2.31)\begin{gather} q^{-}_0=h_0^{3} \left(1-h_0'\right) \quad \text{for}\ \eta<\eta_{L0}, \end{gather}

(2.32)\begin{gather}q^{+}_0={\mathcal{M}}^{{-}1}H_0^{3} \left(1-H_0'\right) \quad \text{for}\ \eta_{L0}<\eta<0, \end{gather}

where the prime $^{\prime }$ denotes differentiation with respect to $\eta$. After a further integration in $\eta$ and simplification, the governing equations in the two regions reduce to

(2.33)\begin{gather} h_0^{2} \left(1-h_0'\right)=1 \quad \text{for}\ \eta<\eta_{L0}, \end{gather}

(2.34)\begin{gather}{\mathcal{M}}^{{-}1}H_0^{2} \left(1-H_0'\right)=1 \quad \text{for}\ \eta_{L0}<\eta<0, \end{gather}

subject to the boundary conditions

(2.35)\begin{gather} q^{-}_0 \to 1 \quad \text{as}\ \eta\to-\infty, \end{gather}

(2.36)\begin{gather}h_0=H_0 \quad \text{at}\ \eta=\eta_{L0}, \end{gather}

(2.37)\begin{gather}q^{-}_0=q^{+}_0 \quad \text{at}\ \eta=\eta_{L0}, \end{gather}

(2.38)\begin{gather}H_0=0 \quad \text{at}\ \eta=0 \end{gather}

and the global mass conservation condition

(2.39)\begin{equation} \int_{\eta_{L0}}^{0} H_0 \,\mathrm{d}\eta=\mathcal{V}. \end{equation}

The travelling-wave solution retains the inherent cube-root singularity near the leading front, of the form

(2.40)\begin{equation} H_0 \sim \left(- 3\mathcal{M} \eta \right)^{1/3} \quad \text{as}\ \eta\to 0^{-}. \end{equation}

After a further integration of the governing equations in $\eta$, it can be seen that there is an analytic solution, given implicitly by

(2.41)\begin{gather} h_0+\frac{1}{2} \log \left|\frac{ h_0-1 }{ h_0+1 }\right| = \eta -\eta_a, \quad \text{for}\ \eta<\eta_{L0}, \end{gather}

(2.42)\begin{gather}H_0+\frac{\sqrt{\mathcal{M}}}{2} \log \left|\frac{ H_0-\sqrt{\mathcal{M}} }{ H_0+\sqrt{\mathcal{M}} }\right| = \eta -\eta_b, \quad \text{for}\ \eta_{L0}<\eta<0. \end{gather}

By applying the vanishing thickness condition (2.38), we find $\eta _b$=0. There are no explicit formulas for the constants $\eta _a$ and $\eta _{L0}$. These can be determined implicitly by applying the global mass balance (2.39) and the thickness continuity condition (2.36).

Two limiting cases are worth noting: the limits of $\mathcal {V}=0$ and $\mathcal {V}\to \infty$. Both of these describe the flow of a single layer of fluid of the corresponding viscosity. In the limit of $\mathcal {V}=0$, the governing equations and boundary conditions reduce to

(2.43)$$\begin{gather} h_0^{2} \left(1-h_0'\right)=1 \quad \text{for}\ \eta<0, \end{gather}$$

(2.44)$$\begin{gather}h_0^{3} \left(1-h_0'\right)\to 1\quad \text{as}\ \eta\to-\infty,\quad h_0=0 \quad\text{at}\ \eta=0, \end{gather}$$

with an analytic solution given implicitly by

(2.45)\begin{equation} h_0+\frac{1}{2} \log \left|\frac{ h_0-1 }{ h_0+1 }\right| = \eta, \quad \text{for}\ \eta<0. \end{equation}

In the limit of $\mathcal {V}\to \infty$, the governing equations and boundary conditions reduce to

(2.46)\begin{gather} \mathcal{M}^{{-}1}H_0^{2} \left(1-H_0'\right)=1 \quad \text{for}\ \eta<0, \end{gather}

(2.47)\begin{gather}\mathcal{M}^{{-}1}H_0^{3} \left(1-H_0'\right)\to 1\quad \text{as}\ \eta\to-\infty,\quad H_0=0 \quad\text{at } \eta=0, \end{gather}

with an analytic solution given implicitly by

(2.48)\begin{equation} H_0+\frac{\sqrt{\mathcal{M}}}{2} \log \left|\frac{ H_0-\sqrt{\mathcal{M}} }{ H_0+\sqrt{\mathcal{M}} }\right| = \eta, \quad \text{for}\ \eta<0. \end{equation}

Illustrative solutions for sample parameter values are shown in figure 2, along with the $\mathcal {V}=0$ and $\mathcal {V}\to \infty$ limiting solutions, and the solution for $-\eta _{L0}$ is shown in figure 3. The solutions are enveloped by the limiting solutions at fixed $\mathcal {M}$. The solution for which ${\mathcal {M}=1}$ corresponds to a classical viscous gravity current of uniform viscosity down an inclined plane. In such a scenario, the viscous fluids to the left and to the right of the intrusion front are identical. The scenario in which $\mathcal {M}>1$ corresponds to a less viscous fluid intruding into a more viscous fluid. In this scenario, the upper-surface slope is discontinuous at the intrusion front. The greater the viscosity ratio, the greater the jump in slope. Such a change in slope corresponds to a change in hydrostatic pressure gradient, as can be deduced from the pressure fields (2.1) and (2.2). The mechanism of instability is directly related to such a jump in hydrostatic pressure gradient, as explained in § 4.

Figure 2. (a) Base-flow solution for a range of values of the viscosity ratio $\mathcal {M}$ and $\mathcal {V}=1$. The jump in upper-surface slope (hence, in the pressure gradient) at the intrusion front increases for increasing viscosity ratios. (b) Base-flow solution for a range of values of $\mathcal {V}$ for $\mathcal {M}=5$. The $\mathcal {V}=0$ limit is shown as a dashed curve and the $\mathcal {V}\to \infty$ limit is shown as a dotted curve.

Figure 3. The extent $-\eta _{L0}$ as a function of the viscosity ratio $\mathcal {M}$ for $\mathcal {V}=1$ (a) and as a function of the volume $\mathcal {V}$ for $\mathcal {M}=2$ (b) on logarithmically scaled axes.

It is worth noting that although the slope discontinuity at the intrusion front may become large in the limits of $\mathcal {M}\to \infty$ or $\mathcal {V}\to \infty$, the upstream upper-surface slope is never positive with respect to the horizontal. This can be seen by solving (2.33) for $h_0'$, which gives $h_0' = 1-1/h_0^{2}\leq 1$. That is, whatever the value of the parameters, the upper-surface slope is smaller than 0 at all points (measuring the slope in coordinates aligned with the horizontal and vertical). This means that the slope does not become positive (with respect to the horizontal), so the depth-integrated flow, and hence the intrusion front, does not travel upwards. Alternatively, this conclusion can be reached by considering the upstream flux as given by (2.31), which is always positive.

3. Perturbations

We examine the stability of the problem by perturbing about the travelling-wave solution of § 2.4, and searching for normal modes of the form

(3.1)\begin{equation} \boldsymbol{\varPsi}(\eta,\zeta,\tau) = \boldsymbol{\varPsi}_0(\eta) + \epsilon \boldsymbol{\varPsi}_1(\eta)\exp(\sigma \tau + \textrm{i} k \zeta),\end{equation}

where

(3.2)\begin{gather} \boldsymbol{\varPsi}(\eta,\zeta,\tau) = \left(h(\eta,\zeta,\tau), H(\eta,\zeta,\tau),\eta_{L}(\zeta,\tau),x_{N}(\zeta,\tau)\right), \end{gather}

(3.3)\begin{gather}\boldsymbol{\varPsi}_0(\eta) = \left(h_0(\eta), H_0(\eta),\eta_{L0},x_{N0}\right), \end{gather}

(3.4)\begin{gather}\boldsymbol{\varPsi}_1(\eta) = \left(h_1(\eta), H_1(\eta),\eta_{L1},x_{N1}\right). \end{gather}

Note that the $\tau$ and $\zeta$ dependence of the perturbations are factored into the $\exp (\sigma \tau +\textrm {i}k\zeta )$ term of (3.1). Here, $\epsilon$ is a small parameter denoting the order of magnitude of the perturbed quantities. We wish to expand in $\epsilon$ to obtain the governing equations, and the boundary, integral and kinematic conditions for the perturbations at first order in $\epsilon$. These will form a system of linear differential equations in $\eta$, which admit the normal mode solutions of the form (3.1).

At first order in $\epsilon$, the perturbations to the fluxes reduce to

(3.5)\begin{equation} \boldsymbol{q_1^{-}} = \left(3h_0^{2}h_1-h_0^{3} h_1'-3h_0^{2}h_1h_0'\right) \boldsymbol{e_x} + \textrm{i} k h_0^{3}\left( x_{N1}h_0' -h_1 \right)\boldsymbol{e_y}, \end{equation}

for $\eta <\eta _{L0}$, and

(3.6)\begin{equation} \boldsymbol{q_1^{+}} = \mathcal{M}^{{-}1}\left[ \left(3H_0^{2}H_1-H_0^{3} H_1'-3H_0^{2}H_1H_0' \right) \boldsymbol{e_x} + \textrm{i} k H_0^{3}\left( x_{N1}H_0' -H_1 \right)\boldsymbol{e_y} \right], \end{equation}

for $\eta _{L0}<\eta <0$. The governing equations become

(3.7)\begin{equation} \sigma h_1 - h_1' - \sigma x_{N1} h_0' =\left( h_0^{3}h_1' + 3h_0^{2}h_1h_0'-3h_0^{2}h_1\right)' - k^{2} h_0^{3}\left(h_1 - x_{N1}h_0'\right), \end{equation}

for $\eta <\eta _{L0}$, and

(3.8)\begin{align} & \sigma H_1 - H_1' - \sigma x_{N1} H_0' \nonumber\\ &\quad =\mathcal{M}^{{-}1}\left[\left( H_0^{3} H_1' + 3 H_0^{2} H_1 H_0'-3H_0^{2} H_1\right)' - k^{2} H_0^{3}\left(H_1 - x_{N1}H_0'\right)\right], \end{align}

for $\eta _{L0}<\eta <0$. The boundary conditions for the perturbations reduce to

(3.9)\begin{gather} q_{1x}^{-} \to 0 \quad \text{as}\ \eta\to-\infty, \end{gather}

(3.10)\begin{gather}h_1 +\eta_{L1} h_0' = H_1 +\eta_{L1} H_0' \quad \text{at}\ \eta=\eta_{L0}, \end{gather}

(3.11)\begin{gather}q_{1x}^{-} + \eta_{L1} \left(q_{0x}^{-}\right)' = q_{1x}^{+} + \eta_{L1} \left(q_{0x}^{+}\right)' \quad \text{at}\ \eta=\eta_{L0}, \end{gather}

(3.12)\begin{gather}H_1=0 \quad \text{at}\ \eta=0. \end{gather}

Note that it is exponential growth in $t$, rather than in $\eta$, that determines stability. As such, although the condition that the perturbation flux vanishes as $\eta \to -\infty$ results in eigenmodes that vanish in the far field, there may well be growth near the intrusion front. That is, the boundary condition that the flux vanishes as $\eta \to -\infty$ does not preclude exponential growth in $t$ near the intrusion front. See Kowal & Worster (Reference Kowal and Worster2019a) for a similar phenomenon.

In contrast to viscous fingering in a porous medium, which is most conveniently explained in terms of associated potentials, the fingering instability examined here is described in terms of depth-averaged fluxes, which depend on the deflection of the upper surface. These, in turn, need to be determined as part of the system of nonlinear differential equations (3.5)–(3.12) describing a local mass balance. Owing to the $\mathcal {M}^{-1}$ prefactor associated with the depth-integrated upstream flux (3.6), higher upstream hydrostatic pressure gradients are necessary to drive the upstream flow if $\mathcal {M}>1$ to preserve continuity of flux (3.11). As such, if a parcel of intruding fluid is perturbed ahead of the initially planar intrusion front, it becomes subject to higher hydrostatic pressure gradients, advancing further downstream.

The global mass balance for the perturbations becomes

(3.13)\begin{equation} \left.\int_{\eta_{L0}}^{0} H_1 \,\mathrm{d}\eta = \eta_{L1} H_0\right|_{\eta=\eta_{L0}}, \end{equation}

the kinematic condition for perturbations to the frontal position reduces to

(3.14)\begin{equation} \sigma x_{N1} = \lim_{\eta\to 0} \left( \frac{q_{1x}^{+}}{H_0} - \frac{q_{0x}^{+}H_1}{H_0^{2}}\right), \end{equation}

and the asymptotic solution near the front becomes

(3.15)\begin{equation} H_1 \sim \tfrac{1}{3}\sigma x_{N1} \left(- 3\mathcal{M} \eta \right)^{1/3} \quad \text{as}\ \eta\to 0^{-}, \end{equation}

which does not diverge as $\eta \to 0^{-}$.

3.1. Numerical solutions

The system of equations specified in the previous section forms a complete set of governing equations, and boundary, integral and kinematic conditions to specify the perturbations as the solution to a nonlinear eigenvalue problem, for which the growth rate $\sigma$ acts as the eigenvalue. We note that by the linearity of this system of differential equations, any multiple of a solution $(h_1, H_1, \eta _{L1}, x_{N1})$ is also a solution. We, therefore, set $x_{N1}=1$, without loss of generality. This condition forms a nonlinear equation

(3.16)\begin{equation} f(\sigma,k,\mathcal{M},\mathcal{V})=0, \end{equation}

which implicitly determines the growth rate $\sigma$ as a function of the wavenumber $k$, viscosity ratio $\mathcal {M}$ and volume $\mathcal {V}$.

Numerically, these equations are solved on the interval $(-L,\eta _{L0})$, approximating the semi-infinite interval $(-\infty ,\eta _{L0})$ and corresponding to the intruding fluid, and on the interval $(\eta _{L0},\delta )$, approximating the full interval $(\eta _{L0},0)$ including the singular nose and corresponding to the viscous fluid ahead of the intrusion front. This has been achieved using Mathematica's built-in numerical solver for differential equations. Here, $L$ is taken to be a large, finite and positive number, and $\delta$ is taken to be small positive number.

Introducing the non-zero parameter $\delta$ serves as a way to avoid numerical instabilities near the singularity at the nose. The equations are discretised on the interval $(\eta _{L0},\delta )$, with boundary conditions applied at $\eta =\delta$ using the asymptotic solution (3.15).

Another source of numerical error originates from the far-field behaviour of the general solution for the perturbations in the limit $\eta \to -\infty$. Specifically, although the desired solution satisfying the far-field condition decays as $\eta \to -\infty$, we find that the general solution for the perturbations consists of a linear combination of terms, one of which grows exponentially for growth rates of interest, specifically when $\sigma >-k^{2}$, as

(3.17)\begin{equation} h_1 \sim c \exp\left(\left(1-\sqrt{1+k^{2}+\sigma}\right)\eta\right) \quad \text{as}\ \eta\to-\infty \end{equation}

for some constant $c$. To ensure numerical stability we solve for the transformed variable $g_1$ instead of $h_1$, where

(3.18)\begin{equation} g_1(\eta) \equiv h_1(\eta) \exp\left(-\left(1-\sqrt{1+k^{2}+\sigma}\right)\eta\right) \end{equation}

and apply the far-field boundary condition in terms of $g_1$. In this way, $g_1$ remains finite as $\eta \to -\infty$.

4. Mechanism of instability

In contrast to widely known viscous fingering instabilities in porous media, which are driven by a change in dynamic pressure gradient across the intrusion front, viscous fingering instabilities may also occur in unconfined environments, such as in the free-surface flow considered in this paper. These instabilities are instead driven by a change in hydrostatic pressure gradient across the intrusion front, and are stabilised by buoyancy forces at large wavenumbers. A necessary, though insufficient, condition for instability is that the intruding fluid be less viscous than the fluid ahead of the intrusion front, that is, $\mathcal {M} > 1$.

This can be understood by considering the dimensionless longitudinal pressure gradient

(4.1)\begin{gather} G_0^{-}= (p_0^{-})' = h_0' \quad \text{for}\ \eta<\eta_L, \end{gather}

(4.2)\begin{gather}G_0^{+}= (p_0^{+})' = H_0' \quad \text{for}\ \eta_L<\eta<\eta_N \end{gather}

of the base flow, made dimensionless by scaling the pressure using the scale $\rho g h_\infty \cos \theta$. As a consequence of the conditions of continuity of thickness (2.36) and flux (2.37), the jump in hydrostatic pressure gradient,

(4.3)\begin{equation} [G_0]^{+}_-= \left(\mathcal{M}-1\right) \left(G_0^{-} -1\right), \end{equation}

changes sign when $\mathcal {M}=1$. Noting that $G_0^{-} <1$ (as can be seen by solving (2.33) for $h_0'$), it follows that when the intruding fluid, upstream of the intrusion front, is less viscous than the viscous fluid ahead of the intrusion front ($\mathcal {M}>1$), there is a jump in the hydrostatic pressure gradient across the intrusion front so that $-[G_0]^{+}_->0$. That is, the negative pressure gradient downstream of the intrusion front exceeds its upstream value. As a parcel of less viscous fluid gets perturbed ahead of the intrusion front, it becomes subject to a steeper hydrostatic pressure gradient, which advances it further away from the intrusion front. The jump $-[G_0]^{+}_-$ versus $\mathcal {M}$ is shown in figure 4, showing that $-[G_0]^{+}_->0$ for $\mathcal {M}>1$.

Figure 4. (a) The jump in pressure gradient $-[G_0]^{+}_-$ across the intrusion front versus the viscosity ratio for various values of $\mathcal {V}$. (b) The rate of work done in deforming the upper surface just to the left (dashed) and right (solid) of the intrusion front as a function of $\mathcal {M}$ for $\mathcal {V}=1$.

If, on the other hand, the intruding fluid is more viscous than the viscous fluid ahead of the intrusion front ($\mathcal {M}<1$) and the upstream pressure force exceeds buoyancy forces ($G_0^{-} >1$), then the pressure gradient downstream of the intrusion front is smaller than the upstream pressure gradient. As a parcel of more viscous fluid gets perturbed ahead of the intrusion front, it becomes subject to a lower hydrostatic pressure gradient and does not advance further away from the intrusion front.

It is also illuminating to think of the instability mechanism from an energy perspective. Out of the various contributions to the right-hand side of the equation of conservation of energy, the ones of interest for this flow include the rate of work done in deforming the upper surface, and the local stress working, owing to viscous dissipation. These two contributions are dominant also for other viscosity-induced interfacial instabilities of thin-film flows, as outlined by Boomkamp & Miesen (Reference Boomkamp and Miesen1997). The former, which are given by $\varPhi ^{-}=h_0h_1(\frac {1}{2}h_0h_1'+h_1(h_0'-1))(h_0'-1)$ for the intruding layer and $\varPhi ^{+}=({1}/{\mathcal {M}})H_0H_1(\frac {1}{2}H_0H_1'+H_1(H_0'-1))(H_0'-1)$ ahead of the intrusion front, per unit area, are largest near the intrusion front and its values at the intrusion front increase with the viscosity ratio $\mathcal {M}$, as seen in figure 4. Both of these are proportional to the upper-surface slope minus the slope of the inclined substrate, over which the two fluids propagate. Equivalently, these are proportional to the upstream and downstream pressure forces, in excess of downslope gravitational forces. As already discussed, a greater imbalance of these across the intrusion front leads to a greater instability. This is consistent with numerical results, where greater values of $\mathcal {M}$ lead to greater instabilities.

The mechanism is similar to the instability found in the experiments of Kowal & Worster (Reference Kowal and Worster2015) of lubricated viscous gravity currents, explained in the stability analyses of Kowal & Worster (Reference Kowal and Worster2019a,Reference Kowal and Worsterb) in various limits. In the context of lubricated viscous gravity currents, the instability is also driven by a jump in hydrostatic pressure gradient; however, no additional longitudinal buoyancy forces act as a stabilising mechanism. As pointed out earlier, the mechanism of instability is also similar to that of classical viscous fingering instabilities in porous media, but is driven by a hydrostatic, rather than dynamic, pressure gradient. This is what highlights the main difference between the instability considered here and classical viscous fingering instabilities.