Two-layer fluid flows on inclined surfaces

Abstract We present a theoretical and experimental study of the dynamics of two-layer viscous fluid flows on inclined surfaces, motivated by natural and industrial phenomena involving the interactions between two fluid layers. A general model describing the evolution of two fluids on an inclined substrate is developed and explored to reveal a rich variety of flow regimes for different modes of release. The asymptotic reduction of this problem due to the dominance of the along-slope component of gravity is shown to yield considerable analytical inroads compared with previous studies of multi-layer flow configurations, which have focused exclusively on the case of horizontal beds. For the canonical example in which two fluids are introduced at a constant flux, the flow forms two regions: an upstream region containing both fluids, and a downstream region comprised purely of the lighter fluid, with a sharp intervening jump in thicknesses between the two. By constructing similarity solutions, we establish a full regime diagram of the possible configurations over all asymptotic limits of the viscosity, flux and density ratios. For the release of two fixed volumes of fluid, the layers separate completely into two disjoint but connected regions, contrasting in essential structure from the constant flux case. Even a small volume of the heavier fluid is able to significantly accelerate the propagation of the lighter fluid in front of it. Excellent agreement is found between our theoretical predictions and the results of a series of laboratory experiments.


Introduction
Fluid systems involving the interactions between two layers driven by gravity arise widely in the natural world.Examples include the lubrication of glaciers by a subglacial hydrological system (Fowler 1982), the flow of stratified layers of glacial material with differing properties (Loewenherz, Lawrence & Weaver 1989), deformation of till underlying glaciers (e.g.Clarke 1987;MacAyeal 1989;Kowal & Worster 2020), the flow of lava where layers form due to differential heating (Griffiths 2000) and subsurface hydrological flows involving two fluids (e.g.Woods & Mason 2000).Related industrial and biological examples of two-layer systems include Leidenfrost droplets (e.g.Biance, Clanet & Quéré 2003), viscous non-wetting droplets on lubricated inclines (e.g.Smith et al. 2013), gravity-driven droplets impinging on a surface underlain by a thin layer of ambient fluid (Hodges, Jensen & Rallison 2004) and bacterial motion at the base of active matter drops driving flow over a substrate (Ramos, Cordero & Soto 2020).The rich variety of behaviour possible in each of these examples arises from the two-way dynamical coupling between the fluid layers.For example, Fowler (1982) applies equations governing the flow of glaciers lubricated by a hydrological system that are based on the coupling of two fluid flows: the viscous flow of the glacier and the hydrological flow of underlying meltwater in the till.In common with two-layer fluid flows, the dynamics produce a two-layer coupled system of kinematic wave equations.In order to understand some general aspects of the control and structure of such systems, we present here the first theoretical study of gravity-driven two-layer fluid systems on inclined substrates.
Thin-layer fluid systems over rigid surfaces are canonically studied as gravity currents, or intrusions of one fluid into another of a different density.Early work on viscous gravity currents at low Reynolds number used lubrication theory to develop a nonlinear diffusion equation governing the evolution of the layer thickness (Mei 1966;Smith 1969).Smith (1969) determined similarity solutions to this equation which describe the release of a fixed volume of fluid on a horizontal substrate in both two-dimensional and axisymmetric geometries, showing that the frontal positions grow as t 1/5 and t 1/8 , respectively.The solutions form broadly convex shapes with the maximum thickness at the input position and large interfacial gradients at the front.Huppert (1982b) generalized these analyses to other release conditions including the case of an input of constant flux, obtaining similarity solutions that grow as t 4/5 and t 1/2 in two-dimensional and axisymmetric geometries, respectively.For an inclined substrate, Huppert (1982a) determined a different similarity solution describing the release of a fixed volume of fluid, showing that the front position instead grows as t 1/3 .In this case, the layer thickens towards a finite thickness at its front (forming a shock), differing qualitatively in essential form from a gravity current on a horizontal substrate.By including also the effect of gravitational spreading due to thickness gradients, Lister (1992) presented a smooth solution for this frontal zone that is steady in the frame of the front.
Studies of two-layer gravity currents have focused to date on the configuration of a horizontal substrate.In particular, Woods & Mason (2000) considered the case of a two-layer gravity current in the context of a semi-infinite porous medium with a horizontal substrate, motivated by subsurface flows.These authors generalized the governing equations of the flow of a gravity current in a porous medium to allow for a second fluid layer and calculated similarity solutions describing the coupled evolution of the two fluid layers.Depending on the viscosity ratio, it was shown that a variety of forms of self-similar solutions are possible for fixed volume releases and constant flux inputs.For example, it is possible for either flow front to outpace the other, and for the lower layer to remain attached to its source position or separate away from it.
For the situation where a two-layer viscous gravity current propagates through an ambient fluid, Kowal & Worster (2015) considered the case of a horizontal substrate where each layer is fed by an input of constant flux.An important distinction between two-layer flows in a porous medium versus a free domain is that the latter introduces the possibility for lubrication by the lower fluid.The study determined similarity solutions Two-layer fluid flows on inclined surfaces for both two-dimensional and axisymmetric configurations and compared the predictions alongside experiments.It was found that the lubricating fluid separates the system into a broad upstream region of relatively mild slope connected to a downstream region comprised purely of the lighter fluid, with lubrication having the potential to increase the overall propagation rate considerably.
Other related studies of gravity currents in two-layer systems have considered situations where the upper fluid layer is confined between two horizontal boundaries (Gunn & Woods 2011;Guo, Zhang & Shi 2014;Pegler, Huppert & Neufeld 2014;Zheng et al. 2015a;Zheng, Rongy & Stone 2015b), and the release of a gravity current into a moving fluid layer (Eames, Gilbertson & Landeryou 2005;Pegler et al. 2017).These studies likewise illustrate the coupling between gravity-driven fluid layers, and the mutual control of the layers by the relative viscosity.Multi-layered and continuously stratified gravity currents have also been shown to exhibit significant variations in the overall shapes of gravity currents compared with the case of a single layer (Pegler, Huppert & Neufeld 2016).The analysis of two-layer axisymmetric gravity currents with equal densities over horizontal substrates has also been considered, with the finding that it is possible for such flows to form interfacial shocks (Dauck et al. 2019).Another group of related studies have analysed viscous instabilities of inclined two-layer films with a free surface (Loewenherz & Lawrence 1989;Loewenherz et al. 1989;Balmforth, Craster & Toniolo 2003) and of viscous stratified flows on slight inclines (Kliakhandler & Sivashinsky 1997).Loewenherz & Lawrence (1989) and Loewenherz et al. (1989) specify the thickness of the two layers as a basic state for their two-dimensional stability analysis, revealing transverse instabilities if the upper layer is less viscous.The dependence of the thicknesses on input flux, viscosity and density, remains an open question.
In this study we investigate for the first time the dynamics of two-layer gravity currents on inclined substrates.We conduct a complete theoretical exploration of the possible flow regimes and structures assuming thin-film flow at low Reynolds number.In solving the model system, we demonstrate the effectiveness of a finite-volume numerical scheme to predict the evolution of a multi-layer fluid system.We focus first on the canonical flow arising from the introduction of two fluids onto a flat slope of uniform inclination, exploring the solutions to the governing equations using a combination of analytical and numerical methods.Finally, we address the simultaneous release of two fluids with fixed volumes.Owing to simplifications arising under the dominant effect of the along-slope component of gravity, our study reveals considerable new analytical inroads for the analysis of two-layer gravity currents compared with previous studies that assume a horizontal substrate.In order to test our theoretical predictions and highlight physical effects not included in our model, we also present a series of laboratory experiments, comparing the observations against our theoretical predictions.
We begin in § 2 by developing from first principles our theoretical model for two-layer gravity currents for general topography.In § 3 we explore the predictions of the model in the case of a constant flux input.Section 4 addresses the case of a fixed volume release.In § 5 we present our laboratory study.We end in § 6 by summarising our key conclusions.

Theoretical model development
Consider a two-dimensional flow formed of two fluid layers flowing on a rigid bed z = b(x), where x and z are the horizontal and vertical coordinates (figure 1).We assume the fluids are incompressible and Newtonian, and allow their viscosities and densities to differ.We neglect surface tension and assume that no mixing between the fluid x l (t) layers occurs.We refer to variables representing the upper layer using the subscript u and the lower layer using subscript l.For layer i, let u i = (u i , w i ) denote the flow velocity, h i (x, t) the thickness, p i (x, z, t) the pressure, ρ i the density and μ i the dynamic viscosity, where t is time.The layers are each assumed to terminate at a time-dependent front with position x = x i (t).
We consider vertical-shear-dominated flow.The x-and z-components of the lubrication equations are respectively.The free-stress condition on the upper surface is and the conditions of continuity on the stress and velocity at the interface read as respectively.The no-slip condition at the base is The model above is based on an approximation of thin layers.It should be noted that in situations where the lower fluid is much less viscous than the upper layer, we also require the viscosity ratio (μ u /μ l ) to be sufficiently small that the extensional stresses in the upper layer are negligible compared with the vertical shear stress.Integrating (2.1b) and applying the continuity of pressure, we determine the pressure fields for the upper and lower fluids, where p a is a constant reference pressure.On substituting (2.3) into (2.1a) and integrating the resulting differential equations subject to (2.2), the velocity profiles of each layer are Both layers are driven by the component of gravity resolved along the bed, represented by the terms proportional to db/dx.The upper layer is also driven by the gradient in its surface slope, as represented by the sum ∂(h u + h l )/∂x in the first set of grouped terms, i.e. the first two lines of (2.4a), and moderated by interfacial stresses exerted by the lower layer, represented by the final line of (2.4a).The lower layer is likewise driven by the gradient in the weight of the fluid columns above it, represented by the terms involving ∂h i /∂x in the first line of (2.4b), and couples to the upper layer through the terms involving h u and ∂h u /∂x.The depth-integrated continuity equations for the two layers are where q i denote the volume fluxes per unit width of the upper and lower layers, namely, respectively.Substituting (2.4) into this integral, the flux expressions are (2.5d) Each term can be interpreted as the ratio of a driving stress due to hydrostatic pressure gradients to a viscous resistance.The first line of (2.5c) describes the balance between the gradient in hydrostatic pressure caused by the gradient in the height of the upper interface, multiplied by (i) the variation in shear rate (the term involving h 3 u ) and (ii) the 'basal drag force' due to the resistance to shear induced in the lower layer (involving h 2 u h l ).The second line represents the contribution to the flow of the upper layer induced by the gravity-driven flow of the lower layer.Equation (2.5d) is similarly structured.The first line represents the pressure gradient caused by the gradient in weight of the fluid columns arising from the variations in both the upper surface and the interface and the resistance due to the shear stress via the viscosity of the upper layer (represented by the term containing μ u h 3 l /μ l ).In (2.5d) there is no analogue of the term representing a 'drag force' from the other layer (represented by the term proportional to h 2 u h l in (2.5c)).This is because, unlike at the substrate, there is no stress exerted at the upper surface, and, hence, the upper layer experiences no resistance to moving the lower layer.The second line of (2.5d) represents the contribution to the flow of the lower layer induced by the gravity-driven flow of the upper layer.The flux expressions (2.5c,d) generalize those applying to the case of a horizontal substrate to allow for a general substrate shape b(x) (note that we exactly recover the equations of Kowal & Worster (2015) by setting b = 0 and by adopting their definitions for layer thicknesses: the total thickness as H(x, t) = h u + h l and the lower layer thickness as h(x, t)).The governing system of equations (2.5) allows for the existence of both two-layer regions (in which both h l > 0 and h u > 0) and single-layer regions in which either h l = 0 and h u > 0 or h u = 0 and h l > 0. For example, if h l = 0, (2.5c) reduces to which recovers the equation describing a single layer on a slope (Lister 1992).The governing equations (2.5) describe the evolution of a two-layer fluid system subject to general variations in the height of the substrate, b(x), and any given initial condition.The two-way coupling between the layers represented by the presence of both h u and h l in both flux expressions (2.5c,d) produces a rich mathematical system to be explored.
The governing equations form a parabolic initial-value problem that can be solved subject to suitable boundary conditions on the thicknesses and fluxes of the layers.At the two flow fronts, we impose representing the independent conditions of vanishing thickness and flux at the flow fronts.Continuity conditions on flux and thickness at the transitions between single-and two-layer regions are imposed automatically by our numerical solver.The other boundary conditions will be specific to the configuration considered, e.g. an input condition on the flux or initial volume conditions, to be introduced in the relevant sections.As we will demonstrate, for flow down a slope, the system can be reduced further under the neglect of the gradients in layer thicknesses (in analogy with corresponding simplifications made in the analysis of single-layer gravity currents on inclined substrates; Mei 1966;Fowler & Larson 1978;Huppert 1982a;Lister 1992).This will reduce (2.5) to a hyperbolic system formed of coupled nonlinear kinematic wave equations (qualitatively similar to those proposed for coupling a glacier and its hydrological system; Fowler 1982), which we will detail later in § 3.3.In the special situation where the two layers have equal density (R = 1), the order of the flux expression (3.6c) reduces because it does not independently depend on the gradient in lower layer thickness.In our reduced asymptotic formulation our flow fronts form shocks and the imposition of zero thickness at the flow fronts is likewise abandoned.Therefore, our theory applies even in the R = 1 situation.Two-layer fluid flows on inclined surfaces

Constant flux input of two layers
We begin by considering the situation where two fluids are introduced onto a slope at constant rates.For this, we impose the flux conditions q l (0, t) = q l0 , q u (0, t) = q u0 , (3.1a,b) where q l0 and q u0 represent the fluxes into the lower and upper layer, respectively.As a consequence of these conditions, the following integral constraints on the total volume are satisfied: These expressions do not add additional constraints to this configuration (because we already have an input flux condition) but will be utilized in our later analysis.As a simple case to demonstrate the general dynamics, we assume a constant slope, where α is a positive constant.
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, , H = αL, (3.4a-c) respectively.These represent the scales on which the slope of the substrate, α = −db/dx, becomes important.We define non-dimensional (hatted) variables by On dropping the hats, (2.5a) becomes for the upper and lower layers, respectively, and the flux expressions become where M = μ u /μ l and R = ρ l /ρ u .The boundary conditions (2.5 f ) become The input conditions (3.1a,b) become The integral constraints (3.2a,b) become The dimensionless model above depends on three dimensionless parameters, 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 (1982) applies to describe the two-way interplay between a glacier and its hydrological system; the main difference being that the flux expressions that Fowler (1982) 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.

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 1992).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 2002) 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 ).
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 2 f ), 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.

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 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.
Integrating the steady-state forms of (3.6a-c) subject to (3.6e), we obtain 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 3 l .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.Layer thickness 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, for M → 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.
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 → ∞, we can neglect all terms in (3.7) except those multiplied by M. Solving the resulting simplified 10 -2 10 -1 10 0 10 1 10 2 10 -1 10 0 10 1 10 2 system, we obtain the asymptotic predictions for M → ∞, where the functions 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 φ(QR) and ψ(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, ψ and φ 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.We note that the M → ∞ limit in (3.9) encompasses limits of small and large input flux ratio Q.In the limit QR → 0, (3.9) reduces to 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 ψ and φ, respectively (figure 5).We name this situation regime A. We can likewise reduce (3.9) in the limit QR → ∞ which asymptotes to 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 → ∞, 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 = MRh 2 l /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 → 0, (3.7) yields the asymptotic predictions 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 → 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∂u u /∂z, where ∂u u /∂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 = u l dz = 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 → ∞, 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 → ∞ 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 1/6Q and Q 3/4R and represents the region in which the numerical solutions fall within 10 % of the asymptotic prediction for the simultaneous limits M → ∞ and Q → 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 M 1/6Q, describing the 'rigid-lid regime' in which Q → 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 M 8/27Q 2 R and describes the region in which the numerical solutions lie within 10 % of the asymptotic prediction for the M → 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 3/4R and M 8/27Q 2 R where the numerical solutions fall within 10 % of the asymptotic predictions for the Q → ∞ 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 10 -6 10 -4 10 -2 10 0 10 2 10 4 10 6 Regime A . 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 .

Regimes
Thickness predictions in the two-layer region 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.
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.

Two-layer fluid flows on inclined surfaces
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 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 (2015).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.

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 ∂h i /∂x become progressively smaller relative to db/dx 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 db/dx) provides the dominant contribution.In dimensionless forms, these expressions reduce to 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 ∼ h 3 /x from (3.13) and hx ∼ t from (3.6e) which are consistent with previous studies of single-layer flows on slopes (e.g.Huppert 1982a;Lister 1992).Combining these, we obtain x ∼ t and h ∼ 1.These scalings differ from those arising for constant flux inputs over horizontal substrates, where the propagation rates go as t 1/5 (Huppert 1982b;Kowal & Worster 2015).
The scaling analysis above indicates the existence of similarity solutions of the form h = h(η), where η = t −1 x, (3.14) with front positions x u = η u t and x l = η l t, where η u and η l are constants.The similarity scaling implies that the flow becomes increasingly slender at long times, ∂h/∂x ∼ 1/t, thus confirming that the approximation of neglecting ∂h/∂x is self-consistent.Recasting (3.6) in terms of the similarity coordinate (3.14), the equations reduce simply to dq i /dη = 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 dq i /dη = 0 subject to the input conditions (3.6e) yields 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 ∼ 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 < η < η 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, η l < η < η u , we set h l = 0 in (3.15a) to obtain The heights in the two regions can thus be determined from (3.15) before the extents of these regions η l and η u are known.The integral constraints (3.6 f ) yield 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 giving explicit formulae for the positions of the flow fronts η l and η 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 Layer thickness 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 η l can then be evaluated using (3.17a).The volume of the lighter fluid in the two-layer region can be determined as η l h u , leaving a residual volume of 1 − η 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 (η u − η l ) = 1 − η 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 1992).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 = η u t and x l = η 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.

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, η u and η 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, η l − η u , is shown in figure 9(a).As M is increased, η u and η 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 M Q Front position  results of the two-layer region given by (3.8) into (3.17) to give the frontal positions as M → 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 η = 0 and η 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 .For M → ∞, the length of the two-layer region becomes comparable to the length of the full system.Substituting (3.9) into (3.17),we obtain 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 ζ(QR) that is independent of M. The value of ζ(QR) across the entire range of QR is strictly bounded within the narrow range 2 < 3 1/3 ζ < 3. The extent of this band of values is indicated by the grey shaded rectangle in figure 9(a).Since η u − η l must also asymptote to the universal value 3 −1/3 as M → 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 → 0, we substitute (3.12) into (3.17) to give 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 → ∞, we substitute (3.11) into (3.17) to give 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 η u − η l ∼ 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

Regimes
Thickness and front position predictions 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.
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.

Fixed volume release of two layers
We consider now the situation where two fluid volumes are released together onto an inclined surface.In this case, we assume the layers are initialised with profiles h u (x, 0) and h l (x, 0), with volumes prescribed by the integral constraints, (4.1a,b)where v u0 and v l0 are the volumes of the upper and lower layers, respectively.We also assume that no fluid enters the domain at x = 0 by imposing q u (0, t) = 0, q l (0, t) = 0, Intrinsic scales and dimensionless model system We non-dimensionalise the system using the intrinsic time, length and height scales, We define dimensionless (hatted) variables according to 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)where V = v l0 /v u0 .The dimensionless system above depends on three dimensionless parameters the ratios of viscosity, density and fluid volumes, respectively.

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.

Similarity solutions
We note that (3.6a-c) and (4.5) yield the scalings h/t ∼ h 3 /x and hx ∼ 1, respectively.Eliminating h, we obtain x ∼ t 1/3 and, hence, h ∼ t −1/3 , which are consistent with the scalings arising in the case of a single layer (Huppert 1982a).Thus, we define the relevant similarity variables as The front positions of the two layers therefore evolve as x l = η l t 1/3 and x u = η 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 where q u and q l are the flux expressions (3.13) in the limit where the along-slope contribution to the flow dominates.
Since these equations depend on η, 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 (η).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 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 2 l in (4.9b) and simplifying, we obtain 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 < η < η l , and that the lower thickness is zero after it terminates, i.e.H l = 0 for η l < η < η u .The volume constraints (4.5) become In order to apply continuity conditions across the shock front at η 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 The thicknesses of each layer therefore increase as η 1/2 towards their respective flow fronts.Substituting these expressions into (4.12),we determine the front positions as and their ratio as 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.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 1982a).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 = η u t 1/3 and x l = η l t 1/3 as black dashed lines, showing excellent agreement with our numerical simulations of the full governing equations at long times.

Laboratory study
In order to test our theoretical results, we conducted a series of laboratory experiments and compared the data with our predictions.Our experimental apparatus comprised of a sloped acrylic channel of dimensions 0.81 m long and 0.14 m wide with 0.05 m high sidewalls (figure 12).A lock gate was formed by placing two sealed acrylic barriers across the sloped channel, approximately 8 cm apart.We conducted two different kinds of experiments: one in which the lighter fluid was more viscous than the heavier fluid (M > 1) and the other in which the lighter fluid was less viscous than the heavier fluid (M < 1).In each, we used pure Karo corn syrup as one of the two fluids.For M > 1, we used a solution of potassium carbonate dissolved in the corn syrup for the more dense, less viscous heavier fluid.For M < 1, we used a solution of the corn syrup dissolved in water for the lighter, less viscous fluid.The heavier and lighter fluids were dyed blue and red, respectively.
The dynamic viscosity of pure corn syrup was determined as a function of temperature using a rheometer (with an error of ±5 × 10 −4 Pa s) at the outset of our experimental investigation.The viscosity of the pure corn syrup, used as one of the layers in each experiment, was determined before each run by measuring its temperature using a thermometer and applying the viscosity vs temperature relationship determined previously from the rheometry.The viscosity of the second fluid used in each experiment was determined by measuring the fall speed of ball bearings through both the corn syrup (whose viscosity was known) and the second fluid, and multiplying the ratio of the fall speeds with the known viscosity of the corn syrup determined previously using the rheometer.The 4 mm ball bearings were dropped into a beaker with a diameter of 15 cm.The densities of the fluids were measured using a hydrometer.The angle of the slope was measured using a digital inclinometer (0.9 • ± 0.05 • ).To prepare the experiment, we first poured the heavier fluid and then the lighter fluid sequentially into the lock.The volume of each fluid released during the experiment was calculated by measuring the Our study has demonstrated considerable analytical inroads for the analysis and classification of regimes applying to multi-layer gravity currents, as compared with previous studies, all of which have focused on horizontal substrates (Woods & Mason 2000;Kowal & Worster 2015;Pegler et al. 2016).Our classification of the regimes for inclined two-layer gravity currents fed by a constant flux input show that the classification can be reduced entirely to the consideration of a cubic equation which independently partitions the M, Q parameter space into four regimes.Remarkably the partitioning curves between our regimes and those in Kowal & Worster (2015) have the same M and Q scalings but represent fundamentally different regimes, specifically, ours is a steady-state solution while theirs is a time-dependent similarity solution.Our framework provides a clearer explanation for the regimes and our asymptotic reduction enables exact descriptions for the separatrices.
The analysis of this paper provides a foundation for the theoretical treatment of two-layer fluid flows on sloped surfaces widely seen in geophysical, environmental and industrial applications.While we have focused here on the idealised case of two-dimensional Newtonian fluid layers, the phenomena revealed are likely to underlie more complex generalized examples of layered flows.Important additional effects we have neglected include those of a porous domain, of surface tension at the interface, of non-Newtonian rheology, of three-dimensionality, and of thermodynamic controls on viscosity and density.Each of these effects presents directions for new research for which the model and analytical inroads revealed here provide a basis for generalization.Since the cells are the same size, a first-order Taylor series is sufficient for stability (LeVeque 1992).
In order to prevent the introduction of spurious maxima and minima in the numerical solution (particularly near any shock fronts), we implement a flux limiter.The limiter used for each layer thickness is with a similar expression for h L .The same limiter is applied to ∂h/∂x.This allows the numerical scheme to maintain a shock front if required.Our finite-volume method implements a flux interpolation scheme between the values to the left and right of each interface according to where the value to the right of the interface (denoted by R+) is the value to the left of the next cell.The value to the left of the interface (denoted by L+) is the value to the right of the cell that is located immediately to the left of cell i.The same construction has been applied to calculate the values of the derivatives at each cell interface.These thickness and derivative values are used to compute the values of q u and q l using expressions (3.6b,c) to the left and right of each cell.A standard Godunov upwind scheme is used to interpolate fluxes (LeVeque 1992).This has been chosen both because it is particularly suited to identifying discontinuities and is inherently mass conservative.We then integrate (3.6a) in time using a forward Runge-Kutta scheme to the next time step and repeat.
The conditions for global volume conservation and continuity of flux and velocity at the lower layer front are automatically satisfied by this numerical method.There is no need for special treatment of a cell that spans a flow front because, like every other point on the numerical grid, a condition of flux continuity is imposed across every cell.It should be noted that the flow fronts predicted by the model generally involve steep gradients in layer thicknesses.Since our flux limiter prevents spurious maxima or minima, our solver is able to resolve these steep gradients sufficiently (without the need to specify a minimum thickness, for example).To increase the spatial resolution at which the flow fronts are extracted from our solutions, we consider the three nodes upstream of x i at which the thickness vanishes, fit a line through these three points and extrapolate the line until it intersects the x-axis.The numerical solver was benchmarked in three ways: first, we checked that the total mass of the numerical solution was conserved.Next, we checked that our numerical solver is stable for very low spatial resolutions.Finally, we checked the agreement with the steady-state asymptotic solutions for both single-layer and two-layer systems (figure 4).

Figure 1 .
Figure 1.Schematic illustrating the general configuration of two fluid layers flowing along an inclined rigid substrate b(x) with an upper fluid (red) of thickness h u (x, t), viscosity μ u and density ρ u with a lower fluid (blue) of thickness h l (x, t), viscosity μ l and density ρ l .The front positions for the upper and lower layers are marked as x u (t) and x l (t), respectively.

917Figure 2 .
Figure 2. Illustrative numerical solutions to the full dimensionless model (3.6) for the situation where two fluids are injected continuously.In example A (panelsa-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.

Figure 3 .
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.

Figure 4 .
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 = μ u /μ 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).

Figure 5 .
Figure 5.The functions ψ(ξ) and φ(ξ) 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.
.13b) 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 917 A54-15 https://doi.org/10.1017/jfm.2021.

Figure 9 .
Figure 9.The length of the frontal (single-layer) region η u − η 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 ζ(QR).The inset in panel (a) plots the function ζ(QR) defined by (3.19c), representing the dimensionless length in the limit M → ∞.

917Figure 10 .
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.

Figure 11 .
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, η u and η l , given by (4.14) plotted as a function of the dimensionless number V(MR) 1/2 .
Figure 12.Cross-section and overhead view of the experimental set-up.