Lubricated viscous gravity currents of power-law fluids. Part 1. Self-similar flow regimes

Abstract

We examine the gravity-driven flow of thin films of viscous fluid of power-law rheology, lubricated by another power-law viscous fluid from below. Such flows are relevant to a range of geophysical and industrial settings, including the flow of ice sheets and fast-flowing ice streams lubricated from below by a layer of subglacial till. We model both layers using lubrication theory in two-dimensional and axisymmetric settings, in the limit in which vertical shear provides the dominant resistance to the flow. The flow is self-similar if the power-law exponents, describing the rheology of the two layers of fluid, are equal. We examine the similarity solutions in both geometries and describe the flow in terms of four distinct flow regimes ranging from thin films of viscous fluid coating a more viscous fluid from above, to thin layers of fluid lubricating a more viscous fluid from below. In contrast to the former scenario, a thin film of a low-viscosity fluid strikingly alters the dynamics of a more viscous fluid when it lubricates it from below: the overlying layer thins, the upper surface gradients lessen and most of the shear is confined to the lower layer. Such features amplify, and this flow regime becomes increasingly dominant when the viscous fluids are shear thinning, like the deformation of glacial ice on the large scale. This flow regime is most relevant to the flow of lubricated ice sheets, which thin and accelerate, forming fast-flowing ice streams, when they are well lubricated from below.

1. Introduction

Interest in gravity-driven flows of viscous fluids spans the fields of geophysics, geology and industrial engineering, involving a range of environmental and industrial phenomena (Simpson 1967, 1982; Benjamin 1968; Hoult 1972; Linden 1999; Huppert 2006). Such flows feature striking dynamics when they involve the superposition of one viscous fluid over another, particularly when the underlying fluid is less viscous (Kowal & Worster 2015; Kumar et al. 2021). These are relevant to the flow of large-scale, glacial ice sheets over till, or a mixture of clay and subglacial sediments (Weertman 1957, 1964; Nye 1969; Kamb 1970; Fowler 1981; Engelhardt et al. 1990; Kyrke-Smith, Katz & Fowler 2014), the flow of layers of lava (Balmforth & Craster 2000) and layered flows in porous media (Woods & Mason 2000), for example.

The glacial analogue involves the flow of a power-law fluid representing ice, the rheology of which is commonly referred to as Glen's flow law in the glaciological literature (Glen 1955). It is more difficult to quantify the rheology of subglacial till, though frequently adopted rheologies for till include a plastic rheology with a Mohr–Coulomb yield stress, relevant on the small scales (Kamb 1970; Tulaczyk, Kamb & Engelhardt 2000; Kamb 2001), or a viscous power-law rheology, relevant to large scales of kilometres or more (Boulton & Hindmarsh 1987; Hindmarsh 1997). However, recent work by Tulaczyk (2006) suggests that the difference between the rheology at small and large scales is not at all clear, as observations of the subglacial till beneath the well-lubricated Whillans Ice Stream, West Antarctica, indicate that its rheology is independent of scale. These observations have been made possible following short speed-up events on the ice plain of Whillans Ice Stream, which provided the opportunity to examine the in situ rheology of till beneath Whillans Ice Stream, integrated over roughly 10–100 km, and compare it to the rheology of small laboratory samples of till from beneath the same ice stream.

We investigate the gravity-driven flow of two superposed thin films of viscous fluids of power-law rheology in two-dimensional and axisymmetric geometries. The flow reduces to that of Kowal & Worster (2015) in the Newtonian limit. We model the flow using lubrication theory by balancing viscous and buoyancy forces, and we assume that inertial effects and the effects of mixing and surface tension at the interface between the layers are negligible. Although it is possible to generalize our work to a mathematical model for which the rheology of the two viscous fluids involves power-law exponents that may, in general, be unequal, we focus on the case in which the two power-law exponents are equal. The flow is self-similar in this scenario in both geometries.

Four distinct flow regimes emerge, from thin films of viscous fluid lubricating a more viscous fluid from below, akin to flows over slippery substrates such as the flow of glacial ice over till, to thin films of fluid coating a more viscous fluid from above, and thin films of fluid forming a more viscous crust over the main current, akin to viscous-gravity currents with a solidifying crust (Fink & Griffiths 1990; Griffiths & Fink 1993; Fink & Griffiths 1998; Balmforth & Craster 2000). In contrast to thin films coating a more viscous fluid from above, the dynamics of the flow are strikingly altered when the less viscous fluid lubricates the more viscous layer from below. This difference becomes greater when the two layers are shear-thinning, as we demonstrate in this paper.

Experiments of similar lubricated viscous gravity currents have been conducted by Kumar et al. (2021), in which the upper layer consisted of a viscous fluid of power-law rheology, and the lower layer consisted of a Newtonian viscous fluid, and by Kowal & Worster (2015), in which both layers were purely Newtonian. In contrast to the Newtonian limit, experiments show that the propagation of a shear-thinning overlying layer follows a different scaling law, in terms of both the thickness and the extent of the two layers of viscous fluid (Kumar et al. 2021).

Experiments show that both Newtonian and non-Newtonian lubricated viscous gravity currents are prone to a viscous fingering instability (Kowal & Worster 2015; Kumar et al. 2021), investigated theoretically in the Newtonian limit by Kowal & Worster (2019a,b) and Kowal (2021), along with the effects of various mechanisms that stabilize the flow, particularly for large wavenumbers. These fall into a class of instabilities referred to as non-porous viscous fingering instabilities (Kowal 2021). We investigate the non-Newtonian analogue, namely the stability of lubricated viscous gravity currents of power-law fluids, in a companion paper (Leung & Kowal 2022), henceforth referred to as Part 2.

We present our theoretical development in § 2, we formulate the governing equations in similarity variables in § 3, we explore various flow regimes in 4 and summarise with concluding remarks in § 5.

2. Theoretical development

We consider the gravity-driven flow of two superposed thin films of viscous fluid of different physical properties, spreading over a horizontal substrate as illustrated in the schematic of figure 1. We assume that the two viscous fluids are immiscible, and that the effects of inertia and surface tension between the two layers are negligible. We use lubrication theory to develop a mathematical model of the flow by balancing viscous and buoyancy forces, assuming that vertical shear stresses provide the dominant resistance to the flow. The flow consists of a lubricated region, in which one layer of fluid intrudes beneath another, and a no-slip region, beyond the intrusion front of the underlying layer. The position of the intrusion front, or lubrication front, is denoted by $\boldsymbol {x}=\boldsymbol {x_L}$, while the position of the leading edge of the no-slip current is denoted by $\boldsymbol {x}=\boldsymbol {x_N}$. The upper surface is denoted by $z=H(\boldsymbol {x},t)$ for the upper layer and $z=h(\boldsymbol {x},t)$ for the lower layer. The velocity field within the upper and lower layers are denoted by $\boldsymbol {u}$ and $\boldsymbol {u_l}$, respectively.

Figure 1. Schematic of the flow of two superposed thin films of power-law viscous fluids spreading under gravity over a horizontal substrate in a two-dimensional geometry. Schematic adapted from Kowal & Worster (2015).

The two layers are of different viscosities $\mu$ and $\mu _l$ and densities $\rho$ and $\rho _l$, where the subscript $_l$ denotes quantities related to the lower layer. We assume that both films of fluid follow a non-Newtonian, power-law rheology, also known as the Ostwald-de Waele power law, so that

(2.1a,b)\begin{equation} \mu = \tilde \mu \left|\frac{\partial \boldsymbol u}{\partial z}\right|^{{1}/{n}-1},\quad \mu_l = \tilde \mu_l \left|\frac{\partial \boldsymbol u_l}{\partial z}\right|^{{1}/{n_l}-1}, \end{equation}

within the limits of lubrication theory, where $\tilde \mu$ and $\tilde \mu _l$ are constant. We derive governing equations for equal power-law exponents ($n= n_l$) and note that unequal exponents ($n\neq n_l$) set an intrinsic length scale for the flow, which precludes self-similarity. This can be seen by considering the flow of a single-layer viscous gravity current of a power-law fluid, in two configurations, each of which involves a different viscous fluid of a different power-law exponent. The flows in each of these configurations are self-similar (Sayag & Worster 2013), and the leading edges of these flows propagate at a rate given by a time-dependent power law of different exponents. Balancing these (time-dependent) similarity length scales of these two flows gives rise to a characteristic length and time scale at which the spreading of the two fluids balances had they been superposed. Exact expressions for this length scale are presented in § 3 in two geometries. However, when $n=n_l$, no such length scale emerges (it can be seen that the said length scale diverges) and the systems tend towards a self-similar flow at late times starting from various initial conditions, similar to what has been demonstrated for single-layer viscous gravity currents (Ball & Huppert 2019). We focus on self-similar flows in the analysis of this paper and the companion paper. These correspond to the flows of Kowal & Worster (2015) in the long-time limit for $n=1$.

We examine the flow in two geometrical settings: (i) a two-dimensional geometry, in which both layers of fluid spread unidirectionally with no cross-flow variations (so that the flow is uniform in the cross-flow variable $y$); and (ii) an axisymmetric geometry, in which both layers of fluid spread radially outward from a specified origin, with no azimuthal variations (so that the flow is uniform in the azimuthal variable $\theta$). Both layers of fluid are supplied at constant line flux at $x=0$, denoted by $q_0$ for the upper layer and $q_{l0}$ for the lower layer in the two-dimensional setting, and at constant point flux at the origin, denoted by $Q_0$ for the upper layer and $Q_{l0}$ for the lower layer in the axisymmetric setting.

In terms of applications to the geophysical scenario of ice lubricated by meltwater-saturated till, we note that the source flux conditions are merely an idealisation for mathematical convenience. Realistically, the accumulation of both layers is actually distributed throughout the domain, with ice accumulation resulting from net snowfall minus melting over the whole domain, and till accumulation resulting from traction and bed failure over the whole domain, with locally distributed water sources from drainage of supraglacial meltwater down crevasses and other water outlets to the bottom of the ice, along with a distributed water source from basal melting throughout the domain. We have idealised both distributed processes as a point source, which is an idealisation necessary for the emergence of similarity solutions, rather than one that directly reflects actual geophysical settings.

For brevity, we derive the governing equations for both geometries simultaneously. To do so, we denote the partial derivates $\partial h/\partial x$ and $\partial H/\partial x$, in the two-dimensional setting, and the partial derivatives $\partial h/\partial r$ and $\partial H/\partial r$, in the axisymmetric setting, collectively by $\delta h$ and $\delta H$, for compactness. We define the horizontal velocities, $\boldsymbol u$ and $\boldsymbol {u_l}$, and depth-integrated fluxes, $\boldsymbol q$ and $\boldsymbol {q_l}$, of the upper and lower layers, respectively, and the position of the lubrication front $\boldsymbol x=\boldsymbol {x_L}$ and leading edge $\boldsymbol x=\boldsymbol {x_N}$ as

(2.2a–f)\begin{equation} \boldsymbol u = u \boldsymbol{e_x}, \quad \boldsymbol{u_l} = {u_l} \boldsymbol{e_x}, \quad \boldsymbol q = q \boldsymbol{e_x}, \quad \boldsymbol{q_l} = {q_l} \boldsymbol{e_x}, \quad \boldsymbol{x_L} = x_L \boldsymbol{e_x}, \quad \boldsymbol{x_N} = x_N \boldsymbol{e_x}, \end{equation}

in the two-dimensional setting, where $\boldsymbol {e_x}$ is the unit basis vector defining the positive $x$ direction. Similarly, in the axisymmetric setting, we define these quantities as

(2.3a–f)\begin{equation} \boldsymbol u = u \boldsymbol{e_r}, \quad \boldsymbol{u_l} = {u_l} \boldsymbol{e_r}, \quad \boldsymbol q = q \boldsymbol{e_r}, \quad \boldsymbol{q_l} = {q_l} \boldsymbol{e_r}, \quad \boldsymbol{x_L} = r_L \boldsymbol{e_r}, \quad \boldsymbol{x_N} = r_N \boldsymbol{e_r}, \end{equation}

where $\boldsymbol {e_r}$ is the unit basis vector defining the radially outward direction.

2.1. Lubricated region

A horizontal force balance under the approximations of lubrication theory gives

(2.4)$$\begin{gather} \frac{\partial}{\partial z}\left(\mu \frac{\partial \boldsymbol u}{\partial z}\right) = \rho g\boldsymbol{\nabla} H, \quad h< z< H, \end{gather}$$

(2.5)$$\begin{gather}\frac{\partial}{\partial z}\left(\mu_l \frac{\partial \boldsymbol u_l}{\partial z}\right) = \rho g(\mathcal{D}\boldsymbol{\nabla} h+\boldsymbol{\nabla} H), \quad 0< z< h, \end{gather}$$

for the two thin films of fluid, where the pressure is hydrostatic to leading order,

(2.6)$$\begin{gather} p= p_0 +\rho g(H-z) \quad \text{for}\ h< z< H, \end{gather}$$

(2.7)$$\begin{gather}p_l = p_0 +\rho g(H-h)+\rho_lg(h-z) \quad \text{for}\ 0< z< h. \end{gather}$$

Apart from the power-law exponents $n$ and $n_l$, which we assume to be equal, the governing equations depend on the dimensionless density difference $\mathcal {D}$ and consistency ratio $\mathcal {M}$, and flux ratio $\mathcal {Q}$ defined by

(2.8a–c)\begin{equation} \mathcal{D} = \frac{\rho_l-\rho}{\rho}, \quad \mathcal{M} = \frac{\tilde \mu}{\tilde \mu_l}, \quad \mathcal{Q} = \frac{q_{l0}}{q_{0}}. \end{equation}

The above definition of $\mathcal {Q}$ holds for the two-dimensional geometry. In the axisymmetric geometry, it is given by $\mathcal {Q} = {Q_{l0}}/{Q_{0}}$.

Integrating (2.4)–(2.5) describing the horizontal force balance, subject to the no-slip boundary condition at $z=0$, continuity of velocity and shear stress at $z=h$, and the no-stress condition at the upper surface $z=H$, yields the following expressions for the horizontal velocities of the upper and lower layers,

(2.9)\begin{align} u &= \frac{1}{n+1}\left(\frac{\rho g}{\tilde \mu_l}\right)^{n} \frac{1}{\mathcal{D}\delta h+\delta H}[\left|(H-h)\delta H\right|^{n+1} -\left|h(\mathcal{D}\delta h+\delta H)+(H-h)\delta H\right|^{n+1}] \nonumber\\ &\quad +\frac{1}{n+1}\left(\frac{\rho g}{\tilde \mu}\right)^n [(H-z)^{n+1} - (H-h)^{n+1}] |\delta H|^{n-1} \delta H, \end{align}

and

(2.10)\begin{align} u_l&=\frac{1}{n+1}\left(\frac{\rho g}{\tilde \mu_l}\right)^{n} \frac{1}{\mathcal{D}\delta h+\delta H} [\left|(h-z)(\mathcal{D}\delta h+\delta H)+(H-h)\delta H\right|^{n+1}\nonumber\\ &\quad -\left|h(\mathcal{D}\delta h+\delta H)+(H-h)\delta H\right|^{n+1}], \end{align}

respectively.

Integrating across the depth of the upper layer of the lubricated region gives the following expression for the depth-integrated flux of fluid in the upper layer,

(2.11)\begin{align} q &= \frac{1}{n+1}\left(\frac{\rho g \mathcal{M}}{\tilde \mu}\right)^{n} \frac{H-h}{\mathcal{D}\delta h+\delta H} \left[\left|(H-h)\delta H\Big|^{n+1} - \Big|h(\mathcal{D}\delta h+\delta H) +(H-h)\delta H\right|^{n+1}\right] \nonumber\\ &\quad -\frac{1}{n+2} \left(\frac{\rho g}{\tilde \mu}\right)^n (H-h)^{n+2} |\delta H|^{n-1}\delta H, \end{align}

comprising Couette terms arising from viscous coupling with the lower layer and Poiseuille terms arising from gravitational spreading under its own weight. In the lower layer, the depth-integrated flux is given by

(2.12)\begin{align} q_l &= \frac{1}{n+1}\left(\frac{\rho g\mathcal{M}}{\tilde \mu}\right)^{n} \frac{1}{\mathcal{D}\delta h+\delta H}\left[-\frac{1}{n+2}\frac{1}{\mathcal{D}\delta h+\delta H} \left[\left|(H-h)\delta H\right|^{n+1} (H-h)\delta H \right.\right.\nonumber\\ &\quad \left.-\left|h(\mathcal{D}\delta h+\delta H)+(H-h)\delta H\right|^{n+1} (h(\mathcal{D}\delta h+\delta H)+(H-h)\delta H)\right] \nonumber\\ &\quad \left.\vphantom{\frac{1}{n+2}\frac{1}{\mathcal{D}\delta h+\delta H}} -h\left|h(\mathcal{D}\delta h+\delta H)+(H-h)\delta H\right|^{n+1}\right], \end{align}

comprising viscous coupling terms as well as terms arising from the spreading of the lower layer under its own weight and the weight of the overlying layer of fluid.

The evolution of the thickness of the two films of fluid is determined by the mass conservation equations

(2.13a,b)\begin{equation} \frac{\partial (H-h)}{\partial t} ={-}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol q, \quad \frac{\partial h}{\partial t} ={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{q_l}, \end{equation}

in the lubricated region.

2.2. No-slip region

In the no-slip region, the horizontal velocity and depth-integrated flux reduce to

(2.14)\begin{equation} u=\frac{1}{n+1}\left(\frac{\rho g}{\tilde \mu}\right)^n [(H-z)^{n+1} - (H)^{n+1}] |\delta H|^{n-1} \delta H, \end{equation}

and

(2.15)\begin{equation} q={-}\frac{1}{n+2} \left(\frac{\rho g}{\tilde \mu}\right)^n H^{n+2} |\delta H|^{n-1}\delta H. \end{equation}

The thickness of the thin film of fluid in the no-slip region is determined from the mass conservation equation

(2.16)\begin{equation} \frac{\partial H}{\partial t} ={-}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol q. \end{equation}

These may be obtained by setting $h=0$ in the expressions relevant to the lubricated region, serving as a way to double check the correctness of the equations derived for the lubricated region. These agree with the governing equations describing classical viscous gravity currents of power-law fluids, spreading axisymmetrically (Sayag & Worster 2013).

2.3. Boundary conditions

The governing equations are supplemented by the condition of continuity of flux

(2.17)\begin{equation} [\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{n}_{\boldsymbol{L}} + \boldsymbol{q}_{\boldsymbol{l}}\boldsymbol{\cdot} \boldsymbol{n}_{\boldsymbol{L}}]^- =[\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{n}_{\boldsymbol{L}}]^+ \quad \text{at}\ \boldsymbol x = \boldsymbol{x_L} \end{equation}

and thickness

(2.18)\begin{equation} [H]^- =[H]^+ \quad \text{at}\ \boldsymbol x = \boldsymbol{x_L} \end{equation}

across the lubrication front, where $\boldsymbol x=\boldsymbol {x_L}$ is a vector describing the position of the lubrication front and $\boldsymbol {n_L}$ is the outward normal at the lubrication front.

Both fronts evolve kinematically,

(2.19a,b)\begin{equation} \boldsymbol{n}_{\boldsymbol{L}} \boldsymbol{\cdot} \dot{\boldsymbol{x}}_{\boldsymbol{L}} = \lim_{\boldsymbol x\to \boldsymbol{x}_{\boldsymbol{L}}} \frac{\boldsymbol{n}_{\boldsymbol{L}} \boldsymbol{\cdot} \boldsymbol{q}_{\boldsymbol{l}}}{h}, \quad \boldsymbol{n}_{\boldsymbol{N}} \boldsymbol{\cdot} \dot{\boldsymbol{x}}_{\boldsymbol{N}} = \lim_{\boldsymbol x\to \boldsymbol{x}_{\boldsymbol{N}}} \frac{\boldsymbol{n}_{\boldsymbol{N}} \boldsymbol{\cdot} \boldsymbol{q}}{H}. \end{equation}

The normal component of the flux of the upper layer vanishes at the leading edge of the no-slip current, so that

(2.20)\begin{equation} \boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{n}_{\boldsymbol{N}} = 0 \quad \text{at}\ \boldsymbol x = \boldsymbol{x_N}. \end{equation}

In contrast, the analogous condition

(2.21)\begin{equation} \boldsymbol{q}_{\boldsymbol{l}}\boldsymbol{\cdot} \boldsymbol{n}_{\boldsymbol{L}} = 0 \quad \text{at}\ \boldsymbol x = \boldsymbol{x_L}, \end{equation}

at the lubrication front holds if and only if $\mathcal {D}\neq 0$. In the special case of the equal density limit ($\mathcal {D}=0$), the stress singularity at the lubrication front reduces to a jump discontinuity, at which the lubrication front ends abruptly with a non-zero thickness (see Kowal & Worster 2015).

What remains to specify are the boundary conditions at the source, at which we assume a constant flux. This type of boundary condition depends on the geometry. In the two-dimensional setting, the two thin films of fluid are supplied at constant flux at the line source $x=0$,

(2.22a,b)\begin{equation} q=q_0, \quad q_l=q_{l0}, \quad\text{at}\ x=0, \end{equation}

while in the axisymmetric setting, the two fluids are supplied at constant flux at the point source $r=0$,

(2.23a,b)\begin{equation} \lim_{r\to0} 2{\rm \pi} r q=Q_0, \quad \lim_{r\to0} 2{\rm \pi} r q_l=Q_{l0}. \end{equation}

It is the dimensional differences in these boundary conditions that leads to different scaling laws for the flow in the two geometrical settings, determined in § 3.

The governing equations and boundary conditions specified here are sufficient to specify the flow completely. In the Newtonian limit ($n=n_l=1$), these reduce to the governing equations and boundary conditions obtained by Kowal & Worster (2015).

3. Similarity solutions

We lead the theoretical development by considering axisymmetric flows, and refer the reader to Appendix A for two-dimensional flows. We consider small perturbations to axisymmetric flows in Part 2.

3.1. Axisymmetric flows

At late times, a radially spreading lubricated viscous gravity current, for which $n=n_l$, tends towards a self-similar flow, for which the thicknesses of the two layers scale as

(3.1)\begin{equation} \left(h(r,t), H(r,t)\right)=\left(\frac{\rho g}{\tilde \mu}\right)^at^bQ_0^c\cdot(f(\xi),F(\xi)) \end{equation}

in terms of the similarity variable $\xi$, which we define in terms of the relationships

(3.2)$$\begin{gather} r = \left(\frac{\rho g}{\tilde \mu}\right)^\alpha t^\beta Q_0^\gamma \xi \xi_L \quad \text{for}\ 0< r< r_L, \end{gather}$$

(3.3)$$\begin{gather}r = \left(\frac{\rho g}{\tilde \mu}\right)^\alpha t^\beta Q_0^\gamma [\xi_L+(\xi-1)(\xi_N-\xi_L)] \quad \text{for}\ r_L< r< r_N, \end{gather}$$

where $\xi _L$ and $\xi _N$ specify the lubrication front $r=r_L$ and leading edge $r=r_N$, respectively. That is,

(3.4a,b)\begin{equation} r_L = \left(\frac{\rho g}{\tilde \mu}\right)^\alpha t^\beta Q_0^\gamma \xi_L \quad \text{and} \quad r_N = \left(\frac{\rho g}{\tilde \mu}\right)^\alpha t^\beta Q_0^\gamma \xi_N. \end{equation}

Such a rescaling maps the lubricated region to $\xi \in (0,1)$ and the no-slip region to $\xi \in (1,2)$. Here,

(3.5a–c)$$\begin{gather} a={-}\frac{2n}{5n+3},\quad b=\frac{n-1}{5n+3},\quad c=\frac{n+1}{5n+3}, \end{gather}$$

(3.6a–c)$$\begin{gather}\alpha =\frac{n}{5n+3},\quad \beta = \frac{2n+2}{5n+3},\quad \gamma=\frac{2n+1}{5n+3}. \end{gather}$$

These scalings are identical to those of single-layer viscous gravity currents of power-law fluids in an axisymmetric geometry (Sayag & Worster 2013). Similar similarity solutions are obtainable in the glaciological context, of isothermal ice sheets with a net accumulation (Bueler et al. 2005).

For unequal exponents, $n\neq n_l$, the front of each viscous fluid propagates at a power law of different exponent, meaning that one of the viscous layers may overtake (or ‘outstrip’) the other. There is a characteristic length scale after which this may happen. This can be found by balancing the above scalings for viscous fluids of differing fluid properties and unequal power-law exponents, giving rise to a radial length scale given by

(3.7)\begin{equation} r\sim \left(L^{N}L_l^{{-}N_l}\right)^{1/(N-N_l)}, \end{equation}

where

(3.8a,b)$$\begin{gather} L=\left(\frac{\rho g}{\tilde \mu}\right)^{n/(5n+3)} Q_0^{(2n+1)/(5n+3)},\quad L_l=\left(\frac{\rho_l g}{\tilde \mu_l}\right)^{n_l/(5n_l+3)} Q_{l0}^{(2n_l+1)/(5n_l+3)}, \end{gather}$$

(3.9a,b)$$\begin{gather}N= \frac{3n+2}{5n+3}, \quad N_l= \frac{3n_l+2}{5n_l+3}. \end{gather}$$

The exponent $1/(N-N_l)$ defining this length scale diverges in the limit of equal power-law exponents ($n=n_l$). That is, no intrinsic length scale exists when $n=n_l$, enabling a self-similar flow regime to exist in this limit.

Under these scaling laws, the governing equations in the lubricated region, $0<\xi <1$, reduce to

(3.10)\begin{gather} q={-}\frac{\left(F-f\right){}^{n+2} \left| F'\right|^{n-1}F' }{(n+2)\xi _{L}^{n}} +\frac{\mathcal{M}^n\left(F-f\right) \left(\left| \left(F-f\right) F'\right| {}^{n+1}- \left|\mathcal{D} f f'+F F'\right| {}^{n+1}\right)}{(n+1) \xi _{L}^{n}\left(\mathcal{D} f'+F'\right)}, \end{gather}

(3.11)\begin{gather} \hspace{-2.3pc} q_{l} ={-}\frac{\mathcal{M}^n}{(n+1) (n+2)} \frac{1}{\xi_{L}^n}\frac{1}{\left(\mathcal{D}f'+F'\right)^2} [\left(F-f\right)^{n+2} \left| F'\right|^{n+1}F' \nonumber\\ +\left|\mathcal{D} f f'+F F'\right|^{n+1} ((n+1)f(\mathcal{D} f'+F')-(F-f) F')], \end{gather}

(3.12)$$\begin{gather} \frac{(n-1)}{5 n+3}f - \frac{2 (n+1)}{5 n+3} \xi f'={-} \frac{\left(\xi q_{l}\right)'}{\xi \xi_{L}} \end{gather}$$

(3.13)$$\begin{gather}\frac{n-1}{5 n+3}(F-f) -\frac{2 (n+1) }{5 n+3}\xi (F' - f') ={-}\frac{\left(\xi q\right)'}{\xi \xi_{L}}. \end{gather}$$

The first term of the expression (3.10) for $q$ is a Poiseuille-like contribution originating from the gravitational spreading of the upper layer under its own weight, while the second term comprises a Couette-like viscous coupling contribution, induced by the motion of the lower layer. Similarly, the last term of the expression (3.11) for $q_l$ reflects a Poiseuille-like flow driven by the spreading of the film under its own weight and the weight of the fluid above it. The pressure gradient, driving the flow of the lower layer, depends on both surface heights in the presence of a non-zero density difference between the two layers. The first term is a Couette-like contribution induced by motion at the interface between the two layers.

In the no-slip region, $1<\xi <2$,

(3.14)\begin{equation} q={-}\frac{F^{n+2}}{(n+2)}\frac{F'}{\left(\xi_{N}-\xi_{L}\right)} \left|\frac{F'}{\xi_{N}-\xi_{L}}\right|^{n-1}, \end{equation}

and

(3.15)\begin{align} & \frac{(n-1)}{5 n+3} F -\frac{2 (n+1) }{5 n+3} \left(\xi-1+\frac{\xi_L}{\xi_N-\xi_L}\right) F' \nonumber\\ &\quad =\frac{q }{ (\xi -2) \xi_{L}-(\xi -1) \xi_{N} } + \frac{ q' }{\xi_{L}-\xi_{N} }. \end{align}

The boundary conditions in similarity coordinates reduce to

(3.16a,b)$$\begin{gather} \lim_{\xi\rightarrow0}2{\rm \pi}\xi_{L}\xi q = 1 \quad \text{and} \quad \lim_{\xi\rightarrow0}2{\rm \pi}\xi_{L}\xi q_{l} = \mathcal{Q}, \end{gather}$$

(3.17a,b)$$\begin{gather}\left[ F \right]^+_- = 0 \quad \text{and}\quad \left[ q \right]^+_- = 0 \quad \text{at}\ \xi=1, \end{gather}$$

(3.18)$$\begin{gather}q_{l} = 0 \quad \text{at}\ \xi=1 \quad\text{if}\ \mathcal{D}\neq0, \end{gather}$$

(3.19)$$\begin{gather}q = 0 \quad \text{at}\ \xi=2, \end{gather}$$

and the kinematic conditions become

(3.20a,b)\begin{equation} \frac{2n+2}{5n+3}\xi_{L} = \lim_{\xi\rightarrow1}\frac{q_{l}}{f} \quad \text{and} \quad \frac{2n+2}{5n+3}\xi_{N} = \lim_{\xi\rightarrow2}\frac{q}{F}. \end{equation}

3.2. Asymptotic solutions

The solutions involve singularities at the boundaries of the lubricated and no-slip regions. The singularities at the two fronts occur as a consequence of diverging stress at the nose of each current. These singularities occur in both geometries. A secondary singularity at the origin occurs for radially spreading flows as a consequence of a finite amount of fluid being supplied at a single point.

The presence of these singularities prompt the need for asymptotic solutions at the end points in order to obtain robust numerical solutions within the whole domain.

We present asymptotic results in the axisymmetric configuration in this section, and refer the reader to Appendix A for asymptotic solutions in the two-dimensional configuration.

3.2.1. Asymptotic solution near the fronts

An asymptotic analysis near the lubrication front in the axisymmetric geometry reveals the asymptotic solution

(3.21)\begin{equation} f \sim\left[\frac{2(n+1) (n+2) \xi_{L}}{5 n+3} \left(\frac{(2 n+1) \xi_{L}}{\mathcal{M} \mathcal{D} n}\right)^n \right]^{{1}/{(2n+1})}(1-\xi )^{n/(2 n+1)} \quad \text{as}\ \xi\to1, \end{equation}

in the lubricated region near the intrusion front and the asymptotic solution

(3.22)\begin{align} F \sim \left[\frac{ 2\left(n+1\right)\left(n+2\right) \xi_{N} }{ 5 n+3 } \left(\frac{(2n+1)(\xi_{N}-\xi_{L})}{n}\right)^{n} \right]^{{1}/{(2 n+1)}}(2-\xi )^{{n}/{(2 n+1)}} \quad \text{as}\ \xi\to2, \end{align}

near the leading edge. These satisfy the kinematic conditions and the governing equations near the two fronts.

These asymptotic solutions are shown in figure 2 in both geometrical settings (see Appendix A for asymptotic solutions in the two-dimensional geometry). A lesser agreement between the asymptotic solution and the full similarity solution is seen near the lubrication front, as forces other than the buoyancy forces associated with the lower layer are driving the flow of the lower layer, away from the lubrication front. However, the agreement becomes better as the lubrication front is approached.

Figure 2. Asymptotic solutions (3.21)–(3.22) and (A21)–(A22) (dotted curves) near the two fronts in similarity coordinates in comparison to the full similarity solution (solid curves) for the profile thicknesses $f$ and $F$ of the two layers in the two-dimensional (a) and axisymmetric (b,c) settings. Asymptotic solutions (3.24) for $n<1$ (dashed curves, b) and (3.23) for $n>1$ (dashed curves, c) depicting the singularity near the origin for the axisymmetric problem are shown near the origin. Parameters used: $\mathcal {M}=20$, $\mathcal {D}=2$, $\mathcal {Q}=0.1$ (a–c), $n=1/2$ (a,b) and $n=2$ (c).

The region of validity, in which the asymptotic solution remains valid, in $\eta$, is largest for small values of $n$, indicating that buoyancy forces associated with the lower layer are most dominant for small values of $n$. As $n$ increases, other forces become increasingly dominant.

The asymptotic solution shows that the frontal singularities become more pronounced the more shear thickening the flow, as shown illustratively in figure 3 for the asymptotic solution near the lubrication front. Gradients near the lubrication front become increasingly steep the lower the value of $n$, leading to a near-vertical front for low $n$.

Figure 3. Asymptotic solution (A21) in similarity coordinates for the lower layer in the two-dimensional geometry for various values of $n$. The frontal singularity becomes more pronounced for shear-thickening flows. Parameters used: $\mathcal {M}=20$, $\mathcal {D}=2$, $\mathcal {Q}=0.1$ and $n=1/10$, $1$, and $10$.

3.2.2. Asymptotic solution near the origin

The solutions to the governing equations become asymptotically large near the origin for radially spreading flows and the character of the singularity depends on the value of $n$. In the Newtonian case, the solutions involve a weak logarithmic singularity near the origin as seen in Huppert (1982) and Kowal & Worster (2015).

The singularity sharpens for shear-thinning flows ($n>1$), for which the slopes of the two films of fluid diverge as

(3.23)\begin{equation} (F,f) \sim (\mathcal{A},\mathcal{a}) (1+v \xi ^{1-{1}/{n}})^{{n}/{(2 n+2)}}\quad \text{as}\ \xi\to0 \end{equation}

to leading order, where $\mathcal {A}$, $\mathcal {a}$ and $v$ are constants.

For shear-thickening flows ($n<1$), the thicknesses of the two layers diverge as

(3.24)\begin{equation} (F, f) \sim (\mathcal{B},\mathcal{b}) \xi^{-({(1-n)}/{2 (n+1)})} \quad \text{as}\ \xi\to0, \end{equation}

where $\mathcal {B}$ and $\mathcal {b}$ are constants that diverge as $(1-n)^{-n/(2n+2)}$ as $n\to 1^-$. The coefficients describing these asymptotic solutions satisfy the system of equations outlined in Appendix B. There is no closed-form solution for this nonlinear system, in general.

The asymptotic solutions (3.23) for $n>1$ and (3.24) for $n<1$ are plotted against the full numerical solution in figure 2.

3.3. Numerical solutions

We use the approach adopted by Kowal & Worster (2015) to integrate the governing equations numerically. In particular, the governing equations have been solved over the interval $(0,\ 1-\delta )$ for the lubricated region, and $(1,\ 2-\delta )$ for the no-slip region, where $\delta \ll 1$ is a small parameter, in the two-dimensional configuration. In the axisymmetric configuration, the equations modelling the lubricated region have instead been solved over the smaller interval $(\varDelta,\ 1-\delta )$, owing to the weak singularity at the origin, where $\varDelta \ll 1$.

The asymptotic solutions (3.21)–(3.22) and (A21)–(A22) have been used in the small intervals $(1-\delta, 1]$ and $(2-\delta, 2]$, and as boundary conditions specifying the layer thickness and slope at $\xi =1-\delta$ and $\xi =2-\delta$ in the axisymmetric configuration and at $\eta =1-\delta$ and $\eta =2-\delta$ in the two-dimensional configuration. We do so by shooting backwards from the leading edge $\xi =2$ in the axisymmetric geometry and $\eta =2$ in the two-dimensional geometry, with $\xi _L$ and $\xi _N$ (respectively, $\eta _L$ and $\eta _N$) being shooting parameters. That is, $\xi _L$ and $\xi _N$ (respectively, $\eta _L$ and $\eta _N$) are found as part of the numerical solution. The numerical integration has been implemented in Mathematica and the numerical solutions have been validated by confirming agreement with those of Kowal & Worster (2015) for $n=1$, and with the asymptotic solutions (3.21)–(3.22) and (A21)–(A22) for a range of values of $\delta$ and $\varDelta$.

Figure 2 displays both the full numerical solution over the interval $[0,1-\delta ]$, and the frontal asymptotic solutions over the interval $[0.7,1]$, which extends beyond its region of validity, for completeness, in the two geometries.

4. Flow regimes

The similarity solutions derived in the prior sections depend upon four dimensionless parameters, namely: the consistency ratio $\mathcal {M}$, the density difference $\mathcal {D}$, the flux ratio $\mathcal {Q}$ and the power-law exponent $n$. The behaviour of the flow, and the mechanisms driving the flow, are sensitive to the values of these parameters. In what follows, we assume that $\mathcal {D}$ is a fixed value bounded away from zero and we investigate the emerging flow regimes in $(\mathcal {Q}, \mathcal {M})$ parameter space.

Akin to the results of Kowal & Worster (2015), the flow can be categorised into four distinct flow regimes, which we label using roman numerals as described in table 1. Regime diagrams depicting the four regimes in the two-dimensional and axisymmetric geometries are shown in figure 4 for various values of $n$, which we refer to in the discussion, below. The curves depicting the boundaries between the regimes are obtained numerically based on the criteria summarised in table 1, as there are no explicit formulae for these in general. The results of this section reduce to those obtained by Kowal & Worster (2015) in the limit of $n=1$, though the correspondence is less clear for the boundary between regimes II and III. It should be noted that the definitions of the regime boundaries in Kowal & Worster (2015) are indicative order-of-magnitude estimates only, rather than quantitative as is the case here.

Figure 4. Regime diagram in $(\mathcal {Q},\mathcal {M})$-space for $\mathcal {D}=1$ and $n=1/2,1,2$ in the two-dimensional (a) and axisymmetric (b) settings. The red curves are contours of $\zeta _L=\zeta _N-\zeta _L$ and the blue curves are contours of $\bar f=\bar F-\bar f$, separating regimes I–IV as outlined in table 1.

Table 1. Summary of parameter regimes.

We lead the discussion by example of axisymmetric flows and present abbreviated results for two-dimensional flows in Appendix A. In what follows, we summarise the leading-order equations for both geometries simultaneously. For brevity, we denote $\eta _L$ and $\eta _N$, for two-dimensional flows, and $\xi _L$ and $\xi _N$, for axisymmetric flows, collectively by $\zeta _L$ and $\zeta _N$, respectively.

4.1. Regime I: classical spreading of the overlying layer

Regime I is one in which the upper layer spreads gravitationally under its own weight, largely unaffected by the underlying layer, and the motion generated at the interface induces a shear flow within the underlying layer. Such a balance of forces results in the following leading-order contributions

(4.1)$$\begin{gather} q_l={-}\frac{\mathcal{M}^n}{2\zeta_L^{n}} f^2 F^n \left| F'\right|^{n-1}F', \end{gather}$$

(4.2)$$\begin{gather}q={-}\frac{1}{(n+2) \zeta_L^n}F^{n+2} F' \left| F'\right| ^{n-1} \end{gather}$$

to the fluxes of the two layers of fluid. These can be deduced by reference to the physical description of the terms comprising the fluxes of the two layers, described in the paragraph following (3.11), and by noting that $F\gg f$. In this case, the lubricated region is relatively thin with a short extent, and its effect on the upper layer is small. The flow of the upper layer is similar to that of a classical viscous gravity current of a single layer of fluid spreading over a rigid, no-slip surface.

For the purposes of the illustration of figure 4, we define this regime quantitatively as one in which

(4.3a,b)\begin{equation} \bar f < \bar F-\bar f, \quad \zeta_L<\zeta_N-\zeta_L, \end{equation}

where $\bar f$ and $\bar F$ are the average lower and upper layer thicknesses, respectively, over the lubricated region. These are defined by

(4.4)\begin{equation} \bar f = \int_0^{1} f(\eta)\,{\rm d}\eta \end{equation}

in the two-dimensional configuration and

(4.5)\begin{equation} \bar f = \int_0^{1} 2 f(\xi)\xi\,{\rm d}\xi \end{equation}

in the axisymmetric configuration, and similarly for $\bar F$.

The remaining part of parameter space is covered by three further parameter regimes, defined in terms of similar inequalities discussed in the subsequent subsections. To compute the boundaries of these regions of parameter space, we replace the two inequalities of (4.3a,b) by corresponding equalities and search for the parameter $\mathcal {M}$ for which the first equality is satisfied (as a function of the remaining parameters) to obtain one of the curves of figure 4 defining the boundary, and repeat for the second equality to obtain the other curve. Regions in between these curves correspond to the four parameter regimes discussed in this section.

Flows corresponding to this regime occur for low values of $\mathcal {Q}$ and intermediate values of $\mathcal {M}$ as shown in the regime diagram of figure 4. It can be seen from the regime diagram that the region of parameter space corresponding to regime I reduces with increasing values of $n$ – both the upper and lower boundaries of this region contract.

Typical profile thicknesses of axisymmetric lubricated viscous gravity currents under regime I are depicted in figure 7. As the lubricating layer becomes very thin in this regime, we have taken $\mathcal {M}=1$ for visibility. The lubricating layer thins and extends and the overlying layer flattens as $n$ increases. The shape of the overlying current is largely unchanged by the presence of the lower layer. The average thickness and frontal positions are displayed in figure 5, transitioning from regime I for small $\mathcal {Q}$ towards regimes III and II for large $\mathcal {Q}$, for various values of $n$, and in figure 6 as a function of $n$. Increasing $n$ promotes the thinning and elongating of both layers. This may be a consequence of shear-thinning currents being thinner and longer than shear-thickening currents, in general.

Figure 5. (a) The average surface heights of the upper and lower layers (solid and dashed curves, respectively) in similarity coordinates in the axisymmetric setting vs $\mathcal {Q}$. (b) The frontal positions $\xi =\xi _N$ and $\xi _L$ (solid and dashed curves, respectively) in the axisymmetric setting vs $\mathcal {Q}$. Parameters used: $\mathcal {M} = 50$, $\mathcal {D} =1$ and $n=1/2,1,2$.

Figure 6. (a) The average upper surface height of the upper and lower layers (solid and dashed curves, respectively) in similarity coordinates in the axisymmetric setting vs $n$. (b) The frontal positions $\xi =\xi _N$ and $\xi _L$ (solid and dashed curves, respectively) in the axisymmetric setting vs $n$. Parameters used: $\mathcal {M} = 5$, $\mathcal {D} = 1$, $\mathcal {Q}= 0.1$.

Figure 7. Profile thicknesses $f$ and $F$ of axisymmetric lubricated viscous gravity currents falling under regime I for $\mathcal {M}=1$, $\mathcal {D}=1$, $\mathcal {Q}=0.05$ and $n=1/2$ (a), $n=1$ (b) and $n=2$ (c). The thicknesses of the overlying and underlying layers are given by the solid and dashed curves, respectively.

4.2. Regime II: flow driven by the underlying layer

Under regime II, the motion is driven by the lower layer, which is spreading gravitationally under its own weight and the weight of the fluid above it. The flow of the overlying layer is driven by interfacial stresses between the two layers. Effectively, the lower layer pulls the upper layer along with it, which behaves like a plug-like flow. Such a balance results in the leading-order contributions

(4.6)$$\begin{gather} q_l={-}\frac{\mathcal{M}^n(\mathcal{D}+1)^n }{(n+2) \zeta_L^n}f^{n+2} \left| f'\right| ^{n-1}f', \end{gather}$$

(4.7)$$\begin{gather}q={-}\frac{\mathcal{M}^n(\mathcal{D}+1)^n }{(n+1) \zeta_L^n}(F-f) f^{n+1} \left| f'\right| ^{n-1}f', \end{gather}$$

to the fluxes of fluid in the two layers. Both of these involve gradients in the lower layer thickness only. These can be deduced by reference to the physical description of the terms comprising the fluxes of the two layers, described in the paragraph following (3.11), and by noting that $f\gg F-f$.

A consequence is that gradients in the upper layer surface height are not influential to leading order and slope reversals in the upper surface height can occur, as shown in figure 8. Figure 8 shows typical profile thicknesses for $n<1$, $n=1$ and $n>1$. A slope reversal near the lubrication front occurs in all three cases, although in reduced form for $n>1$. Despite the reversed slope of the upper surface height near the lubrication front, which gives rise to a negative contribution to the flux of fluid in the upper layer, the dominant contribution to the flow of the upper layer originates from the motion of the lower layer, which is downstream.

Figure 8. Profile thicknesses $f$ and $F$ of axisymmetric lubricated viscous gravity currents falling under regime II for $\mathcal {M}=100$, $\mathcal {D}=1$, $\mathcal {Q}=2$ and $n=1/2$ (a), $n=1$ (b) and $n=2$ (c). The thicknesses of the overlying and underlying layers are given by the solid and dashed curves, respectively.

These slope reversals are most prominent for $n<1$. Increasing $n$ thins and elongates both layers, as seen for the other regimes, and as seen in the plots of the average thickness and frontal positions in figure 5. As mentioned earlier, this may be a consequence of shear-thinning currents being thinner and longer than shear-thickening currents, in general. A prominent feature that occurs for increasing values of $n$ is the increasing thinning of the overlying layer, which behaves like a viscous coating over the underlying layer of fluid, which is driving the flow.

We define this regime quantitatively as one in which

(4.8a,b)\begin{equation} \bar f > \bar F-\bar f, \quad \zeta_L>\zeta_N-\zeta_L. \end{equation}

As seen in the regime diagrams of figure 4, the flows of regime II occur for large values of the flux ratio $\mathcal {Q}$ and moderate values of the consistency ratio $\mathcal {M}$. This means that the overlying viscous coating may be either more or less viscous than the driving viscous gravity current. This regime occurs more frequently for small values of $n$. Increasing values of $n$ lead to a diminishing region of parameter space describing region II.

4.3. Regime III: thin lubricating layer

Regime III characterises itself by a thin lubricating layer, or a shear layer, to which most of the shear of the system is confined to. Despite its thinness, its presence significantly alters the dynamics of the overlying layer, which, in this regime, behaves like a plug flow. The thin lubricating shear layer deforms under hydrostatic pressure gradients associated with the upper layer. Even though these gradients are weak, the consistency ratio $\mathcal {M}$ is large enough that the underlying, much less viscous layer is driven by these upper layer gradients.

This regime is unique in that the flow of neither layer is similar to that of a single-layer classical viscous gravity current consisting of a fluid with the properties of one or the other layer. The flow is significantly altered in both layers.

Such a balance results in the leading-order contributions

(4.9)$$\begin{gather} q_l={-}\frac{\mathcal{M}^n}{2 \zeta_L^n}f^2 F^n \left| F'\right| ^{n-1}F', \end{gather}$$

(4.10)$$\begin{gather}q={-}\frac{\mathcal{M}^n }{\zeta_L^n} f F^{n+1} \left| F'\right| ^{n-1}F', \end{gather}$$

to the fluxes of the two layers of fluid. Both of these involve gradients in the upper surface height of the overlying layer of fluid. These can be deduced by reference to the physical description of the terms comprising the fluxes of the two layers, described in the paragraph following (3.11), and by noting that $f\ll F-f$.

We define this regime quantitatively as one in which

(4.11a,b)\begin{equation} \bar f < \bar F-\bar f, \quad \zeta_L>\zeta_N-\zeta_L. \end{equation}

As seen in the regime diagrams of figure 4, the flows of regime III occur for large values of the consistency ratio $\mathcal {M}$ and moderate values of the flux ratio $\mathcal {Q}$. Regime III becomes increasingly prominent for increasing values of $n$.

Typical profile thicknesses of the two layers of fluid are shown in figure 9 for $n<1$, $n=1$ and $n>1$. There is a visible flattening of the upper surface of the overlying layer of fluid, which becomes more apparent for $n>1$. There is a departure to this observation near the origin, at which there is a weak singularity. A similar observation can be made about the lubricating layer, apart from the vicinity of the lubrication front, at which the underlying layer spreads under its own weight owing to large thickness gradients of the underlying layer in this region.

Figure 9. Profile thicknesses $f$ and $F$ of axisymmetric lubricated viscous gravity currents falling under regime III for $\mathcal {M}=100$, $\mathcal {D}=1$, $\mathcal {Q}=0.1$ and $n=1/2$ (a), $n=1$ (b) and $n=2$ (c). The thicknesses of the overlying and underlying layers are given by the solid and dashed curves, respectively.

Apart from this flattening, increasing values of $n$ thin and elongate both layers of fluid, as seen in figure 9 and the plots of the average thickness and frontal positions in figure 5, corresponding to regime III. As mentioned earlier, this may follow from shear-thinning currents being thinner and longer than shear-thickening currents, in general.

4.4. Regime IV: low-viscosity thin film coating a more viscous current

Regime IV consists of flows driven mainly by the gravitational spreading of the underlying layer, which spreads mainly under its own weight with negligible influence of the overlying layer. The overlying layer behaves as a thin coating film of fluid, which exerts negligible shear stress on the underlying layer and provides little resistance to the flow. Such a balance results in the leading-order contributions

(4.12)$$\begin{gather} q_l={-}\frac{\mathcal{M}^n(\mathcal{D}+1)^n }{(n+2) \zeta_L^n}f^{n+2} \left| f'\right| ^{n-1}f', \end{gather}$$

(4.13)$$\begin{gather}q={-}\frac{1}{(n+2) \zeta_L^n} (F-f)^{n+2} \left| F'\right| ^{n-1}F', \end{gather}$$

to the fluxes of the two layers of fluid. These can be deduced by reference to the physical description of the terms comprising the fluxes of the two layers, described in the paragraph following (3.11), and by noting that $f\gg F-f$.

In contrast to regime II, flows under regime IV do not involve upper surface slope reversals near the lubrication front, as the flow of the upper layer is driven by hydrostatic pressure gradients established within the overlying layer, rather than the underlying layer. The coupling between the two layers is purely a consequence of geometry, rather than dynamics.

Both this regime and regime III involve thin films of fluid – a thin coating layer, overlying the driving viscous gravity current in the case of regime IV, and a thin lubricating layer, underlying the driving viscous layer in the case of regime III. In contrast to regime III, the thin coating film exerts negligible shear stress on the driving layer as it is a coating covering the driving layer from above. A simple change in location of the film produces significant dynamical changes. That is, the flow changes significantly when the thin film is a lubricating film, underlying the driving layer from below as in regime III.

We define this regime quantitatively as one in which

(4.14a,b)\begin{equation} \bar f > \bar F-\bar f, \quad \zeta_L<\zeta_N-\zeta_L. \end{equation}

Flows under regime IV occur for low consistency ratios and moderate flux ratios, as shown in figure 4. Such flows become increasingly common for increasing values of $n$ – similarly to regime III.

Illustratively, such flows involve a low-viscosity thin film of fluid coating a more viscous layer. Plots of typical profile thicknesses of the two layers of fluid are shown in figure 10 for various values of $n$. Slope reversals do not occur near the lubrication front and the overlying, less viscous layer freely runs out ahead of the front. The thin coating film becomes less noticeable, and runs ahead further, for increasing values of $n$.

Figure 10. Profile thicknesses $f$ and $F$ of axisymmetric lubricated viscous gravity currents falling under regime IV for $\mathcal {M}=0.005$, $\mathcal {D}=1$, $\mathcal {Q}=0.4$ and $n=1/2$ (a), $n=1$ (b) and $n=2$ (c). The thicknesses of the overlying and underlying layers are given by the solid and dashed curves, respectively.

5. Conclusions

We investigated the effect of shear-thinning and shear-thickening rheology on the flow of lubricated viscous gravity currents in self-similar form, both in a two-dimensional and in an axisymmetric geometry. Similarity solutions exist in both geometries as long as the power-law exponents, describing the rheology of the viscous fluids, are identical in both layers. Four distinct flow regimes emerge depending on the consistency ratio $\mathcal {M}$, the density difference $\mathcal {D}$ and the flux ratio $\mathcal {Q}$. The flow regimes involve either (i) classical spreading of the overlying layer, largely unaffected by the underlying layer, in which the motion generated at the interface induces a shear flow within the underlying layer (regime I); (ii) a flow driven by the underlying layer, in which the interfacial stresses between the two layers drive a plug-like flow of the overlying layer (regime II); (iii) the flow of a thin lubricating layer, driven by hydrostatic pressure gradients associated with the upper layer, which, unlike its thinness, generates a plug flow within the overlying layer (regime III); and (iv) the flow of a low-viscosity thin film coating a more viscous current from above (regime IV).

Similar to purely Newtonian flows, a stress singularity is found at the nose of both layers of fluid for all power-law exponents. The frontal singularity becomes more pronounced the more shear thickening the rheology. Shear thickening leads to increasingly steep gradients near the two fronts, leading to the fronts becoming near vertical.

The flow dynamics is similar for shear thinning ($n>1$), purely Newtonian ($n=1$), and shear-thickening flows ($n<1$) within each regime. However, the scalings, frontal positions, thicknesses, the structure of the singularity at the nose of both layers and the regime boundaries change intrinsically depending upon the rheology. Shear thinning, in general, promotes the elongating and thinning of both layers except when the overlying layer forms a low-viscosity thin film coating a more viscous current. This may be a consequence of shear-thinning currents being thinner and longer than shear-thickening currents, in general. As such, regimes III and IV become increasingly prominent for increasing values of the power-law exponent. The former regime, in which the main current is shear thinning and lubricated by a much less viscous thin layer of fluid from below, is most relevant to the flow of lubricated ice sheets and fast-flowing ice streams.

Figure 11. Profile thicknesses $f$ and $F$ of two-dimensional lubricated viscous gravity currents falling under regime I for $\mathcal {M}=1$, $\mathcal {D}=1$, $\mathcal {Q}=0.05$ and $n=1/2$ (a), $n=1$ (b) and $n=2$ (c). The thicknesses of the overlying and underlying layers are given by the solid and dashed curves, respectively.

Figure 12. Profile thicknesses $f$ and $F$ of two-dimensional lubricated viscous gravity currents falling under regime II for $\mathcal {M}=100$, $\mathcal {D}=1$, $\mathcal {Q}=2$ and $n=1/2$ (a), $n=1$ (b) and $n=2$ (c). The thicknesses of the overlying and underlying layers are given by the solid and dashed curves, respectively.

Figure 13. Profile thicknesses $f$ and $F$ of two-dimensional lubricated viscous gravity currents falling under regime III for $\mathcal {M}=100$, $\mathcal {D}=1$, $\mathcal {Q}=0.1$ and $n=1/2$ (a), $n=1$ (b) and $n=2$ (c). The thicknesses of the overlying and underlying layers are given by the solid and dashed curves, respectively.

Figure 14. Profile thicknesses $f$ and $F$ of two-dimensional lubricated viscous gravity currents falling under regime IV for $\mathcal {M}=0.005$, $\mathcal {D}=1$, $\mathcal {Q}=0.4$ and $n=1/2$ (a), $n=1$ (b) and $n=2$ (c). The thicknesses of the overlying and underlying layers are given by the solid and dashed curves, respectively.