2.1. No-slip region

Integrating the horizontal force balance

(2.10)\begin{equation} \frac{\partial}{\partial z}\left(\mu\,\frac{\partial \boldsymbol u}{\partial z}\right) = \rho g\, \boldsymbol{\nabla} H \end{equation}

in the no-slip region $r_L< r< r_N$ results in the velocity field

(2.11)\begin{equation} \boldsymbol{u} =\frac{1}{n+1} \left(\frac{\rho g}{\tilde \mu}\right)^n (H^{n+1}-(H-z)^{n+1}) |\boldsymbol{\nabla} H|^{n-1} (-\boldsymbol{\nabla} H), \end{equation}

and corresponding depth-integrated flux

(2.12)\begin{equation} \boldsymbol{q} =\frac{1}{n+2} \left(\frac{\rho g}{\tilde \mu}\right)^n H^{n+2} |\boldsymbol{\nabla} H|^{n-1} (-\boldsymbol{\nabla} H), \end{equation}

which are of the same functional form as for axisymmetric flows, including the non-axisymmetric contributions. These agree with Kowal & Worster (Reference Kowal and Worster2015). Linearising gives rise to the components

(2.13)\begin{gather} q_{r1} ={-}\frac{1}{n+2} \left(\frac{\rho g}{\tilde\mu }\right)^n H_0^{n+1} \left| \frac{\partial H_0}{\partial r}\right|^{n-1}\left(nH_0\,\frac{\partial H_1}{\partial r}+ (n+2)\,\frac{\partial H_0}{\partial r}\,H_1\right), \end{gather}

(2.14)\begin{gather}q_{\theta1} ={-}\frac{1}{n+2}\left(\frac{\rho g}{\tilde\mu }\right)^n H_0^{n+2} \left| \frac{\partial H_0}{\partial r}\right|^{n-1}\frac{1}{r}\,\frac{\partial H_1}{\partial\theta}, \end{gather}

of the perturbations to the flux.

Mass conservation, at first order, is described by

(2.15)\begin{equation} \frac{\partial H_1}{\partial t} ={-} \frac{1}{r}\,\frac{\partial (rq_{r1})}{\partial r} - \frac{1}{r}\,\frac{\partial q_{\theta1}}{\partial \theta}, \end{equation}

within the no-slip region $r_L< r< r_N$. We note that additional terms are required when transforming to similarity variables (2.42) to capture terms involving the base state flow owing to perturbations in the frontal position.

2.2. Lubricated region

Unlike single-layer flows for any value of $n$, and lubricated flows for $n=1$, there are no closed-form expressions for the velocity and flux, which include non-axisymmetric contributions, unless linearised.

We proceed by starting from the horizontal force balance

(2.16)\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.17)\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}

in the upper and lower layers, supplemented by the stress-free boundary condition at $z=H$, continuity of velocity and shear stress at $z=h$, and the no-slip boundary condition at $z=0$.

For the upper layer, this can be integrated directly so that

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

where $\boldsymbol u_I$ is the interfacial velocity, to be determined by matching with the velocity of the lower layer. Linearising gives rise to the perturbed velocity

(2.19) \begin{align} \boldsymbol{u}_1 &={-}\left(\frac{\rho g}{\tilde \mu}\right)^n \frac{1}{n+1} [[(H_0-z)^{n+1}- (H_0-h_0)^{n+1}]|\boldsymbol{\nabla} H_0|^{n-1} (-\boldsymbol{\nabla} H_1)\nonumber\\ &\quad + (n-1)[(H_0-z)^{n+1}- (H_0-h_0)^{n+1}]|\boldsymbol{\nabla} H_0|^{n-3}(\boldsymbol{\nabla} H_1\boldsymbol{\cdot}\boldsymbol{\nabla} H_0) (-\boldsymbol{\nabla} H_0) \nonumber\\ &\quad{}+(n+1)[H_1(H_0-z)^n - (H_1-h_1)(H_0-h_0)^n|\boldsymbol{\nabla} H_0|^{n-1} (-\boldsymbol{\nabla} H_0)]] + \boldsymbol u_{I1}, \end{align}

where $\boldsymbol u_{I1}$ is the perturbed part of the interfacial velocity $\boldsymbol u_{I}$.

For the lower layer, we obtain

(2.20)\begin{equation} \frac{\partial\boldsymbol{u_l}}{\partial z} = \lvert\boldsymbol{a}-z\boldsymbol{c}\rvert^{n-1}(\boldsymbol{a}-z\boldsymbol{c}), \end{equation}

where

(2.21)\begin{gather} \boldsymbol{a} ={-}\frac{\rho g}{\tilde{\mu}}\,\mathcal{M}\left(H\boldsymbol{\nabla} H+\mathcal{D}h\boldsymbol{\nabla} h\right), \end{gather}

(2.22)\begin{gather} \boldsymbol{c} ={-}\frac{\rho g}{\tilde{\mu}}\,\mathcal{M}\left(\boldsymbol{\nabla} H+\mathcal{D}\boldsymbol{\nabla} h\right). \end{gather}

Linearising in $\epsilon$ and integrating the linearised expressions yields

(2.23) \begin{align} u_{lr1} &= \frac{1}{(n+1) c_{{r0}}^2} [(c_{{r1}} \left(a_{{r0}}+n z c_{{r0}}\right)-(n+1) a_{{r1}} c_{{r0}}) \left| a_{{r0}}-z c_{{r0}}\right| {}^{n-1} \left(a_{{r0}}-z c_{{r0}}\right)\nonumber\\ &\quad+\left((n+1) a_{{r1}} c_{{r0}}-a_{{r0}} c_{{r1}}\right) \left| a_{{r0}}\right| {}^{n-1}a_{{r0}} ], \end{align}

(2.24) \begin{align} u_{l\theta1} &=\frac{1}{n (n+1) c_{{r0}}^2} [ \left(c_{{\theta 1}} \left(a_{{r0}}+n z c_{{r0}}\right)-(n+1) a_{{\theta 1}} c_{{r0}}\right) \left| a_{{r0}}-z c_{{r0}}\right| {}^{n-1} \left(a_{{r0}}-z c_{{r0}}\right)\nonumber\\ &\quad + \left((n+1) a_{{\theta 1}} c_{{r0}}-a_{{r0}} c_{{\theta 1}}\right) \left| a_{{r0}}\right| {}^{n-1} a_{{r0}} ], \end{align}

from which the interfacial velocity $\boldsymbol {u_I}$ can be deduced. Explicitly,

(2.25)\begin{equation} \boldsymbol{u_I}= u_{I0}\boldsymbol{e_r}+\epsilon( u_{Ir1}\boldsymbol{e_r}+ u_{I\theta1}\boldsymbol{e_\theta}), \end{equation}

where

(2.26)\begin{gather} u_{I0}= \frac{1}{n+1}\left(\frac{\rho g}{\tilde \mu_l}\right)^{n} \frac{1}{\mathcal{D}\,\partial h_{0}/\partial r+\partial H_0/\partial r}\left[ \left|(H_0-h_0)\,\frac{\partial H_0}{\partial r}\right|^{n+1} \right.\nonumber\\ \hspace{12pt}\left. - \left|h_0\left(\mathcal{D}\,\frac{\partial h_0}{\partial r}+\frac{\partial H_0}{\partial r}\right)+(H_0-h_0)\,\frac{\partial H_0}{\partial r}\right|^{n+1} \right], \end{gather}

(2.27)\begin{gather} u_{Ir1}= \left[h_1\,\frac{\partial u_{lr0}}{\partial z}+u_{lr1}\right]_{z=h_0}, \end{gather}

(2.28)\begin{gather} u_{I\theta1}= \left[h_1\,\frac{\partial u_{l\theta0}}{\partial z}+u_{l\theta1}\right]_{z=h_0}. \end{gather}

Note that since the basic state is axisymmetric, it follows that $a_{\theta 0}=c_{\theta 0}=0$. Expressions for $a_{r0}, a_{r1}, a_{\theta 1}, c_{r0}, c_{r1}$ and $c_{\theta 1}$ are specified explicitly in Appendix B.

Further integration yields the following expressions for the $r$-components,

(2.29)\begin{gather} q_{lr1} = A_1h_1+A_2a_{r1}+A_3c_{r1}, \end{gather}

(2.30)\begin{gather} q_{r1} = A_4\,\frac{\partial H_1}{\partial r}+ A_5H_1+ A_6h_1+A_7a_{r1}+A_8c_{r1}, \end{gather}

and the $\theta$-components,

(2.31)\begin{gather} q_{l\theta1} = A_9a_{\theta1}+A_{10}c_{\theta1}, \end{gather}

(2.32)\begin{gather} q_{\theta1} = A_{11}\,\frac{1}{r}\,\frac{\partial H_1}{\partial\theta}+ A_{12}a_{\theta1}+A_{13}c_{\theta1}, \end{gather}

of the perturbations to the fluxes of the two layers in the lubricated region, where the $A_i$ are specified in Appendix B. These expressions reduce to those of Kowal & Worster (Reference Kowal and Worster2019b) for $n=1$.

Mass conservation, at first order in $\epsilon$, is described by

(2.33)\begin{equation} \frac{\partial h_1}{\partial t} ={-} \frac{1}{r}\,\frac{\partial (rq_{lr1})}{\partial r} - \frac{1}{r}\,\frac{\partial q_{l\theta1}}{\partial \theta} \end{equation}

for the lower layer, and

(2.34)\begin{equation} \frac{\partial (H_1-h_1)}{\partial t} ={-} \frac{1}{r}\,\frac{\partial (rq_{r1})}{\partial r} - \frac{1}{r}\,\frac{\partial q_{\theta1}}{\partial \theta} \end{equation}

for the upper layer within the lubricated region $0< r< r_L$. Similarly to the no-slip region, additional terms are required when transforming to similarity variables (2.41) to capture terms involving the base state flow owing to perturbations in the frontal position.

2.3. Boundary conditions

We apply the source flux conditions

(2.35a,b)\begin{equation} \lim_{r\to0} 2{\rm \pi} r q_{lr} = Q_{l0}, \quad \lim_{r\to0} 2{\rm \pi} r q_{r} = Q_{0}, \end{equation}

the thickness and height continuity conditions

(2.36a,b)\begin{equation} [H]^+_-= 0 \quad \text{and} \quad [(\boldsymbol{q}+\boldsymbol{q}_l)\boldsymbol\cdot \boldsymbol{n_L}]^+= [\boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{n_L}]^- \quad (r=r_L), \end{equation}

where $\boldsymbol {n_L}=\boldsymbol {e}_r -\boldsymbol {e}_\theta ({1}/{r_L})\,{\partial r_L}/{\partial \theta }+ {O}(\epsilon ^2)$ is an outward normal vector at the lubrication front, and the kinematic conditions

(2.37)\begin{equation} \dot r_L = \lim_{r\to r_L} \left.\left[ q_{lr} - q_{l\theta}\,\frac{1}{r_L}\,\frac{\partial r_L}{\partial \theta}\right] \right/h, \end{equation}

for the lubrication front and

(2.38)\begin{equation} \dot r_N = \lim_{r\to r_N} \left.\left[ q_{r} - q_{\theta}\,\frac{1}{r_N}\,\frac{\partial r_N}{\partial \theta}\right] \right/H, \end{equation}

for the leading edge. We also apply the zero-flux condition

(2.39)\begin{equation} \boldsymbol{q_l}\boldsymbol{\cdot} \boldsymbol{n_L} = 0 \quad (r=r_L) \end{equation}

at the lubrication front for $\mathcal {D}\neq 0$, and

(2.40)\begin{equation} \boldsymbol{q}\boldsymbol{\cdot} \boldsymbol{n_N}=0\quad (r=r_N) \end{equation}

at the leading edge, where $\boldsymbol {n_N}=\boldsymbol {e}_r -\boldsymbol {e}_\theta ({1}/{r_N})\,{\partial r_N}/{\partial \theta }+ {O}(\epsilon ^2)$ is an outward normal vector at the leading edge.

2.4. Similarity coordinates

To conduct a linear stability analysis about the self-similar axisymmetric flow of Part 1, we revert to the similarity coordinates $(\xi, \phi, \tau )$ defined by

(2.41)\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}

(2.42)\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}

(2.43a,b)\begin{gather}\tau=\log t, \quad \phi=\theta. \end{gather}

where $0<\xi <1$ corresponds to the lubricated region $0< r< r_L$, and $1<\xi <2$ corresponds to the no-slip region $r_L< r< r_N$. The constants $\alpha$, $\beta$ and $\gamma$ are given by

(2.44a–c)\begin{equation} \alpha =\frac{n}{5n+3}, \quad \beta = \frac{2n+2}{5n+3}, \quad \gamma=\frac{2n+1}{5n+3}, \end{equation}

as specified in Part 1.

The lubricated region is therefore mapped to the interval $(0,1)$, and the no-slip region is mapped to the interval $(1,2)$. Perturbations to the two fronts can be read from

(2.45a,b)\begin{align} \xi_L(\phi,\tau) = \xi_{L0}+\epsilon \xi_{L1}\exp({\sigma \tau+{\rm i}k\phi}), \quad \xi_N(\phi,\tau) = \xi_{N0}+\epsilon \xi_{N1}\exp({\sigma \tau+{\rm i}k\phi}), \end{align}

in similarity coordinates. Here, $\xi _{L0}$ and $\xi _{N0}$ correspond to the unperturbed positions of the intrusion front and leading edge, respectively. Both $\xi _{L0}$ and $\xi _{N0}$ are constants. We are searching for normal mode solutions of growth rate $\sigma$ and azimuthal wavenumber $k$, which exist under the change of variables (2.41)–(2.43a,b). Under this transformation, contributions owing to the perturbations to the two frontal positions are reflected through appropriate terms in the governing equations, rather than through the boundary conditions. Such an approach eliminates difficulties associated with the stress singularities at the two fronts.

The zeroth- and first-order surface heights are transformed as

(2.46) \begin{equation} \left(\begin{array}{@{}c@{}} h_0(r,t)\\H_0(r,t)\\h_1(r,\theta,t)\\H_1(r,\theta,t) \end{array}\right) = \left(\frac{\rho g}{\tilde \mu}\right)^at^bQ_0^c\boldsymbol{\cdot} \left(\begin{array}{@{}c@{}}f_0(\xi)\\F_0(\xi)\\f_1(\xi) \exp({\sigma \tau+{\rm i}k\phi})\\F_1(\xi) \exp({\sigma \tau+{\rm i}k\phi})\end{array}\right), \end{equation}

and the components of the flux of the two layers are transformed as

(2.47)\begin{equation} \left( \begin{array}{@{}c@{}} q_{lr0}(r,t)\\ q_{r0}(r,t) \end{array} \right) = \left(\frac{\rho g}{\tilde \mu}\right)^{-\alpha}t^{-\beta}Q_0^{1-\gamma} \left( \begin{array}{@{}c@{}} \tilde q_{lr0}(\xi)\\ \tilde q_{r0}(\xi) \end{array} \right) \end{equation}

at zeroth order, and

(2.48)\begin{equation} \left(\begin{array}{@{}c@{}} q_{lr1}(r,\theta,t)\\ q_{l\theta1}(r,\theta,t)\\ q_{r1}(r,\theta,t)\\ q_{\theta1}(r,\theta,t) \end{array} \right)= \left(\frac{\rho g}{\tilde \mu}\right)^{-\alpha}t^{-\beta}Q_0^{1-\gamma} \exp({\sigma \tau+{\rm i}k\phi}) \left( \begin{array}{@{}c@{}} \tilde q_{lr1}(\xi)\\ \tilde q_{l\theta1}(\xi)\\ \tilde q_{r1}(\xi)\\ \tilde q_{\theta1}(\xi) \end{array} \right) \end{equation}

at first order, where the constants $a$, $b$ and $c$ are given by

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

as functions of $n$.

Correspondingly, after dropping tildes for convenience, the components of the flux perturbations are given by the expressions

(2.50)\begin{gather} q_{lr1}=B_1 f_1'+B_2 F_1'+B_3 f_1+B_4 F_1+B_5 \xi _{{L1}}, \end{gather}

(2.51)\begin{gather} q_{l\theta 1}={\rm i}k (B_6 f_1 + B_7 F_1), \end{gather}

for the lower layer, and

(2.52)\begin{gather} q_{r1} = B_8 f_1'+B_9 F_1'+B_{10} f_1+B_{11} F_1+B_{12} \xi _{{L1}}, \end{gather}

(2.53)\begin{gather} q_{\theta 1} = {\rm i}k(B_{13} f_1+B_{14} F_1) -{\rm i}k\xi\,\frac{\xi_{L1}}{\xi_{L0}}\,(B_{13} f_0'+B_{14} F_0'), \end{gather}

for the upper layer. In the no-slip region, the components become

(2.54)\begin{gather} q_{r1} = B_{15}F_1'+B_{16}F_1 + B_{17} \left(\xi _{{N1}}-\xi _{{L1}}\right), \end{gather}

(2.55)\begin{gather} q_{\theta 1} = {\rm i}kB_{18}\left((\xi_{L0}-\xi_{N0}) F_1 + \xi_{{N1}}(\xi-1) F_0' - \xi_{{L1}}(\xi-2)F_0' \right), \end{gather}

where the $B_i$ are specified in Appendix B. These expressions reduce to those of Kowal & Worster (Reference Kowal and Worster2019b) for $n=1$.

The mass conservation equations become

(2.56)\begin{equation} \left(\sigma+\frac{n -1}{5 n+3}\right)f_1 -\frac{2 (n+1)}{5 n+3}\,\xi f_1' -\frac{ \sigma \xi _{L1}}{\xi _{{L0}}}\,\xi f_0' ={-}\frac{ (\xi q_{{lr1}})' +{\rm i} k q_{{l\theta1}} }{\xi \xi _{{L0}}} + \frac{ \xi _{L1} \left(\xi q_{{lr0}}\right)'}{\xi \xi _{{L0}}^2} \end{equation}

for the lower layer of the lubricated region, and

(2.57)$$\begin{gather} \left(\sigma+\frac{n -1}{5 n+3}\right)(F_1-f_1) -\frac{2 (n+1)}{5 n+3}\,\xi (F_1'-f_1') -\frac{ \sigma \xi _{L1}}{\xi _{{L0}}}\,\xi (F_0'-f_0')\nonumber\\ ={-}\frac{ (\xi q_{{r1}})' +{\rm i} k q_{{\theta1}} }{\xi \xi _{{L0}}} + \frac{ \xi _{L1} \left(\xi q_{{r0}}\right)'}{\xi \xi _{{L0}}^2} \end{gather}$$

for the upper layer of the lubricated region. These include contributions owing to the perturbations to the frontal positions. The mass conservation equation in the no-slip region becomes

(2.58)$$\begin{gather} \left(\sigma+\frac{n -1}{5 n+3}\right) C_{1} F_1 - \frac{2(n+1)}{5 n+3}\,C_{2} F_1' - \frac{2(n+1)}{5 n+3}\,C_{3} F_0' + \sigma C_{4}F_0' \nonumber\\ =C_{5}(q_{{r1}}-{\rm i} k q_{{\theta1}}) +q_{{r1}}' +C_{6}q_{{r0}}' + C_{7}q_{r0}, \end{gather}$$

where the $C_i$ are specified in Appendix B.

The source flux boundary conditions reduce to

(2.59a,b)\begin{equation} \lim_{\xi\rightarrow0}2{\rm \pi}\xi(\xi_{L0} q_{lr1}+\xi_{L1} q_{lr0}) = 0, \quad \lim_{\xi\rightarrow0}2{\rm \pi}\xi(\xi_{L0} q_{r1}+\xi_{L1} q_{r0}) = 0, \end{equation}

and the matching conditions at the lubrication front reduce to

(2.60)\begin{gather} \left[ F_1 \right]^+_-= 0 \quad (\xi=1), \end{gather}

(2.61)\begin{gather} \left[ q_{r1} \right]^+_-= 0 \quad (\xi=1). \end{gather}

Note that contributions owing to the perturbations to the frontal positions do not appear in these matching conditions as they are built into the governing equations instead. The remaining boundary conditions are the zero flux conditions

(2.62)\begin{equation} q_{lr1} = 0 \quad (\xi=1) \end{equation}

at the lubrication front, and

(2.63)\begin{equation} q_{r1} = 0 \quad (\xi=2) \end{equation}

at the leading edge.

Note that the fronts are given by $\xi =1$ and $\xi =2$ by the definition (2.41)–(2.42) of the scaled similarity coordinate, as $\xi _L$ and $\xi _N$ are scaled out. The perturbations to the front (from linearising $\xi _L=\xi _{L0}+\epsilon \xi _{L1}$ and $\xi _N=\xi _{N0}+\epsilon \xi _{N1}$) are factored into the governing equations, rather than the radial coordinate by scaling $\xi _L$ and $\xi _N$ out as in (2.41)–(2.42).

The kinematic conditions become

(2.64)\begin{equation} \frac{2n+2}{5n+3}\,(1+\sigma)\xi_{L1} = \lim_{\xi\to1} \left[ \frac{q_{lr1}}{f_0} - \frac{ q_{lr0}f_1}{f_0^2} \right] \end{equation}

at the lubrication front, and

(2.65)\begin{equation} \frac{2n+2}{5n+3}\,(1+\sigma)\xi_{N1} = \lim_{\xi\to2} \left[ \frac{q_{r1}}{F_0} - \frac{ q_{r0}F_1}{F_0^2} \right] \end{equation}

at the leading edge, which lead to the asymptotic solutions described in the following subsection.

3.1. Asymptotic solutions near the two fronts

An asymptotic analysis near the two fronts, in which the governing equations (2.56) and (2.58) are solved in an inner region by rescaling $f_1=\delta ^p \hat f_1$, $F_1=\delta ^p \hat F_1$, $\xi =1-\delta X$ (near the intrusion front) and $\xi =2-\delta X$ (near the leading edge), and balancing dominant terms in the limit $\delta \ll 1$, gives rise to $p={n}/{(2n+1)}$ and the asymptotic solutions

(3.1)\begin{equation} f_1\sim \frac{(5 n\!+\!3) \sigma +2 (n\!+\!1)^2 }{(n+1) (2 n+1) } \left[\frac{(n+1) (n+2)}{5 n+3} \left(\frac{2 n+1}{4 n\mathcal{MD}\xi_{L0}}\right)^{n} \right]^{{1}/({2 n+1})} \xi_{L1}(1\!-\!\xi)^{{n}/({2n+1})}, \end{equation}

as $\xi \to 1^-$, near the lubrication front, and

(3.2)\begin{equation} F_1\sim \mathcal{A}(2-\xi)^{{n}/({2n+1})}, \end{equation}

as $\xi \to 2^-$, near the leading edge, where

(3.3)\begin{align} \mathcal{A}=\left(\frac{2 (n\!+\!1) (n\!+\!2)\xi_{N0} (\xi_{N0}\!-\!\xi_{L0})^n}{(5 n\!+\!3) (2 n\!+\!1)^{n+1}n^n }\right)^{{1}/({2 n\!+\!1})} \left[\frac{n (\xi_{N1}\!-\!\xi_{L1})}{\xi_{N0}\!-\!\xi_{L0}}+\frac{\xi_{N1}}{\xi_{N0}} \left(\frac{(5 n+3) \sigma }{2 (n+1)}+1\right)\right]. \end{align}

These asymptotic solutions are of the same spatial structure as those of the basic state, with prefactors proportional to a linear combination of the perturbations to the frontal positions. These reduce to the asymptotic solutions of Kowal & Worster (Reference Kowal and Worster2019b) in the limit $n=1$. The asymptotic solutions are used to alleviate difficulties associated with the stress singularities that occur at the two fronts, when solving for the solutions numerically.

3.2. Transformation near the origin

An artefact of radially spreading lubricated viscous gravity currents, supplied at constant flux at the origin, is that the thickness of both layers of fluid approaches a point singularity at the origin, as a finite amount of fluid is being supplied from a single point. The form of the solutions, towards which the surface heights approach at zeroth order in $\epsilon$, are specified in Part 1. The asymptotic behaviour is of different character depending on the value of $n$, specifically, depending on whether $n<1$, $n=1$ or $n>1$. A similar phenomenon occurs at first order, which we examine by rescaling $\xi =\delta X$, $f_1=\delta ^p \hat f_1$, $F_1=\delta ^p \hat F_1$ and balancing dominant terms of (2.56)–(2.57) in the limit $\delta \ll 1$.

For $n<1$, the general solution for the perturbations to the surface heights $f_1$ and $F_1$ approaches the functional form $\xi ^\lambda$, where

(3.4)\begin{equation} \lambda=\lambda_{{\pm}} = \frac{1-n\pm(n+1) \sqrt{4 k^2 n +(n-1)^2}}{2 n (n+1)}. \end{equation}

For $n>1$, the exponent is, instead, given by

(3.5)\begin{equation} \lambda=\lambda_{{\pm}} = \frac{n-1\pm\sqrt{4 k^2 n+(n-1)^2}}{2 n}. \end{equation}

The dominant term as $\xi \to 0$ corresponds to $\lambda =\lambda _{-}$. In the limit $n\to 1$, approaching from either the left or the right, the power-law dependence of $f_1(\xi )$ and $F_1(\xi )$ is of the form $\xi ^{-k}$.

These exponents become large in magnitude for large $k$, for any $n$. Therefore, to resolve this singularity at the origin for all wavenumbers and to ensure numerical stability, we reformulate the problem in terms of $g_1(\xi )=\xi ^{-\lambda _-}\,f_1(\xi )$ and $G_1(\xi )=\xi ^{-\lambda _-}\,F_1(\xi )$, instead of $f_1(\xi )$ and $F_1(\xi )$, and revert back to $f_1(\xi )$ and $F_1(\xi )$ through a change of variables after the governing equations have been solved numerically. Although it does not provide a formal asymptotic solution, this is useful in regularising numerical computations by providing a convenient choice for a scaling factor.

As described in Kowal & Worster (Reference Kowal and Worster2019b), for $n=1$ we instead solve for

(3.6)\begin{equation} ( g_1 , G_1 ) = \xi^{k} ( - \log \xi)^{3/4} (f_1(\xi), F_1(\xi)). \end{equation}

The prefactor, similarly, involves an exponent that grows with $k$.