Lubricated axisymmetric gravity currents of power-law fluids

The motion of glaciers over their bedrock or drops of fluid along a solid surface can become unstable when these substrates are lubricated. Previous studies modeled such systems as coupled gravity currents (GCs) consisting of one fluid that lubricates the flow of another fluid, and having two propagating fronts. When both fluid are Newtonian and discharged at constant flux, global similarity solutions were found. However, when the top fluid is strain-rate softening experiments have shown that each fluid front evolved with a different exponent. Here we explore theoretically and numerically such lubricated GCs in a model that describes the axisymmetric spreading of a power-law fluid on top of a Newtonian fluid, where the discharge of both fluids is power law in time. We find that the model admits similarity solutions only in specific cases, including the purely Newtonian case, for a certain discharge exponent, at asymptotic limits of the fluids viscosity ratio, and at the vicinity of the fluid fronts. Generally, each fluid front has a power-law time evolution with a similar exponent as a non-lubricated GC of the corresponding fluid, and intercepts that depend on both fluid properties. Consequently, we identify two mechanisms by which the inner lubricating fluid front outstrips the outer fluid front. Many aspects of our theory are found consistent with recent laboratory experiments. Discrepancies suggest that hydrofracturing or wall slip may be important near the fronts. This theory may help to understand the dynamics of various systems, including surges and ice streams.


Introduction
Gravity-driven flows of one fluid over another can involve complex interactions between the two fluids, which can lead to a rich dynamical behavior. Such flows occur in a wide range of natural and human-made systems, as in lava flow over less viscous lava (Griffiths 2000;Balmforth et al. 2000), spreading of the lithosphere over the mid-mantle boundary (Lister & Kerr 1989;Dauck et al. 2019), ice flow over an ocean (DeConto & Pollard 2016;Kivelson et al. 2000) and over bedrock consisting of sediments and water (Stokes et al. 2007;Fowler 1987), flows in permeable rocks (Woods & Mason 2000), liquid drops deforming over lubricated surfaces (Daniel et al. 2017), and droplets motion on liquid-infused surfaces (Keiser et al. 2017).
The flow of GCs in circular geometry has been studied with a range of boundary conditions. In the absence of a lubricating layer, a common boundary condition along the base of a sole GC is no slip. Such GCs of Newtonian fluids that are discharged at a rate proportional to , where is time and is a non-negative scalar, has a similarity solution in which the front evolves proportionally to (3 +1)/8 (Huppert 1982). Similar axisymmetric gravity currents of power-law fluids having exponent , where = 1 represents a Newtonian fluid and > 1 represents a strain-rate softening fluid, also have similarity solutions in which the front propagation is proportional to [ (2 +1)+1]/(5 +3) (Sayag & Worster 2013).
On the other extreme, the presence of a lower fluid layer can significantly reduce friction at the base of the top fluid, resulting in extensionally dominated GCs. This is the case, for example, for ice shelves, which deform over the relatively inviscid oceans with weak friction along their interface. The late-time front evolution of such axisymmetric GCs of Newtonian fluids is proportional to and is believed to be stable (Pegler & Worster 2012). However, when the top fluid is strain-rate softening, an initially axisymmetric front can destabilise and develop fingering patterns that consist of tongues separated by rifts (Sayag & Worster 2019).
In the more general case friction along the base of GCs can vary spatiotemporally, as their stress field evolves. For example, the interface of an ice sheet with its underlying bed rock can include distributed melt water and sediments, which impose nonuniform and time-dependent friction along the ice base, and evolve spatiotemporally under the stresses imposed by the ice layer (Schoof & Hewitt 2013;Fowler 1981). Consequently, the coupled ice-lubricant system may evolve various flow patterns such as ice streams (Stokes et al. 2007;Kyrke-Smith et al. 2013) and glacier surges (Fowler 1987). Such flows were also modelled as two coupled GCs of Newtonian fluids spreading axisymmetrically one on top of the other (Kowal & Worster 2015). The early stage of these flows follows a self-similar evolution, in which the fronts of the two fluids evolve like 1/2 , as in non-lubricated (no-slip) GCs (Huppert 1982), but they can have a radially non-monotonic thickness. Furthermore, laboratory experiments were found to be consistent with the similarity solutions after an initial transient state, but at a later stage they became unstable and developed fingering patterns (Kowal & Worster 2015). It has been suggested that such instabilities appear when the jump in hydrostatic pressure gradient across the lubrication front is negative (Kowal & Worster 2019a,b).
Despite the wide range of natural lubricated GCs that involve non-Newtonian fluids, the radial flow of a non-Newtonian fluid over a lubricated layer of Newtonian fluid has just recently been explored experimentally (Kumar et al. 2021). Motivated by glacier flow over lubricated bedrock, the experimental setup consisted of a GC of a strain-rate-softening fluid (Xanthan-gum solution) that has a power-law viscous deformation similar to ice (Glen 1952), and a lubricating GC of a less viscous Newtonian fluid (diluted sugar solution). The pattern of both fluids in those constant-flux ( = 1) experiments remained axisymmetric throughout the flow (Figure 1), in contrast to the fingering pattern that emerged in the purely Newtonian experiments (Kowal & Worster 2015), as long as the flux ratio of the lubricating fluid to the non-Newtonian fluid was lower than ∼ 0.06. The fronts of the two fluids appeared to have a power-law time evolution with different exponents. Particularly, the front of the top non-Newtonian fluid evolved with the same exponent (2 + 2)/(5 + 3) as a non-lubricated power-law fluid (Sayag & Worster 2013), whereas the front of the Newtonian lubricating fluid evolved with an exponent 1/2 similar to a Newtonian non-lubricated GC (Huppert 1982). Despite the similarity of the exponents, the fronts of those lubricated GCs evolved faster than the corresponding non-lubricated GCs owing to larger intercepts. In addition, in contrast with the monotonically declining thickness of non-lubricated GCs, the thickness of the lubricated, non-Newtonian fluid was found to be nearly uniform in the lubricated part of the flow, while that of the lubricating fluid was non-monotonic with localised spikes.
Following up the experimental study of Kumar et al. (2021), here we develop the theory for lubricated axisymmetric GCs of power-law fluids, and explore its major consequences. Specifically, we develop a mathematical model for axisymmetric flow of a viscous gravity current of a power-law fluid lubricated by a Newtonian gravity current, considering a general input flux of the form , and show that in the general case the flow has no global similarity solution of the 1st kind ( §2). We then describe the numerical solver ( §3) and investigate / = 1 / = 1.13 / = 1.33 several special cases that have similarity solutions ( §4). We then explore the possibility that the inner lubricating front outstrip the outer front ( §5). Finally, we investigate the case of constant flux discharge ( §6), and compare our theoretical predictions with the laboratory experiments of Kumar et al. (2021) ( §7).

Mathematical Model
Consider a power-law fluid having viscosity and density that spreads axisymmetrically under its own weight over a horizontal rigid surface. Simultaneously, a lubricating film of Newtonian fluid of viscosity ℓ and density ℓ spreads axisymmetrically over the substrate and below the power-law fluid ( Figure 2). Both fluids are discharged at the origin of a cylindrical coordinate system in which is the radial coordinate and the vertical coordinate. The upper and lower fluid fronts are denoted by ( ) and ( ) respectively, and the corresponding fluid thicknesses are denoted by ( , ) − ℎ( , ) and ℎ( , ) respectively, where is the time variable. Beginning the discharge of the lubricating fluid at a time delay with respect to the upper fluid, creates two regions in the flow: an inner, lubricated region ( ) in which the power-law fluid flows at a finite velocity along the interface with the lubricating fluid, and an outer, non-lubricated region ( < ) in which the power-law fluid meets the substrate and has zero velocity along it.
Assuming that the radial extent of the flow is much greater than its thickness in both fluid layers, we can apply the lubrication approximation, which implies that the flow is primarily horizontal (a GC) and that the dominant strain rate in each fluid layer is ( , , ) is the radial velocity component. Therefore, the axisymmetric Cauchy equations for each fluid layer, simplified for flows of low Reynolds numbers and lubrication approximations are to leading order, where is the vertical velocity component, is the pressure field, and is the gravitational acceleration. The pressure is determined hydrostatically in equation 2.1b because we assume a large Bond number, implying that the effects of surface tension are negligible compared to gravity. We note that equations 2.1 satisfy the lubrication approximation, so this model does not accurately describe the problem at early times, before the fluid radial extent is sufficiently greater than its thickness.
The viscosity of the power-law fluid is determined by where e is the strain-rate tensor, is the consistency coefficient, and is the power-law exponent, which determines the fluid's response to stress. Specifically, = 1 represents a Newtonian fluid with dynamical viscosity , > 1 a shear-thinning fluid, and < 1 a shear-thickening fluid. In the lubrication limit e : e ≈ ( / ) 2 to leading order, so that the viscosity of the power-law fluid simplifies to We assume that the total volumes of the power-law fluid and the Newtonian fluid evolve following a power law in time given by respectively, where is a constant exponent, and and ℓ are constant coefficients representing, for instance, the discharge flux of each fluid when = 1.

The Non-Lubricated Region
In the non-lubricated region, the power-law fluid meets the substrate and the lubricating fluid is absent. Integrating 2.1b, the pressure distribution is where 0 is the ambient pressure over the top free surface of the power-law fluid. Integration of the radial force balance (2.1a) across the depth of the fluid layer together with no-slip boundary conditions along the substrate and no stress along the free surface gives the radial velocity field (2.7) Similar integration of the continuity equation (2.1c), accounting for a free surface at the upper boundary = ( , ), and assuming no normal flow through the substrate, gives the Reynolds equation which should satisfy the boundary conditions where the local flux is determined using eq. 2.7 to give (2.10)

The Lubricated Region
As in the non-lubricated region and assuming continuity of pressure at the fluid-fluid interface = ℎ, the pressure in the lubricated region is determined hydrostatically by The radial force balances simplify to 2 1 −1 together with the boundary conditions, which represent, respectively, no-slip along the solid substrate, continuous flow velocity and shear stress at the fluid-fluid interface, and no shear stress along the free surface of the power-law fluid. Integrating the radial force balance (equations 2.12) across the thickness of each fluid layer, and using the boundary conditions (equations 2.13) we get the radial velocity fields in each fluid layer (2.14b) The Reynolds equations corresponding to each fluid layer have a similar form as that of the non-lubricated region and should satisfy the boundary conditions 16a) and where equation 2.16b signifies flux and height continuity across . The local fluxes result from integrating equations 2.14 to get (2.17b) We note that without a lubricant, this set of PDEs is reduced to the non-lubricated region, and that this PDE set is consistent with the non-lubricated power-law GC model (Sayag & Worster 2013), which converges to the Newtonian non-lubricated GC when = 1 and the fluid dynamic viscosity is = (Huppert 1982). In addition, the full PDE set, for the Newtonian case ( = 1) and constant source fluxes ( = 1) is consistent with the model for lubricated, Newtonian GCs (Kowal & Worster 2015).

Dimensionless Equations
We non-dimensionalize equations 2.4, 2.8-2.10 and 2.15-2.17 with the following time, length and height scales ≡ˆ, (2.18a) where hats denote dimensionless quantities. The resulted dimensionless model, dropping hats, in the non-lubricated region is (2.20) and in the lubricated region 0 it is and where the boundary conditions are (2.24) represent respectively the discharge flux ratio, the relative density difference of the lubricating and power-law fluids, the power-law fluid exponent, and the dynamic viscosity ratio. The viscosity ratio 2.24(iv) implies that the timescale in (2.18) is scaled out of the equations in two special cases = 1 and = 5. This implies that asymptotically in time ( / 1) our PDE set admits a global similarity solution of the first kind in these two special cases, whereas for any other value of or , including the constant flux case ( = 1), such global similarity solutions do not exist. Nevertheless, as we show in the section that follows ( §4), we find several additional asymptotic limits in which part of the flow evolves in a self-similar manner.
In the general case ≠ 1 and ≠ 5 there are enough relations to determine the time scale (2.18) independently of the height and radial scales. This is primarily because the flux in the top fluid layer (2.22a) consists of two contributions that are not necessarily of the same scale. Requiring that the scales of those two contributions balance M ( ) = N leads to an independent constraint for the time scale , which upon substitution in (2.18) leads to a set of independent scales .
(2.25c) Therefore, in the more general case ( ≠ 1 and ≠ 5) there is no global similarity solution of the 1st kind, and the above scales may represent the time, radius and thickness.

Numerical Solution
We solve the dimensionless model (eqs. 2.19 -2.24) numerically explicitly using the Matlab PDEPE solver, with an open-ended, non-uniform, and time-dependent adaptive spatial mesh. We find that lower spatial resolution can be used in most of the domain if the mesh is logarithmically spaced, and is denser around the fluid fronts, where the fluid heights drop sharply and the slopes / and ℎ/ become singular. Therefore, we use a non-uniform spatial mesh that consists of two Gaussian distributions centered at the fronts and , and uniform distributions between the origin and , the two fronts, and beyond ( Figure 3). The instantaneous front positions ( ) and ( ) are determined where ( , ) 10 −6 and ℎ( , ) 10 −6 respectively. Since the fronts evolve, we enlarge the flow domain progressively and re-mesh the computational grid points every several time steps in order to keep the higher mesh density centered at the front positions. For the re-mesh we predict the front positions in the next time step using a discretization of the front evolution equations (3.1) The size of the time step needed in order to solve the PDE is determined by the solver, and we scale with , so thatˆ= 1. We validated our numerical solver using several known asymptotic solutions. In the limit Q = 0 (no lubricating fluid) our model converges to a GC that propagates under no-slip condition along the substrate, which has a similarity solution ( ) ∝ [ (2 +1)+1]/(5 +3) (Huppert 1982; Sayag & Worster 2013). We find that our numerical solutions for the fluid heights and for the leading front are consistent with those theoretical predictions (Appendix A). Particularly, discrepancy between the predicted and computed exponents is less than 5 · 10 −3 and can be minimized further with a higher spatial resolution. In the limit = 1 and = 1 our model converges to a purely Newtonian lubricated GC released at constant flux, which has a similarity solution , ∝ 1/2 (Kowal & Worster 2015). We find that our numerical solutions are consistent with both the front and thickness predictions for a wide range of M and Q values (Appendix A).

Similarity solutions
The model described by eqs. 2.19-2.24 admits similarity solutions in several special cases. Specifically, global similarity solutions exist when the top fluid is also Newtonian ( = 1), or when the mass-discharge exponent is = 5. In addition, part of the flow evolves in a selfsimilar manner in the asymptotic limits of the viscosity ratio M . Finally, for any parameter combination the solutions at the vicinity of each front evolve in a self-similar manner.
4.1. The similarity solution for Newtonian fluids = 1 When the top fluid is Newtonian ( = 1) the viscosity ratio M is independent of , and the resulted purely Newtonian coupled flow admits a global similarity solution. This case, restricted to a constant flux ( = 1), has been thoroughly explored (Kowal & Worster 2015). For the general case of 0, the PDE set can be reduced to an ODE set with a similarity variable and a solution of the form where and are dimensionless functions of . Therefore, the fronts of both the lubricated and lubricating fluids evolve like where and are numerical coefficients of order 1. These solutions have identical structure as the classical solution of Newtonian GCs (Huppert 1982). Upon substitution, the Reynolds equations (2.8 and 2.15) become, for the non-lubricated region ( ) where prime denotes a derivative with respect to , and for the lubricated region (0 ) The corresponding boundary conditions (2.23) become where the last condition implies convergence to a similarity solution long time after the initiation of the fluid discharge, . For = 1 the model converges to that of Kowal & Worster (2015), and the similarity solutions we predict for general are consistent with our full numerical solution (Figures 4,5).

The similarity solution for mass-discharge exponent = 5
When the top fluid is non-Newtonian ≠ 1 the viscosity ratio M becomes independent of the time scale for a discharge exponent = 5, and the resulted flow also admits a similarity solution with a similarity variable where and are dimensionless functions of . Therefore, the fronts of both the lubricated and lubricating fluids evolve like where and are numerical coefficients of order 1. Upon substitution, the Reynolds equations (2.8 and 2.15) become, for the non-lubricated region ( ) and for the lubricated region (0 ) The corresponding boundary conditions (2.23) become (4.12c) Figure 4 shows numerical solution of the full PDE set for = 5 for varying power-law exponent values. We find that the numerical solutions are consistent with the predicted similarity solution 4.8. This consistency further validates the numerical simulation for the combination of lubricated gravity currents with lubricated power-law fluid.

Similarity solutions for the asymplotic limits in M :
The upper-fluid solid and liquid limits In the general case of ≠ 5 and ≠ 1 our model does not admit global similarity solutions of the first kind. However, a similarity solution arises in part of the flow domain in each of the asymptotic limits of the viscosity ratio, in which the upper fluid layer is either relatively more viscous (M 1, the "solid" limit) or relatively less viscous (M 1, the "liquid" limit) compared with the lower fluid layer. 4.3.1. The top-layer solid limit, M 1 The limit M 1 represents the case where the top fluid is much more viscous than the lubricating fluid. One motivation to study this limit is the geophysical setting of ice sheets creeping over less-viscous lubricated beds. In this case, the leading-order scalling of the power-law fluid local flux in the lubricated region (2.22a) is which is identical to the scaling of the Newtonian lubricating fluid flux (2.22b). In such a case, the three scales , and in the lubricated region cannot be determined independently despite the fact that the upper fluid is a power-law fluid. Therefore, both fluids in the lubricated region have a similarity solution in which the dimensionless front in the lubricated region evolves with an -independent exponent where the intercept may depend on the system dimensionless numbers. The exponent in (4.14) is consistent with the predictions made for the special cases discussed above. Specifically, for = 1 it is (3 + 1)/8 as we predict in equation (4.2b), and for = 5 it is 2 independently of as we predict in equation (4.8b). This similarity solution does not reveal the evolution of the fluid front in the non-lubricated region, which depends on the fluid exponent .
4.3.2. The top-layer "liquid" limit, M 1 The second asymptotic limit, M 1, in which the top fluid is much less viscous than the lubricating fluid, also admits a similarity solution, but of a different sub set of the model. In this case, the leading-order scaling of the power-law fluid local flux in the lubricated region (2.22a) is which is identical to the scaling of the power-law fluid in the non-lubricated region (2.20). Therefore, the evolution of the power-law fluid in the entire domain converges when ℎ to the similarity solution of a non-lubricated power-law GC, in which case the upper fluid front evolves with an -dependent exponent where the intercept may depend on the system dimensionless numbers. This similarity solution does not hold for the lubricating fluid layer since the local flux of the Newtonian lubricating fluid sets a different scaling. Therefore, it does not predict the evolution of the lubricating-fluid front . For − ℎ > 1 the coupling between the upper fluid layer and the lower one grows stronger the more shear-thinning the fluid is ( → ∞), as the exponent in (4.15) grows like . Therefore, (4.16) is expected to be less accurate the more shear-thinning the upper fluid is when − ℎ > 1, and more accurate when − ℎ < 1.

Similarity solutions at the vicinity of the fluid fronts
The thicknesses of the two fluid layers at the vicinity of the fronts have also got self-similar forms, even at the absence of a global similarity solution. These self-similar forms arise because the dynamics at the vicinity of both fronts is dominated by the same dynamics that governs non-lubricated GCs, which are known to have self-similar solutions for any and (Sayag & Worster 2013;Huppert 1982).
The dominance of a non-lubricated-GC dynamics is naturally expected at the vicinity of the front , which is the edge of the non-lubricated region in the flow. In this case we expect consistency of the fluid thickness with Sayag & Worster (2013), which implies a similarity solution of the form (4.17b) for the model described by equations (2.19, 2.20), where is the similarity variable and ≡ ( = ). Therefore, the asymptotic solution for the function near the front in the similarity space is (Sayag & Worster 2013) Near the front the dominance of non-lubricated-GC dynamics is also valid, but requires more careful supporting argument. Specifically, we assert that the free surface slope / is finite at , whereas the slope of the lubricating fluid is singular. Therefore, provided the density difference is non zero, the leading-order flux ℓ (2.22b) is Consequently, the model (2.21a) for the lubricating fluid near is decoupled from , and describes a modified non-lubricated GC of a Newtonian fluid (Huppert 1982). In this case the thickness ℎ has a similarity solution of the form we obtain the closed-form solutions for the fronts intercept (4.23b) The asymptotic solutions 4.18 and 4.21 that we obtain based on the non-lubricated theory appear to provide a precise prediction to the fluid thicknesses of the lubricated GC at the vicinity of the fronts when using the numerically calculated values for , of the lubricated GC ( Figure 5). The values for and (4.23) based on the non-lubricated GC theory provide a rough estimate for those of the lubricated theory. Discrepancies between the two seem to decline the more shear-thickening the upper fluid is (Figure 6a).

Front outstripping
In the absence of similarity solutions the fronts propagate at different rates, implying that under certain conditions the inner lubricating front can outstrip the outer front . We identify two outstripping mechanisms, one driven by the solution intercepts and the other by the exponents.

Intercept outstripping
The intercept difference Δ ≡ − predicted by the non-lubricated theory (4.23) diminishes with M and becomes negative when M grows beyond a threshold value that we denote by M (Figure 6b). This implies that across that threshold in the viscosity ratio > and the lubricating front outstrips the upper fluid front. The threshold viscosity ratio M is found by setting Δ = 0 and solving for M M ≡ 2 Q 9 28 3 9(3 + 1) 8D Therefore, in the purely Newtonian case the threshold viscosity ratio is and is independent of . In the case of constant flux = 1 the threshold is This behaviour reasoned based on the non-lubricated theory is consistent with the trend implied from the numerical solution of the full lubricated theory (Figure 6b), though we do not explicitly observe an outstripping of the outer front by the inner one since the present simulation assumes < .

Exponent outstripping
Another mechanism of outstripping can arise due to the different time exponents of the two fronts. This difference can be evaluated from the similarity solutions at the vicinity of each front. Specifically, denoting the front exponents by ≡ [ (1 + 2 ) + 1]/(5 + 3) (4.17a) and ≡ (3 + 1)/8 (4.20a) then .

(5.4)
When the flow has a global similarity solution ( = 5 or = 1) the exponents of the two fluid fronts are identical Δ = 0, implying that asymptotically in time the gap between the two fronts evolves with the same exponent and the ratio of their position / is constant. When ≠ 1 and ≠ 5 there is no global similarity solution, implying that ≠ and that the fronts ratio evolves proportionally to Δ . Consequently, the gap between the fronts closes down in time and front outstripping occurs when Δ > 0, which is when < 5 & > 1 or  (Figures 8, 9). When the upper fluid is non-Newtonian ( ≠ 1), no global similarity solution arises for a constant-flux discharge ( §2.3). Nevertheless, we find numerically that the characteristic distributions of the fluid thickness in each of the four regimes are qualitatively preserved (Figure 9). Moreover, the thickness of the lubricating fluid in the non-Newtonian case varies very weakly from that in the purely Newtonian case, and the fronts in both cases appear to follow a similar evolution pattern (Figure 9). In contrast, the propagation rate of the front varies strongly with the non-Newtonian properties of the upper fluid, becoming faster than the purely Newtonian case for shear-thickening fluids and slower than the purely Newtonian case for shear-thinning fluids (Figure 9). This pattern in all of the four regimes can be physically rationalised through a dimension-based argument: ( / ℓ ) 1/2 /(2 ) declines radially. Consequently, a shear-thinning fluid ( > 1) becomes increasingly more viscous with radius than a shear-thickening fluid, leading to its relatively slower propagation.
The slower front speed of a shear-thinning fluid compared to a Newtonian fluid implies that outstripping by the front of the Newtonian lubricating fluid may occur after some time ( §5). We investigate this more carefully by considering the evolution of the two fronts within a range of viscosity ratios 10 −2 M 10 4 , a range of fluid exponent = 10, 1, 0.5 and for fixed flux ratio Q = 0.2 and density ratio D = 0.1. Even though no global similarity solution of the first kind exists when , ≠ 1, we find that each front has a power-law evolution in time of the form ( ) ≈ , ( ) ≈ ( − 1) (6.1) in dimensionless units, where we compute the intercepts , and the exponents , , which may depend on the dimensionless groups M , Q, D and , through regression to the numerical solution (Figures 10, 11, and 11). Our numerical solver does not admit front outstripping since it constrains . Therefore, we infer that front outstripping is in progress when the interval − is closing in time. In the upper-fluid "liquid" limit, M 1, we find that the numerical solution for the front of the upper fluid evolves consistently with (4.16), in which = (2 + 2 )/(5 + 3) and = N 1/(5 +3) . For example, when M = 0.01 the computed exponents for = 0.5, 1, 10 differ from the corresponding theoretical values 6/11, 1/2, 22/53 by less than 2% (Figure 10a). The larger discrepancy among those is for the shear-thinning fluid, which is expected due to the stronger coupling in this case between the upper and the lower fluid layers, as mentioned in §4.3.2. The corresponding intercepts differ by less than 5% from the theoretical values. The lubricating fluid front evolves with exponent ≈ 1/2 for = 1, larger than 1/2 for a shear-thinning fluid (by 10% for = 10), and smaller than 1/2 for shear-thickening fluid (by 1.2% for = 1/2) (Figure 10a). The corresponding intercepts have small variation with , being in the range 0.17 0.21, which is much smaller than the intercept in non-lubricated Newtonian GCs (1/3) 1/8 ≈ 0.62.
In the upper-fluid "solid" limit, 1, we find that the front of the lubricating fluid evolves consistently with (4.14), in which → 1/2 is -independent, and the solution for also tends to 1/2 (Figure 10b). For intermediate M values we find that the front exponent varies weakly with M initially (M 10), remaining close to the M 1 asymptotic values (Figure 11a), and the intercept grows with M and varies weakly with ( Figure  11d). The exponent of the lubricating front grows with M over 1/2 for shear-thickening fluids and less than 1/2 for shear-thinning fluids (Figure 11b), whereas the intercept grows monotonically with M for all s while varying very weekly with (Figure 11e).
Throughout the range 10 −2 M 10 4 we find that > for shear-thickening ( < 1), < for shear-thinning fluids ( > 1), and = = 1/2 ± 0.0008 when = 1 (Figure 11c), consistently with our predictions of the similarity solutions near the front ( §5.2). This implies that in the long-time limit the front outstrips when the upper fluid is shear thinning, independently of M , while no exponent outstripping is expected when the upper fluid is shear thickening. However, the exponent differences |Δ | ≡ | − | declines with M so that when M = 10 4 > M it is marginally zero (Figure 11c), implying termination of the exponent-driven outstripping. Simultaneously, the intercept difference − is positive throughout the range of M but approaches zero at the vicinity of M = 10 4 (Figure 11f). This may reflect the approach to an intercept-driven outstripping ( − < 0) across a critical viscosity ratio M > M , given by (5.3), while both exponents converge to 1/2. We note again that numerically we do not resolve explicit outstripping, so the behaviour as M approaches M and grows beyond it reflects mostly the qualitative trend.

Comparison with experimental evidence for = 1 and > 1
The theory we have developed can be validated with the recent laboratory experiments of lubricated GCs at constant-flux ( = 1) that consisted of a strain-rate softening fluid that was lubricated by a denser Newtonian fluid (Kumar et al. 2021). The experimental findings have shown that for flux ratio Q 0.06 both fronts were highly axisymmetric, in experiments that lasted over 10 in some cases, implying that a comparison with our axisymmetric theory should be valid.

Proparation of the fronts
Experimentally it was found that the two fronts evolved faster than those of non-lubricated GCs of the corresponding fluids. Nevertheless, each front had a power-law time evolution with a similar exponent as non-lubricated GCs of the corresponding fluids. In particular, the front evolved with exponent (2 + 2)/(5 + 3) before the lubrication fluid was introduced ) as well as when 5 , while the exponent was larger during the transition between the two intervals. At the same time, the front evolved with exponent 1/2.
To compare our theoretical predictions to the experimental measurements we compute numerical solutions to each of the 15 experiments described in Kumar et al. (2021 , Table 2c and eq. 2.2) (Figures 13, 14) based on the measured quantities Q, M , D and and without any fitting parameter. The viscosity ratio in each of these experiments was in the range 1 M < M , implying that the experiments were all in the exponent-outstripping regime ( §5.2) as well as in the solid-limit regime ( §4.3.1). Consequently, we expect the lubrication front to evolve like ∝ 1/2 . Comparing the time exponents, we find that the numerical solutions and the theoretical predictions are consistent with those measured experimentally for both (Figure 13a, b) and (Figure 14a, b). Specifically, during < 5 the solutions to the front evolve with larger exponent than non-lubricated current (Figure 13a,b), and when 5 the exponent diminishes and converges back to that of non-lubricated GCs (Figure 13b). Similarly, the solutions to the front evolve with an exponent 1/2, consistently with the experimental measurements (Figure 14a, b). In spite of the exponent consistency between the experiment and the theory, the experimental fronts advance faster than the predicted ones (Figure 13c,  d and 14c, d). This discrepancy is larger in experiments with higher polymer concentration (Figure 13d and 14d). It implies that the numerically predicted intercepts are lower than those that were measured experimentally, suggesting that additional physical processes that contribute to a faster front propagation are not accounted for by our theoretical model. We elaborate on potential additional mechanisms in §8.

Thickness evolution
Experimentally, the thickness of the lubricated, non-Newtonian fluid was found to be nearly uniform in the lubricated region. The layer of the lubricating fluid was also largely uniform with localised spikes, and its average thickness was approximated through mass conservation to be ℎ ℓ = ℓ ( − )/ ℓ 2 (Kumar et al. 2021). This nearly uniform pattern differs substantially from the monotonically diminishing thickness of a non-lubricated GC under similar conditions. Our numerical solutions, which do not involve any fitting parameter, follow a similar pattern as in the experiments (Figure 15). The theoretical and experimental patterns are highly consistent during most of the flow (e.g., Figure 15a-c), but discrepancy between the two grows progressively near the fronts , and in the non-lubricated region (e.g., Figure 15d). In particular, the fronts in the experiment evolve faster than those computed numerically, and the thickness distribution of the top fluid at the vicinity of the lubricating front changes more sharply in the computed solutions than in the measurements. This growing discrepancy may imply the action of additional mechanisms, as we elaborate in §8.

Discussion
The lubricated GCs that we consider involve five dimensionless parameters Q, D, M , and , associated with significant qualitative transitions in the structure of the solution, in the relative motion of the fluid fronts and their stability, and in the thickness distribution of the two fluids.
The fluid exponent and the discharge exponent have a dramatic qualitative impact on the solutions. Specifically, similarity solutions exist only when = 1, in which both fluids follow a similar constitutive law, and when = 5. In those cases the two fronts evolve with the same power law in time and as a result the ratio / is constant. In all other cases ( ≠ 1, ≠ 5) the fronts also appear to follow a power law evolution, but each with a different exponent, implying that asymptotically in time the ratio / evolves following a power law in time with exponent Δ (5.4). Therefore, there are solutions where / declines in time resulting in the lubrication front outstripping the front of the upper fluid ( > 1 & < 5, and < 1 & > 5), and otherwise / grows in time (Figure 7).
The flux ratio Q can lead to the emergence of two significantly different patterns, and may have a critical impact on the front stability. When Q < 1 the flux of the lubricating fluid is lower than the top fluid layer. Consequently, the thickness of the lubricating layer is significantly smaller and the propagation of the front is affected by the relatively larger pressure imposed by the thicker top layer ( Figure 9III, and Figure 15). The opposite occurs when Q > 1 -the lubricating fluid is discharged at a larger flux and its thickness is significantly larger than the top fluid layer ( Figure 9II). Q may also have a crucial impact on the stability of the axisymmetric fronts. Preliminary experimental evidence indicate that when the flux is constant ( = 1), the top fluid layer is strain-rate softening ( = 6) and M 1, the initially axisymmetric fronts become unstable when Q 0.1 and develop fingering patterns after an initial axisymmetric spreading (Kumar et al. 2021). Similar symmetry breaking also emerges in the purely Newtonian case when 0.14 Q 0.44 (Kowal & Worster 2015). The viscosity ratio M affects the relative motion of the fronts, and the relative thickness of the fluid layers. At a high viscosity ratio (M 1) the more viscous top fluid is effectively solid-like compared with the less viscous lower fluid. Consequently, the flow in the lubricated region is independent of the fluid exponent , and both fluid layers in that region follow the same similarity solution, in which the front of the lubrication fluid evolves with a time exponent (3 + 1)/8, same as a non-lubricated Newtonian GC (4.14). In the low viscosity ratio (M 1) the top fluid is significantly more mobile than the lower fluid layer, which does not provide an effective lubrication. In this case the two fluid layers in the lubricated region do not exhibit a global similarity solution, but the top fluid layer along the whole domain does. Consequently, a self-similar solution exist in the top-fluid layer, in which the front evolves with a time exponent [ (2 + 1) + 1]/(5 + 3), same as a non-lubricated GCs (4.16). The impact of the viscosity ratio on the fluid thickness distributions can be appreciated through the constant flux case ( = 1), in which the free surface of the top fluid is substantially flatter in the M 1 case than in the M 1 case (Figure 9). Independently of the value of the dimensionless parameters, the solutions at the vicinity of the fronts are also self similar, with exponents consistent with those of non-lubricated GCs of power-law fluids (Huppert 1982;Sayag & Worster 2013). Particularly, both fronts evolve with an exponent [ (2 + 1) + 1]/(5 + 3), which simplifies to (3 + 1)/8 for the Newtonian lubricating fluid. The intercepts of the fronts depend on the different dimensionless numbers of the system. Together, the exponents and the intercepts provide insights into the interaction between the two fronts. One important consequence of that interaction is the outstripping of the upper fluid front by the lower lubricating fluid front, which can occur either through the intercept difference Δ or through the exponent difference Δ . The condition for an interceptdriven outstripping can be formalised in terms of the critical viscosity ratio M (Q, D, , ) (5.2), so that outstripping occurs when M > M . Physically this implies that when Q 1 then M ∝ 1/Q 3 1 and the top fluid should be significantly more viscous than the lower fluid for outstripping to occur, in which case the top fluid deforms substantially slower making it easier for the lubricating fluid front to outstrip. Alternatively, when Q 1 then for shear-thinning fluids M ( → ∞) ∝ Q 1/5 1, whereas for shear-thickening fluids M ( → 0) ∝ 1/Q 1/3 1. In either case the superiority of the lubricating fluid flux leads to front outstripping. When M < M then Δ > 0 and outstripping can only be driven by the exponent difference. As discussed above, the condition for that mechanism depends on the values of and (Figure 7). For example, shear-thinning fluids at relatively low discharge exponent ( < 5) become increasingly more viscous as they expand radially, resulting in slower front velocity than the lubricating fluid front. Moreover, their thickness and correspondingly the pressure they apply on the lubricating fluid is relatively larger and contribute further to the radial spreading of the lubricating fluid. We find that the interceptdriven outstripping can occur significantly faster than exponent-driven outstripping. For example, as in the constant flux ( = 1) case intercept-driven outstripping occurs at roughly / 10, which is much faster than the / 10 3 in the exponent-driven case ( Figure 10). It is important to note that the global similarity solution that we find for a discharge exponent = 5 arises in additional axisymmetric GCs of different settings, which a priori appear remotely related. This includes for example, isothermal lava domes that are modeled as axisymmetric GCs of visco-plastic fluids (Balmforth et al. 2000). The structure of the similarity solution in that case is identical to the lubricating GCs that we consider, in which the front evolves like ∝ 2 and the fluid thickness evolves like ℎ ∝ independently of the fluid exponent . Another system with a similarity solution at = 5 is the axisymmetric viscous GCs flowing over a porous medium (Spannuth et al. 2009). Such a similarity among a broad range of physical systems may not be coincidental and could imply a more general symmetry associated with the circular geometry.
Many aspects of the theory were found consistent with experiments performed for the dimensionless parameters D ≈ 0.15, Q < 0.06, 1 M < M , = 1, and > 1 (Kumar et al. 2021). In particular, the time evolution of both fluid fronts predicted by the theory is consistent with the power law measured in the experiments. In addition, the thickness distribution we predict for the top fluid layer is in good agreement with the experimental measurements. However, some discrepancies that arise may imply that the theory is not entirely complete. Specifically, the theoretical predictions for the intercepts do not accurately capture the measured ones, particularly in the case of the lubricating front, which evolves faster than the theoretical predictions. Several potential physical mechanisms that the present theory does not account for may explain these discrepancies. One possible mechanism is that the lubrication front advances as a hydrofracture in between the substrate and the relatively solid viscous fluid layer (Ball & Neufeld 2018), or as a shock in the fluid-fluid interface at (Dauck et al. 2019). In addition, discrepancy between the experimental measurements of before introducing the lubrication fluid ( / < 1) and the theoretical prediction of a nonlubricated GC, particularly for the 2% polymer concentration (Kumar et al. 2021), may imply that the power-law constitutive equation that we use is incomplete. Specifically, the time for the viscosity to adjust to the evolving strain rates may not be instantaneous as we assume, but finite. In addition polymer entanglements may arise in higher polymer concentrations that potentially drive wall slip at the fluid-solid interface through adhesive failure of the polymer chains at the solid surface or through cohesive failure owing to disantanglement of chains in the bulk from chains adsorbed at the wall (Brochard & Gennes 1992). The implications of these potential physical mechanisms will be addressed in future studies.

Conclusions
Lubricated gravity currents are controlled by complex interactions between two fluid layers. The lower lubricating layer modifies the friction between the substrate and the top layer, which in turn applies stresses that affect the distribution of the lubricating layer. The resulting flow can vary dramatically from non-lubricated gravity currents.
Unlike previous axisymmetric gravity current models that involve a single fluid layer (Huppert 1982; Sayag & Worster 2013), or two coupled layers that have the same constitutive structure (Kowal & Worster 2015), the flow we consider does not in general admit a global self-similar solution. Exceptional cases, in which the model has a global similarity solutions are the purely Newtonian case and the case of a discharge exponent = 5. Several other situations admit a similarity solution in part of the domain. This includes the asymptotic limits of the viscosity ratio, corresponding to the top layer solid (M 1) and liquid (M 1) limits, and the solution at the vicinity of the fluid fronts. The latter implies that the time evolution of each front is proportional to that of the non-lubricated front of the corresponding fluid. This implies that generally the ratio of the two fronts positions / either diverge in time or converge to zero, and that the difference of the front intercepts can change sign. Consequently, there are situations where the lubricating fluid front can outstrip the outer-fluid front. This situation can arise when the top fluid is non-Newtonian, having a different exponent than the lubricating fluid, but also when there is global similarity and the viscosity ratio is larger than a critical value M that depends on the fluids flux and density ratios, on the discharge exponent and on the power-law fluid exponent. In the canonical case of constant flux the flow we consider has no global similarity when the top fluid is non-Newtonian. The evolution of the lubricating front does not differ substantially from the purely Newtonian case, but the evolution of the top fluid front differs substantially. Our model solutions are found consistent with laboratory experiments (Kumar et al. 2021), particularly in predicting the time exponents of the front evolution and the thickness fields. Discrepancies that we find in the intercepts predictions, particularly that of the lubricating fluid front, suggest that additional physical mechanisms may contribute to the front evolution, such as hydrofracturing or wall-slip along the substrate. Exploration of these mechanisms will be the topic of future studies. Our results have implications to the understanding of spatiotemporal distribution of lubrication networks beneath ice sheets and to understanding their stability. AG was partially supported by VATAT High-TEC fellowship for excellent women in science. This research was supported by the GERMAN-ISRAELI FOUNDATION (grant No. I240430182015). Declaration of Interests. The authors report no conflict of interest. for constant flux = 1. We use this solution to validate our numerical solution in the Q = 0 limit. Specifically, we solve the dimensionless equation set ( §2.3) with Q = 0, zero initial thickness ( , 0) = 0, and logarithmically spaced spatial mesh with 1200 points. We find our solutions for the fluid height and for the leading front consistent with the theoretical predictions ( Figure 16). Repeating the same computation for varying spatial resolutions, we find that the convergence accuracy of the front exponent to the predicted theoretical value grows with a the number of spatial grid points (Figure 16d, inset).

A.2. Lubricated Newtonian gravity currents
In the limit = 1, = 1 our numerical model converges to the purely Newtonian lubricated GC discharged at constant flux (Kowal & Worster 2015). Considering first the specific case where M = 10000, Q = 0.2 and D = 0.1, we find that both the fronts , and the upper fluid height at the lubricant front ( ) converge to the theoretical values , and ( ), respectively ( Figure 17). Second, we find that our solutions for the coefficients and is consistent with the theoretical predictions for a wide range of M values (Figures 11. As shown in Appendix A.1, small discrepancies from the predicted values are due to low spatial resolution. Lastly, the solutions to the specific regime discussed in §6 (Figures 9) is another evidence for the consistency between our numerical results and those of (Kowal & Worster 2015, Figure 13).