3.1. Intrinsic scales and dimensionless model system

In order to reduce the number of parameters in the problem, we non-dimensionalise the system using the intrinsic time, length and height scales,

(3.4a–c)\begin{equation} \mathcal{T} = \left( \frac{\mu_u^2}{\alpha^5 \rho_u^2 g^2 q_{u0}} \right)^{{1}/{3}}, \quad \mathcal{L} = \left( \frac{q_{u0} \mu_u}{\alpha^4 \rho_u g} \right)^{{1}/{3}}, \quad \mathcal{H} = \alpha \mathcal{L} , \end{equation}

respectively. These represent the scales on which the slope of the substrate, $\alpha = - \textrm {d} b / \textrm {d} x$, becomes important. We define non-dimensional (hatted) variables by

(3.5a–c)\begin{equation} t = \mathcal{T} \hat{t}, \quad x = \mathcal{L} \hat{x}, \quad h_i = \mathcal{H} \hat{h}_i. \end{equation}

On dropping the hats, (2.5a) becomes

(3.6a)\begin{equation} \frac{\partial h_u}{\partial t} + \frac{\partial q_u}{\partial x} = 0 \quad\text{and}\quad \frac{\partial h_l}{\partial t} + \frac{\partial q_l}{\partial x} = 0, \end{equation}

for the upper and lower layers, respectively, and the flux expressions become

(3.6b)\begin{gather} q_u = \left( \frac{h_u^3}{3} + M h_u^2 h_l \right) \left( 1 - \frac{\partial h_u}{\partial x} - \frac{\partial h_l}{\partial x} \right) + \frac{M}{2} h_u h_l^2 \left( R - \frac{\partial h_u}{\partial x} - R \frac{\partial h_l}{\partial x} \right) , \end{gather}

(3.6c)\begin{gather}q_l = \frac{M}{3} h_l^3 \left( R - \frac{\partial h_u}{\partial x} - R \frac{\partial h_l}{\partial x} \right) + \frac{M}{2} h_u h_l^2 \left( 1 - \frac{\partial h_u}{\partial x} - \frac{\partial h_l}{\partial x} \right), \end{gather}

where $M = \mu _u/\mu _l$ and $R = \rho _l/\rho _u$. The boundary conditions (2.5f) become

(3.6d)\begin{equation} h_u(x_u,t) = 0, \quad h_l(x_l,t) = 0, \quad q_u(x_u,t) = 0, \quad q_l(x_l,t) = 0. \end{equation}

The input conditions (3.1a,b) become

(3.6e)\begin{equation} q_{l0}(0,t) = Q, \quad q_{u0}(0,t) = 1. \end{equation}

The integral constraints (3.2a,b) become

(3.6f)\begin{equation} \int_0^{x_l} h_l \,\textrm{d} x = Q t, \quad \int_0^{x_u} h_u \,\textrm{d} x = t. \end{equation}

The dimensionless model above depends on three dimensionless parameters,

(3.6g)\begin{equation} M = \frac{\mu_u}{\mu_l}, \quad R = \frac{\rho_l}{\rho_u}, \quad Q = \frac{q_{l0}}{q_{u0}}, \end{equation}

representing the ratios of viscosity, density and input flux, respectively.

Cases of $M > 1$ represent configurations where the lighter fluid is more viscous than the heavier fluid. In such cases, (3.6) is structurally similar to the coupled kinematic wave equations that Fowler (Reference Fowler1982) applies to describe the two-way interplay between a glacier and its hydrological system; the main difference being that the flux expressions that Fowler (Reference Fowler1982) develop depends on interstitial water flux while ours depend on layer thicknesses. Additionally, our study offers a preliminary framework for understanding more complicated rheologies for the lower layer, such as a granular medium (like till) beneath a power-law viscous glacier.

3.2. Illustration of phenomena

To illustrate the general dynamics, we begin by presenting a series of numerical solutions to the dimensionless model (3.6). For this, we used a finite-volume numerical scheme based on converting the conservative form of the partial differential equations to surface integrals and solving them by evaluating the flux into and out of ‘cells’ (small volumes) at each spatial node (e.g. LeVeque Reference LeVeque1992). The details of the scheme are provided in Appendix A. Note that the global volume constraint and continuity of flux and velocity at the lower layer flow front are automatically satisfied by this numerical method. The finite-volume method is particularly suited to the system (2.5) because it can resolve steep gradients, is automatically mass conservative and remains stable during the development of any shock fronts (such as configurations where $R=1$). Despite these advantages, few studies of thin-film dynamics use finite-volume schemes (one example is Grün, Lenz & Rumpf Reference Grün, Lenz and Rumpf2002) and our application here demonstrates in particular its effectiveness for multi-layer fluid flows.

Two illustrative examples are shown in figure 2. For example A, we set $M=10$ (corresponding to the upper fluid being ten times more viscous than the lower fluid) and for example B, we set $M=0.1$ (corresponding to the lower fluid being ten times more viscous than the upper fluid). In each, $R=2$ and $Q=0.5$. Time slices of each solution are shown at a progression of times, $t=1$ and 10, in figure 2(a,b,d,e). The evolutions of the flow fronts are illustrated in figure 2(c,f).

Figure 2. Illustrative numerical solutions to the full dimensionless model (3.6) for the situation where two fluids are injected continuously. In example A (panels a–c) the viscosity ratio is $M = 10$, corresponding to the upper layer being ten times more viscous than the lower layer (forming a lubricating film). In example B (panels d–f) $M = 0.1$, corresponding to the upper layer being ten times less viscous than the lower layer. In both, the flux ratio is $Q = 0.5$ and the density ratio is $R = 2$. Profiles show the solutions at $t=1$ (a,d), $t=10$ (b,e). Panels (c,f) show the evolution of the two flow fronts, $x_u(t)$ and $x_l(t)$, towards linear growth. The inset in (c) shows the convergence of the numerical solutions to the asymptotic theory $900<t<1000$.

The solution for example A illustrates the development of two distinct flow regions. Upstream, the flow forms a region containing both fluid layers. Beyond a critical position, the lower layer terminates and the flow forms a downstream region comprised purely of the lighter fluid. In this example, the thickness of the downstream single-layer region is larger than that of the two-layer region, with a shock-like transition between the two. The plot of the frontal evolutions in figure 2(c) indicate that the fronts propagate at a constant speed at late times, such that the length of the two-layer region approaches a fixed proportion of the length of the total flow. For example B, the upper layer is instead ten times less viscous than the lower layer. The flow likewise forms two regions comprising a region of lighter fluid extending ahead. A key difference compared with example A is that the thickness of the downstream single-layer region is now thinner than the total thickness of the upstream two-layer region. The front positions in case B again transition towards linear growth (figure 2f), showing that the distance between the fronts continues to grow to late times.

The solutions have illustrated a number of universal features. One is that a two-layer release always forms a two-region structure in which the lighter fluid extends ahead of the heavier fluid. Another is that the flow fronts approach linear growth at long times. The solutions also illustrate that it is possible for the frontal thickness to be either greater than or less than the total thickness of the two-layer region. Our analysis will determine solutions describing the long-term asymptotic flow. To this end, we split the analysis into two components. First, we investigate the independent question of how the thickness of the layers in the two-layer region are controlled. With this analysis in hand, we consider the full system comprising both the two-layer region and the frontal region.

3.3. Two-layer region

Our numerical solutions of (3.6) indicated that the solutions approach steady-state thicknesses over time. To confirm this, in figure 3 we plot the thicknesses of both fluid layers as a function of time at the position $x=10$. The lighter fluid front passes first through this point, followed by the two-layer region. By $t=10$ in figure 3(a) and by $t=50$ in figure 3(b), the layer thicknesses in the two-layer region converge to steady-state values. We compare the time-dependent layer thicknesses from the numerical solutions in figure 2 at $x=10$ to the theoretical thickness values (obtained under a steady-state approximation of (3.6), whose calculation we describe in this section) and find agreement. This confirms the approach of the upstream two-layer region to a steady state, whilst also corroborating the validity of our numerical solver and the validity of the steady-state approximation.

Figure 3. Thicknesses of the individual fluid layers at the fixed position $x=10$, showing the convergence of the numerical solutions in the upstream two-layer region towards steady values. Thicknesses plotted here are from our full numerical simulations shown in figure 2 for the following dimensionless parameters: $R=2$, $Q=0.5$ and (a) $M=10$, (b) $M=0.1$. The numerical lower layer thickness is represented by the thin blue line and the numerical upper layer thickness by the thin red line. The steady-state thicknesses predicted by the theory in (3.7) are overlaid in thick black lines.

Integrating the steady-state forms of (3.6a–c) subject to (3.6e), we obtain

(3.7a)\begin{align} 1 &= \frac{h_u^{3}}{3} + M h_u^{2} h_l + \frac{MR}{2} h_l^2 h_u, \end{align}

(3.7b)\begin{align}Q &= M \left( \frac{R}{3} h_l^3 + \frac{1}{2} h_u h_l^2 \right). \end{align}

Since these equations are independent of $x$, it follows that the thicknesses of the layers, $h_u$ and $h_l$, are spatially uniform, confirming the form indicated by our numerical solutions (figures 2 and 3). It should be noted that, in deriving the algebraic equations above, we did not impose the condition that the layers are uniform, only the weaker assumption that the gradient in thicknesses of the layers tends to zero. Since the equations are purely algebraic on neglect of these gradients, it follows that the long-term asymptotic state of the layers is spatially uniform. To solve the equations above, we could use (3.7b) to eliminate $h_u$ in (3.7a) to obtain a cubic equation for $h_l^3$. Since the order of the equation is odd, it is guaranteed to have at least one real root. In order to conduct a complete parameter sweep over values of $M$ and $Q$ using parameter continuation, we opt to solve for $h_u$ and $h_l$ in (3.7) numerically using a Newton–Raphson iterative solver. In doing this, the closest real, positive solution found previously by the solver is used as an initial guess for the next value. There is only one root for which both $h_u$ and $h_l$ are real and positive (and, therefore, physically relevant). As a general point of contrast, we note that previous studies of multi-layered gravity currents have required considerably more detailed analysis of differential systems with unknown free-parameters (cf. Woods & Mason Reference Woods and Mason2000; Kowal & Worster Reference Kowal and Worster2015; Pegler et al. Reference Pegler, Huppert and Neufeld2016). The reduction of our governing equations (3.6) to a simple coupled algebraic system demonstrates both the analytical tractability of the present problem and, as we will show below, will establish a new basis for regime classification of such flows based on the dynamics of the two-layer region alone.

We begin by considering the effect of the viscosity ratio $M$. In our dimensionless model $M$ represents the dimensionless viscosity of the upper layer, whilst the viscosity of the lower layer is fixed at unity. Figure 4(a) shows the thicknesses of the two layers, $h_l$ and $h_u$, as functions of $M$, with $Q=0.5$ and $R=2$ each held fixed. At small values of $M$ (for which the upper layer has a very small viscosity), both layer thicknesses asymptote to the thicknesses that would apply if the layers were considered in isolation, namely,

(3.8a,b)\begin{equation} h_l \sim \left( \frac{3Q}{MR} \right)^{1/3}, \quad h_u \sim 3^{1/3}, \end{equation}

for $M \to 0$. The approach of the upper layer towards its single-layer thickness (3.8b) can be understood by noting that, in view of its small viscosity, the upper layer is much thinner than the lower layer and, therefore, exerts a negligible stress on it. The upper layer therefore has no effect on the lower layer. Conversely, the small viscosity of the upper layer means that the lower layer is effectively static relative to the speed of the upper layer, and, thus, acts as a rigid surface for the upper layer. In this limit, the lighter fluid therefore flows effectively as a single layer over the top of the more viscous heavier fluid. We refer to this situation as regime C.

Figure 4. Thicknesses of the two layers, $h_u$ and $h_l$, with fixed $R=2$ for the case of varying (a) viscosity ratio $M=\mu _u/\mu _l$ and (b) input flux ratio $Q=q_{l0}/q_{u0}$. Here, $h_i$ denotes the steady-state thicknesses for layer $i$. The fixed parameters are: (a) $Q = 0.5$; (b) $M = 1$. The numerical solutions of (3.7) are represented by the thin blue (lower layer) and red (upper layer) solid lines. The asymptotes are overlaid in thick black dashed lines. Layer thicknesses are equal when $M=1$ in panel (a) and when $Q = 1/2$ in panel (b).

As $M$ is increased, the thicknesses of the two layers approach comparable values for $M=O(1)$. For large $M$, both layers thin together as $M^{-1/3}$ (figure 4a). For the lower layer, this decrease follows the trend established from the small $M$ limit (albeit with a different prefactor) and is consistent with the layer flowing faster as its viscosity is reduced. The concurrent decrease in the thickness of the upper layer as $M^{-1/3}$ instead differs from its small-$M$ trend. The decrease is an effect of the lubrication provided by the lower layer: the reduction in underlying shear stress caused by lubrication increases the flow rate of the upper layer and, to maintain the same flux, results in a thinner layer. For $M \to \infty$, we can neglect all terms in (3.7) except those multiplied by $M$. Solving the resulting simplified system, we obtain the asymptotic predictions

(3.9a,b)\begin{equation} h_l \sim \frac{\psi(QR)}{(MR^2)^{1/3}} , \quad h_u \sim \frac{\phi(QR)}{(M/R)^{1/3}}, \end{equation}

for $M \to \infty$, where the functions

(3.9c)\begin{align} \psi(\xi) &\equiv \left[ \frac{3}{2} \left( 5 \xi - \left[ (3 \xi + 5)^2 - 16 \right]^{1/2} + 3 \right) \right]^{1/3}, \end{align}

(3.9d)\begin{align}\phi(\xi) &\equiv \psi(\xi) \left( \frac{2\xi}{\psi(\xi)^3} - \frac{2}{3} \right). \end{align}

The results of (3.9) are shown as dashed black lines in figure 4(a) and confirm the mutual $M^{-1/3}$ thinning trends. Via the functions of $\phi (QR)$ and $\psi (QR)$, the large-$M$ limit retains a relatively complex dependence of the layer thicknesses on the flux ratio $Q$. The parameter $Q$ can be interpreted as the dimensionless flux inputted into the lower layer, whilst the dimensionless flux inputted into the upper layer is fixed at unity. As shown in figure 5, $\psi$ and $\phi$ are increasing and decreasing functions, respectively. The former is consistent with the anticipated thickening of the lower layer as the flux $Q$ introduced into it is increased. The latter implies that the thickness of the upper layer decreases as the flux inputted into the lower layer is increased. This applies because increasing $Q$ results in more lubrication, reducing the stress on the upper layer and causing it to flow faster. Since the input flux is fixed, the faster flow rate results in a thinner upper layer.

Figure 5. The functions $\psi (\xi )$ and $\phi (\xi )$ defined by (3.9c,d) arising in the large-$M$ asymptotic theory for the thicknesses of two fluid layers on a slope, $h_u$ and $h_l$. Asymptotes are overlaid as black dashed lines.

We note that the $M\to \infty$ limit in (3.9) encompasses limits of small and large input flux ratio $Q$. In the limit $QR \to 0$, (3.9) reduces to

(3.10a,b)\begin{equation} h_l \sim \left( \frac{4 Q^2}{M} \right)^{{1}/{3}}, \quad h_u \sim \frac{1}{(2QM)^{{1}/{3}}}, \end{equation}

confirming that the upper layer thickens as $Q^{2/3}$ and the lower layer thins as $Q^{-1/3}$ as $Q$ is increased, both of which are indicated by the increasing and decreasing functions of $\psi$ and $\phi$, respectively (figure 5). We name this situation regime A. We can likewise reduce (3.9) in the limit $QR \to \infty$ which asymptotes to

(3.11a,b)\begin{equation} h_l \sim \left( \frac{3Q}{MR} \right)^{1/3}, \quad h_u \sim \frac{2}{(9 M Q^2 R)^{1/3}}. \end{equation}

In this limit the lower layer grows in proportion to $Q^{1/3}$ whilst the upper decreases in proportion to $Q^{-2/3}$. As $Q \to \infty$, the upper fluid layer is very thin and has negligible impact on the lower layer. The lower layer therefore approaches the thickness it would have in isolation, which is consistent with (3.11a). The upper layer simply translates with the top of the lower layer at the velocity $u_u = MR h_l^2/2$, where $h_l$ is given by (3.11a). The combination of this expression with the flux constraint, $h_u u_u = 1$, indeed produces (3.11b). We refer to this situation as regime D.

To illustrate the effect of the flux ratio $Q$, we plot $h_u$ and $h_l$ as functions of $Q$ in figure 4(b) for $M=1$ and $R=2$. For $Q \to 0$, (3.7) yields the asymptotic predictions

(3.12a,b)\begin{equation} h_l \sim \left( \frac{2Q}{3^{1/3} M} \right)^{1/2}, \quad h_u \sim 3^{1/3}, \end{equation}

which we have plotted as dashed lines in figure 4(b). The thickness of the upper layer asymptotes to the value that would apply if it were in isolation (3.12b), implying a negligible effect of the lower layer on the upper in this limit. Equation (3.12a) shows that the lower layer assumes a small thickness proportional to $Q^{1/2}$. To understand this, we note that, for $Q \to 0$, the lower layer forms a thin Couette flow with a shear rate controlled by the stress applied at the base of the upper layer. More specifically, the stress continuity condition (2.2b) implies that the shear rate is $M \partial u_u / \partial {z}$, where $\partial u_u / \partial {z} = 3^{1/3}$ is the shear rate of the upper layer near the base, resulting in a leading-order velocity profile of $u_l = 3^{1/3}Mz$ in the lower layer. Integrating this shear profile over the depth of the lower layer and applying the flux constraint $q_l = \int u_l \,\textrm {d} z = Q$, we obtain (3.12a). This situation is referred to as regime B.

As $Q$ is increased to values of order unity, the two layers approach comparable thicknesses. For $Q \to \infty$, we recover the asymptotes (3.11) in regime D which are shown to agree with our numerical predictions in figure 4(b). That there are two routes to deriving (3.11), one involving the reduction of (3.7) and one involving the general expression for the $M \to \infty$ limit, indicates that (3.11) describes a limit that encompasses all values of $M$ for sufficiently large $Q$, and this will be illustrated in § 3.5.

To visualise all the configurations of a two-layer gravity current on a single parameter–regime diagram, we present in figure 6 a map of the $(M,Q)$ parameter space divided into four regions by four curves. These regimes are determined by evaluating when the Newton–Raphson numerical solutions for the layer thicknesses in (3.7) fall within 10 % of the theoretical predictions in table 1. The unshaded regions in the $(M,Q)$ parameter space represent parameters over which the system transitions between these regimes. The asymptotes describing the partitioning curves between these regimes are annotated. Regime A (in red) falls within the curves $M \gg 1/6Q$ and $Q \ll 3/4R$ and represents the region in which the numerical solutions fall within 10 % of the asymptotic prediction for the simultaneous limits $M\to \infty$ and $Q\to 0$ (3.10). In regime A the layer thicknesses in the two-layer region thin as $M^{-1/3}$. Regime B is represented by the purple region $8QR^2/27 \ll M \ll 1/6Q$, describing the ‘rigid-lid regime’ in which $Q\to 0$. In this regime the upper layer thickness asymptotes to the value that would apply if it were in isolation and the lower layer forms a thin Couette flow at the base of the upper layer. Regime C (in green) falls within $8QR^2/27 \ll M \ll 8/27Q^2R$ and describes the region in which the numerical solutions lie within 10 % of the asymptotic prediction for the $M \to 0$ limit (3.8). In this regime both layer thicknesses in the two-layer region asymptote to the value that would apply if the layers were released in isolation. Regime D falls in the remaining blue region $Q \gg 3/4R$ and $M \gg 8/27Q^2R$ where the numerical solutions fall within 10 % of the asymptotic predictions for the $Q\to \infty$ limit (3.11). In this regime the lower layer acts as a ‘conveyor belt’ on top of which the upper layer translates. This regime covers a significant region of the parameter space, encompassing asymptotic limits involving both large and small $M$. The structure produces a thick region of heavier fluid in the two-layer region.

Figure 6. Universal regime diagram illustrating the configurations of a two-layer fluid flow introduced at a constant flux onto an inclined plane for $R=2$. The diagram is constructed by determining where the Newton–Raphson numerical solutions for the layer thicknesses from (3.7) fall within 10 % of the asymptotic thickness predictions in table 1. The parameter space of $(M, Q)$ is partitioned into four regimes. Regime A (red region) describes the region in which layer thicknesses mutually thin as $M^{-1/3}$. Regime B (purple region) describes the limit in which the upper layer thickness asymptotes to the value that would apply if the layer was released in isolation. Regime C (green region) describes the region in which the thicknesses of both fluid layers asymptote to the values that would apply if the layers were released in isolation. Regime D (blue region) describes a region that encompasses a wide range of $M$ and $Q$. The blue line marks the demarcating boundary $h_l=h_u$.

Table 1. Predictions for layer thicknesses in the two-layer region (§ 3.3) in each of the four regimes A, B, C and D that are identified in figure 6.

To visualise the control of the relative thicknesses of the two layers in this parameter–regime diagram, we also plot in blue the curve along which the thicknesses in the two-layer region are equal to each other and by setting $h_u = h_l$ in (3.7), we determine the exact seperatrix between these two configurations to be $M = 2Q/(3+2R-6Q-3QR)$. The diagram shows that for $M \gtrsim 1$, the ratio of thicknesses is controlled almost independently by the flux and density ratios alone. An interesting feature is that if the flux ratio is sufficiently large, $Q > (3+2R)/(6 + 2R)$, the lower layer is certain to be thicker than the upper layer, irrespective of the value of the viscosity ratio $M$. Regimes A and B have a thicker upper layer and regimes C and D have a thicker lower layer. We establish a simple demarcation between these pairs of regimes based on layer thicknesses.

Remarkably, the results above show that the dynamics of inclined two-layer gravity currents can be classified entirely from the thicknesses in the two-layer region (3.7). These regimes are not necessarily associated with finite fluid layers. Notably, we emphasise the equivalence of the $M,Q$ scalings that define the demarcations of the regimes in figure 6 and the regimes in Kowal & Worster (Reference Kowal and Worster2015). Unlike previous work, we have derived full asymptotic conditions, including prefactors and dependence on the density ratio $R$. The second crucial distinctions between the regimes in these two studies arises from the scalings in the evolution equations for layer thicknesses and front positions: the horizontal substrate thicknesses scale as $1/5$-powers in $M$ and $Q$, while the inclined substrate thicknesses scale as $1/3$-powers. That said, their regime I (describing situations in which the upper layer has a higher viscosity and a higher input flux than the lower layer and spreads under its own weight) shares similarities with our regime B, in which the lower layer forms a Couette flow due to the stress applied at the base of the upper layer. The theoretical predictions we obtain for regime B (table 1) are identical to the scalings for thickness in regime I of their study. We note that despite the similarities between their regime I and our regime B, the analytical reductions differ. In view of the fundamental difference between the two problems, we show that our coupled cubic algebraic equations (3.7) likely underlie the regimes of all two-layer flows.

3.4. Initiation of two-layer fluids into an empty domain

Having understood the independent dynamics of the two-layer region, we are in a position to construct solutions to the problem in which two finite layers are released and form independently propagating flow fronts.

At long times, the flow becomes increasingly slender as it extends down-slope. Therefore, the thickness gradients $\partial h_i / \partial {x}$ become progressively smaller relative to $\mathrm {d} b/\mathrm {d} x$ in the governing flux expressions (3.6b,c) as $t$ increases. Thus, at long times, we assume that the along-slope component of gravity (represented by $\textrm {d} b / \textrm {d} x$) provides the dominant contribution. In dimensionless forms, these expressions reduce to

(3.13a)\begin{align} \frac{\partial h_u}{\partial t} + \frac{\partial q_u}{\partial x} &= 0, \quad q_u = \frac{h_u^3}{3} + M h_u^2 h_l + \frac{MR}{2} h_u h_l^2 , \end{align}

(3.13b)\begin{align}\frac{\partial h_l}{\partial t} + \frac{\partial q_l}{\partial x} &= 0, \quad\, q_l = \frac{MR}{3} h_l^3 + \frac{M}{2} h_u h_l^2. \end{align}

Equation (3.13a) describes the flux of the upper layer where each term represents contributions from shear variations of the upper layer itself, from the basal stress acting on the upper layer and from the drive induced by the motion by the lower layer, respectively. Equation (3.13b) represents the flux of the lower layer with contributions from flow induced by the motion of the upper layer and from the shear variations in the lower layer, respectively.

By conducting a scaling analysis of the system, we obtain the scalings of $h/t \sim h^3/x$ from (3.13) and $hx \sim t$ from (3.6e) which are consistent with previous studies of single-layer flows on slopes (e.g. Huppert Reference Huppert1982a; Lister Reference Lister1992). Combining these, we obtain $x \sim t$ and $h \sim 1$. These scalings differ from those arising for constant flux inputs over horizontal substrates, where the propagation rates go as $t^{1/5}$ (Huppert Reference Huppert1982b; Kowal & Worster Reference Kowal and Worster2015).

The scaling analysis above indicates the existence of similarity solutions of the form $h =h(\eta )$, where

(3.14)\begin{equation} \eta = t^{{-}1}x, \end{equation}

with front positions $x_u = \eta _u t$ and $x_l = \eta _l t$, where $\eta _u$ and $\eta _l$ are constants. The similarity scaling implies that the flow becomes increasingly slender at long times, $\partial h / \partial {x} \sim 1/t$, thus confirming that the approximation of neglecting $\partial h / \partial {x}$ is self-consistent. Recasting (3.6) in terms of the similarity coordinate (3.14), the equations reduce simply to $\mathrm {d} q_i/\mathrm {d} \eta = 0$, where $q_i$ are given by (3.13). Thus, the fluxes through both layers are uniform through the length of the similarity solutions. Integration of $\mathrm {d} q_i/\mathrm {d} \eta = 0$ subject to the input conditions (3.6e) yields

(3.15a)\begin{gather} \frac{h_u^{3}}{3} + M h_u^{2} h_l + \frac{MR}{2} h_l^2 h_u = 1, \end{gather}

(3.15b)\begin{gather}M \left( \frac{R}{3} h_l^3 + \frac{1}{2} h_u h_l^2 \right) = Q. \end{gather}

These expressions are equivalent to those arising in our analysis of steady two-layer flows in § 3.3. The equivalence reflects a property of the similarity solutions considered here, namely, that the height profiles that conform to the similarity scalings are steady (inherent in the scaling $h \sim 1$). The time dependence of the flow therefore arises entirely from the evolution of the moving flow fronts, $x_l(t)$ and $x_u(t)$.

Within the upstream region comprising both fluid layers, $0 < \eta < \eta _l$, we apply (3.15) to calculate the thicknesses of the two layers. This forms a component of the similarity solution that is equivalent to the problem explored in § 3.3. To determine the thickness in the downstream region comprising the lighter fluid alone, $\eta _l < \eta < \eta _u$, we set $h_l = 0$ in (3.15a) to obtain

(3.15c)\begin{equation} h_u = 3^{1/3}. \end{equation}

The heights in the two regions can thus be determined from (3.15) before the extents of these regions $\eta _l$ and $\eta _u$ are known.

The integral constraints (3.6f) yield

(3.16a,b)\begin{equation} \int_0^{\eta_l} h_l \, \mathrm{d} \eta = Q, \quad \int_0^{\eta_l} h_u \, \mathrm{d} \eta + \int_{\eta_l}^{\eta_u} 3^{1/3} \, \mathrm{d} \eta = 1, \end{equation}

where we have substituted for the thickness in the single-layer region using (3.15c). Since $h_l$ and $h_u$ are uniform, these expressions reduce to

(3.17a,b)\begin{equation} \eta_l = \frac{Q}{h_l}, \quad \eta_u = \eta_l + 3^{{-}1/3} \left(1 - \frac{Q h_u}{h_l} \right), \end{equation}

giving explicit formulae for the positions of the flow fronts $\eta _l$ and $\eta _u$ in terms of the known thicknesses $h_l$ and $h_u$ given by (3.15). It is evident from the similarity solutions (3.17) that the upper fluid always extends ahead of the lower layer to form a single-layer region in front of a two-layer region, confirming our observation from the numerical solutions in figure 2. This essential structure can be understood by noting that, in the two-layer region, the lighter fluid must always flow at least as fast as the heavier fluid below it (its speed is given by the speed of the top surface of the lower layer plus an additional contribution due to the slope of its upper surface). Therefore, the lighter upper fluid must eventually flow ahead of the heavier lower fluid to form a secondary, independent single-layer region.

The similarity solutions can be constructed systematically as follows. First, we determine the thicknesses in the two-layer region, $h_l$ and $h_u$, by solving (3.15a,b) using the Newton–Raphson solver applied in § 3.3. The length of the upstream region $\eta _l$ can then be evaluated using (3.17a). The volume of the lighter fluid in the two-layer region can be determined as $\eta _l h_u$, leaving a residual volume of $1-\eta _l h_u$ to extend ahead and form the single-layer region. Equating the cross-sectional area of the downstream region with this residual, $3^{1/3} (\eta _u - \eta _l) = 1 - \eta _l h_u$ and rearranging, we obtain (3.17b).

Two illustrative solutions are shown in figure 7. For these, we apply the same parameter settings used for our two numerical examples shown earlier in figure 2. We note that while the similarity solutions predict a discontinuity in thickness at the lower layer flow front, in the full unsimplified model (3.6), the higher-order diffusive contributions intervene to smooth the discontinuities (cf. a single-layer flow on a slope Lister Reference Lister1992). The solutions exhibit differences in the thicknesses and lengths of the two regions. While the thickness of the single-layer region is universal, the total thickness of the two-layer region can either be greater than or less than the thickness of the single-layer region. The predictions for the frontal positions, $x_u = \eta _u t$ and $x_l = \eta _l t$, are shown as dashed lines in figure 2(c,f), where they are shown to successfully match our numerically determined predictions at long times.

Figure 7. Similarity solutions constructed using (3.15) and (3.17) for two examples (a) $M=10$ and (b) $M=0.1$ for fixed $Q=0.5$ and $R=2$.

3.5. The control of the propagation rate

In order to understand the general parametric control of the flow structure predicted by the system of (3.15) and (3.17), we have plotted the positions of the layer fronts, $\eta _u$ and $\eta _l$, as functions of the viscosity ratio $M$ for fixed values of $Q=0.5$ and $R=2$ in figure 8(a). The corresponding length of the single-layer region, $\eta _l - \eta _u$, is shown in figure 9(a). As $M$ is increased, $\eta _u$ and $\eta _l$ both increase, which is consistent with both layers propagating faster as the viscosity of the lower layer is reduced. For small $M$, we substitute the asymptotic results of the two-layer region given by (3.8) into (3.17) to give the frontal positions

(3.18a,b)\begin{equation} \eta_l \sim \left( \frac{MRQ^2}{3} \right)^{1/3}, \quad \eta_u \sim \eta_l + 3^{{-}1/3}, \end{equation}

as $M \to 0$. In this limit, the two-layer region is considerably thicker and shorter than the single-layer region. In effect, the heavier fluid acts as a near-rigid, narrow topographic obstruction for the lighter fluid between $\eta = 0$ and $\eta _l$. Once the lighter fluid flows over this, its extent is then effectively equivalent to that which would apply if it were in isolation, forming a dimensionless length of $3^{-1/3}$.

Figure 8. Numerical solutions of the front positions in (3.17) (thin unbroken lines) overlaid by their asymptotes (thick dashed lines) for fixed $R=2$ over a wide range of the (a) viscosity ratio $M$ and (b) the flux ratio $Q$. The fixed parameters are (a) $Q=0.5$ and (b) $M=1$. In panel (a) the function $\zeta (QR) = 3^{-1/3}[1 - QR \phi (QR)/\psi (QR)]$, as defined by (3.19c).

Figure 9. The length of the frontal (single-layer) region $\eta _u-\eta _l$ overlaid with its asymptotes for fixed $R=2$ over a wide range of (a) the viscosity ratio $M$ (fixed $Q=0.5$) and (b) the flux ratio $Q$ (fixed $M=1$). The plots illustrate the approach of the length towards constants in each limit. In panel (a) the grey rectangle indicates the range of possible values spanned by $\zeta (QR)$. The inset in panel (a) plots the function $\zeta (QR)$ defined by (3.19c), representing the dimensionless length in the limit $M \to \infty$.

For $M \to \infty$, the length of the two-layer region becomes comparable to the length of the full system. Substituting (3.9) into (3.17), we obtain

(3.19a,b)\begin{equation} \eta_l \sim \frac{Q (MR^2)^{1/3}}{\psi(QR)}, \quad \eta_u \sim \eta_l + \zeta(QR), \end{equation}

for $M \to \infty$, where

(3.19c)\begin{equation} \zeta(QR) \equiv 3^{{-}1/3}\left( 1 - QR \frac{\phi(QR)}{\psi(QR)} \right) \end{equation}

is a function of $QR$ alone, representing the length of the single-layer region. The front position of the two-layer region thus grows as $M^{1/3}$, representing the acceleration of the rate of propagation by lubrication. By contrast, the length of the single-layer region approaches a constant $\zeta (QR)$ that is independent of $M$. The value of $\zeta (QR)$ across the entire range of $QR$ is strictly bounded within the narrow range $2 < 3^{1/3} \zeta < 3$. The extent of this band of values is indicated by the grey shaded rectangle in figure 9(a). Since $\eta _u - \eta _l$ must also asymptote to the universal value $3^{-1/3}$ as $M \to 0$, the dimensionless length of the frontal region is an almost universal function of the viscosity ratio $M$ alone.

As shown in figure 9(b), both layers grow as the flux ratio $Q$ is increased. For $Q \to 0$, we substitute (3.12) into (3.17) to give

(3.20a,b)\begin{equation} \eta_l \sim \left( \frac{3^{1/3}}{2} MQ \right)^{1/2}, \quad \eta_u \sim 3^{{-}1/3}. \end{equation}

The two-layer region thus occupies a short region of length $Q^{1/2}$ in front of the source, while the flow is primarily comprised of the lighter fluid. As $Q$ increases, the two-layer region develops and eventually dominates the flow domain. The upper layer front position is asymptotically described only by the propagation of this lower layer frontal region. For $Q \to \infty$, we substitute (3.11) into (3.17) to give

(3.21a,b)\begin{equation} \eta_l \sim \left( \frac{MRQ^2}{3} \right)^{1/3}, \quad \eta_u \sim \eta_l + 3^{{-}4/3}, \end{equation}

showing that the extent of the two-layer region (and, thus, of the entire flow) grows as $Q^{2/3}$. As illustrated in figure 9(b), the length of the single-layer region approaches the universal asymptotic value $\eta _u-\eta _l \sim 3^{-4/3}$, equal to one third of the length predicted by single-layer theory, e.g. (3.20b). Remarkably, for large input flux of the heavier fluid, the lighter fluid will always partition its volume such that $2/3$ of it lies in the two-layer region while the remaining $1/3$ lies in the frontal region. This property applies universally for all $M$ and $R$, assuming large $Q$. As a result, the thickness and length of the single-layer region is almost entirely insensitive to the specific value of all the dimensionless parameters.

We compare the convergence of the front positions of the numerical solutions towards the asymptotic predictions determined here (figure 2c,f). In both cases, the numerical solutions converge to the theoretical prediction. There is exact agreement between the prediction for both layers’ front position in the case of $M=0.1$ and between the prediction for the lower layer front position in the case of $M=10$. The result for $M=10$ shows that the front of the upper layer takes relatively longer to approach the asymptotic prediction of the reduced theory, as compared with the lower layer front (as indicated by the inset of figure 2c). This is likely because of the considerably larger viscosity of the upper layer, which acts to increase the time it takes for the layer to adjust to its corresponding asymptotic prediction. Conversely, when the lighter fluid is much less viscous than the heavier fluid ($M=0.1$), the convergence of the front position of the lighter fluid is relatively faster.

Our description of the four regimes in figure 6 is now enriched with information gained from our finite layer analysis. A summary of the front position predictions in each regime are presented in table 2, building on the summary of partitioning curves and layer thickness predictions in table 1. In regime A the structure involves a region of lubrication connected to a thick region of lighter fluid forming the nose. In this case, the lubricating region acts as a ‘bulldozer’, pushing the relatively larger quantity of more viscous fluid considerably further than it would spread in isolation. (In this instance, a horizontal pressure gradient drives the system forward, transferring the lubricating lower layer velocity to the single-layer region.) Regime B describes a configuration in which the lighter fluid layer occupies a considerable proportion of the length of the two-layer system. In contrast to the configuration in regime A, the structure in regime C produces a thin region of lighter fluid that extends downstream of the heavier fluid front. In regime D the structure produces a thick region of heavier fluid upstream of a thin single-layer region. Single-layer frontal regions that are thicker and thinner than the upstream two-layer region fall within this region.

Table 2. Predictions for layer thicknesses (§ 3.3) and front positions (§ 3.5) in each of the four regimes A, B, C and D identified in figure 6.

4.1. Intrinsic scales and dimensionless model system

We non-dimensionalise the system using the intrinsic time, length and height scales,

(4.3a–c)\begin{equation} \mathcal{T} = \left( \frac{\mu_u^2}{\alpha^5 \rho_u^2 g^2 v_{u0}} \right)^{1/2}, \quad \mathcal{L} = \left( \frac{\alpha \rho_u g \mathcal{T} v_{u0}^2}{\mu_u} \right)^{1/3}, \quad \mathcal{H} = \alpha {\mathcal{L}}. \end{equation}

We define dimensionless (hatted) variables according to

(4.4a–c)\begin{equation} t = \mathcal{T} \hat{t}, \quad x = \mathcal{L} \hat{x}, \quad h_i = \mathcal{H} \hat{h}_i, \end{equation}

The dimensionless governing equations are identical to (3.6a–d). The only difference is that we now impose the volume constraints (4.1) and the flux conditions (4.2), the former taking the dimensionless forms

(4.5a,b)\begin{equation} \int_{0}^{x_u} h_{u} \,\textrm{d} x = 1, \quad \int_{0}^{x_l} h_{l} \,\textrm{d} x = V, \end{equation}

where $V = v_{l0}/v_{u0}$. The dimensionless system above depends on three dimensionless parameters

(4.6a–c)\begin{equation} M = \frac{\mu_u}{\mu_l}, \quad R = \frac{\rho_l}{\rho_u}, \quad V = \frac{v_{l0}}{v_{u0}}, \end{equation}

the ratios of viscosity, density and fluid volumes, respectively.

4.2. Illustration of phenomena

We illustrate the general behaviour using our numerical solver in figure 10. To initialise the layers, we introduce them into an initially empty domain at constant dimensionless fluxes ($q_u = 1$ and $q_l = V$) over the short interval $-1 < t < 0$. The input was terminated at $t=0$, producing fixed volumes for all subsequent times, $t > 0$, in accordance with (4.5). We show two example viscosity ratios given by (A) $M=10$ and (B) $M=0.1$, with $R=2$ and $V=0.5$. Each evolution is shown as a progression of time slices at $t=1$ (figure 10a,d) and $t=10$ (figure 10b,e). In both examples, the thickness of the lighter fluid above the heavier fluid thins progressively, to the extent that the two fluids largely occupy independent layers at long times. This can be understood by noting that the fluid on top of the lower layer always flows faster than the lower layer and, hence, without a resupply of lighter fluid, the upper layer progressively ‘spills’ ahead the heavier fluid. The lighter fluid always advances ahead of the heavier fluid because its speed is at least as fast as the top surface of the heavier fluid below it. In contrast with example A, in case B the lighter fluid front is consistently thinner than the heavier fluid upstream of it. There is a drop, rather than a jump, in thickness at the front of the heavier fluid at large times. A second key difference between the examples is the extent of each fluid. The heavier fluid is longer in case A than it is in case B and the lighter fluid is longer in case B than in case A. This suggests that, despite the separation of the flows, the less viscous layer plays a key role in propagating the frontal region of lighter fluid. In contrast to the constant flux case, the front positions both grow with a sublinear trend, but likewise move progressively further apart in time.

Figure 10. Numerical solutions of the dimensionless model (3.6a–d), (4.5) and (3.6) from an initial injection of a fixed volume release of both layers, with $V = 0.5$, $R = 2$ and (a,b) $M = 10$, (d,e) $M = 0.1$. Both fluids are injected into an initially empty domain at $-1<t<0$. Profiles at $t=1$ (a,d) and $t=10$ (b,e) are illustrated. (c,f) The evolution of the front positions from the numerical output are overlaid with their theoretical predictions.

4.3. Similarity solutions

We note that (3.6a–c) and (4.5) yield the scalings $h/t \sim h^3/x$ and $hx \sim 1$, respectively. Eliminating $h$, we obtain $x \sim t^{1/3}$ and, hence, $h \sim t^{-1/3}$, which are consistent with the scalings arising in the case of a single layer (Huppert Reference Huppert1982a). Thus, we define the relevant similarity variables as

(4.7a,b)\begin{equation} \eta = t^{{-}1/3} x , \quad h = t^{{-}1/3} H(\eta). \end{equation}

The front positions of the two layers therefore evolve as $x_l = \eta _l t^{1/3}$ and $x_u = \eta _u t^{1/3}$, confirming the trends indicated earlier by our numerical solution of figure 10(c,f). Recasting (3.6a–c) in terms of the similarity variables (4.7), we obtain

(4.8a)\begin{align} \frac{1}{3} \frac{\mathrm{d}}{\mathrm{d}\eta} (H_u \eta) &= \frac{\mathrm{d}}{\mathrm{d}\eta} \left( \frac{H_u^3}{3} + M H_u^2 H_l + \frac{MR}{2} H_l^2 H_u \right) \equiv \frac{\textrm{d} q_u}{\textrm{d} \eta}, \end{align}

(4.8b)\begin{align}\frac{1}{3} \frac{\mathrm{d}}{\mathrm{d}\eta} (H_l\eta) &= \frac{\mathrm{d}}{\mathrm{d}\eta} \left( \frac{MR}{3} H_l^3 + \frac{M}{2} H_u H_l^2 \right) \equiv \frac{\textrm{d} q_l}{\textrm{d} \eta}, \end{align}

where $q_u$ and $q_l$ are the flux expressions (3.13) in the limit where the along-slope contribution to the flow dominates.

Integrating (4.8) subject to the flux conditions (4.2), and simplifying, we obtain

(4.9a)\begin{align} \frac{\eta}{3} &= \frac{H_u^2}{3} + M H_u H_l + \frac{MR}{2} H_l^2 , \end{align}

(4.9b)\begin{align}\frac{\eta}{3} &= \frac{MR}{3} H_l^2 + \frac{M}{2} H_u H_l. \end{align}

Since these equations depend on $\eta$, it follows that the thickness profiles of the fluid layers both vary along the flow, differing qualitatively from the asymptotic form of layers produced by a constant flux (§ 3). In principle, these equations could be solved by using (4.9b) to eliminate $H_l$ in (4.9a) and attempting to solve the resulting polynomial equation for $H_u(\eta )$. As noted in § 4.2, however, the lighter fluid spills progressively ahead of the heavier fluid. This indicates that the fluids should partition at long times into two distinct single-layer regions. In other words, we can anticipate that there is no part of the final similarity solution in which both $H_l$ and $H_u$ are positive simultaneously.

In order to prove this mathematically, we proceed by contradiction and assume that both $H_l$ and $H_u$ are positive. Using (4.9b) to eliminate $H_u$ in (4.9a) and applying the quadratic formula on the resulting quadratic, we obtain

(4.10)\begin{equation} \frac{MR H_l^2}{\eta} = 1 \pm \left(1- \frac{8R}{9M} \right)^{{-}1/2}. \end{equation}

Since the second term on the right-hand side is greater than unity, we must take the positive root in order for $H_l$ to be real (we also require $8R < 9M$; if this inequality is not satisfied then $H_l$ has a complex solution which is unphysical). Using the resulting expression to eliminate $H_l^2$ in (4.9b) and simplifying, we obtain

(4.11)\begin{equation} -\frac{M}{2} H_u H_l = \frac{\eta}{3} \left(1- \frac{8R}{9M} \right)^{{-}1/2}. \end{equation}

Since $H_u>0$ and the right-hand side is positive, it follows that $H_l$ is negative, leading to a contradiction of the requirement that the layer thicknesses are positive. Therefore, the initial assumption that $H_l$ and $H_u$ can be simultaneously positive is false. Consequently, there is no location along the similarity solution in which the two layers exist together.

Thus, we can proceed under the assumption that the upper layer thickness is zero in the region where the lower layer flows, i.e. $H_u=0$ for $0 < \eta < \eta _l$, and that the lower thickness is zero after it terminates, i.e. $H_l=0$ for $\eta _l < \eta < \eta _u$. The volume constraints (4.5) become

(4.12a,b)\begin{equation} \int_0^{\eta_l} H_l \,\textrm{d} \eta = V, \quad \int_{\eta_l}^{\eta_u} H_u \,\textrm{d} \eta = 1. \end{equation}

In order to apply continuity conditions across the shock front at $\eta _l$, it is required that the velocity of the shock at the front of the heavier fluid moves at the same velocity as the back of the lighter fluid. In accordance with mass conservation, the fronts move at velocities $q_u / h_u$ and $q_l / h_l$. Hence, we impose the continuity condition $q_u / h_u = q_l / h_l$.

Integrating (4.8) subject to (4.2) and the continuity condition above, we obtain

(4.13) \begin{equation} \left.\begin{aligned} H_l &= \left( \dfrac{\eta}{MR} \right)^{1/2} &\quad (0 < \eta < \eta_l), \\ H_u &= \eta^{1/2} &\quad\, (\eta_l < \eta < \eta_u). \end{aligned}\right\} \end{equation}

The thicknesses of each layer therefore increase as $\eta ^{1/2}$ towards their respective flow fronts. Substituting these expressions into (4.12), we determine the front positions as

(4.14a,b)\begin{equation} \eta_l = \left( \frac{3}{2} V(MR)^{1/2} \right)^{2/3}, \quad \eta_u = \left( \frac{3}{2} ( 1 + V(MR)^{1/2}) \right)^{2/3}, \end{equation}

and their ratio as

(4.14c)\begin{equation} \frac{\eta_u}{\eta_l} = \frac{1}{V(MR)^{{1}/{2}}} + 1. \end{equation}

This expression is a dimensionless ratio that is a decreasing function of $V(MR)^{1/2}$, representing the shortening length of the lighter fluid as the volume and density of the heavier fluid is increased and the viscosity of the heavier fluid is reduced relative to that of the lighter fluid. Two illustrative solutions given by (4.13) and (4.14) are shown in figure 11 for (a) $M=10$ and (b) $M=0.1$, each for $V = 0.5$ and $R=2$. The solutions illustrate the intervening shock layer, the curved shape of the profiles, and the potential for the thickness to increase or decrease across the shock layer.

Figure 11. Similarity solutions describing the release of two fluids for (a) $M=10$ and (b) $M=0.1$, both with $R=2$, $V=0.5$. The similarity solutions are given by (4.13) and (4.14). (c) The front positions, $\eta _u$ and $\eta _l$, given by (4.14) plotted as a function of the dimensionless number $V(MR)^{1/2}$.

The self-similar shape of the heavier fluid predicted by (4.13) is exactly the same as that which would apply if just the heavier fluid were released in isolation (corresponding to the similarity solution describing a single layer on a slope Huppert Reference Huppert1982a). Therefore, the lighter fluid has no effect on the dynamics of the heavier fluid at long times. By contrast, the frontal region comprising the lighter fluid differs from the single-layer prediction because it must reside in front of the upstream region. In effect it is ‘pushed’ forward by the heavier fluid. The added distance travelled by the lighter fluid compared with the prediction that would apply if it were in isolation is represented by the $V(MR)^{1/2}$ term in (4.14b). As shown in figure 11(c), the front positions of the two layers are given by a universal function of the grouping $V(MR)^{1/2}$. For the situation where the heavier fluid is much less viscous than the lighter fluid (large $M$), the control of the propagation of the lighter fluid by the heavier fluid can be considerable even for a relatively small volume of the heavier fluid $V$. This dependence is encapsulated in the $M^{1/3}$ factor in (4.14). Figure 10(c,f) plots the predicted front position $x_u = \eta _u t^{1/3}$ and $x_l = \eta _l t^{1/3}$ as black dashed lines, showing excellent agreement with our numerical simulations of the full governing equations at long times.