1. Introduction
Viscous, or Saffman–Taylor, fingering is one of the canonical fluid-mechanical instabilities, which occurs when a lower-viscosity fluid is driven into a restricted environment occupied by fluid with a (sufficiently) higher viscosity. In the context of flow in a porous medium (e.g. Homsy Reference Homsy1987), it has huge economic impact by significantly reducing the efficiency of oil extraction from reservoirs by water injection (Lake Reference Lake1989). In the context of flow in a Hele-Shaw cell, it has generated significant scientific interest following the seminal paper by Saffman & Taylor (Reference Saffman and Taylor1958) as a prototypical example of pattern formation and selection, particularly in the limit of small or vanishing surface tension (see e.g. Couder Reference Couder, G.K., Moffatt and Worster2000; Bischofberger, Ramachandran & Nagel Reference Bischofberger, Ramachandran and Nagel2014; Andersen et al. Reference Andersen, Lustri, McCue and Trinh2024). Other recent work considers the possible suppression or modification of Saffman–Taylor fingering by varying the geometry of the Hele-Shaw cell boundaries (e.g. Pihler-Puzovic et al. Reference Pihler-Puzovic, Illien, Heil and Juel2012; Al-Housseiny, Tsai & Stone Reference Al-Housseiny, Tsai and Stone2012; Zheng, Kim & Stone Reference Zheng, Kim and Stone2015; Peng & Lister Reference Peng and Lister2019), viscous fingering in two-layer viscous gravity currents (e.g. Kowal & Worster Reference Kowal and Worster2019a,Reference Kowal and Worsterb; Dauck Reference Dauck2020) and bubble compressibility (Cuttle, Morrow & MacMinn Reference Cuttle, Morrow and MacMinn2023).
The linear instability of immiscible displacement in a Hele-Shaw cell or porous medium has been variously analysed for both unidirectional and radial motion (e.g. Hill Reference Hill1952; Saffman & Taylor Reference Saffman and Taylor1958; Chouke, van Meurs & van der Pol Reference Chouke, van Meurs and van der Pol1959; Wilson Reference Wilson1975; Paterson Reference Paterson1981). For simplicity, it is assumed that the intruding fluid displaces all of the ambient fluid, though several authors comment that the analysis is easily adapted to the case where a constant-thickness layer of ambient fluid is left behind (cf. Park & Homsy Reference Park and Homsy1984) and coats the cell walls. In the simple case, the jump in the viscous pressure gradient at the interface drives growth of interfacial perturbations at a rate proportional to
$(\mu _a-\mu _i)kV/(\mu _a+\mu _i)$
, where
$\mu _a$
is the ambient viscosity,
$\mu _i$
is the intruding fluid viscosity,
$k$
is the transverse wavenumber, and
$V$
is the displacement velocity. Gravity, surface tension and a radial geometry may provide stabilising effects, but if the viscosity ratio
$m\equiv \mu _a/\mu _i$
exceeds
$1$
, and
$V$
is sufficiently large, then the flow will be unstable to what has become known as Saffman–Taylor fingering. Surface tension does stabilise the short wavelengths, leading to a most-unstable wavelength proportional to, and a growth rate inversely proportional to, the square root of the surface tension. Hence the limit of zero surface tension appears singular, but it can be regularised (e.g. Paterson Reference Paterson1985; Dias & Miranda Reference Dias and Miranda2013; Nagel & Gallaire Reference Nagel and Gallaire2013) by re-including finite-aspect-ratio effects, which lead to a most-unstable wavelength comparable to the cell thickness in accord with observations.
The case of zero surface tension also arises naturally when considering miscible displacement in a Hele-Shaw cell. Wooding (Reference Wooding1969) first described experimentally viscous fingering in a Hele-Shaw cell with miscible fluids similar to that seen with immiscible fluids. Paterson (Reference Paterson1985) presented a stability analysis for an inviscid intrusion spreading from a point source into a miscible viscous ambient with diffusion assumed negligible. Significantly, Paterson neglected radial variations in the thickness of the intrusion by assuming that the ambient leaves at most a thin and immobile film of constant thickness behind the front. His model is then effectively equivalent to the case of complete displacement with immiscible fluids in the zero-surface-tension limit. However, if the intrusion fluid is viscous, then miscible displacement results in an intruding tongue of fluid along the centre of the channel, whose thickness varies with radial position (Petitjeans & Maxworthy Reference Petitjeans and Maxworthy1996; Chen & Meiburg Reference Chen and Meiburg1996; Rakotomalala, Salin & Watzky Reference Rakotomalala, Salin and Watzky1997; Yang & Yortsos Reference Yang and Yortsos1997). The resultant viscosity variation across the cell causes the velocity profile to differ from the simple Poiseuille profile of immiscible displacement (see figure 1). Moreover, at low or moderate Péclet number, the radial and vertical viscosity structure of miscible displacement is also affected by cross-flow diffusion and radial dispersion (e.g. Tan & Homsy Reference Tan and Homsy1987; Goyal & Meiburg Reference Goyal and Meiburg2006; Nijjer, Hewitt & Neufeld Reference Nijjer, Hewitt and Neufeld2018; Videbæk & Nagel Reference Videbæk and Nagel2019; Sharma et al. Reference Sharma, Nand, Pramanik, Chen and Mishra2020), and even the unperturbed base state is time-dependent and must be determined numerically. At large Péclet number, however, it is reasonable to neglect diffusion until the dimensionless radius is comparably large, and it is possible to make more progress analytically.
Figure 1. A radial cross-section of an axisymmetric intrusion with constant influx
$2Q$
into a Hele-Shaw cell with gap thickness
$2h_0$
. The shape of the intrusion is described by the intruding fluid fraction
$\lambda (r,\theta ,t)$
and its radial extent
$r_*(\theta ,t)$
. The viscosities
$\mu _i$
and
$\mu _a$
of the intruding and ambient fluids give rise to a viscosity ratio
$m=\mu _a/\mu _i$
. The densities
$\rho$
are equal. The velocity profile is piecewise parabolic and given by (2.3).
Yang & Yortsos (Reference Yang and Yortsos1997) analysed unidirectional miscible displacement with negligible diffusion and obtained a kinematic-wave equation for the height of the intruding tongue of fluid. For viscosity ratios
$m\lt {3/2}$
, they found a smooth similarity solution (with no shocks) as a function of
$x/t$
. However, for
$m\gt {3/2}$
, the kinematic-wave equation necessarily forms a frontal shock of a height that they recognised, in principle, might require a fully two-dimensional Stokes flow calculation near the nose to determine (cf. Goyal & Meiburg Reference Goyal and Meiburg2006). As discussed further in § 3.1, they instead presented a classical tangent construction of a so-called ‘contact’ shock from the flux function of the intruding fluid. In terms of the intruding fluid fraction
$\lambda _*$
at the nose, the contact-shock height is given by
$\lambda _*=\lambda _c(m)$
, where
\begin{align} \lambda _c= 2\left (\tfrac {2}{3}m-1\right )^{-1/2}\sinh \left [\tfrac {1}{3}\sinh ^{-1}\left \{(m-1)^{-1}\left (\tfrac {2}{3}m-1\right )^{3/2}\right \}\right ] \quad \big(m\gt \tfrac {3}{2} \big) \end{align}
is the real root of a certain cubic polynomial. They note that experiments and numerical models (Petitjeans & Maxworthy Reference Petitjeans and Maxworthy1996; Chen & Meiburg Reference Chen and Meiburg1996; Rakotomalala et al. Reference Rakotomalala, Salin and Watzky1997) suggest that (1.1) underestimates the shock height, particularly for
$m\gt 5$
, and more recent experiments (Bischofberger et al. Reference Bischofberger, Ramachandran and Nagel2014; Videbæk Reference Videbæk2020) support this. Limited data for smaller
$m$
are approximately consistent with (1.1). To the best of our knowledge, the exact nature of these shocks and how their height is determined are not yet fully understood, perhaps because they are experimentally unstable to fingering in the transverse direction.
Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala and Salin1997, Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999, Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos2001) and later Bischofberger et al. (Reference Bischofberger, Ramachandran and Nagel2014) conducted experiments in Hele-Shaw cells with miscible fluids and negligible diffusion, and observed a viscous fingering instability at the front of the intrusion if
$m$
was sufficiently large. The minimum viscosity ratio for instability to be seen was 2–3, a little larger than the critical value
$m={3/2}$
derived by Yang & Yortsos (Reference Yang and Yortsos1997), and clearly larger than the critical value
$m=1$
for the case of immiscible displacement. Both Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala and Salin1997, Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999, Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos2001) and Bischofberger et al. (Reference Bischofberger, Ramachandran and Nagel2014) suggested that the existence of a flat shock front and the associated jump in pressure gradient across the front are crucial to the development of a fingering instability similar to the classical Saffman–Taylor instability. Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos2001) approximated the intrusion as being of uniform thickness equal to the shock height, and by adapting the Saffman and Taylor analysis appropriately, obtained a good prediction of the instability threshold for vertical displacement. Recently, Videbæk (Reference Videbæk2020) adopted the same approach to obtain a good prediction for the instability onset radius for radial flow. (Videbæk also provides an interesting synthesis of experimental observations of immiscible and miscible displacements in linear and radial geometries.)
From the preceding work, it seems to be accepted that
$m={3/2}$
is the predicted stability boundary for miscible intrusion without diffusion. However, to the best of our knowledge, it has not been demonstrated theoretically that the flow is stable for
$m\lt {3/2}$
, and it has been difficult to demonstrate experimentally that the flow is unstable for
$m$
less than approximately 2–3 (Lajeunesse et al. Reference Lajeunesse, Martin, Rakotomalala and Salin1997; Bischofberger et al. Reference Bischofberger, Ramachandran and Nagel2014; Videbæk Reference Videbæk2020). Indeed, Bischofberger et al. (Reference Bischofberger, Ramachandran and Nagel2014) note: ‘It is important to point out that both the connection between the shock-front formation and the onset of the lateral instability, and the suppression of any instability (for example, of the kind from the original Saffman–Taylor analysis) for
$0.67\lt \eta _{in}/\eta _{out}\lt 1$
[i.e.
$1\lt m\lt {3/2}$
] remain to be explained.’ It is our intention to provide some explanation in this paper.
We consider the linear stability of miscible displacement with negligible diffusion from a point source in Hele-Shaw flow, which is parametrised by the viscosity ratio
$m$
between the ambient and intruding fluids. The set-up of the mathematical model, its governing equations, assumptions and non-dimensionalisation are described in § 2. An analytic solution to the initial-value problem is found in § 3 for the special case of axisymmetric flow using the method of characteristics. In the absence of any non-axisymmetric perturbations, this kinematic-wave solution tends towards a simple axisymmetric similarity solution like
$t^{-1}$
. The central point of the paper is a linear stability analysis of this similarity solution in § 4. Working in similarity space, in § 4.1 we derive coupled ordinary differential equations for the radial structure of eigenmodes of specified azimuthal wavenumber, and in § 4.2 we present numerical results for their growth rates as functions of wavenumber and viscosity ratio. Further insight into the structure of the problem and some good asymptotic approximations to the numerical results are obtained in § 4.3 and in Appendix B by analysing the various modes using the Wentzel–Kramers–Brillouin (WKB) method for large azimuthal wavenumber. We confirm stability for
$m\lt {3/2}$
, and explain why instability is hard to observe for
$m$
only somewhat larger than
$3/2$
. We discuss our conclusions in § 5.
2. Model description
2.1. Governing equations
Consider radial flow from a point source in a Hele-Shaw cell consisting of infinite parallel rigid plates separated by a constant distance
$2h_0$
. The cell is initially filled with ambient fluid of viscosity
$\mu _a$
and density
$\rho$
. For
$t\gt 0$
, fluid of viscosity
$\mu _i$
and (equal) density
$\rho$
is injected at the origin at a constant volume flux
$2Q$
. As shown in figure 1, we use cylindrical polar coordinates
$(r,\theta ,z)$
to describe the horizontal extent of the current
$r_{*}(\theta ,t)$
and the vertical thickness of the intrusion
$2h_0\,\lambda (r,\theta ,t)$
, where
$\lambda (r,\theta ,t)$
is the local fluid fraction of injected fluid. Surface tension, diffusion and inertia are all assumed to be negligible.
After an initial transient, the horizontal extent of the intrusion is much greater than its vertical extent,
$r_*\gg h_0$
. In this limit, the vertical velocity is negligible, and the horizontal velocity
$\textbf {u}(r,z,t)$
is related to the horizontal pressure gradient
$\boldsymbol {\nabla }\tilde {p}$
by the lubrication approximation
\begin{align} \mu \frac {\partial {^2\textbf {u}}}{\partial {z^2}}=\boldsymbol {\nabla }\tilde {p} \end{align}
subject to boundary conditions
\begin{align} \textbf {u}=\textbf {0}\text { at }z=\pm h_0,\quad \left [\textbf {u}\right ]_-^+=\textbf {0}\quad\mbox{and}\quad \left [\mu \frac {\partial {\textbf {u}}}{\partial {z}}\right ]_-^+=\textbf {0}\text { at }z=\pm \lambda h_0, \end{align}
which impose no-slip on the boundaries, and velocity and stress continuity at the interfaces between the fluids. Solution of (2.1) and (2.2) yields the velocity profile
\begin{align} \textbf {u}&=\frac {\boldsymbol {\nabla }\tilde {p} }{2\mu _a}\left (mz^2-h_0^2\{1+(m-1)\lambda ^2\}\right )\quad \text {for}\; |z|\lt \lambda h_0, \end{align}
\begin{align} \textbf {u} &=\frac {\boldsymbol {\nabla }\tilde {p} }{2\mu _a}\left (z^2-h_0^2\right ) \;\qquad \qquad\qquad\qquad\ \text {for}\; \lambda h_0\lt |z|\lt h_0, \end{align}
whose shape depends on the viscosity ratio
$m=\mu _a/\mu _i$
and the intruding fluid fraction
$\lambda$
. Integration of (2.3a) between
$\pm \lambda h_0$
gives the horizontal flux
$2\tilde {\textbf {q}}_{i}$
of intruding fluid, while integration of (2.3) between
$\pm h_0$
gives the total flux
$2\tilde {\textbf {q}}$
. We obtain
\begin{align} 2\tilde {\textbf {q}}_i&= -\frac {h_0^3}{3\mu _a} \boldsymbol {\nabla }\tilde {p} \left \{3\lambda +(2m-3)\lambda ^3\right \}, \end{align}
\begin{align} 2\tilde {\textbf {q}}&= -\frac {2h_0^3}{3\mu _a}\boldsymbol {\nabla }\tilde {p} \left \{1+(m-1)\lambda ^3\right \}. \\[-6pt] \nonumber \end{align}
Using these fluxes, we can straightforwardly obtain two local mass-conservation equations for
$r\lt r_*(\theta ,t)$
:
\begin{align} \frac {\partial {\lambda }}{\partial {t}}= \frac {h_0^2}{6\mu _a}\boldsymbol {\nabla }\boldsymbol {\cdot }\left (\boldsymbol {\nabla }\tilde {p} \left \{3\lambda +(2m-3)\lambda ^3\right \}\right ), \end{align}
\begin{align} \boldsymbol {\nabla }\boldsymbol {\cdot }\left (\boldsymbol {\nabla }\tilde {p} \left \{1+(m-1)\lambda ^3\right \}\right )=0. \\[6pt] \nonumber \end{align}
Equation (2.5a) determines the evolution of
$\lambda$
from conservation of intruding fluid, while (2.5b) determines the pressure gradient
$\boldsymbol {\nabla }\tilde {p}$
from a divergence-free constraint on the total flux due to the fixed cell boundaries. Ahead of the intrusion, in
$r\gt r_*(\theta ,t)$
, we have
$\lambda =0$
, thus (2.5a) is not relevant, and (2.5b) reduces to
\begin{align} \nabla ^2\tilde {p} =0. \end{align}
We assumethat there is no imposed far-field pressure gradient, i.e.
$\boldsymbol {\nabla }\tilde {p} \to 0$
as
$r\to \infty$
, and in the absence of surface tension, there is no capillary pressure jump at the nose of the intrusion.
The injection of intruding fluid (only) at the origin at a constant flux
$2Q$
corresponds to the boundary conditions that
$\lambda =1$
at
$r=0$
and
\begin{align} -\frac {mh_0^3}{3\mu _a}\lim _{r\to 0}\left (2\pi r\,\textbf {e}_r\boldsymbol {\cdot }\boldsymbol {\nabla }\tilde {p} \right ) = Q, \end{align}
where
$\textbf {e}_r$
is the radially outward unit vector.
At
$r=r_*(\theta ,t)$
, continuity of pressure and of the total flux normal to the nose yields
\begin{align} \left [\tilde {p} \right ]_-^+&=0, \quad \left [\frac {h_0^2\left (\textbf {n}\boldsymbol {\cdot }\boldsymbol {\nabla }\tilde {p} \right )}{3\mu _a}\left \{1+(m-1)\lambda ^3\right \}\right ]_-^+=0 \quad \text {at}\ r=r_*, \end{align}
where
$\textbf {n}$
is the normal to the perimeter
$r=r_*(\theta ,t)$
of the intrusion. The flux of intruding fluid normal to the nose also has to be consistent with the normal velocity of the nose, which gives
\begin{align} \left (\textbf {n}\boldsymbol {\cdot }\textbf {e}_r\right )\frac {\partial {r_*}}{\partial {t}}=-\frac {h_0^2\left (\textbf {n}\boldsymbol {\cdot }\boldsymbol {\nabla }\tilde {p} \right )}{6\mu _a}\left \{3+(2m-3)\lambda _*^2\right \}, \end{align}
where
$\lambda _*$
is the limiting value of
$\lambda$
as
$r\to r_{*-}$
. For the case of a rounded nose with
$\lambda _*=0$
, (2.8) is just the kinematic condition that the nose moves with the centreline velocity. For the case of a frontal shock with
$\lambda _*\gt 0$
, (2.8) is the condition of mass conservation of intruding fluid across the shock.
Equations (2.5)–(2.8), together with a condition such as (1.1) on any frontal shock height, describe the evolution of the spreading intrusion in terms of the dimensional pressure
$\tilde {p}$
, the intruding fluid fraction
$\lambda$
, and the shape of the perimeter
$r_*$
. These equations for non-axisymmetric flow are equivalent to those of Yang & Yortsos (Reference Yang and Yortsos1997) for unidirectional flow.
2.2. Non-dimensionalisation and similarity variables
The intrusion volume suggests an approximate scaling
$r_*^2h_0\sim Qt$
, and (2.8) suggests
$r_*/t\sim h_0^2\,\tilde {p} /r_* \mu _a$
. More detailed scaling of (2.5a) and (2.6) provides numerical factors and motivates definition of a radial similarity variable
$\xi$
and a dimensionless pressure
$p$
by
\begin{align} \xi =\left (\frac {2\pi h_0}{Q}\right )^{1/2}\frac {r}{t^{1/2}} \quad \text {and}\quad p(\xi ,\theta ,t)=\left (\frac {2\pi h_0^3}{3\mu _aQ}\right )\,\tilde {p} (r,\theta ,t). \end{align}
We also define a mobility function
$\mathcal {M}$
and a flux function
$\mathcal {F}$
by
\begin{align} \mathcal {M}(\lambda ;m)=1+(m-1)\lambda ^3 \quad \text {and}\quad \mathcal {F}(\lambda ;m)=\frac {3\lambda +(2m-3)\lambda ^3}{2+2(m-1)\lambda ^3}, \end{align}
which give the relative mobility for the total flux and the flux fraction of intruding fluid, respectively. As usual for description of evolution towards self-similarity (e.g. Witelski & Bernoff Reference Witelski and Bernoff1999; Leppinen & Lister Reference Leppinen and Lister2003; Mathunjwa & Hogg Reference Mathunjwa and Hogg2006; Peng & Lister Reference Peng and Lister2014), we define a dimensionless time variable by
$\tau =\ln (t/\hat {t}\,)$
, where
$\hat {t}$
is a reference time scale such as
$h_0^3/Q$
.
In terms of the new dimensionless variables, the local mass-conservation equations (2.5) can be written as the coupled partial differential equations
\begin{align} \textbf {q}=-\mathcal {M}\,\boldsymbol {\nabla } p, \quad \boldsymbol {\nabla }\boldsymbol {\cdot }\textbf {q}=0, \quad \frac {\partial {\lambda }}{\partial {\tau }}-\frac {\xi }{2}\,\frac {\partial {\lambda }}{\partial {\xi }}= -\boldsymbol {\nabla }\boldsymbol {\cdot }\left (\mathcal {F}\textbf {q}\right )\quad \text {for}\ \xi \lt \xi _*, \end{align}
\begin{align} \textbf {q}=-\boldsymbol {\nabla } p,\quad \boldsymbol {\nabla }^2p=0 \quad \text {for}\ \xi \gt \xi _*, \\[6pt] \nonumber \end{align}
where
$2\textbf {q}(\xi ,\theta ,\tau )$
is the total flux, and
$\boldsymbol {\nabla }=\textbf {e}_\xi \,\partial /\partial \xi +\textbf {e}_\theta \, \xi ^{-1}\,\partial /\partial \theta$
now denotes the horizontal gradient operator in similarity space. The boundary conditions (2.6)–(2.8) can be written as
\begin{align} \textbf {q}\to \xi ^{-1}\textbf {e}_\xi \ \text {as}\ \xi \to 0, \quad \textbf {q}\to \textbf {0} \ \text {as}\ \xi \to \infty , \end{align}
\begin{align} [p]^+_-=0, \quad \left [\textbf {n}\boldsymbol {\cdot }\textbf {q}\right ]^+_-=0 \quad \text {and}\quad \frac {\partial {\xi _*}}{\partial {\tau }}=\frac {\textbf {n}\boldsymbol {\cdot }\textbf {q}}{\textbf {n} \boldsymbol {\cdot }\textbf {e}_\xi }\,\frac {\mathcal {F}_*}{\lambda _*}-\frac {\xi }{2} \quad \text {at}\ \xi =\xi _*, \\[-6pt] \nonumber \end{align}
where
$\partial \xi _*/\partial \tau$
is the dimensionless speed of the nose in similarity space.
If the evolution of the system becomes self-similar and independent of
$\tau$
at late times, then (2.11) reduces to a system of coupled ordinary differential equations.
3. Axisymmetric flows and similarity solutions
If the flow is axisymmetric, i.e.
$\partial /\partial \theta =0$
, then
$q_\theta =0$
, so
$\boldsymbol {\nabla }\boldsymbol {\cdot }\textbf {q}=0$
gives
$\textbf {q}=\xi ^{-1}\textbf {e}_\xi$
everywhere. Therefore, (2.11c) becomes
\begin{align} \frac {\partial {\lambda }}{\partial {\tau }}=\left (\frac {\xi }{2}-\frac {\mathcal {F}^{\prime}(\lambda )}{\xi }\right )\frac {\partial {\lambda }}{\partial {\xi }}, \end{align}
where
$\mathcal {F}^{\prime}(\lambda )$
denotes
$\textrm {d}\mathcal {F}/\textrm {d}\lambda$
. This equation is a simple quasi-linear hyperbolic equation, which, in the absence of shocks (discontinuities in
$\lambda$
), can be solved analytically by the method of characteristics. (A similar kinematic-wave construction is given by Yang & Yortsos (Reference Yang and Yortsos1997) and Lajeunesse et al. (Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos1999) for related unidirectional flows.)
Equation (3.1) implies that
$\lambda$
is constant along characteristic curves
$\xi (\tau )$
defined by
\begin{align} \frac {\mathrm {d}{\xi }}{\mathrm {d}{\tau }}=\frac {\mathcal {F}^{\prime}(\lambda )}{\xi }-\frac {\xi }{2} \quad \Longleftrightarrow \quad \left (\xi ^2-2\mathcal {F}^{\prime}\right )\mathrm {e}^{\tau }=\mathrm {const}. \end{align}
Thus
$\lambda$
maintains its initial value
$\lambda (\xi _{{init}},0)$
on the characteristic that passes through
$\xi =\xi _{{init}}$
at
$\tau =0$
. Solving (3.2) for
$\xi _{{init}}(\xi ,\tau )$
thus leads to a solution of (3.1) in the form
\begin{align} \lambda (\xi ,\tau )=\lambda \left (\left \{\xi ^2\mathrm {e}^\tau +2\mathcal {F}^{\prime}\left (1-\mathrm {e}^\tau \right )\right \}^{1/2},0\right ). \end{align}
This is implicit in
$\lambda$
as the value of
$\xi _{{init}}(\xi ,\tau )$
depends on
$\lambda$
through
$\mathcal {F}^{\prime}$
, and in general it is not possible to solve (3.3) for
$\lambda (\xi ,\tau )$
explicitly. If
$\lambda$
varies monotonically with
$\xi$
in some region, then it is, however, possible to exploit a change of variable from
$\lambda (\xi ,\tau )$
to
$\xi (\lambda ,\tau )$
to obtain an explicit solution for
$\xi (\lambda ,\tau )$
:
\begin{align} \xi (\lambda ,\tau )=\left \{\xi _{{init}}^2(\lambda )\,\mathrm {e}^{-\tau }+2\mathcal {F}^{\prime}\left (1-\mathrm {e}^{-\tau }\right )\right \}^{1/2}. \end{align}
We will see later that the fluid fraction
$\lambda$
is often a more convenient independent variable than the radial distance
$\xi$
.
Provided that the characteristics do not cross (no shocks form), (3.4) shows that
\begin{align} \xi (\lambda ,\tau )\to X_0(\lambda )\equiv \big \{2\mathcal {F}^{\prime}(\lambda )\big \}^{1/2} \quad \text {as}\ \tau \to \infty . \end{align}
For
$m\leqslant {3/2}$
,
$\mathcal {F}^{\prime}(\lambda )$
is a monotonically decreasing function for
$\lambda \in [0,1]$
, and (3.5) describes the shape
$X_0(\lambda )$
of a long-time similarity solution in which
$\lambda$
varies smoothly from
$\lambda =1$
at
$\xi =0$
to
$\lambda =0$
at
$\xi =X_{0*}=\{2\mathcal {F}^{\prime}(0)\}^{1/2}$
. For
$m\gt {3/2}$
,
$\mathcal {F}^{\prime}(\lambda )$
is an increasing function for
$\lambda$
in a certain range
$[0,\lambda _m]$
, where
$\lambda _m\gt 0$
and
$\mathcal {F}^{\prime\prime}(\lambda _m)=0$
; hence a frontal shock must form by some characteristics for some
$\lambda \in (0,\lambda _m)$
overtaking the characteristic for
$\lambda =0$
, as discussed further below. Neverthless, (3.5) still gives the shape of a long-time similarity solution for
$\lambda \in [\lambda _*,1]$
, where
$\lambda _*$
is the frontal shock height. Similar results were obtained by Yang & Yortsos (Reference Yang and Yortsos1997) for unidirectional flow.
Importantly, as we are interested in the linear stability of radial intrusions into a Hele-Shaw cell, we can expand (3.4) as
$\tau \to \infty$
, to obtain
\begin{align} \xi \sim X_0+\frac {\xi _{{init}}^2-X_0^2}{2X_0}\,\mathrm {e}^{-\tau }+\cdots . \end{align}
A key implication of (3.6) is that any axisymmetric perturbations left over from the initial conditions decay as
$O (\mathrm {e}^{-\tau } )$
, or equivalently,
$O(t^{-1})$
. The decay of all axisymmetric perturbations at the same rate in this problem may be contrasted with perturbations from self-similarity in other problems (see e.g. Witelski & Bernoff Reference Witelski and Bernoff1999; Leppinen & Lister Reference Leppinen and Lister2003; Mathunjwa & Hogg Reference Mathunjwa and Hogg2006) where there is a discrete spectrum of distinct eigenmodes with different decay rates.
3.1. Solutions with shocks
For
$m\gt {3/2}$
,
$\mathcal {F}^{\prime}(\lambda )$
does not vary monotonically with
$\lambda$
, as noted above, so some characteristics overtake others, which leads to interfacial steepening and shock formation – in physical terms, we interpret shocks as regions where
$\lambda (r)$
varies significantly on the short length scale
$h_0$
rather than the long length scale
$r_*$
. (An interpretation that the interface folds over to give a multivalued solution for
$\lambda$
is unphysical as the flow profile (2.3) decreases monotonically away from the maximum velocity on the centreline.) Depending on initial conditions, the shock may initially form in the interior of the flow (Dauck et al. Reference Dauck, Box, Gell, Neufeld and Lister2019), but it will eventually overtake the front and become a frontal shock of some height
$\lambda _*$
. For axisymmetric flow, the frontal condition (2.12e) reduces to
\begin{align} \frac {\mathrm {d}{\xi _*}}{\mathrm {d}{\tau }}=\frac {1}{\xi _*}\frac {\mathcal {F}(\lambda _*)}{\lambda _*}-\frac {\xi _*}{2}. \end{align}
If
$\mathcal {F}(\lambda _*)/\lambda _*\lt \mathcal {F}^{\prime}(\lambda _*)$
, then characteristics continue to overtake the front, and the shock height increases until
$\mathcal {F}(\lambda _*)/\lambda _*\geqslant \mathcal {F}^{\prime}(\lambda _*)$
. There are two cases to consider (see figure 2), which previous work suggests may be relevant for smaller and larger values of
$m$
, respectively.
Figure 2. Two possible axisymmetric similarity solutions with different frontal shock heights for an intrusion with viscosity ratio
$m=10$
. The curved profile for
$\lambda \gt \lambda _*$
is given by (3.5) in both cases. The minimal shock height is
$\lambda _*=\lambda _c\approx 0.34$
for a contact shock (solid line). Also shown is a possible undercompressive shock of height
$\lambda _*=0.55$
(dashed), which travels faster than the characteristics with
$\lambda \gt 0.55$
, but slower than a contact shock.
If the shock height is determined by the information reaching it along characteristics, then
$\lambda _*$
will tend towards the equilibrium height
$\lambda _c$
of a so-called ‘contact’ shock, where
$\mathcal {F}(\lambda _c)/\lambda _c=\mathcal {F}^{\prime}(\lambda _c)$
, with
$\lambda _c$
given by (1.1). A consequence of this condition is that
$\{\mathcal {F}(\lambda _c)/\lambda _c\}^{\prime}=0$
, so the shock speed
$\mathcal {F}(\lambda _*)/\lambda _*$
for small perturbations differs from the equilibrium value only at
$O ((\lambda _*-\lambda _c)^2 )$
. This does not affect the linear behaviour, hence the shock position
$\xi _*$
tends towards its equilibrium position as
$\mathrm {e}^{-\tau }$
just like the characteristics (compare (3.2) and (3.7)).
Alternatively, as seems to be the case for at least
$m\gt 5$
, the shock height is determined by local dynamics on the length scale
$h_0$
of two-dimensional Stokes flow around the front of the intrusion (Yang & Yortsos Reference Yang and Yortsos1997). In this case, we have
$\lambda _*\gt \lambda _c$
,
$\xi _*\lt \xi _c$
and
$\mathcal {F}(\lambda _*)/\lambda _*\gt \mathcal {F}^{\prime}(\lambda _*)$
, and the so-called ‘undercompressive’ shock (see Bertozzi, Münch & Shearer Reference Bertozzi, Münch and Shearer1999) outpaces the characteristics to leave a flat region behind it where
$\lambda =\lambda _*$
(dashed line in figure 2). There is currently no theory for
$\lambda _*$
, but prior work suggests that
$\lambda _*$
increases from about 0.45 to about
$0.6$
as
$m$
increases from 5 to
$\infty$
(Reinelt & Saffman Reference Reinelt and Saffman1985; Rakotomalala et al. Reference Rakotomalala, Salin and Watzky1997; Videbæk Reference Videbæk2020). Since the shock height is determined by local dynamics, it will become constant as the front moves some
$O(1)$
multiple of
$h_0$
, which is on a much shorter time scale than the evolution of the whole flow.
To summarise, we have shown in this section that radial intrusions into a Hele-Shaw cell with or without a shock are stable to axisymmetric perturbations, with all perturbations decaying like
$\mathrm {e}^{-\tau }=\hat {t}/t$
. We have in (3.5) the shape
$X_0(\lambda )$
of a steady axisymmetric similarity solution. We now proceed to the central point of the paper, a linear stability analysis of this base state to determine the growth rate of possible fingering instabilities.
4. Linear stability analysis
4.1. Formulation of the equations
We wish to consider small non-axisymmetric perturbations to the axisymmetric similarity solution of § 3. For simplicity, we will assume that any frontal shock for
$m\gt {3/2}$
is a contact shock, and note that this gives the smallest shock height and plausibly the smallest tendency to instability. We return to the case of undercompressive shocks in § 4.4.
Introducing the function
$\Phi =\xi q_\xi$
for convenience, we can write (2.11) in the form
\begin{align} \Phi =-\mathcal {M}\xi\, {\partial p\over \partial \xi },\quad \xi\, {\partial \Phi \over \partial \xi }=\mathcal {M}\,{\partial ^2 p\over \partial \theta ^2}+O(2), \end{align}
\begin{align} \qquad {\partial \lambda \over \partial \tau }+{2\mathcal {F}^{\prime}\Phi -\xi ^2\over 2\xi }\,{\partial \lambda \over \partial \xi }=O(2)\quad \text {for}\ \xi \lt \xi _*, \\[-6pt] \nonumber \end{align}
where
$O(2)$
denotes terms proportional to
$(\partial \mathcal {M}/ \partial \theta )(\partial p/ \partial \theta )$
in (4.1b) and
$q_\theta (\partial \lambda / \partial \theta )$
in (4.1c) that are both quadratically small in the perturbation and can thus be neglected in a linear analysis. Neglecting the
$O(2)$
term, (4.1c) implies that
$\lambda$
is constant to linear order along radial characteristics defined by
\begin{align} \frac {\mathrm {d}{\xi ^2}}{\mathrm {d}{\tau }}+\xi ^2=2\mathcal {F}^{\prime}\Phi ,\quad \frac {\mathrm {d}{\theta }}{\mathrm {d}{\tau }}=0. \end{align}
The axisymmetric base state is given by
$\xi =X_0(\lambda )$
,
$p=P_0(\lambda )$
and
$\Phi =\Phi _0=1$
, where
\begin{align} X_0=\left \{2\mathcal {F}^{\prime}(\lambda )\right \}^{1/2},\quad \mathcal {M} P^{\prime}_0=-\mathcal {X}(\lambda )\quad \text {and}\quad \mathcal {X}(\lambda )\equiv \frac {X^{\prime}_0}{X_0}=\frac {\mathcal {F}^{\prime\prime}}{2\mathcal {F}^{\prime}}. \end{align}
The base state and the functions
$\mathcal {M}$
,
$\mathcal {F}$
and
$\mathcal {X}$
are all given as functions of
$\lambda$
. Hence it is again convenient to use
$\lambda$
as the independent radial variable in place of
$\xi$
, and to pose the perturbation expansion in
$\xi \lt \xi _*(\theta ,\tau )$
in the form
\begin{align} \xi (\lambda ,\theta ,\tau )&=X_0(\lambda )+X_1(\lambda )\,\mathrm {e}^{{\mathrm {i}}k\theta +\sigma \tau }+\cdots , \end{align}
\begin{align} p(\lambda ,\theta ,\tau )&=P_0(\lambda )+P_1(\lambda )\,\mathrm {e}^{{\mathrm {i}}k\theta +\sigma \tau }+\cdots , \end{align}
\begin{align} \Phi (\lambda ,\theta ,\tau )&=1+\Phi _1(\lambda )\,\mathrm {e}^{{\mathrm {i}}k\theta +\sigma \tau }+\cdots , \\[6pt] \nonumber \end{align}
where
$\sigma$
is the growth rate. The azimuthal wavenumber
$k$
takes integer values for
$2\pi$
-periodicity, but can be treated as a continuous variable for convenience without loss of generality. We neglect terms that are quadratic or higher in the perturbation quantities
$X_1$
,
$P_1$
and
$\Phi _1$
.
Applying the chain rule to the transformation from
$(\xi ,\theta ,\tau )$
to
$(\lambda ,\theta ,\tau )$
, we transform the derivatives in (4.1a,b
) using
\begin{align} \xi\, {\partial \over \partial \xi }=\frac {\xi }{\partial \xi /\partial \lambda }\,{\partial \over \partial \lambda }\quad \text {and}\quad \left ({\partial \over \partial \theta }\right )_\xi =\left ({\partial \over \partial \theta }\right )_\lambda - \frac {\partial \xi /\partial \theta }{\partial \xi /\partial \lambda }\,{\partial \over \partial \lambda }\,. \end{align}
We then substitute the expansion (4.4) into (4.1a,b ) and (4.2), linearise the result, and use (4.3) to simplify the equations further. After some algebra, we obtain
\begin{align} \mathcal {X}\Phi _1+\left ({X_1\over X_0}\right )^{\prime} =-\mathcal {M} P^{\prime}_1,\quad {\Phi^{\prime}_1}=-k^2 \mathcal {X}\left (\mathcal {M} P_1+{X_1\over X_0} \right ) \end{align}
and
\begin{align} 2(\sigma +1){X_1\over X_0}=\Phi _1\,. \end{align}
The special case
$\sigma =-1$
provides stable perturbations, such as the axisymmetric perturbations of § 3, and it will not be considered further. If
$\sigma \neq -1$
, then we can use (4.7) to eliminate
$X_1/X_0$
from (4.6) and obtain the coupled system
\begin{align} \begin{pmatrix} \mathcal {M} P^{\prime}_1 \\ \Phi^{\prime}_1 \end{pmatrix} = \mathcal {X} \begin{pmatrix} \dfrac {k^2}{2(1+\sigma )} & \dfrac {k^2}{4(1+\sigma )^2}-1 \\ -k^2 & -\dfrac {k^2}{2(1+\sigma )} \end{pmatrix} \boldsymbol {\cdot } \begin{pmatrix} \mathcal {M} P_1 \\ \Phi _1 \end{pmatrix}. \end{align}
As shown in Appendix A, the general boundary conditions (2.12) reduce to the boundary conditions
\begin{align} \frac {P_1}{\Phi _1}=\frac {1}{k}-\frac {1}{2(1+\sigma )} \ \text {at}\ \lambda =\lambda _* ,\quad \Phi _1\to 0\ \text {as}\ \lambda \to 1 \, \end{align}
on (4.8a). Equations (4.8) form a linear homogeneous system, which constitutes an eigenvalue problem to determine the growth rates
$\sigma (k;m)$
of perturbations with radial structure given by eigenfunctions
$P_1$
and
$\Phi _1$
. It can be solved numerically in this form.
Alternatively, we can eliminate
$P_1$
to obtain a second-order equation for
$\Phi _1$
:
\begin{align} \mathcal {X}\mathcal {M}\left ({\Phi ^{\prime}_1\over \mathcal {X}\mathcal {M}}\right )^{\prime}= k^2 \mathcal {X}^2\left ( 1+{\mathcal {N}\over \sigma +1} \right )\Phi _1 ,\quad \text {where }\mathcal {N}(\lambda )\equiv {\mathcal {M}^{\prime}\over 2\mathcal {M}\mathcal {X}} , \end{align}
with boundary conditions
\begin{align} \frac {\Phi^{\prime}_1}{k^2\mathcal {X}_*}=\left (\frac {\mathcal {M}_*-1}{2(1+\sigma )}-\frac {\mathcal {M}_*}{k}\right )\Phi _1 \ \text {at}\ \lambda =\lambda _*,\quad \Phi _1\to 0 \ \text {as}\ \lambda \to 1. \end{align}
This second-order form is convenient for WKB analysis of the limit
$k\to \infty$
.
Though it is slightly unusual for the eigenvalue
$\sigma$
to appear in the boundary condition (4.9b) as well as the differential equation (4.9a), this second-order form is sufficiently close to a standard Sturm–Liouville eigenvalue problem to expect, as proves to be the case (see Appendix C), that there is a discrete spectrum of eigenmodes for each
$k$
and
$m$
. We label these modes by an integer
$n\geqslant 0$
, which is equal to the number of zeros of the eigenfunction away from the zero boundary condition at
$\lambda =1$
, or equivalently, at
$\xi =0$
(see figure 5, for example).
4.2. Numerical solution and results
We solved the boundary-value problem (4.8) numerically using continuation methods implemented with the software package AUTO-07P (freely available at http://indy.cs.concordia.ca/auto). The strategy for obtaining an eigenmode is analogous to that detailed in Ribe, Lister & Chiu-Webster (Reference Ribe, Lister and Chiu-Webster2006): start with a guess for
$\sigma$
and a non-zero solution of (4.8a) that satisfies one of the homogeneous boundary conditions (4.8b), but will not, in general, satisfy the other; then use continuation to slowly impose the other boundary condition, keeping the solution non-zero and allowing
$\sigma$
to vary slowly until it reaches an eigenvalue at the point where the second boundary condition is satisfied. Having obtained the eigenmode for one set of parameters, continuation can again be used to track its variation with
$k$
and
$m$
. We present results for the first three eigenmodes,
$n=0,1,2$
, but, given the discrete nature of the spectrum, it is easy to find starting values of
$\sigma$
that give the higher eigenmodes.
Figure 3 shows the analytical base-state profiles (4.3) for various values of
$m$
. For
$m\lt {3/ 2}$
, the profile has a rounded nose with the tip position at
$\xi _*=\sqrt 3$
, as determined by a combination of the centreline velocity for
$\lambda =0$
and radial spreading. For
$m\ll 1$
, the more viscous intruding fluid is lubricated by the less viscous ambient fluid near the origin, and the intrusion there is wider and closer to plug flow than for
$m=1$
. For
$1\lt m\lt {3/2}$
, the lower viscosity intrusion is narrower near the origin, and wider near the rounded nose than for
$m=1$
. For
$m\gt {3/2}$
, there is a frontal shock, which we are assuming has height given by (1.1).
Figure 3. The axisymmetric base-state profiles
$\xi =X_0(\lambda )$
, as given by (4.3), of self-similar solutions for intrusions with viscosity ratios
$m\in \{0.15,0.5,1.25,5\}$
. For
$m=5$
, there is a frontal shock at
$\xi _*=1.815$
of height
$\lambda _*=0.354$
, rather than the unphysical non-monotonic profile (dashed) that would be predicted by ignoring the crossing of characteristics.
In the results below, we will use viscosity ratios
$m\in \{0.15,1.25,5\}$
as illustrative examples of the stability behaviours found in the three distinct cases: a more viscous intrusion (
$m\lt 1$
), a less viscous intrusion without a shock (
$1\lt m\lt {3/2}$
), and a less viscous intrusion with a shock (
$m\gt {3/2}$
). We observe briefly that
$m=1$
(equal viscosities) is a very special case as the lack of any viscosity differences means that the flow is always radial with flux
$\textbf {q}=\xi ^{-1}\textbf {e}_\xi$
and a parabolic (Poiseuille) profile. Hence there is no perturbation flow (
$P_1=\Phi _1=0$
), the interface is simply a passive tracer in the base-state flow, and any perturbations to
$X_0(\lambda )$
decay purely kinematically with
$\sigma =-1$
.
Figure 4. The growth rates
$\sigma$
corresponding to the first three eigenmodes
$n\in \{0,1,2\}$
with viscosity ratios
$m\in \{0.15,1.25,5\}$
as functions of the azimuthal wavenumber
$k$
. For
$m=5$
, the fundamental mode
$n=0$
is unstable if
$k$
exceeds a critical value
$\approx 18$
where
$\sigma =0$
(blue dot).
Figure 4 shows the growth rates
$\sigma$
of the first three eigenmodes
$n\in \{0,1,2\}$
for the three illustrative viscosity ratios as functions of the azimuthal wavenumber
$k$
. Some general observations can be understood physically. First, for all the eigenmodes,
$\sigma \to -1$
as
$k\to 0$
, which reflects the result
$\sigma =-1$
for axisymmetric perturbations in § 3. (Recall that we are treating
$k$
as a continuous variable for convenience, rather than imposing
$2\pi$
-periodicity.) A second, related observation is that in each graph,
$\sigma$
becomes closer to
$-1$
as
$n$
increases. Increasing
$n$
corresponds to increasing the number of zeros in the eigenfunctions and hence to increasing the amount of radial structure. We can reasonably expect that as
$n\to \infty$
for fixed
$k$
, the radial structure dominates the azimuthal variation, and therefore the growth rate again approaches the
$\sigma =-1$
result for axisymmetric perturbations. Third,
$\sigma \lt -1$
for
$m\lt 1$
, and
$\sigma \gt -1$
for
$m\gt 1$
. This is consistent with the observation that
$\sigma =-1$
for all perturbations when
$m=1$
, and is also consistent with intuition derived from the Saffman–Taylor instability mechanism that pushing a more viscous fluid into a less viscous fluid tends to be stable, whereas pushing a less viscous fluid into a more viscous fluid tends to promote instability. Nevertheless,
$m=1$
is definitely stable (
$\sigma =-1$
), so
$m\gt 1$
is not sufficient to produce instability, as can been seen, for example, in figure 4 for
$m=1.25$
.
For
$m=0.15$
, all modes are stable with
$\sigma \lt -1$
, the fundamental mode
$n=0$
is the most stable, and as
$k$
increases the perturbations become more stable. For
$m=1.25$
, all modes are again stable, but with
$-1\lt \sigma \lt 0$
, the fundamental mode
$n=0$
is the least stable, and though the perturbations become less stable as
$k$
increases, the growth rates appear to tend to a limit that is still negative as
$k\to \infty$
. For
$m=5$
, which has a base state with a shock, modes
$n=1,2$
are again stable with
$-1\lt \sigma \lt 0$
. The crucial difference for
$m=5$
is that the fundamental mode
$n=0$
becomes unstable at
$k\approx 18$
, and the growth rate increases rapidly as
$k\to \infty$
. In § 4.3, we show that the instability mechanism is essentially the Saffman–Taylor mechanism acting on the jump in mobility at the frontal shock.
Figure 5 shows the radial structure of the first three eigenmodes of the pressure perturbation
$P_1$
and flux perturbation
$\Phi _1$
as functions of the radial similarity variable
$\xi$
. The plots show solutions for the three illustrative viscosity ratios and for four azimuthal wavenumbers
$k\in \{1,5,25,100\}$
. As expected, the number of zeros (additional to
$\xi =0$
) increases with mode number
$n$
. We note that
$\Phi _1$
is approximately in phase with
$P_1$
for
$m=0.15$
, but has approximately the opposite phase (sign) for
$m=1.25$
and
$m=5$
. From the form of the matrix in (4.8a), it can be inferred that this is largely a consequence of the sign of the factor
$1/(\sigma +1)$
, which is different for
$m\gt 1$
and
$m\lt 1$
.
Figure 5. Numerical solutions for the first three eigenmodes
$n\in \{0,1,2\}$
in terms of the perturbation pressure
$P_1$
(solid) and perturbation flux
$\Phi _1$
(dotted) for wavenumbers
$k\in \{1,5,25,100\}$
and viscosity ratios
$m\in \{0.15,1.25,5\}$
. The nose position is at
$\xi _*=\sqrt {3}$
for
$m\in \{0.15,1.25\}$
, and at
$\xi _*\approx 1.81$
for
$m=5$
.
For
$k=1$
(sideways displacement, perturbations
$\propto \cos \theta$
) the figure shows that the eigenmodes giving relaxation back to axisymmetry have a comparable length scale to the full extent of the intrusion. As
$k$
increases, the eigenmodes become increasingly localised radially. In Appendix B, we show that as
$k\to \infty$
for
$m\lt {3/2}$
, the eigenmodes become localised about an interior position between the origin and the nose; for
$m\gt {3/2}$
, the eigenmodes become localised near the frontal shock.
Figure 4 showed that for
$m=5$
, the fundamental mode is unstable for
$k\gtrsim 18$
. Figure 6 extends this result by showing the regions in the
$(k,m)$
-plane where
$\sigma \gt 0$
or
$\sigma \lt 0$
, and the curve of marginal stability where
$\sigma =0$
. For
$m\gt {3/ 2}$
, the fundamental mode is always unstable for sufficiently large
$k$
, while for
$m\lt {3/2}$
, it is stable for all
$k$
. In particular, for
$1\lt m\lt {3/2}$
, the flow is stable to all perturbations, which agrees with experimental observations of stability in this regime, but contrasts with instability in the classical Saffman–Taylor problem for
$m\gt 1$
.
Figure 6. The marginal-stability contour
$\sigma =0$
for the fundamental mode
$n=0$
in the
$(m,k)$
-plane compared to the asymptotic result (4.16) for
$k\gg 1$
.
As
$m$
decreases towards
$3/2$
, these calculations show that the flow is only unstable to very large
$k$
, for example
$k\gt 10^3$
for
$m\lt 1.75$
. However, if
$k$
is too large (
$k\gg r_*/h_0$
), then the horizontal length scale
$r_*/k$
of perturbations near the front is less than the thickness
$2h_0$
of the Hele-Shaw cell, and the horizontal viscous stresses, which are neglected in the lubrication model (2.1), will stabilise the flow and provide a large wavenumber cut-off (cf. Paterson Reference Paterson1985). It follows that as
$m$
decreases towards
$3/ 2$
, the instability would manifest only at very large
$r_*/h_0$
, and would thus be difficult to observe experimentally in practice.
Figure 6 also shows asymptotic results for the large-wavenumber limit
$k\to \infty$
, which are derived in the next subsection. The excellent agreement with the numerical results supports the accuracy of the calculations.
4.3. Analysis of the limit
$k\gg 1$
for
$m\gt {3/2}$
and
$n=0$
Numerically, instability occurs only for
$n=0$
, sufficiently large
$k$
and
$m\gt {3/ 2}$
(figures 4 and 6). We now pursue an analysis of the limit
$k\to \infty$
to confirm these results and examine the nature of the instability as
$m\to {3/2}$
. Numerically, it appears that
$\sigma \sim k$
as
$k\to \infty$
, and we note that this is also true for the Saffman–Taylor instability in the limit of vanishing surface tension.
For convenience of notation, we define
\begin{align} T={k\over 2(1+\sigma )}\,, \end{align}
which we anticipate will be an
$O(1)$
quantity if
$\sigma \sim k$
. Equations (4.9a,b) can then be written in the form
\begin{align} \Phi^{\prime\prime}_1- k^2 \mathcal {X}^2\Phi _1 = \left ({\mathcal {X}^{\prime}\over \mathcal {X}}+{\mathcal {M}^{\prime}\over \mathcal {M}}\right )\Phi^{\prime}_1 +k\mathcal {X}\, {T\mathcal {M}^{\prime}\over \mathcal {M}}\, \Phi _1 , \end{align}
\begin{align} \frac {\Phi^{\prime}_1}{k\mathcal {X}_*}=\left (T\mathcal {M}_*-T-\mathcal {M}_*\right )\Phi _1 \ \text {at}\ \lambda =\lambda _*,\quad \Phi _1\to 0 \ \text {as}\ \lambda \to 1. \\[-6pt] \nonumber \end{align}
We make the usual WKB assumption that
$\Phi _1=\exp [S(\lambda )]$
with
$S=kS_0+S_1+O(k^{-1})$
, and substitute into (4.11a) to obtain
\begin{align} k^2S^{\prime}_0{}^2- k^2 \mathcal {X}^2+ kS^{\prime\prime}_0 +2kS^{\prime}_0S^{\prime}_1=kS^{\prime}_0\left ({\mathcal {X}^{\prime}\over \mathcal {X}}+{\mathcal {M}^{\prime}\over \mathcal {M}}\right ) + k\mathcal {X}\, {T\mathcal {M}^{\prime}\over \mathcal {M}} +O(1). \end{align}
At
$O(k^2)$
, we have the decaying solution
$S^{\prime}_0=+\mathcal {X}$
. (Note that
$\mathcal {X}\lt 0$
, so this solution has
$\Phi _1\to 0$
as
$\lambda \to 1$
.) This also implies that
$kS^{\prime\prime}_0 =kS^{\prime}_0\mathcal {X}^{\prime}/\mathcal {X}$
in (4.12). At
$O(k)$
, the remaining terms simplify to
\begin{align} 2S^{\prime}_1= (1+T) {\mathcal {M}^{\prime}\over \mathcal {M}} \quad \Rightarrow \quad S_1= {1+T\over 2} \ln \mathcal {M} +c\,. \end{align}
We note that we can substitute the solutions for
$S_0$
and
$S_1$
into the WKB ansatz to obtain a uniformly asymptotic expression
$\Phi _1=A\mathcal {M}^{(1+T)/2} (X_0/X_{0*})^k$
for the structure of the fundamental eigenmode, which includes the effect of radial mobility variations.
The boundary condition (4.11b) at
$\lambda =\lambda _*$
yields
\begin{align} \frac {S^{\prime}_0}{\mathcal {X}_*}+{S^{\prime}_1\over k\mathcal {X}_*}+O(k^{-2})=T\mathcal {M}_*-T-\mathcal {M}_* \,, \end{align}
where
$S^{\prime}_0$
and
$S^{\prime}_1$
are now known. This equation can be rearranged using (4.10) to give the desired asymptotic result for the growth rate:
\begin{align} \sigma (k;m)\sim {k\over 2}{\mathcal {M}_*-1\over \mathcal {M}_*+1}-1+ {1\over 2|\mathcal {X}_*|}{\mathcal {M}^{\prime}_*\over (\mathcal {M}_*+1)^2}+O(k^{-1}) \,. \end{align}
As shown in figure 7, (4.15) gives a very good approximation to the full numerical result, even at moderate
$k$
. The first term agrees with a simple analysis of Saffman–Taylor instability at the front of an intrusion of uniform thickness and hence uniform mobility
$\mathcal {M}_*\gt 1$
(see e.g. Saffman & Taylor Reference Saffman and Taylor1958; Lajeunesse et al. Reference Lajeunesse, Martin, Rakotomalala, Salin and Yortsos2001; Videbæk Reference Videbæk2020). The second term,
$-1$
, reflects the stabilising effect of radial geometry (cf. Wilson Reference Wilson1975; Paterson Reference Paterson1981, Reference Paterson1985).
Figure 7. Comparison of the numerically computed growth rate
$\sigma (k)$
for
$n\in \{0,1,2\}$
and
$m=5$
(solid lines) with the corresponding WKB solution (dashed lines). Note that
$\sigma$
changes sign for
$n=0$
at
$k\approx 18$
.
The third term, involving
$\mathcal {M}^{\prime}_*$
, is new and describes the effects of the base-state variation of the intrusion thickness away from the front. The term is positive as the mobility
$\mathcal {M}$
is an increasing function of
$\lambda$
for
$m\gt 1$
, so
$\mathcal {M}^{\prime}_*\gt 0$
. Hence it represents an additional destabilisation relative to the Saffman–Taylor result for an intrusion of uniform thickness in radial geometry. This may be understood physically as the effect of greater mobility (
$\mathcal {M}\gt \mathcal {M}_*$
) behind the frontal shock, which facilitates the growth of perturbations. The size of the third term varies only slowly with
$m$
from
$1/4$
at
$m={3/2}$
to
$3/{16}$
as
$m\to \infty$
.
By setting
$\sigma =0$
in (4.15), we can also obtain an asymptotic estimate for the curve of marginal stability. We substitute for
$\mathcal {M}_*$
and
$\mathcal {X}_*$
from (2.10) and (4.3), and rearrange (4.15) to obtain
\begin{align} k(m)\approx \frac {3}{(m-1)\lambda _*^3}+\frac {3+2(m-1)\lambda _*^3}{2+(m-1)\lambda _*^3}\,. \end{align}
This result is asymptotic as
$m\to ({3/ 2})_+$
, since this gives
$\lambda _*\to 0$
and
$k\to \infty$
, and it agrees remarkably well with the numerically calculated curve of marginal stability (figure 6), even for relatively small wavenumber
$k$
.
The WKB analysis for
$n=0$
with
$m\gt {3/2}$
has thus provided excellent confirmation of the numerical results, and gives asymptotic expressions as
$k\to \infty$
for the unstable eigenmode, growth rate and marginal stability curve that include the leading-order effects of the radial variation in the mobility
$\mathcal {M}$
of the base state. In Appendix B, we analyse the limit
$k\to \infty$
for
$n\gt 0$
with
$m\gt {3/2}$
, and for
$n\geqslant 0$
with
$m\lt {3/2}$
, which again provides good confirmation of the numerical results. For
$1\lt m\lt {3/2}$
, we find that
$-1\lt \sigma \lt -{3/4}$
.
4.4. Undercompressive shocks
A formal stability analysis of undercompressive shocks is somewhat complicated by the need to split the intrusion into a central region where
$\lambda \gt \lambda _*$
that is reached by characteristics from the origin, and an annular region where
$\lambda =\lambda _*$
that lies between the shock and the central region. The equations in the annular region are similar to those ahead of the shock, but with mobility
$\mathcal {M}_*$
. As noted earlier, there is currently no theory for
$\lambda _*$
. Fortunately, we can, nevertheless, derive a simple result from the large-
$k$
analysis in the previous section.
Recalling that
$\mathcal {X}=X^{\prime}_0/X_0$
, we can rewrite
$\mathcal {M}^{\prime}/\mathcal {X}$
in (4.15) as
$X_0\,\textrm {d}\mathcal {M}/\textrm {d}X_0$
. For an undercompressive shock, this derivative is zero in the annular region of constant
$\lambda$
, thus
\begin{align} \sigma (k;m)\sim {k\over 2}{\mathcal {M}_*-1\over \mathcal {M}_*+1}-1 \,. \end{align}
This result is the same as for an entirely uniform thickness intrusion (Paterson Reference Paterson1985; Videbæk Reference Videbæk2020). It appears here as the asymptotic result, despite the non-uniform central region, because the asymptotic radial eigenfunction
$\Phi _0\propto (X_0/X_{0*})^k$
is concentrated near the front, and lies in the uniform annular region (cf. the final panel of figure 5).
The loss of the third term from (4.15) makes the disturbance more stable, but the fact that
$\lambda _*\gt \lambda _c$
for an undercompressive shock makes the disturbance less stable. The net effect depends on
$\lambda _*$
. However, provided that
$\lambda _*\to 0$
as
$m\to {3/2}$
, the marginal stability estimate from (4.17) also has
$k\to \infty$
.
5. Discussion and conclusions
Using lubrication theory, we have derived the equations governing the shape of a fluid tongue intruding from a point source into a Hele-Shaw cell filled with another fluid of the same density but differing viscosity, neglecting both diffusion and surface tension. For the case of perfectly axisymmetric flow, changing variables to the relative fluid fraction
$\lambda$
allowed an explicit solution from initial conditions, from which we could show that the initial-value problem approaches a late-time similarity solution with all relative perturbations decaying like
$t^{-1}$
. This similarity solution is the axisymmetric equivalent of that found by Yang & Yortsos (Reference Yang and Yortsos1997) for rectilinear flow, and it grows like
$t^{1/2}$
instead of
$t$
. As in Yang & Yortsos (Reference Yang and Yortsos1997), if the viscosity ratio
$m$
exceeds
$3/2$
, then there must be a frontal shock of height
$\lambda _*\geqslant \lambda _c$
.
Our linear stability analysis of non-axisymmetric perturbations to this self-similar base state has provided clear theoretical confirmation of the previously stated hypothesis that the frontal shock is crucial for the development of instability in miscible Hele-Shaw displacement with negligible diffusion: axisymmetric flows with viscosity ratio
$m\lt {3/2}$
, and hence without a shock, were shown to be stable; axisymmetric flows with viscosity ratio
$m\gt {3/2}$
, and hence with a shock, were shown to be unstable. In particular, this means that intrusions of less viscous fluid into a more viscous ambient can be stable, provided that the viscosity ratio does not exceed
$3/2$
.
The analysis gives an indication why the critical viscosity ratio is
$m={3/2}$
and not
$m=1$
as for immiscible displacement. For miscible displacement with
$m=1$
, the fluid interface is just a passive tracer, which is advected by the base-state radial Poiseuille flow. The kinematics are the same as described in § 3, and give stable behaviour with
$\sigma =-1$
. Crucially, the radial gradient of the mobility
$\mathcal {M}$
remains finite for
$1\leqslant m\lt {3/2}$
, so
$\sigma$
varies continuously from
$\sigma =-1$
as
$m$
varies from 1. As might be expected physically, the flow does becomes less stable as
$m$
increases, but as shown in Appendix B,
$\sigma$
actually remains less than
$-{3/4}$
for any
$m\lt {3/2}$
, hence the flow remains stable. (By contrast, immiscible displacement has a sharp jump in mobility at the front for
$m\gt 1$
, and this causes a significant change in the behaviour of
$\sigma$
at large
$k$
, which drives the Saffman–Taylor instability.)
For
$m\gt {3/2}$
, there must be a frontal shock and a mobility jump, and we find the expected Saffman–Taylor-like instability. Our analysis puts this expectation on a firm theoretical footing by including the effects of the radial variation of the base state. There is a discrete spectrum of radial eigenmodes with differing numbers of zeros, and only the fundamental mode
$n=0$
can be unstable. It is unstable for sufficiently large azimuthal wavenumber
$k$
and, within the lubrication model, is most unstable as
$k\to \infty$
. This divergence for very large
$k$
could be regularised (see e.g. Paterson Reference Paterson1985; Dias & Miranda Reference Dias and Miranda2013; Nagel & Gallaire Reference Nagel and Gallaire2013) by re-including neglected horizontal stresses, which would stabilise wavelengths less than the cell thickness and predict the most unstable wavelength to be a few times the cell thickness.
Comparison of wavelengths to the cell thickness is also essential for prediction of the onset radius at which the instability becomes manifest in experiments, particularly as
$m$
decreases towards
$3/2$
. As shown in figure 6, the sufficiently large
$k$
for instability is predicted to increase rapidly as
$m$
decreases, and it diverges as
$m\to {3/2}$
. For example,
$k\gt 10^3$
is necessary for instability with
$m=1.75$
. For the unstable wavelengths to exceed the cell thickness, one also needs
$r_*\gt h_0 k$
, making the instability more difficult to observe in a given experimental cell. The predicted divergence in the onset radius is consistent with experimental data in Videbæk (Reference Videbæk2020) and with the difficulty in observing instability for
$m$
less than approximately 2–3 (Lajeunesse et al. Reference Lajeunesse, Martin, Rakotomalala and Salin1997; Bischofberger et al. Reference Bischofberger, Ramachandran and Nagel2014). (A further contributing factor may be the time take to form the frontal shock structure: for
$m=2$
, the velocity of even the fastest characteristic exceeds the centreline velocity by less than 1 %.)
In addition to the numerical results, we also found asymptotic solutions to the linear stability problem by using a WKB analysis for large
$k$
. These confirm the numerical results and provide useful analytic expressions for the growth rates, for example, the remarkably good marginal stability curve (4.16) for a contact shock or the asymptotic growth rate (4.17) of an undercompressive shock.
It would be nice to have a better understanding which of the two theoretical shock structures applies for
$m\lt 5$
, and why. We note, however, both from above and from Goyal & Meiburg (Reference Goyal and Meiburg2006), that for
$m\lt 2{-}3$
, the front is likely very slow to become quasi-steady (even in a tube where it is stable), and that a very large, or infinite, Péclet number may be required for diffusion to remain negligible. Full direct numerical simulations or experiments could be challenging. The effects of diffusion (physical or numerical) increase with radial distance in an axisymmetric geometry since the velocity
$\dot r_*$
of the front decreases like
$t^{-1/2}$
. Goyal & Meiburg (Reference Goyal and Meiburg2006) studied diffusive effects on instability numerically in a linear geometry, and found that the most unstable wavelength and the growth rate increase only slightly, by of order 10 %, for
$m\gt 7$
as the Péclet number
$2h_0 \dot x_*/D$
increased from 500to 2000. Similarly, Videbæk & Nagel (Reference Videbæk and Nagel2019) found similarly in radial (and linear) experiments that the most unstable wavelength and also the onset radius showed no discernible trend for
$m=5$
as the Péclet number at onset ranged from 1000 to 20 000, but also that a quite different fingering instability is seen at lower Péclet numbers. The boundary between their low- and high-Péclet-number regimes increases as
$m$
decreases, which might be linked to the increase in the marginally stable wavenumber, and hence onset radius, shown here (figure 6) as
$m$
decreases towards the critical value
$m={3/2}$
. While questions remain over the detailed effects of diffusion, these studies are at least indicative of a high-Péclet-number regime in which diffusion can be neglected. Perhaps the shock structure could be investigated by Stokes flow calculations in this limit.
We have concentrated here on the case of radial flow, in part because of its experimental relevance, and in part because the radial geometry allows fully separable eigenfunctions with a certain radial structure in similarity space, fixed azimuthal wavenumber and power-law time dependence. (Analysis of unidirectional flow is significantly more difficult.) To summarise our main conclusions: the case
$m = 1$
is kinematically stable with algebraic decay like
$t^{-1}$
; intrusions with
$1 \lt m \lt {3/2}$
are stable, though less so; intrusions with
$m \gt {3/ 2}$
are unstable for sufficiently large wavenumber due to the jump in mobility, and hence pressure gradient, at the frontal shock; the instability is hard to observe for
$m$
only a little larger than
$3/2$
.
Acknowledgements.
The authors are grateful to F. Box for the background experimental image in the graphical abstract.
Funding.
T.-F.D. was supported by an Engineering and Physical Sciences Research Council PhD studentship.
Declaration of interests.
The authors report no conflict of interest.
Data availability statement.
The research data supporting the findings of this study are available within the paper.
Appendix A. Boundary conditions on the perturbations
The perturbation equations (4.8a) are a second-order system of linear ordinary differential equations, hence we expect two boundary conditions. Physically, these are given by a matching condition to the flux and pressure distribution ahead of the nose, and a condition on there being no perturbation flux at the origin. Mathematically, we solve the system ahead of the nose analytically to find a matching condition via (2.12b–d), and we use a local expansion near the origin to find an appropriate boundary condition from (2.12a).
Ahead of the nose (
$\xi \gt \xi _*$
), we have
$\lambda \equiv 0$
and must pose a perturbation expansion of the form
\begin{align} p(\xi ,\theta ,\tau )&=p_0(\xi )+p_1(\xi )\,\mathrm {e}^{\mathrm {i} k\theta +\sigma \tau }+\cdots , \end{align}
\begin{align} \Phi (\xi ,\theta ,\tau )&=1+\phi _1(\xi )\,\mathrm {e}^{\mathrm {i} k\theta +\sigma \tau }+\cdots , \\[6pt] \nonumber \end{align}
rather than (4.4). Since
$\lambda =0$
, we have
$\mathcal {M}=1$
,
$\textbf {q}=-\boldsymbol {\nabla } p$
and
$\nabla ^2 p=0$
. The base state
$\Phi =1$
gives
$p_0=-\ln \xi$
to within an unimportant additive constant. Since the perturbation term in (A1a) also satisfies Laplace’s equation, we must have
$p_1\propto \xi ^{\pm k}$
. Assuming
$k\gt 0$
for definiteness, the condition (2.12a) of decaying flow in the farfield rules out
$\xi ^{+k}$
, hence we obtain
\begin{align} p_1(\xi )=A \xi ^{-k}\quad \text {and}\quad \phi _1(\xi )=A k\xi ^{-k} \,, \end{align}
where
$A$
is the (small) perturbation amplitude ahead of the nose.
Since
$q_\theta n_\theta$
is quadratically small in the perturbation, the frontal boundary conditions (2.12c) and (2.12 d) reduce to
$\Phi$
and
$p$
being continuous across
$\xi =\xi _*$
. As discussed in § 3.1, contact shocks tend to their equilibrium height like
$\mathrm {e}^{-\tau }$
. Hence for
$\sigma \ne -1$
, we can assume that the frontal height
$\lambda _*$
is not perturbed, so the frontal position from the expansion (4.4a) is
$\xi _*=X_0(\lambda _*)+X_1(\lambda _*)\,\mathrm {e}^{\mathrm {i} k\theta +\sigma \tau }$
. We substitute this position into (A1) and equate the linearised results to the expansions (4.4b,c) to obtain
\begin{align} P_1=p_1(X_0)+X_1\,{\textrm {d}p_0\over \textrm {d}\xi }=p_1(X_0)- {X_1\over X_0} \quad \text {and}\quad \Phi _1=\phi _1(X_0) \quad \text {at}\ \lambda =\lambda _*. \end{align}
We use (4.7) to replace
$X_1/X_0$
by
$\Phi _1/2(1+\sigma )$
. We then divide (A3a) by (A3b) to obtain the desired matching condition
\begin{align} \frac {P_1}{\Phi _1}=\frac {1}{k}-\frac {1}{2(1+\sigma )} \quad \text {at}\ \lambda =\lambda _*, \end{align}
which ensures that the solution to (4.8a) in
$x\lt x_*$
can be matched to the decaying solution (A2) in
$x\gt x_*$
.
At the origin
$\xi =0$
, we have
$\lambda =1$
. Expanding (3.5) about this point yields
$\lambda \sim 1-{1\over 6}m\xi ^2$
as
$\xi \to 0$
, thus
$\mathcal {M}=m+O(\xi ^2)$
. We can thus approximate (4.1) by
\begin{align} \Phi \sim -m\xi\, {\partial p\over \partial \xi }\quad \text {and}\quad \xi\, {\partial \Phi \over \partial \xi }\sim m\,{\partial ^2 p\over \partial \theta ^2}\quad \text {as}\ \xi \to 0. \end{align}
Hence
$p$
again satisfies Laplace’s equation at leading order, and substitution of the expansion (A1) again leads to
$p_1\propto \xi ^{\pm k}$
. This time it is the singular solution
$\xi ^{-k}$
that is ruled out by the origin condition (2.12a), hence we obtain
\begin{align} p_1(\xi )\sim B \xi ^{k}\quad \text {and}\quad \phi _1(\xi )\sim -B m k\xi ^{k} \quad \text {as}\ \xi \to 0, \end{align}
where
$B$
is an amplitude. Again by equating the expansions (4.4) and (A1), this time as
$\xi \to 0$
or
$\lambda \to 1$
, we obtain
\begin{align} P_1=p_1(X_0)- {X_1\over m X_0} \quad \text {and}\quad \Phi _1=\phi _1(X_0) \quad \text {as}\ \lambda \to 1. \end{align}
The condition
$\Phi _1(1)=0$
is the simplest way of imposing the boundary condition (2.12a) on (4.9a). Alternatively, we again eliminate
$X_1/X_0$
and divide the equations to obtain the equivalent condition
\begin{align} \lim _{\lambda \to 1}\frac {P_1}{\Phi _1}=-\frac {1}{mk}-\frac {1}{2m(1+\sigma )}\,, \end{align}
which is more convenient numerically for (4.8a). Either form of boundary condition ensures regularity of the perturbation as
$\xi \to 0$
.
Appendix B. Analysis of the limit
$k\to \infty$
for the stable modes
We wish to solve (4.9) as
$k\to \infty$
for the cases
$n\geqslant 1$
with
$m\gt {3/2}$
and
$n\geqslant 0$
with
$m\lt {3/ 2}$
. From figure 4, we expect that
$\sigma \to 0$
in the first case, and that
$\sigma$
tends to a finite negative limit (distinct from
$-1$
) in the second.
From the form of (4.9a), we expect that
$\Phi _1(\lambda )$
varies rapidly on a scale that is
$O ((k\mathcal {X})^{-1} )$
except perhaps where
$\sigma \approx -\mathcal {N}-1$
. On the other hand, the functions
$\mathcal {X}(\lambda )$
,
$\mathcal {M}(\lambda )$
and
$\mathcal {N}(\lambda )$
vary relatively slowly over the
$O(1)$
interval
$\lambda _*\lt \lambda \lt 1$
. Hence we expect
$\Phi ^{\prime\prime}\gg \Phi^{\prime}\mathcal {X}^{\prime}/\mathcal {X}$
,
$\Phi ^{\prime\prime}\gg \Phi^{\prime}\mathcal {M}^{\prime}/\mathcal {M}$
, and we approximate (4.9) asymptotically by
\begin{align} \Phi^{\prime\prime}_1= k^2 \mathcal {X}^2\left ( 1+{\mathcal {N}\over \sigma +1} \right )\Phi _1 ,\quad \text {where }\mathcal {N}(\lambda )={\mathcal {M}^{\prime}\over 2\mathcal {M}\mathcal {X}} = {\mathcal {M}^{\prime}\mathcal {F}^{\prime}\over \mathcal {M}\mathcal {F}^{\prime\prime}}, \end{align}
\begin{align} \Phi _1=0 \ \text {at}\ \lambda =\lambda _*,\quad \Phi _1\to 0 \ \text {as}\ \lambda \to 1. \\[-6pt] \nonumber \end{align}
For
$m\gt {3/2}$
, the asymptotic approximation of (4.9b) by (B1b) as
$k\to \infty$
follows from
$\mathcal {M}_*\ne 1$
,
$\sigma =O(1)$
and
$\Phi^{\prime}_1/k^2=O(k^{-1})$
; for
$m\lt {3/2}$
, the approximation will be justified in § B.2.
To analyse (B1), we again make the WKB assumption that
$\Phi _1=\exp [kS_0+S_1+O(k^{-1})]$
, and find straightforwardly that
\begin{align} S^{\prime}_0= \pm |\mathcal {X}|\left ( 1+{\mathcal {N}\over \sigma +1} \right )^{1/2} \,. \end{align}
If
$1+\mathcal {N}/(\sigma +1)\gt 0$
everywhere, then
$S^{\prime}_0$
is real, the WKB solutions are exponential in character, and it is impossible to satisfy the boundary conditions (B1b). Hence to satisfy the boundary conditions, there must be a region where
$1+\mathcal {N}/(\sigma +1)\lt 0$
and the WKB solutions are oscillatory. For
$n=O(1)$
, this region must be relatively small since the oscillations are rapid.
The full solutions to (B1) are constructed by matching the oscillatory solution in the region where
$1+\mathcal {N}/(\sigma +1)\lt 0$
through the ‘turning point(s)’ where
$1+\mathcal {N}/(\sigma +1)=0$
to exponentially decaying behaviour in the region(s) where
$1+\mathcal {N}/(\sigma +1)\gt 0$
. The value of
$\sigma$
is determined by the criterion that matching through the turning point results in only the decaying, and not the growing, solution. There are two cases to consider.
B.1. Case
$m\gt {3/2}$
and
$n\geqslant 1$
For
$m\gt {3/2}$
, it can be shown that
$\mathcal {N}(\lambda )$
varies monotonically from
$\mathcal {N}(1)=0$
to
$\mathcal {N}(\lambda _*)=-1$
. Hence if
$\sigma \gt 0$
or
$\sigma \lt -1$
, then there is only exponential behaviour and no solution. We deduce that
$-1\lt \sigma \lt 0$
. Moreover, if the region of oscillatory behaviour is relatively small for
$n=O(1)$
, then it must be near
$\lambda =\lambda _*$
, and we must have
$|\sigma |\ll 1$
.
We expand locally using
$|\sigma |\ll 1$
and
$\mathcal {N}_*=-1$
to obtain
\begin{align} 1+{\mathcal {N}(\lambda )\over 1+\sigma } = 1+\mathcal {N}_*(1-\sigma +\cdots)+(\mathcal {N}-\mathcal {N}_*)(1-\sigma +\cdots)= \sigma +(\lambda -\lambda _*)\mathcal {N}^{\prime}_*+\cdots , \end{align}
where
$\mathcal {N}^{\prime}_*\gt 0$
. At leading order, (B1a) reduces to a shifted form of Airy’s equation:
\begin{align} \Phi^{\prime\prime}_1 = k^2 \mathcal {X}^2\big \{ \sigma +(\lambda -\lambda _*)\mathcal {N}^{\prime}_* \big \}\Phi _1 \,. \end{align}
The condition of matching to a decaying solution as
$\lambda \to 1$
requires
\begin{align} \Phi _1\sim \textrm {Ai}\left ( \big (k^2\mathcal {X}_*^2\mathcal {N}^{\prime}_*\big )^{1/3}\left \{\lambda -\lambda _*+{\sigma \over \mathcal {N}^{\prime}_*}\right \}\right ), \end{align}
where Ai(
$z$
) denotes the Airy function. We note that
$\Phi^{\prime}_1/k^2=O(k^{-4/3})$
, which is consistent with our previous approximation of (4.9b) by (B1b).
The condition
$\Phi _1(\lambda _*)=0$
requires
$(k\,|\mathcal {X}_*|/\mathcal {N}^{\prime}_*)^{2/3}\sigma$
to be one of the roots
$z_n$
of
$\textrm{Ai}(z)=0$
. These roots satisfy
$({2/ 3})(-z_n)^{3/2}\sim (n-({1/4}))\pi$
as
$n\to \infty$
(Abramowitz & Stegun Reference Abramowitz and Stegun1964), and the formula gives even the first root correct to within 1 %. Using this approximation, we obtain the growth rate of the
$n$
th eigenmode as
\begin{align} \sigma \sim -\left (\frac {3}{2}\left \{\left(n-\frac {1}{4}\right)\pi \right \} \frac {\mathcal {N}^{\prime}_*}{|\mathcal {X}_*|}\right )^{2/3}k^{-2/3} \quad \text {for }n\geqslant 1, \end{align}
which, as anticipated, is negative and tends to zero as
$k\to \infty$
. The local expansion can be continued to calculate an
$O(k^{-4/3})$
correction (Dauck Reference Dauck2020). A global WKB approximation can be obtained by integration of the equations for
$S_0$
and
$S_1$
away from the turning point, and the result of this approximation gives excellent agreement with the full numerical results for
$m=5$
and
$n=1,2$
shown in figure 7.
B.2. Case
$m\lt {3/2}$
and
$n\geqslant 1$
For the case
$1\lt m\lt {3/ 2}$
, it can be shown that
$\mathcal {N}(0)=\mathcal {N}(1)=0$
,
$\mathcal {N}(\lambda )\lt 0$
for
$\lambda \in (0,1)$
, and
$\mathcal {N}(\lambda )$
has a unique minimum
$\mathcal {N}_m=\mathcal {N}(\lambda _m)$
, where
$\mathcal {N}_m\in (-{1/4},0)$
and depends on
$m$
. Hence if
$\sigma +1\gt -\mathcal {N}_m$
or
$\sigma +1\lt 0$
, then
$1+\mathcal {N}/(\sigma +1)\gt 0$
everywhere,and there is only exponential behaviour and no solution. We deduce that
$-1\lt \sigma \lt -1-\mathcal {N}_m$
, that the small region of oscillatory behaviour for
$n=O(1)$
must be near
$\lambda _m$
, and that
$\sigma +1$
must be just a little smaller than
$-\mathcal {N}_m$
.
We expand (B1a) near
$\lambda _m$
to obtain
\begin{align} \Phi^{\prime\prime}_1 =k^2 \mathcal {X}_m^2\left ( {\sigma +1+\mathcal {N}_m+{1\over 2}\mathcal {N}^{\prime\prime}_m (\lambda -\lambda _m)^2\over -\mathcal {N}_m} \right )\Phi _1 , \end{align}
where
$\mathcal {N}^{\prime\prime}_m\gt 0$
and we have used
$\sigma +1\sim -\mathcal {N}_m$
in the denominator. Equation (B7) has the same form as the equation of a quantum harmonic oscillator, and the eigenvalues and eigenfunctions are determined similarly by the condition of matching to decaying solutions outside the oscillatory region (i.e. for
$|\lambda -\lambda _m|\gg k^{-1/2}$
). On rescaling the standard results for the harmonic oscillator, we obtain the asymptotic growth rate of the
$n$
th eigenmode as
\begin{align} \sigma \sim -1-\mathcal {N}_m+{(2n+1)\big (-\mathcal {N}^{\prime\prime}_m /2\mathcal {N}_m\big )^{1/2}\mathcal {N}_m\over |\mathcal {X}_m|}\, k^{-1}\quad \text {for }n\geqslant 0, \end{align}
which, as anticipated, is just a little smaller than
$-1-\mathcal {N}_m$
. The value of
$\mathcal {N}_m$
varies monotonically in
$1\lt m\lt {3/2}$
, with
$\mathcal {N}_m\to 0$
as
$m\to 1_+$
, and
$\mathcal {N}_m\to -{1/4}$
as
$m\to ({3/ 2})_-$
. Thus as
$m$
approaches 1, the full range of decay rates
$-1\lt \sigma (k)\lt -1-\mathcal {N}_m$
is confined increasingly closely to the
$t^{-1}$
behaviour obtained at
$m=1$
when the interface is just a passive tracer. As
$m$
approaches
$3/ 2$
, the range of decay rates expands to rates between
$t^{-3/4}$
and
$t^{-1}$
at large and small
$k$
, respectively, but remains well short of instability. (There is a non-uniformity in the double limit
$m\to {3/2}$
as
$k\to \infty$
, which in principle should allow matching to the behaviours found for
$m\gt {3/2}$
in §§ 4.3 and B.1.)
For the case
$m\lt 1$
, we find that
$\mathcal {N}(0)=\mathcal {N}(1)=0$
again, but
$\mathcal {N}(\lambda )\gt 0$
for
$\lambda \in (0,1)$
, and
$\mathcal {N}(\lambda )$
has a unique maximum
$\mathcal {N}_m$
. Hence if
$\sigma +1\lt -\mathcal {N}_m$
or
$\sigma +1\gt 0$
, then
$1+\mathcal {N}/(\sigma +1)\gt 0$
everywhere, and we deduce that
$-1-\mathcal {N}_m\lt \sigma \lt -1$
. The signs of
$\mathcal {N}_m$
,
$\mathcal {N}^{\prime\prime}_m $
and
$\sigma +1$
are the opposite of those for
$1\lt m\lt {3/2}$
, but the same argument can be followed through and leads to the same expression (B8) for the asymptotic behaviour of the growth rates.
The asymptotic analysis again agrees well with the numerical calculations. For example, as
$k\to \infty$
,
$\sigma \to -1-\mathcal {N}_m$
and the eigenfunctions concentrate around
$\xi _m=X_0(\lambda _m)$
. For
$m=0.15$
this gives
$\sigma \to -2.93$
and
$\xi _m=1.10$
, while for
$m=1.25$
it gives
$\sigma \to -0.87$
and
$\xi _m=1.53$
; cf. figures 4 and 5.
Appendix C. Connection with Sturm–Liouville theory
To establish a connection between the linear stability problem of § 4.1 and Sturm–Liouville theory, it is convenient to change variables using
$\xi =\big \{\mathcal {F}(\lambda )\big \}^{1/2}$
and
$\mathcal {X}^{-1}\,\textrm{d}/\textrm{d}\lambda = \xi \, \textrm {d}/\textrm {d}\xi$
to rewrite (4.9) as
\begin{align} \frac {\textrm {d}}{\textrm {d}\xi }\left (\frac {\xi }{\mathcal {M}}\,\frac {\textrm {d}\Phi _1}{\textrm {d}\xi }\right )= \frac {k^2}{\xi \mathcal {M}} \left ( 1+\frac {\mathcal {N}}{\sigma +1} \right )\Phi _1 ,\quad \text{where }\mathcal {N}=\frac {\xi }{2\mathcal {M}}\,\frac {\textrm {d}\mathcal {M}}{\textrm {d}\xi }, \end{align}
with
$\mathcal {M}$
and
$\mathcal {N}$
now implicit functions of
$\xi$
. We compare (C1) with the standard form
$(py^{\prime})^{\prime}-qy+\lambda w y=0$
, where
$\lambda$
is the eigenvalue,
$y$
is the eigenfunction, and
$p$
,
$q$
and
$w$
are coefficient functions. We set
\begin{align} p(\xi )={\xi \over \mathcal {M}},\quad q(\xi )=\frac {k^2}{\xi \mathcal {M}},\quad w(\xi )=\frac {k^2}{ \mathcal {M}^2}\left |\frac {\textrm {d}\mathcal {M}}{\textrm { d}\xi }\right |,\quad \lambda =\frac {\textrm {sgn}(m-1)}{ 2(\sigma +1)}, \end{align}
so
$p$
and
$w$
are, as usual, positive on
$(0,\xi _*)$
; the factor
$\textrm {sgn}(m-1)$
in
$\lambda$
then accounts for
$\textrm {d}\mathcal {M}/ \textrm {d}\xi \lt 0$
if
$m\gt 1$
and
$\textrm {d}\mathcal {M}/ \textrm {d}\xi\gt 0$
if
$m\lt 1$
. Also
$p$
,
$p^{\prime}$
,
$q$
and
$w$
are continuous on
$(0,\xi _*)$
.
Equation (C1) has a singular point at
$\xi =0$
, where
$p=0$
, but it is regular and non-oscillatory with local solutions
$\Phi _1=\xi ^{\pm k}$
. Hence boundedness of
$\Phi _1$
at
$\xi =0$
would be sufficient as a boundary condition to eliminate the singular solutions. (The situation at this point is essentially the same as for Bessel’s equation at the origin.)
For
$m\lt {3/ 2}$
, we have
$\mathcal {M}_*=1$
and the boundary conditions (4.9b) reduce to
\begin{align} \Phi _1\to 0 \ \text {as}\ \xi \to 0,\quad \xi\, \frac {\textrm {d}\Phi _1}{\textrm {d} \xi }+k\Phi _1=0 \ \text {at}\ \xi =\xi _* , \end{align}
which are standard separated boundary conditions. Thus (C1) and (C3) have the form of a Sturm–Liouville problem. From the general theory of such problems (see e.g. Morse & Feshbach Reference Morse and Feshbach1953, pp. 719ff), we deduce that (C1) and (C3) have a discrete spectrum of distinct eigenvalues
$\lambda _n$
with
$\lambda _n\to \infty$
as
$n\to \infty$
. It follows from (C2) that
$\sigma _n\to -1$
from above (below) if
$m\gt 1$
(
$m\lt 1$
) as
$n\to \infty$
for fixed
$k$
, which confirms and extends the trends seen for small
$n$
in figure 4.
For
$m\gt {3/ 2}$
, the argument is complicated slightly by the appearance of the eigenvalue
$\lambda$
in the boundary condition:
\begin{align} \xi\, \frac {\textrm {d}\Phi _1}{\textrm {d} \xi }+\mathcal {M}_*k\Phi _1-\lambda (\mathcal {M}_*-1)k^2\Phi _1=0 \quad \text {at}\ \xi =\xi _* . \end{align}
However, if (C4) is replaced by either of the boundary conditions (i)
$\xi \Phi _{1\xi }+\mathcal {M}_*k\Phi _1=0$
or (ii)
$\Phi _1(\xi _*)=0$
, then in each case, we obtain a standard Sturm–Liouville problem with a discrete unbounded spectrum. It can be argued that the sequence of eigenvalues with condition (C4) is bracketed by the sequences with conditions (i) and (ii), hence that the eigenvalues with (C4) also form a discrete unbounded spectrum; once again, for fixed
$k$
,
$\sigma _n\to -1$
from above as
$n\to \infty$
.
Open access
