# Buoyancy segregation suppresses viscous fingering in horizontal
displacements in a porous layer

Edward M. Hinton1 and Apoorv Jyoti2<EMAIL_ADDRESS>1School of
Mathematics and Statistics, The University of Melbourne, Victoria 3010,
Australia 2Peter Cook Centre for CCS Research, School of Geography, Earth and
Atmospheric Sciences, The University of Melbourne, Victoria 3010, Australia

###### Abstract

We consider the axisymmetric displacement of an ambient fluid by a second
input fluid of lower density and lower viscosity in a horizontal porous layer.
If the two fluids have been vertically segregated by buoyancy, the flow
becomes self-similar with the input fluid preferentially flowing near the
upper boundary. We show that this axisymmetric self-similar flow is stable to
angular-dependent perturbations for any viscosity ratio. The Saffman-Taylor
instability is suppressed due to the buoyancy segregation of the fluids. The
radial extent of the segregated current is inversely proportional to the
viscosity ratio. This horizontal extension of the intrusion eliminates the
discontinuity in the pressure gradient between the fluids associated with the
viscosity contrast. Hence, at late times viscous fingering is shut down even
for arbitrarily small density differences. The stability is confirmed through
numerical integration of a coupled problem for the interface shape and the
pressure gradient, and through complementary asymptotic analysis, which
predicts the decay rate for each mode. The results are extended to anisotropic
and vertically heterogeneous layers. The interface may have steep shock-like
regions but the flow is always stable when the fluids have been segregated by
buoyancy, as in a uniform layer.

##  1 Introduction

Displacement flows in porous media occur in many important environmental
settings including geological CO2 sequestration, underground hydrogen storage,
geothermal energy generation and the infiltration of sea water into fresh
water aquifers. If the displacing fluid is of lesser viscosity than the
ambient or displaced fluid, then the fluid-fluid interface can become
unstable. This phenomenon is known as the ‘Saffman-Taylor instability’ or as
‘viscous fingering’ owing to the fingers of low viscosity input fluid that
penetrate the ambient fluid (Saffman & Taylor, 1958; Paterson, 1981; Homsy,
1987). In many applications, understanding this instability is a key concern.
For example, sequestered carbon dioxide is much less viscous than the host
brine in the aquifer and viscous fingering may significantly reduce the
fraction of the aquifer that can be accessed by the CO2, which in turn reduces
the storage efficiency (Bachu, 2015). Multiple strategies have been developed
to suppress viscous fingering including controlling the flow geometry and
varying the volume flux of the input fluid (Al-Housseiny et al., 2012; Zheng
et al., 2015a). Significant miscibility between the fluids has also been shown
to shut down fingering at later times (Nijjer et al., 2018; Sharma et al.,
2020).

In the case that the two fluids have different densities, as is relevant to
CO2 sequestration, variations in the hydrostatic pressure provide an extra
horizontal force that drives (or opposes) the fluid displacement; such flows
are known as ‘gravity currents’ (Huppert & Woods, 1995). There has been much
research into gravity currents in porous media, particularly in confined
layers in which the ambient fluid must be displaced (Hesse et al., 2007;
Juanes et al., 2010; Zheng & Stone, 2022). Viscous fingering can occur at
early times if the input fluid is of lower viscosity. However, it is generally
assumed that at later times a stable interface between the input and ambient
fluids develops, which seems to be corroborated by laboratory experiments
(Nordbotten & Celia, 2006; Pegler et al., 2014; Guo et al., 2016). It has also
been shown that the combination of buoyancy and mixing can suppress the
Saffman-Taylor instability at long times (Tchelepi, 1994; Riaz & Meiburg,
2003; Nijjer et al., 2022). However, for sharp-interface gravity currents in
confined porous media, stability of the gravity current solution has not yet
been confirmed and the flow physics that suppresses viscous fingering is not
fully understood.

In the present paper, we consider injection of a relatively buoyant fluid into
an aquifer confined above and below by two impermeable layers. The injection
source is on the upper boundary and the aquifer is initially filled with a
second fluid of greater viscosity (see figure 1). We show that if the fluids
have segregated owing to buoyancy, the gravity current solution is stable to
both axisymmetric and angular-dependent perturbations for any values of the
viscosity ratio and the input flux relative to the buoyancy velocity.

The classical Saffman-Taylor instability is associated with the discontinuity
in the pressure gradient across the fluid-fluid interface (Paterson, 1981).
For example, in the case of flow in a porous medium with no vertical
dimension, the base radial displacement has a circular interface. The total
radial flux in each fluid in a circular cross-section is given by

$q_{r}=\frac{-2\pi rk}{\mu}\frac{\partial p}{\partial r},$ (1)

where $k$ is the permeability, $\mu$ is the viscosity, $r$ is the radial
coordinate and $p$ is the pressure. This flux is continuous across the
interface so the pressure gradient either side of the interface is
proportional to the viscosity of the fluid there. If the pressure gradient is
larger in the ambient fluid then perturbations to the interface into the
ambient fluid experience a larger pressure gradient and tend to grow. Hence
viscous fingers can arise if the ambient fluid is of greater viscosity.

In a three-dimensional porous layer with fluids of different densities that
have been vertically segregated by buoyancy, the input fluid will form an
intrusion with large radial extent at later times. The radial flux at the
interface no longer needs to be continuous because the input intrusion
advances into the ambient fluid. Indeed, we will show that the vertical
structure of the flow ensures that there is no destabilising pressure gradient
across the interface. This stabilisation occurs for any density difference
between the two fluids that eventually leads to buoyancy segregation.

The combination of buoyancy segregation and an intrusion of large radial
extent is a fundamentally different stability mechanism to the related problem
of two-layer viscous gravity currents where two viscous fluids co-flow driven
by gravity alone. In this case, buoyancy segregation does not guarantee
stability. Instead, stability is controlled by competing contributions to the
pressure gradient from buoyancy and the viscosity contrast (Kowal & Worster,
2019). The Saffman-Taylor instability is suppressed provided that the density
difference is sufficiently large. A similar phenomena occurs in unconfined
flows with a longitudinal viscosity contrast (Kowal, 2021).

The stability of single-layer unconfined gravity currents (in which the
aquifer is bounded only on one side by a single impermeable layer), has been
investigated extensively (Pattle, 1959; Grundy & McLaughlin, 1982; Newman,
1984; Bernoff & Witelski, 2002; Chertock, 2002; Mathunjwa & Hogg, 2006). In
this case, the motion of the ambient fluid is unimportant and the gravity
current is not resisted by its displacement; the solutions are always stable.
The stability of confined gravity currents is different because the
displacement of the more viscous ambient fluid influences the physics. For
these confined flows, stability has been established for perturbations that do
not depend on the transverse or azimuthal coordinate (Nordbotten & Celia,
2006). The stability analysis of both unconfined gravity currents and
axisymmetric variations to confined gravity currents is significantly
simplified because the driving pressure gradient in the input fluid can be
eliminated from the problem. The present study of three-dimensional
perturbations to confined currents requires the solution to a coupled system
for both the pressure gradient and the interface shape.

The paper is structured as follows. The model for the buoyancy-segregated
displacement flow is formulated in §2. The coupled system that governs the
background pressure and the interface shape is integrated numerically in §3
with the results demonstrating that the axisymmetric solutions are stable to
angular-dependent perturbations for all values of the viscosity ratio and the
density difference. In §4, the physical mechanism that suppresses the Saffman-
Taylor instability is analysed. A linear stability analysis is carried out in
§5 for the case of buoyancy-segregated fluids with an arbitrarily small
density difference. In §6, we show that the stabilising effect of buoyancy
segregation generalises to layers with vertically varying permeability for
which the interface may contain steep shock-like regions. Concluding remarks
are made in §7.

Figure 1: (a) Schematic of the setup. (b) Cross-section of the setup.

##  2 Model

The setup is shown in figure 1. A similar geometry was studied by Nordbotten &
Celia (2006) and Guo et al. (2016). Fluid of density $\rho$ and viscosity
$\mu_{i}$ is injected into a porous medium of thickness $H_{0}$ and porosity
$\phi$. The medium is initially filled with a second fluid of greater density
$\rho+\Delta\rho$ and different viscosity $\mu_{a}$. The top and bottom
boundaries of the layer are impermeable. The input fluid is injected with
volume flux $Q$ from a point source at the upper boundary, which we take to be
the origin of our coordinate system. The $Z$ coordinate is measured downwards
from the upper boundary at which $Z=0$. We use cylindrical polar coordinates
with $R$ denoting the distance from the $Z$ axis and $\theta$ denoting the
angular coordinate. Time is denoted by $T$. We neglect any intermingling of
the fluids and assume that there is a sharp interface at $Z=H(R,\theta,T)$
between the fluids. We also assume that the two fluids are segregated by
buoyancy so that $H$ is a single-valued function.

When the input fluid has spread far from the source, the flow is predominantly
in the radial direction. We may then apply the shallow flow approximation and
the pressure is approximately hydrostatic,

$P=\bar{P}+\rho gZ\qquad\text{for}\quad 0<Z<H,$ (2) $P=\bar{P}-\Delta\rho
gH+(\rho+\Delta\rho)gZ\qquad\text{for}\quad H<Z<H_{0},$ (3)

in the input and ambient fluids, respectively, where $\bar{P}(R,\theta,T)$ is
the pressure at the upper boundary, $Z=0$, which is to be determined. The
Darcy velocities are given by

$\mathbf{U}_{i}=-\frac{k}{\mu_{i}}\nabla\bar{P},$ (4)
$\mathbf{U}_{a}=-\frac{k}{\mu_{a}}\Big{(}\nabla\bar{P}-\Delta\rho g\nabla
H\Big{)},$ (5)

in the input and ambient fluids, respectively. The governing equations of the
flow are

$\phi\frac{\partial H}{\partial T}+\nabla\cdot\big{(}H\mathbf{U}_{i}\big{)}=0$
(6)

and

$\nabla\cdot\big{[}H\mathbf{U}_{i}+(H_{0}-H)\mathbf{U}_{a}\big{]}=0,$ (7)

which represent mass conservation integrated over the thickness of the input
fluid, and mass conservation integrated over the thickness of the layer,
respectively.

At $T=0$, the porous layer is entirely filled with the ambient fluid and
constant-flux injection of the input fluid begins. Global mass conservation of
the input fluid takes the form

$\int_{0}^{2\pi}\int_{0}^{R_{f}(\theta,T)}\phi
RH\,\mathrm{d}R\,\mathrm{d}\theta=QT,$ (8)

where $R_{f}(\theta,T)$ denotes the frontal contact line where $H=0$ (see
figure 1b).

Eventually, the fluid-fluid interface will touch the lower boundary
(Nordbotten & Celia, 2006). It is well-known that the predicted interface
shape of an unconfined axisymmetric gravity current has an unphysical
singularity at the origin, whilst in confined layers, the flow will always
encompass the entire layer near the source, regardless of the thickness of the
layer (Lyle et al., 2005; Guo et al., 2016). Recently, Benham et al. (2022)
carried out a detailed analysis of the flow near the source, incorporating the
pressure-driven vertical flow there. Their results corroborate the conclusion
that a confined layer becomes filled with input fluid at sufficiently late
times (see also Huppert & Pegler, 2022). After the fluid-fluid interface has
touched the lower boundary, there is a second contact line at
$R=R_{b}(\theta,T)$ where $H=H_{0}$. The flow then has three regions: there is
the ‘fully-flooded’ zone in which the input fluid fills the layer ($H=H_{0}$;
$0\leq R<R_{b}(\theta,t)$); next is the partially filled zone in which the
fluid-fluid interface lies between the top and bottom boundaries ($0<H<H_{0}$;
$R_{b}(\theta,t)<R<R_{f}(\theta,t)$); finally, beyond this there is an
uninvaded zone in which the ambient fluid fills the layer ($H=0$;
$R>R_{f}(\theta,t)$); see figure 1b.

When the input fluid has filled the layer near the source, the boundary
conditions are as follows. At the source,

$H\to H_{0},\qquad\frac{\partial\bar{P}}{\partial R}\to\frac{-\mu_{i}Q}{2\pi
kH_{0}R}\qquad\text{as}\quad R\to 0,$ (9)

where the latter arises from volume conservation. In the far field, the
boundary condition is

$H\to 0,\qquad\frac{\partial\bar{P}}{\partial R}\to\frac{-\mu_{a}Q}{2\pi
kH_{0}R}\qquad\text{as}\quad R\to\infty.$ (10)

The dependent variables, $\bar{P}(R,\theta,T)$ and $H(R,\theta,T)$, are
defined on $\theta\in[0,2\pi)$ with periodic boundary conditions:
$\bar{P}(R,0,T)=\bar{P}(R,2\pi,T)$ and $H(R,0,T)=H(R,2\pi,T)$. We also require
that $0\leq H\leq H_{0}$. Equations (6) and (7) combined with global mass
conservation (8) and the boundary conditions (9), (10) define a complete
coupled system for determining the two dependent variables
$\bar{P}(R,\theta,T)$ and $H(R,\theta,T)$.

Although the system derived above is sufficient to solve the problem
numerically (see §3), the stability investigation requires detailed treatment
of the behaviour in the vicinity of the two contact lines where the interface
touches the layer boundaries. Across these lines, the gradient of the
interface and the gradient of the pressure at the top of the layer can be
discontinuous. Indeed, it is well-known that the stability of the fluid-fluid
interface for equally dense fluids is controlled by the jump in the pressure
gradient at the interface (see §1). Therefore, we derive boundary conditions
for the interface shape and the pressure gradient across each contact line.

At the trailing contact line, $R=R_{b}(\theta,t)$, continuity of the
interface, continuity of the pressure and continuity of the flux of the input
fluid take the following forms

$H(R=R_{b}^{+})=H(R=R_{b}^{-})=H_{0},\qquad\big{[}\bar{P}\big{]}^{+}_{-}=0,\qquad\Bigg{[}\frac{\partial\bar{P}}{\partial
R}-\frac{\partial
R_{b}}{\partial\theta}\frac{1}{R}\frac{\partial\bar{P}}{\partial\theta}\Bigg{]}^{+}_{-}=0,$
(11)

respectively, where $[\cdot]^{+}_{-}$ denotes the discontinuity in the
quantity in square brackets at the contact line $R=R_{b}(\theta,t)$. In
addition, we have a kinematic boundary condition; at the trailing contact
line, the radial velocity of the ambient fluid is equal to the radial velocity
of the interface, which takes the following form

$-\frac{k}{\mu_{a}}\Bigg{(}\frac{\partial\bar{P}}{\partial R}-\Delta\rho
g\frac{\partial H}{\partial
R}\Bigg{)}\Bigg{\rvert}_{R=R_{b}^{+}}=\frac{\partial R_{b}}{\partial
T}-\frac{k}{\mu_{a}}\frac{\partial
R_{b}}{\partial\theta}\Bigg{(}\frac{1}{R}\frac{\partial\bar{P}}{\partial\theta}-\frac{\Delta\rho
g}{R}\frac{\partial H}{\partial\theta}\Bigg{)}\Bigg{\rvert}_{R=R_{b}^{+}}.$
(12)

At the leading contact line, $R=R_{f}(\theta,t)$, the corresponding boundary
conditions are

$H(R=R_{f}^{+})=H(R=R_{f}^{-})=H_{0},\qquad\big{[}\bar{P}\big{]}^{+}_{-}=0,$
(13) $\Bigg{[}\frac{\partial(\bar{P}-\Delta\rho gH)}{\partial
R}-\frac{\partial
R_{f}}{\partial\theta}\frac{1}{R}\frac{\partial(\bar{P}-\Delta\rho
gH)}{\partial\theta}\Bigg{]}^{+}_{-}=0,$ (14)
$-\frac{k}{\mu_{i}}\frac{\partial\bar{P}}{\partial
R}\Bigg{\rvert}_{R=R_{f}^{-}}=\frac{\partial R_{f}}{\partial
T}-\frac{k}{\mu_{i}R}\frac{\partial
R_{f}}{\partial\theta}\frac{\partial\bar{P}}{\partial\theta}\Bigg{\rvert}_{R=R_{f}^{-}},$
(15)

where $[\cdot]^{+}_{-}$ denotes the discontinuity in the quantity in square
brackets at the contact line $R=R_{f}(\theta,t)$.

###  2.1 Transformed coordinate system

Mass conservation of the input fluid suggests that the radial extent of the
currents grows with $R^{2}\sim QT/(\phi H_{0})$. This motivates introducing
the following transformed dimensionless coordinates

$\eta=\Bigg{(}\frac{2\pi\phi
H_{0}}{Q}\Bigg{)}^{1/2}\frac{R}{T^{1/2}},\qquad\tau=\log\big{(}T/T_{\mathrm{ref}}\big{)},$
(16)

where $T=T_{\mathrm{ref}}>0$ is a reference time, which corresponds to
$\tau=0$. We also scale the dependent variables to obtain the following
dimensionless quantities

$\mathcal{H}(\eta,\theta,\tau)=\frac{H(R,\theta,T)}{H_{0}},\qquad\mathcal{P}(\eta,\theta,\tau)=\frac{2\pi
kH_{0}\bar{P}(R,\theta,T)}{\mu_{i}Q}.$ (17)

The two mass conservation equations (6, 7) become

$\frac{\partial\mathcal{H}}{\partial\tau}-\frac{\eta}{2}\frac{\partial\mathcal{H}}{\partial\eta}-\tilde{\nabla}\cdot\big{(}\mathcal{H}\tilde{\nabla}\mathcal{P}\big{)}=0,$
(18)
$\tilde{\nabla}\cdot\Big{[}\big{(}M+(1-M)\mathcal{H}\big{)}\tilde{\nabla}\mathcal{P}-GM(1-\mathcal{H})\tilde{\nabla}\mathcal{H}\Big{]}=0,$
(19)

where

$\tilde{\nabla}=\Bigg{(}\frac{\partial}{\partial(\eta\cos\theta)},\,\frac{\partial}{\partial(\eta\sin\theta)}\Bigg{)}$
(20)

is the gradient operator with respect to the transformed coordinates and we
have introduced the following two dimensionless groups

$M=\frac{\mu_{i}}{\mu_{a}},\qquad G=\frac{2\pi k\Delta\rho
gH_{0}^{2}}{\mu_{i}Q},$ (21)

which are the viscosity ratio and gravity number, respectively. The latter
represents the magnitude of pressure gradients associated with buoyancy
relative to pressure gradients associated with injection of fluid. In the case
of equally dense fluids, the gravity number is $G=0$. The boundary conditions
at the source and in the far field are

$\mathcal{H}\to
1,\qquad\frac{\partial\mathcal{P}}{\partial\eta}\to\frac{-1}{\eta}\qquad\text{as}\quad\eta\to
0,$ (22)

and

$\mathcal{H}\to
0,\qquad\frac{\partial\mathcal{P}}{\partial\eta}\to\frac{-M^{-1}}{\eta}\qquad\text{as}\quad\eta\to\infty.$
(23)

The trailing contact line in the transformed coordinates is denoted by
$\eta=\eta_{b}(\theta,\tau)$ and the boundary conditions there (11, 12) are
recast as

$\mathcal{H}(\eta=\eta_{b}^{+})=\mathcal{H}(\eta=\eta_{b}^{-})=1,\qquad\big{[}\mathcal{P}\big{]}^{+}_{-}=0,\qquad\Bigg{[}\frac{\partial\mathcal{P}}{\partial\eta}-\frac{\partial\eta_{b}}{\partial\theta}\frac{1}{\eta}\frac{\partial\mathcal{P}}{\partial\theta}\Bigg{]}^{+}_{-}=0,$
(24)
$-M\Bigg{(}\frac{\partial\mathcal{P}}{\partial\eta}-G\frac{\partial\mathcal{H}}{\partial\eta}\Bigg{)}\Bigg{\rvert}_{\eta=\eta_{b}^{+}}=\frac{\partial\eta_{b}}{\partial\tau}+\frac{1}{2}\eta_{b}-M\frac{\partial\eta_{b}}{\partial\theta}\Bigg{(}\frac{1}{\eta}\frac{\partial\mathcal{P}}{\partial\theta}-\frac{G}{\eta}\frac{\partial\mathcal{H}}{\partial\theta}\Bigg{)}\Bigg{\rvert}_{\eta=\eta_{b}^{+}},$
(25)

whilst at the leading contact line, $\eta=\eta_{f}(\theta,\tau)$,

$\mathcal{H}(\eta=\eta_{f}^{+})=\mathcal{H}(\eta=\eta_{f}^{-})=0,\qquad\big{[}\mathcal{P}\big{]}^{+}_{-}=0,$
(26)
$\Bigg{[}\frac{\partial(\mathcal{P}-G\mathcal{H})}{\partial\eta}-\frac{\partial\eta_{b}}{\partial\theta}\frac{1}{\eta}\frac{\partial(\mathcal{P}-G\mathcal{H})}{\partial\theta}\Bigg{]}^{+}_{-}=0,$
(27)
$-\frac{\partial\mathcal{P}}{\partial\eta}\Bigg{\rvert}_{\eta=\eta_{f}^{-}}=\frac{\partial\eta_{f}}{\partial\tau}+\frac{1}{2}\eta_{f}-\frac{1}{\eta}\frac{\partial\eta_{b}}{\partial\theta}\frac{\partial\mathcal{P}}{\partial\theta}\Bigg{\rvert}_{\eta=\eta_{f}^{-}}.$
(28)

Mass conservation (8) takes the form

$\frac{1}{2\pi}\int_{0}^{2\pi}\int_{0}^{\eta_{f}}\eta\mathcal{H}_{0}\,\mathrm{d}\eta\,\mathrm{d}\theta=1.$
(29)

###  2.2 Self-similar axisymmetric solution

The flow has a self-similar axisymmetric solution with (Nordbotten & Celia,
2006; Guo et al., 2016)

$\mathcal{H}(\eta,\theta,\tau)=\mathcal{H}_{0}(\eta),\qquad\mathcal{P}(\eta,\theta,\tau)=\mathcal{P}_{0}(\eta),$
(30)

and

$\eta_{b}(\theta,\tau)=\eta_{b_{0}},\qquad\eta_{f}(\theta,\tau)=\eta_{f_{0}},$
(31)

so that the contact lines are circles, fixed in similarity space. We are
interested in analysing the stability of and convergence to these solutions in
§§3, 4 and 5. In the present section, we derive a single ordinary differential
equation which governs the axisymmetric self-similar interface shape and we
recall an important analytical solution to this equation from Nordbotten &
Celia (2006) and Guo et al. (2016). For axisymmetric flow, the pressure at the
top boundary, $\mathcal{P}_{0}(\eta)$, can be eliminated from the problem as
follows. We integrate equation (19) with respect to $\eta$ and apply the
source flux boundary condition (22) to obtain

$\frac{\mathrm{d}\mathcal{P}_{0}}{\mathrm{d}\eta}=\frac{-1}{\eta\big{(}M+(1-M)\mathcal{H}_{0}\big{)}}+\frac{GM(1-\mathcal{H}_{0})}{\big{(}M+(1-M)\mathcal{H}_{0}\big{)}}\frac{\mathrm{d}\mathcal{H}_{0}}{\mathrm{d}\eta}.$
(32)

Upon substituting (32) into (18), we obtain the following ordinary
differential equation for $\mathcal{H}_{0}(\eta)$,

$-\frac{\eta}{2}\frac{\mathrm{d}\mathcal{H}_{0}}{\mathrm{d}\eta}+\frac{1}{\eta}\frac{\mathrm{d}}{\mathrm{d}\eta}\Bigg{(}\frac{\mathcal{H}_{0}}{M+(1-M)\mathcal{H}_{0}}\Bigg{)}=\frac{1}{\eta}\frac{\mathrm{d}}{\mathrm{d}\eta}\Bigg{(}\eta\frac{GM\mathcal{H}_{0}(1-\mathcal{H}_{0})}{M+(1-M)\mathcal{H}_{0}}\frac{\mathrm{d}\mathcal{H}_{0}}{\mathrm{d}\eta}\Bigg{)}.$
(33)

Figure 2: Self-similar axisymmetric flow. (a) Interface shape for a viscosity
ratio of $M=0.2$ and $G=0.1$ (blue line) and the small-$G$ approximate
solution (35, black dashed line). (b) Small-$G$ approximate interface shapes
(35) for four values of $M$.

This axisymmetric system can be solved numerically for $\mathcal{H}_{0}(\eta)$
as described in Guo et al. (2016) for a wide range of $M$ and $G$. For $G\gg
1$ (buoyancy-dominated flow), the input fluid forms a thin current near the
top of the layer and the displacement of the ambient fluid is unimportant. The
flow is effectively ‘unconfined’ and single phase. Hence, we do not expect
Saffman-Taylor instabilities in this regime. Instead, we focus on the case of
relatively larger input flux ($G\leq\mathcal{O}(1)$) and a less viscous input
fluid ($M<1$) for which stability is not well-understood. For $G\ll 1$ and
$M<1$, (33) has an approximate analytic solution. We neglect the right-hand
side of (33). Although this removes the highest derivative of
$\mathcal{H}_{0}$, the boundary conditions (28) and (25) become independent of
$\mathcal{H}_{0}$ for $G\ll 1$ and simply give the locations of the contact
points so the problem remains complete. Integrating the left-hand side of (33)
furnishes the following interface shape (Guo et al., 2016)

$\displaystyle\mathcal{H}_{0}$ $\displaystyle=1,$ $\displaystyle 0<$
$\displaystyle\eta<\eta_{b_{0}},$ (34) $\displaystyle\mathcal{H}_{0}$
$\displaystyle=\frac{(2M)^{1/2}}{(1-M)\eta}-\frac{M}{1-M},$
$\displaystyle\eta_{b_{0}}<$ $\displaystyle\eta<\eta_{f_{0}},$ (35)
$\displaystyle\mathcal{H}_{0}$ $\displaystyle=0,$ $\displaystyle\eta_{f_{0}}<$
$\displaystyle\eta,$ (36)

where the contact lines are given by

$\eta_{b_{0}}=\sqrt{2M}\qquad and\qquad\eta_{f_{0}}=\sqrt{2/M}.$ (37)

Note that these provide inner bounds on the locations of the contact lines for
general $G>0$ provided that $M<1$:

$\eta_{b_{0}}\leq\sqrt{2M},\qquad\eta_{f_{0}}\geq\sqrt{2/M},$ (38)

with equality as $G\to 0^{+}$. These inequalities arise because buoyancy acts
to extend the interface (its effect is proportional to $G$).

The corresponding pressure at the top of the layer takes the form

$\displaystyle\mathcal{P}_{0}$
$\displaystyle=-1-\log\Bigg{(}\frac{\eta}{(2M)^{1/2}}\Bigg{)}$ $\displaystyle
0<$ $\displaystyle\eta<\eta_{b_{0}},$ (39) $\displaystyle\mathcal{P}_{0}$
$\displaystyle=-\frac{\eta}{(2M)^{1/2}},$ $\displaystyle\eta_{b_{0}}<$
$\displaystyle\eta<\eta_{f_{0}},$ (40) $\displaystyle\mathcal{P}_{0}$
$\displaystyle=-M^{-1}-M^{-1}\log\Bigg{(}\frac{\eta}{(2/M)^{1/2}}\Bigg{)}$
$\displaystyle\eta_{f_{0}}<$ $\displaystyle\eta,$ (41)

where $\mathcal{P}_{0}(\eta)$ is defined up to addition of an arbitrary
constant.

The numerical solution to (33) for $\mathcal{H}_{0}$ is plotted as a blue line
in figure 2a for $G=0.1$ and $M=0.2$. The small-$G$ analytical solution (35)
is shown as a black dashed line and there is good agreement with the numerical
result. Figure 2b shows the analytical solutions for various values of $M<1$.
When the input fluid is of relatively lower viscosity (smaller $M$), it forms
a finger that intrudes further into the ambient fluid at the top of the layer.

The axisymmetric analytic solution (35) is relevant to the regime in which
pressure gradients associated with buoyancy are weak compared to those
associated with injection ($G\ll 1$). It is worth considering the behaviour in
the limit $G\to 0$, corresponding to equally dense fluids and no buoyancy
force. As $G\to 0$, the axisymmetric analytic solution (35) becomes a more
accurate approximation of the numerical solution, with the input fluid
preferentially flowing near the upper boundary (figure 2). However, $G=0$
corresponds to no buoyancy segregation and so the input fluid should have no
preference for the top of the layer. Indeed, in the case of equally dense
fluids or an arbitrarily thin layer ($G=0$), and $M<1$, one would observe
classical Saffman-Taylor fingering with angular dependence and no vertical
preference (Paterson, 1981). The derivation of the axisymmetric analytical
solution (35) implicitly assumes that the input and ambient fluid have been
segregated by buoyancy, which requires $G>0$ and takes a time proportional to
$G^{-1}$. Hence, the case of $G=0$ has qualitatively different late-time
behaviour to the case of small but non-zero $G$ and the limit $G\to 0$ is
singular. In this paper, we analyse the stability when $G$ is small but
positive; we exclude the case $G=0$ for which buoyancy has no effect.

For $G$ small but non-zero, the interfaces in figure 2b occur provided that
(i) buoyancy segregation has occurred and (ii) the radial extent of the flow
is much greater than the layer thickness so that the shallow approximation
applies. For $M<1$, in dimensional terms, the former requires

$T\gg\frac{\mu_{a}H_{0}}{\Delta\rho gk},$ (42)

whilst the latter requires

$T\gg\frac{2\pi\phi H_{0}^{3}}{Q}.$ (43)

Prior to buoyancy segregation, Saffman-Taylor fingering can occur. The present
paper is concerned with the post-segregation stability of the axisymmetric
solutions to (33) to angular-dependent perturbations with $\Delta\rho>0$. For
$G\ll 1$, the first condition (42) implies the second (43). Notice that for
small density differences, $\Delta\rho$, the time for buoyancy segregation is
very large. It is also important to note that in an anistropic porous medium,
the permeability $k$ that appears in segregation timescale (42) will depend
only on the vertical permeability, whereas the gravity number (and hence the
self-similar solutions) depends only on horizontal permeability (cf. Benham et
al., 2022). The results we derive in this paper concerning flow stability
extend to porous layers that have different vertical and horizontal
permeabilities (as is common in many aquifers); see also §6. The case of more
complicated anisotropy such as different horizontal permeabilities in
different directions is beyond the scope of the present work.

For simplicity, we have modelled the case where the source of input fluid is
located at a point on the upper boundary. However, the analysis and results of
this paper apply for any location of the injection source (for example, at the
lower boundary, or over a vertical line within the layer). This is because the
input fluid will always eventually fill the thickness of the layer near the
source. Subsequently, the flow becomes predominantly radial. There will, of
course, be different early-time transient behaviour for different source
locations but once (42) and (43) apply, the exact source location becomes
unimportant.

##  3 Numerical results

Figure 3: Suppression of mode six fingers for $G=0.1$ and $M=0.2$. (a) Initial
condition at the reference time ($\tau=0$) for the flow thickness,
$\mathcal{H}$, in transformed coordinates. (b,c) The flow thickness at
$\tau=1.25$ and $\tau=4$. (d) The evolution of the contact line, where
$\mathcal{H}=0$. Figure 4: Suppression of mode three fingers for $G=0.5$ and
$M=0.5$. (a) Initial condition ($\tau=0$) for the interface shape. (b) Initial
condition ($\tau=0$) for the pressue gradient (log of its magnitude is shown).
(c,d) Corresponding plots at $\tau=1$. Radial cross-sections of these results
are shown in figure 5. Figure 5: Suppression of mode three fingers for $G=0.5$
and $M=0.5$ (same parameters and initial condition as figure 4). Radial cross-
sections at $\theta=0$ (black lines) and $\theta=\pi/3$ (blue lines) of the
flow thickness, $\mathcal{H}$, and radial pressure gradient at the upper
boundary $\partial\mathcal{P}/\partial\eta$, at three times: (a,d)
$\tau=0.25$, (b,e) $\tau=1$, (c,f) $\tau=4$. The red dotted lines denote the
locations of the frontal contact line, $\eta_{f}$, at $\theta=0$ and
$\theta=\pi/3$.

In order to investigate the response of the flow to perturbations with angular
dependence, we develop a numerical method to solve the coupled system for the
evolution of the fluid-fluid interface, $\mathcal{H}(\eta,\theta,\tau)$, and
the pressure at the upper boundary, $\mathcal{P}(\eta,\theta,\tau)$. This
consists of equations (18) and (19) with boundary conditions (22) and (23),
and an appropriate initial condition. We solve this system numerically using a
finite-difference scheme, details of which are given in Appendix A.

The numerical method is applied to a wide range of axisymmetric and non-
axisymmetric initial conditions over many values of the parameters, $G$ and
$M$. In all cases, the solution converges to the axisymmetric similarity
solution. This demonstrates both the veracity of the numerical method and that
the similarity solution provides the intermediate asymptotics for many initial
conditions (cf. Mathunjwa & Hogg, 2006). Examples are shown in figures 3 and 4
illustrating how mode six and mode three fingers are suppressed even though
the input fluid is of lesser viscosity than the ambient. We find that any
perturbation is suppressed, including radial perturbations, provided that the
fluids remain buoyancy segregated (if the interface becomes nonmonotonic as a
function of vertical coordinate then the model breaks down).

The numerical results also demonstrate that there is a jump in the radial
pressure gradient at the upper boundary across the leading contact line (see
figure 5). This pressure gradient jump is associated with the density
difference between the two fluids and the discontinuity vanishes as $G\to 0$.
The discontinuity has a stabilising effect because the magnitude of the
pressure gradient is smaller in the ambient fluid (see §4). However, this
discontinuity in the hydrostatic pressure is not required for the interface to
be stable. Indeed, in §5 we show that the interface is stable even in the
limit $G\to 0$.

##  4 Stability mechanism

In the present section, we describe the mechanism by which buoyancy
segregation suppresses viscous fingering. In the classical case of radial flow
with no vertical dimension (discussed in §1), the pressure gradient either
side of the fluid-fluid interface is proportional to the fluid viscosity owing
to continuity of the radial flux, which drives the well-known instability.

We now analyse the more complicated three-dimensional case in which the fluids
have different densities and are segregated by buoyancy. Between the contact
lines, the pressure gradient in the input fluid is
$\mathrm{d}\mathcal{P}_{0}/\mathrm{d}\eta$ and the pressure gradient in the
ambient fluid is
$\mathrm{d}(\mathcal{P}_{0}-G\mathcal{H}_{0})/\mathrm{d}\eta$, where the
$-G\mathcal{H}_{0}$ contribution arises from the density difference between
the fluids. Since both $\mathrm{d}\mathcal{P}_{0}/\mathrm{d}\eta$ and
$\mathrm{d}\mathcal{H}_{0}/\mathrm{d}\eta$ are negative, the magnitude of the
pressure gradient is larger in the input fluid (in contrast to the classical
case with no vertical dimension), which implies that between the contact lines
the interface is stable for $G>0$. The stability at the contact lines must be
treated separately because there is a discontinuity in
$\mathrm{d}\mathcal{P}_{0}/\mathrm{d}\eta$ and
$\mathrm{d}\mathcal{H}_{0}/\mathrm{d}\eta$ between the input and ambient
fluids as shown in figure 5. We consider the driving pressure gradients in the
input and ambient fluids either side of the contact lines (these locations are
labelled (i)-(iv) in figure 2a).

At location (i), just ahead of the leading contact line, at
$\eta=\eta_{f_{0}}^{+}$, the interface gradient and pressure gradient at the
top boundary are

$\frac{\mathrm{d}\mathcal{H}_{0}}{\mathrm{d}\eta}=0,\qquad\frac{\mathrm{d}\mathcal{P}_{0}}{\mathrm{d}\eta}=-\frac{1}{M\eta_{f_{0}}}.$
(44)

The latter is calculated from (32). The pressure gradient driving the ambient
fluid at (i) is

$\frac{\mathrm{d}(\mathcal{P}_{0}-G\mathcal{H}_{0})}{\mathrm{d}\eta}=-\frac{1}{M\eta_{f_{0}}}.$
(45)

At location (ii), just behind the leading contact line,
$\eta=\eta_{f_{0}}^{-}$, the interface and pressure gradients are

$\frac{\mathrm{d}\mathcal{H}_{0}}{\mathrm{d}\eta}=\frac{1}{\eta_{f_{0}}MG}-\frac{\eta_{f_{0}}}{2G},\qquad\frac{\mathrm{d}\mathcal{P}_{0}}{\mathrm{d}\eta}=-\frac{\eta_{f_{0}}}{2},$
(46)

The former is calculated by taking $\mathcal{H}\to 0$ in (33). There is a
discontinuity in the pressure gradient at the top boundary,
$\mathrm{d}\mathcal{P}_{0}/\mathrm{d}\eta$ across $\eta=\eta_{f_{0}}$, which
can be observed in the numerical results in figure 5. This discontinuity
vanishes as $G\to 0$ even for fluids with different viscosities; the jump in
the pressure gradient associated with the viscosity contrast is eliminated by
the radial extension of the interface. The pressure gradient driving the input
fluid at (ii) is given by (46)b, and its magnitude is greater than (or equal
to) the magnitude of the pressure gradient in the ambient fluid (45) because
$\eta_{f_{0}}\geq\sqrt{2/M}$ (see 38). Hence small perturbations to the
contact line decay because they experience unfavourable pressure gradients,
which suggests interfacial stability provided that $G>0$. The driving pressure
gradients either side of the contact line are equal in the limit $G\to 0$;
this case is discussed in more detail in §5.

Similar analysis applies at the trailing contact line. At location (iii), just
downstream of the trailing contact line, $\eta=\eta_{b_{0}}^{+}$, the
interface gradient, the pressure gradient at the upper boundary, and the
pressure gradient driving the ambient fluid are given by

$\frac{\mathrm{d}\mathcal{H}_{0}}{\mathrm{d}\eta}=\frac{\eta_{b_{0}}}{2MG}-\frac{1}{G\eta_{b_{0}}},\qquad\frac{\mathrm{d}\mathcal{P}_{0}}{\mathrm{d}\eta}=-\frac{1}{\eta_{b_{0}}},\qquad\frac{\mathrm{d}(\mathcal{P}_{0}-G\mathcal{H}_{0})}{\mathrm{d}\eta}=-\frac{\eta_{b_{0}}}{2M},$
(47)

respectively. At location (iv), just behind the trailing contact line,
$\eta=\eta_{b_{0}}^{-}$, the interface gradient and the pressure gradient at
the upper boundary (which drives the input fluid) are given by

$\frac{\mathrm{d}\mathcal{H}_{0}}{\mathrm{d}\eta}=0,\qquad\frac{\mathrm{d}\mathcal{P}_{0}}{\mathrm{d}\eta}=-\frac{1}{\eta_{b_{0}}}.$
(48)

The magnitude of the pressure gradient driving the input fluid,
$1/\eta_{b_{0}}$, is greater than (or equal to) the magnitude of the pressure
gradient driving the ambient fluid, $\eta_{b_{0}}/(2M)$ since
$\eta_{b_{0}}\leq\sqrt{2M}$ with equality when $G=0$. This suggests that the
trailing contact line is stable to small perturbations for $G>0$.

It is worth noting that for small $G$, at the location where
$\mathcal{H}_{0}=0$, the pressure gradient in the ambient fluid (44)b is
proportional to the relative viscosity of the ambient fluid divided by the
radial distance and similarly at $\mathcal{H}_{0}=1$, the pressure gradient in
the input fluid (48)b is proportional to the relative viscosity of the input
fluid divided by the radial distance. This dependence on the relative
viscosity and radial location is as in the classical case (see 1). However,
the key difference is that, in contrast to the classical instability, the
contact lines are separated so that the radial coordinate is different at
$\mathcal{H}_{0}=0$ and $\mathcal{H}_{0}=1$. The change in the pressure
gradient between location (iv) and location (i) (figure 2a) is controlled by
the viscosity contrast (relatively larger ambient viscosity leads to a larger
increase in the magnitude of the pressure gradient) and how far apart the
contact lines are, which reduces the magnitude of the pressure gradient
downstream. These two effects act in opposition. If the intrusion of input
fluid has a relatively short radial extent then the flow is unstable and the
intrusion grows radially until the distance between the contact lines acts to
stabilise the interface. Hence the self-similar axisymmetric solutions have
larger radial extent at smaller viscosity ratios. In summary, the radially
extensive intrusion of the interface along the upper boundary associated with
buoyancy segregation counteracts the destabilising discontinuity in the
pressure gradient associated with the viscosity contrast. The vertically-
segregated radial intrusion could be thought of as a single axisymmetric
viscous finger, which is stable to angular-dependent fingers.

We have shown that it is not gravity that acts against the pressure gradient
to stabilise the contact lines. Indeed, the radial extent of the interface is
insensitive to the size of $G$ for small $G$. The role of gravity is simply to
segregate the fluids and stability occurs after the segregation has occurred.
In section 5, we show that provided that buoyancy has segregated the fluids,
the small-$G$ axisymmetric self-similar solutions are linearly stable. This
formally confirms the results in the present section.

##  5 Linear stability for small density difference ($G\ll 1$)

In the present section, we demonstrate the linear stability of the $G\ll 1$
axisymmetric self-similar solutions (given by 35). Although buoyancy does not
appear in the form of the solution, some buoyancy is required to segregate the
fluids. Confirming the stability of the $G\ll 1$ solution will demonstrate
stability for $G>0$ as discussed in §4. The linear stability analysis also
gives the rate of decay of each mode and its dependence on the viscosity
ratio.

We consider $\theta$-dependent perturbations to the axisymmetric self-similar
solutions of the form

$\displaystyle\mathcal{H}$
$\displaystyle=\mathcal{H}_{0}(\eta)+\epsilon\mathcal{H}_{1}(\eta)e^{\sigma\tau+in\theta},$
(49) $\displaystyle\eta_{f}$
$\displaystyle=\eta_{f_{0}}+\epsilon\eta_{f_{1}}e^{\sigma\tau+in\theta},$ (50)
$\displaystyle\eta_{b}$
$\displaystyle=\eta_{b_{0}}+\epsilon\eta_{b_{1}}e^{\sigma\tau+in\theta},$ (51)
$\displaystyle\mathcal{P}$
$\displaystyle=\mathcal{P}_{0}(\eta)+\epsilon\mathcal{P}_{1}(\eta)e^{\sigma\tau+in\theta},$
(52)

where $\epsilon\ll 1$ and $n$ is an integer. We seek to determine the
stability of these perturbations. Note that the perturbation corresponds to
$e^{\sigma\tau}\sim T^{\sigma}$ in terms of real time $T$ (see 16).

Often linear stability analyses of viscous gravity currents requires rescaling
the spatial domain owing to the singular behaviour of the interface at the
contact lines (e.g. Mathunjwa & Hogg, 2006; Kowal & Worster, 2019). For the
present porous gravity current, the interface is linear at the contact lines
so such a transformation is not required. Instead, we linearise the boundary
conditions about the leading order locations of the contact lines,
$\eta_{f}=\eta_{f_{0}}$ and $\eta_{b}=\eta_{b_{0}}$. We consider perturbations
with $\theta$-dependence ($n\geq 1$) first as the stability in this case has
not yet been investigated. Note that global mass conservation of the input
fluid (29) is identically satisfied by the form of the perturbations for
$n\geq 1$. The case of axisymmetric perturbations ($n=0$) is included at the
end of this section for completeness.

In the single-phase regions (upstream of the trailing contact line and
downstream of the leading contact line), the pressure satisfies Laplace’s
equation and the hence pressure perturbation is given by the solution to

$\frac{1}{\eta}\frac{\mathrm{d}}{\mathrm{d}\eta}\Bigg{(}\eta\frac{\mathrm{d}\mathcal{P}_{1}}{\mathrm{d}\eta}\Bigg{)}-\frac{n^{2}}{\eta^{2}}\mathcal{P}_{1}=0.$
(53)

Given that the pressure perturbation remains finite as $\eta\to\infty$ and as
$\eta\to 0$, the solution in the single-phase regions is

$\mathcal{P}_{1}=c_{n}\eta^{n}\qquad\text{for}\qquad\eta<\eta_{b_{0}}$ (54)
$\mathcal{P}_{1}=d_{n}\eta^{-n}\qquad\text{for}\qquad\eta>\eta_{f_{0}}$ (55)

where $c_{n}$ and $d_{n}$ are constants.

After some algebra, the linearised governing equations (18), (19) in the
interface region ($\eta_{b_{0}}<\eta<\eta_{f_{0}}$) become

$(M-1)\frac{\mathrm{d}}{\mathrm{d}\eta}\big{(}\eta\mathcal{H}_{1}\big{)}+2M\Bigg{(}\frac{\mathrm{d}^{2}\mathcal{P}_{1}}{\mathrm{d}\eta^{2}}-\frac{n^{2}}{\eta^{2}}\mathcal{P}_{1}\Bigg{)}=0,$
(56)
$(1-M)(\sigma+1/2)\mathcal{H}_{1}=-\frac{M}{\eta}\frac{\mathrm{d}\mathcal{P}_{1}}{\mathrm{d}\eta}.$
(57)

We can then eliminate $\mathcal{H}_{1}$ and obtain a single equation for
$\mathcal{P}_{1}$ in the interface region between the contact lines:

$\Bigg{(}\frac{\sigma+1}{\sigma+1/2}\Bigg{)}\frac{\mathrm{d}^{2}\mathcal{P}_{1}}{\mathrm{d}\eta^{2}}-\frac{n^{2}}{\eta^{2}}\mathcal{P}_{1}=0.$
(58)

Figure 6: The dispersion function, $\mathcal{F}(\sigma,n,M)$ (69) plotted
versus the growth rate $\sigma$ for $M=0.5$. It has been normalised by its
maximum value over the range of $\sigma$, $\mathcal{F}_{\mathrm{max}}$. The
zeros of the curves correspond to admissible decay rates $\sigma$ (68) with
the exception of the largest zero (shown as stars). The next largest zero
(shown as circles) will be the decay rate observed (further details in the
text).

To determine the pressure perturbation in the interface region, we solve (58)
with three boundary conditions at each contact line, which arise from the
linearised version of continuity of the pressure, continuity of the flux, and
the kinematic boundary condition (see 24, 25, 26, 28). At
$\eta=\eta_{b_{0}}^{+}$, these boundary conditions take the form

$\displaystyle\mathcal{P}_{1}\rvert_{\eta_{b_{0}}^{+}}$
$\displaystyle=c_{n}\eta_{b_{0}}^{n},$ (59)
$\displaystyle\frac{\mathrm{d}\mathcal{P}_{1}}{\mathrm{d}\eta}\bigg{\rvert}_{\eta_{b_{0}}^{+}}$
$\displaystyle=\frac{\eta_{b_{1}}}{\eta_{b_{0}}^{2}}+nc_{n}\eta_{b_{0}}^{n-1},$
(60) $\displaystyle(\sigma+1/2)\eta_{b_{1}}$
$\displaystyle=-M\frac{\mathrm{d}\mathcal{P}_{1}}{\mathrm{d}\eta}\bigg{\rvert}_{\eta_{b_{0}}^{+}}=-M\Bigg{(}\frac{\eta_{b_{1}}}{\eta_{b_{0}}^{2}}+nc_{n}\eta_{b_{0}}^{n-1}\Bigg{)},$
(61)

where we have used the upstream behaviour (54). At $\eta=\eta_{f_{0}}^{-}$,
the analogous boundary conditions take the form

$\displaystyle\mathcal{P}_{1}\rvert_{\eta_{f_{0}}^{-}}$
$\displaystyle=d_{n}\eta_{f_{0}}^{-n},$ (62)
$\displaystyle\frac{\mathrm{d}\mathcal{P}_{1}}{\mathrm{d}\eta}\bigg{\rvert}_{\eta_{f_{0}}^{-}}$
$\displaystyle=\frac{\eta_{f_{1}}}{2}-nd_{n}\eta_{f_{0}}^{-n-1},$ (63)
$\displaystyle(\sigma+1/2)\eta_{f_{1}}$
$\displaystyle=-\frac{\mathrm{d}\mathcal{P}_{1}}{\mathrm{d}\eta}\bigg{\rvert}_{\eta_{f_{0}}^{-}}=-\Bigg{(}\frac{\eta_{f_{1}}}{2}-nd_{n}\eta_{f_{0}}^{-n-1}\Bigg{)},$
(64)

where we have used the downstream behaviour (55). The system for
$\mathcal{P}_{1}$ between the contact lines is linear and has six boundary
conditions with six unknown constants: $c_{n}$, $d_{n}$, $\eta_{b_{1}}$,
$\eta_{f_{1}}$ and two constants arising from solving the ordinary
differential equation (58). The system is an eigenvalue problem for the growth
rate $\sigma$. The general solution to equation (58) is characterised by the
value of $\sigma$ relative to the critical value

$\sigma_{c}=\frac{-2n^{2}-1}{4n^{2}+1}.$ (65)

Figure 7: Eigenfunctions for the interface perturbation, $\mathcal{H}_{1}$,
where we have set $\mathcal{H}_{1}(\eta=1)=1$. (a) $M=0.2$ and (b) $M=0.5$.
The eigenfunctions are calculated from (57) and (66). The $x$ axis is
$\eta_{b_{0}}<\eta<\eta_{f_{0}}$, where in (a) $\eta_{b_{0}}=0.63$,
$\eta_{f_{0}}=3.16$ and in (b) $\eta_{b_{0}}=1$, $\eta_{f_{0}}=2$.

This critical value of $\sigma$ corresponds to a repeated root in the
auxiliary equation for the ordinary differential equation (58). There are no
non-trivial solutions that satisfy all the boundary conditions for
$\sigma\geq\sigma_{c}$. For $-1<\sigma<\sigma_{c}$, the solution takes the
form

$\mathcal{P}_{1}=\eta^{1/2}\big{[}a_{1}\cos(\alpha\log\eta)+a_{2}\sin(\alpha\log\eta)\big{]},$
(66)

where

$\alpha=\sqrt{-\frac{1}{4}-\bigg{(}\frac{\sigma+1/2}{\sigma+1}\bigg{)}n^{2}},$
(67)

and $a_{1}$ and $a_{2}$ are constants. We apply the six boundary conditions to
obtain the following dispersion relation for $\sigma$

$\mathcal{F}(\sigma,n,M)=0,$ (68)

where

$\displaystyle\mathcal{F}$ $\displaystyle(\sigma,n,M)=$ (69)
$\displaystyle\Big{[}$
$\displaystyle\big{(}\sigma+\textstyle{\frac{1}{2}}\big{)}ns_{b}-(\sigma+1)\big{(}\textstyle{\frac{1}{2}}s_{b}+\alpha
c_{b}\big{)}\Big{]}\Big{[}\big{(}\sigma+\textstyle{\frac{1}{2}}\big{)}nc_{f}+(\sigma+1)\big{(}\textstyle{\frac{1}{2}}c_{f}-\alpha
s_{f}\big{)}\Big{]}-$ $\displaystyle\Big{[}$
$\displaystyle\big{(}\sigma+\textstyle{\frac{1}{2}}\big{)}ns_{f}+(\sigma+1)\big{(}\textstyle{\frac{1}{2}}s_{f}+\alpha
c_{f}\big{)}\Big{]}\Big{[}\big{(}\sigma+\textstyle{\frac{1}{2}}\big{)}nc_{b}-(\sigma+1)\big{(}\textstyle{\frac{1}{2}}c_{b}-\alpha
s_{b}\big{)}\Big{]},$

and

$\displaystyle s_{b}=$ $\displaystyle\sin(\alpha\log\eta_{b_{0}}),$
$\displaystyle s_{f}=$ $\displaystyle\sin(\alpha\log\eta_{f_{0}}),$ (70)
$\displaystyle c_{b}=$ $\displaystyle\cos(\alpha\log\eta_{b_{0}}),$
$\displaystyle c_{f}=$ $\displaystyle\cos(\alpha\log\eta_{f_{0}}).$ (71)

The locations of the contact lines, $\eta_{b_{0}}$ and $\eta_{f_{0}}$, are
functions of the viscosity ratio $M$ and so the dispersion equation (68) also
depends on $M$.

The dispersion equation (68) can be solved to obtain the growth rate, $\sigma$
for any $n\geq 1$ and $0<M<1$. The corresponding eigenfunctions are given by
(66). Since $\sigma<\sigma_{c}<0$, the growth rate is always negative and
hence the axisymmetric similarity solutions are linearly stable.

The dispersion equation (68) can have multiple solutions. The function
$\mathcal{F}(\sigma,n,M)$ is plotted against $\sigma$ in figure 6 for $M=0.5$
and four values of $n$. For each $n$, the largest zero for $\sigma$ is
$\sigma=\sigma_{c}$ (these solutions are indicated by stars in figure 6).
However, the case $\sigma=\sigma_{c}$ corresponds to a trivial solution for
$\mathcal{P}_{1}$, which is not admissible. The next largest solution $\sigma$
of (68) corresponds to the slowest decaying non-trivial solution and these
values are indicated by circles in figure 6. We expect this to be the decay
rate that is observed and we take this solution as the correct value of
$\sigma$.

The eigenfunctions (66) that correspond to this correct solution $\sigma$ are
defined up to a multiplicative constant. The form of the interfacial
perturbation, $\mathcal{H}_{1}$ for each eigenfunction (66) can be calculated
from (57). These interfacial perturbation eigenfunctions are shown in figure 7
for different values of the viscosity ratio, $M$ and the mode, $n$. The shapes
are more oscillatory at higher values of $n$.

Figure 8: The decay rate, $\sigma$, for
$M=0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9$ and $0\leq n\leq 10$. The crosses show
the decay rate observed from our numerical results (see Appendix A.1).

The growth rate, $\sigma$, calculated for different values of $M$ and $n$
using the dispersion relation (68) is plotted using continuous lines in figure
8. The rate of decay is slower for larger values of $n$ and smaller viscosity
ratios, $M$. The crosses in figure 8 denote the decay rate calculated from the
numerical method for $M=0.5$ and $M=0.8$ for three values of $n$ with
$G=0.01$; details of this calculation are given in Appendix A.1. There is
excellent agreement between the numerically-derived prediction for $\sigma$
and the values obtained from the linear stability analysis.

Finally, we consider the case $n=0$. In $\eta<\eta_{b_{0}}$, the pressure
perturbation is constant, $\mathcal{P}_{1}=c_{0}$, whilst in
$\eta>\eta_{f_{0}}$, $\mathcal{P}_{1}=d_{0}$. In the interface region, the
pressure perturbation is linear; $\mathcal{P}_{1}=a+b\eta$. The kinematic
boundary conditions (61, 64) furnish $\sigma=-1$. The thickness perturbation
takes the form (see 57)

$\mathcal{H}_{1}=\frac{2Mb}{1-M}\eta^{-1}.$ (72)

The numerical results for axisymmetric perturbations to the self-similar
solutions corroborated that $\sigma=-1$ (shown as a cross at $n=0$ in figure
6). Axisymmetric perturbations decay faster than $\theta$-dependent
perturbations with the error decaying proportional to $T^{-1}$ as has been
found previously for unconfined gravity currents (Grundy & McLaughlin, 1982;
Mathunjwa & Hogg, 2006). For a different stability argument in the case of
axisymmetric perturbations, see the Appendix of Nordbotten & Celia (2006).

##  6 Vertically varying permeability

Figure 9: Self-similar axisymmetric interface shapes in a layer with
vertically varying permeability. (a) Small-$G$ analytic solutions (83) for
different linear permeability variations (87). (b) Small-$G$ analytic
solutions in the case that there is a shock-like region.

The results obtained in previous sections generalise to layers in which the
permeability varies vertically with $k=k(Z)$. The flow is axisymmetric and
self-similar after buoyancy segregation. Here, we find the analytical
solutions that arise in the small $G$ regime with a less viscous input fluid
($M<1$). One key difference to the uniform case is that shock-like regions in
which the interface is relatively steep can occur even for $M<1$ (see §6.1).
Despite this, we show that at late times the axisymmetric solutions for any
continuous permeability profile, $k(Z)$, are stable.

We define the dimensionless permeability as

$\mathcal{K}(Z/H_{0})=\frac{k(Z)}{\frac{1}{H_{0}}\int_{0}^{H_{0}}k(Z)\,\mathrm{d}Z},$
(73)

where the denominator is the mean permeability. We also define the
dimensionless depth-integrated permeability (Hinton & Woods, 2018),

$\Phi(\mathcal{H})=\int_{0}^{\mathcal{H}}\mathcal{K}(s)\,\mathrm{d}s.$ (74)

Note that $\Phi(0)=0$ and $\Phi(1)=1$. For a uniform layer, $\mathcal{K}\equiv
1$ and $\Phi(\mathcal{H})=\mathcal{H}$. The two mass conservation equations
(18, 19) generalise to

$\frac{\partial\mathcal{H}}{\partial\tau}-\frac{\eta}{2}\frac{\partial\mathcal{H}}{\partial\eta}-\tilde{\nabla}\cdot\big{(}\Phi(\mathcal{H})\tilde{\nabla}\mathcal{P}\big{)}=0,$
(75)
$\tilde{\nabla}\cdot\Big{[}\big{(}M+(1-M)\Phi(\mathcal{H})\big{)}\tilde{\nabla}\mathcal{P}-GM(1-\Phi(\mathcal{H}))\tilde{\nabla}\mathcal{H}\Big{]}=0.$
(76)

This model applies provided that the fluids have become segregated by
buoyancy. In a uniform layer this applies at times given by (42) with $k$
representing the vertical permeability. In a layer with vertically varying
permeability and less viscous input fluid ($M<1$), the dimensional time for
buoyancy segregation is adjusted to

$T\gg\frac{\mu_{a}}{\Delta\rho
g}\int_{0}^{H_{0}}\frac{1}{k_{v}(Z)}\mathrm{d}Z,$ (77)

where $k_{v}(Z)$ is the vertical permeability. In an isotropic, vertically
heterogeneous layer, $k_{v}(Z)=k(Z)$. The time (77) becomes singular if the
permeability is zero at any height in the layer (and in this case our model
does not apply). The time for the shallow approximation to apply is given by
(43).

As for a uniform layer, there are self-similar axisymmetric solutions,
$\mathcal{H}_{0}=\mathcal{H}_{0}(\eta)$, $\mathcal{P}=\mathcal{P}_{0}(\eta)$
to (75), (76), with the pressure gradient at the top boundary given by

$\frac{\mathrm{d}\mathcal{P}_{0}}{\mathrm{d}\eta}=\frac{-1}{\eta\big{(}M+(1-M)\Phi(\mathcal{H}_{0})\big{)}}+\frac{GM(1-\Phi(\mathcal{H}_{0}))}{\big{(}M+(1-M)\Phi(\mathcal{H}_{0})\big{)}}\frac{\mathrm{d}\mathcal{H}_{0}}{\mathrm{d}\eta},$
(78)

and $\mathcal{H}_{0}$ satisfies the following ordinary differential equation

$-\frac{\eta}{2}\frac{\mathrm{d}\mathcal{H}_{0}}{\mathrm{d}\eta}+\frac{1}{\eta}\frac{\mathrm{d}}{\mathrm{d}\eta}\Bigg{(}\frac{\Phi(\mathcal{H}_{0})}{M+(1-M)\Phi(\mathcal{H}_{0})}\Bigg{)}=\frac{1}{\eta}\frac{\mathrm{d}}{\mathrm{d}\eta}\Bigg{(}\eta\frac{GM\Phi(\mathcal{H}_{0})(1-\Phi(\mathcal{H}_{0}))}{M+(1-M)\Phi(\mathcal{H}_{0})}\frac{\mathrm{d}\mathcal{H}_{0}}{\mathrm{d}\eta}\Bigg{)}.$
(79)

The interface gradient and pressure gradient at the leading contact line are

$\frac{\mathrm{d}\mathcal{H}_{0}}{\partial\eta}=\frac{1}{\eta_{f_{0}}MG}-\frac{\eta_{f}}{2G\mathcal{K}(0)},\qquad\frac{\mathrm{d}\mathcal{P}_{0}}{\partial\eta}=\frac{-\eta_{f_{0}}}{2\mathcal{K}(0)},$
(80)

respectively. Whilst at the trailing contact line

$\frac{\mathrm{d}\mathcal{H}_{0}}{\mathrm{d}\eta}=\frac{\eta_{b_{0}}}{2MG\mathcal{K}(1)}-\frac{1}{\eta_{b_{0}}G},\qquad\frac{\mathrm{d}\mathcal{P}_{0}}{\mathrm{d}\eta}=\frac{-1}{\eta_{b_{0}}}.$
(81)

The small-$G$ analytic solutions are given by neglecting the right-hand side
of (79) from which we obtain the following implicit solution

$\displaystyle\mathcal{H}_{0}$ $\displaystyle=1,$ $\displaystyle 0<$
$\displaystyle\eta<\eta_{b_{0}},$ (82) $\displaystyle\eta$
$\displaystyle=\frac{\sqrt{2M\mathcal{K}(\mathcal{H}_{0})}}{M+(1-M)\Phi(\mathcal{H}_{0})},$
$\displaystyle 0<$ $\displaystyle\mathcal{H}_{0}<1,$ (83)
$\displaystyle\mathcal{H}_{0}$ $\displaystyle=0,$ $\displaystyle\eta_{f_{0}}<$
$\displaystyle\eta,$ (84)

where the contact lines are

$\eta_{f_{0}}=\sqrt{\frac{2\mathcal{K}(0)}{M}},\qquad\eta_{b_{0}}=\sqrt{2M\mathcal{K}(1)}.$
(85)

The associated pressure gradient in the interface region is given by

$\frac{\mathrm{d}\mathcal{P}_{0}}{\mathrm{d}\eta}=-\frac{1}{\sqrt{2M\mathcal{K}(\mathcal{H}_{0})}}.$
(86)

In a uniform layer, the interface (83) is monotonic whenever $M<1$ but in a
layer with a vertical permeability variation, turning points can arise even
when the input fluid is less viscous than the ambient (Hinton & Woods, 2018).
The solution (83) is only valid when the interface is monotonic. Otherwise,
the interface has a turning point with buoyant input fluid lying below denser
ambient fluid. This would invalidate the model, which assumed that the fluids
have segregated owing to buoyancy. For example, for a linear permeability
structure with dimensionless variation, $\Delta k$,

$\mathcal{K}(s)=1+\Delta k(s-1/2),$ (87)

the interface (83) is monotonic if and only if (Hinton & Woods, 2018)

$M<\frac{(2-\Delta k)^{2}}{4-2\Delta k+\Delta k^{2}}.$ (88)

Provided (83) is monotonic, the interface forms an axisymmetric self-similar
intrusion, which extends along the upper boundary in a qualitatively identical
fashion to the case of a uniform layer (see figure 9a). In the case that (83)
has a turning point, a shock must be introduced; this regime is studied in
§6.1.

We apply our numerical method to layers with vertical variations in
permeability for a wide range of initial conditions and find that the self-
similar axisymmetric solutions are stable to both $\theta$-dependent and
axisymmetric perturbations for any $G>0$. Next, we generalise §4 to show how
the pressure gradients suppress instabilities in a layer with vertically
varying permeability (but no shock-like regions; for that case, see §6.1).
Between the contact lines, the magnitude of the pressure gradient in the
ambient fluid is less than or equal to that in the input fluid owing to the
contribution of the term $-G\partial\mathcal{H}_{0}/\partial\eta$ (as in §4).
The interface gradient at the contact lines is discontinuous and the stability
there is treated separately.

We can bound the location of the contact lines of the self-similar
axisymmetric solution. In the case that the small-$G$ analytic solution (83)
is monotonic, the contacts lines for $G>0$ satisfy

$\eta_{f_{0}}\geq\sqrt{\frac{2\mathcal{K}(0)}{M}},\qquad\eta_{b_{0}}\leq\sqrt{2M\mathcal{K}(1)},$
(89)

with equality as $G\to 0^{+}$. These inequalities reflect the action of
buoyancy to extend the interface.

The gradient of the interface and the pressure at the upper boundary at the
contact lines are given by (80) and (81). Thus, the pressure gradient in the
ambient fluid just ahead of the leading contact line is given by

$\frac{\mathrm{d}(\mathcal{P}_{0}-G\mathcal{H}_{0})}{\mathrm{d}\eta}=-\frac{1}{M\eta_{f_{0}}},$
(90)

whilst in the input fluid just upstream of the leading contact line, the
pressure gradient is

$\frac{\mathrm{d}\mathcal{P}_{0}}{\mathrm{d}\eta}=-\frac{\eta_{f_{0}}}{2\mathcal{K}(0)},$
(91)

and the inequality (89)a ensures that the magnitude of the pressure gradient
is larger in the input fluid, which stabilises the interface, with equality as
$G\to 0^{+}$.

The pressure gradient in the ambient fluid just ahead of the trailing contact
line is given by

$\frac{\mathrm{d}(\mathcal{P}_{0}-G\mathcal{H}_{0})}{\mathrm{d}\eta}=-\frac{\eta_{b_{0}}}{2M\mathcal{K}(1)},$
(92)

whilst in the input fluid just upstream of the trailing contact line, the
pressure gradient is

$\frac{\mathrm{d}\mathcal{P}_{0}}{\mathrm{d}\eta}=-\frac{1}{\eta_{b_{0}}},$
(93)

and the inequality (89)b ensures that the magnitude of the pressure gradient
is larger in the input fluid, which stabilises the interface.

###  6.1 Shock-like regions

In the case in which the small-$G$ interface (83) has a turning point, a shock
must be introduced to vertically segregate the fluids. This shock may occupy
part or all of the layer (see figure 9b). By way of an example, we consider a
linear variation in permeability with $\Delta k=1.5$ and viscosity ratio,
$M=0.3$. The axisymmetric self-similar profile has a shock-like region near
the top of the layer, even though $M<1$ (red line in figure 9b). The shock at
the top of the layer is at location $\eta=\eta_{s}$ and is of thickness
$\mathcal{H}_{s}$, which satisfy

$\eta_{s}=\sqrt{2F^{\prime}(\mathcal{H}_{s})},\qquad\eta_{s}^{2}\mathcal{H}_{s}=2F(\mathcal{H}_{s}),\quad\text{where}\quad
F(\mathcal{H})=\frac{\Phi(\mathcal{H})}{M+(1-M)\Phi(\mathcal{H})}.$ (94)

These arise from continuity with the solution (83), which is valid in
$\mathcal{H}>\mathcal{H}_{s}$, and mass conservation, respectively. The shock-
like region may also extend across the entire layer (e.g. black line in figure
9b). For more complicated permeability profiles, there can be multiple shock-
like regions separated by rarefaction-like regions (Hinton & Woods, 2018).

Figure 10: Shutdown of mode five fingers in a layer with vertically varying
permeability. We use a linear variation with $\Delta k=1.5$ and viscosity
ratio, $M=0.3$, and $G=0.05$. Flow thickness at three radial cross-sections at
(a) $\tau=0.5$, (b) $\tau=1.5$ and (c) $\tau=20$ (results calculated
numerically; see §3). In (c), the black dashed line shows the axisymmetric
small-$G$ solution with a shock-like region in the upper portion of the layer.
(d) Location of the trailing contact line $\eta=\eta_{b}(\theta,\tau)$ at unit
time intervals between $\tau=0$ and $\tau=20$ (lighter colours correspond to
later times).

We apply our numerical method to layers with shock-like regions and find that
the self-similar axisymmetric solutions are stable to both $\theta$-dependent
and axisymmetric perturbations. An example is shown in figure 10 with a linear
variation in permeability with $\Delta k=1.5$, viscosity ratio $M=0.3$, and
$G=0.05$, demonstrating that mode five fingers are suppressed. The small-$G$
analytic solution is shown as a dashed line in in figure 10c; there is a
shock-like region in the top part of the layer. For $G>0$, the shock-like
regions are smoothed by buoyancy so that they are not vertical (Hinton &
Woods, 2018). The smoothed shock-like region has radial extent proportional to
$G$ in the similarity coordinate $\eta$. Hence, in dimensional terms, the
aspect ratio of the shock evolves with

$\frac{L}{H_{0}}\sim G\Bigg{(}\frac{QT}{2\pi\phi H_{0}^{3}}\Bigg{)}^{1/2},$
(95)

and eventually the shock-like region of the flow is long and thin (when the
right-hand side becomes large) and the shallow flow approximation applies at
sufficiently late times. This smoothing of the shock-like region also occurs
when the shock encompasses the entire thickness of the layer (e.g. figure 9b).
Thus the fluid-fluid interface is not cylindrical in the case. Instead, due to
buoyancy, the interface becomes gradually shallower over time, which
suppresses the classical instability.

The stability analysis generalises to the case in which there are shock-like
regions. At the top of the layer, the leading contact line for $G>0$ satisfies

$\eta_{f_{0}}\geq\eta_{s}\geq\sqrt{\frac{2\mathcal{K}(0)}{M}}.$ (96)

The first inequality arises from buoyancy smoothing the interface beyond the
shock location (see figure 10c) whilst the second arises because the profile
(83) has a turning point in the shock-like region and the shock conserves mass
so it must be ahead of this profile at $Z=0$. Hence (89)a applies even when
there is a shock-like region at the top of the layer. A similar argument
applies when there is a shock-like region at the bottom of the layer and (89)b
remains correct.

The gradient of the interface and the pressure at the upper boundary at the
contact lines are given by (80) and (81), even when there are shock-like
regions. Hence the stability argument of §6 applies.

In summary, the self-similar axisymmetric solutions in a layer with vertically
varying permeability are stable owing to the same mechanism as in a uniform
layer. The less viscous input fluid intrudes into the ambient fluid
sufficiently far that the destabilising increase in the magnitude of the
pressure gradient associated with the viscosity contrast is smoothed over the
extent of the interface. Buoyancy segregation corresponds to a monotonic
interface and the input fluid may be forced into the low permeability region
(see also Debbabi et al., 2018). We note that there may be significant
intermingling of the fluids in the transient evolution to the buoyancy-
segregated flow (e.g. Huppert et al., 2013). There may also be a competition
between viscous fingering and heterogeneity-driven fingering prior to buoyancy
segregation (De Wit & Homsy, 1997).

##  7 Discussion and conclusion

We have examined the stability of the horizontal flow of one fluid injected
into another fluid of greater viscosity and density. When there is a sharp
interface and the two fluids have segregated owing to buoyancy, the flow
evolves in an axisymmetric self-similar fashion. Whilst Saffman-Taylor fingers
may occur prior to buoyancy segregation, we have demonstrated that the
buoyancy-segregated self-similar flow is stable to both axisymmetric and
angular-dependent perturbations. The input fluid forms an intrusion with large
radial extent neighbouring the upper boundary. This disperses the
destabilising pressure gradient jump associated with the viscosity contrast
between the fluids that drives the classical instability. The mechanism
generalises to layers with a vertical variation in permeability and
anisotropic layers with different horizontal and vertical permeability. After
buoyancy segregation, the flow in a heterogeneous layer is always stable even
though the interface may contain steep shock-like regions owing to slow flow
in the low permeability zone. The time for buoyancy segregation can be long if
the porous layers has zones of very low permeability.

In future work, these ideas could be extended to viscous displacements in
horizontal channels as is relevant to mantle plumes (Schoonman et al., 2017).
The base self-similar flow was found for a line source and a point source by
Zheng et al. (2015b) and Hinton (2020) but the stability has not yet been
investigated. The analysis would be more complicated than the present work
owing to the effect of the no-slip boundaries and the shear flow (Snyder &
Tait, 1998; John et al., 2013).

In the context of CO2 sequestration, viscous fingering is undesirable as it
may reduce storage efficiency. Our study could be usefully extended to
consider strategies to vary the injection rate to avoid viscous fingering.
Given that the interface is stable after buoyancy segregation (even for
relatively high injection rates) one could analyse how to vary the input flux
during the pre-segregation transient to ensure that fingers never occur.

## Acknowledgements

E.M.H. is grateful to the University of Melbourne for the award of a Harcourt-
Doig research fellowship.

## Declaration of interests

The authors report no conflict of interest.

##  Appendix A Numerical method

In this appendix, the numerical method for calculating the evolution of the
fluid-fluid interface and the pressure at the upper boundary is described. The
two dependent variables, $\mathcal{H}(\eta,\theta,\tau)$, and
$\mathcal{P}(\eta,\theta,\tau)$ are calculated from equations (18) and (19)
with boundary conditions (22) and (23) using a finite difference method.

The numerical scheme requires initial data for the interface shape at the
reference time ($\tau=0$). The initial condition was chosen to satisfy the
mass conservation condition (29). The pressure at the upper boundary,
$\mathcal{P}$, and interface shape, $\mathcal{H}$, must satisfy (19) initially
since the fluid is incompressible and the upper and lower boundaries are
impermeable. They must also satisfy the boundary conditions (22) and (23).
Admissible choices for the initial conditions are constructed by first
selecting any $\mathcal{H}$, which satisfies the boundary conditions and then
solving (19) (described below) to obtain an initial form for $\mathcal{P}$.

The independent variables are discretised on an annular domain with
$(\eta,\theta,\tau)\in[\eta_{l},\eta_{r}]\times[0,2\pi)\times[0,\tau_{1}]$,
where $\tau_{1}$ is the time to which the simulation is run, and the inner and
outer radii are chosen such that the domain fully incorporates the contact
lines and $\eta_{l}$ is taken to be small but non-zero so that the boundary
condition (22) may be applied. Similarly $\eta_{r}$ is chosen to be large and
(23) is applied. Typically we used $\eta_{l}=0.01$ and $\eta_{r}=10$ and we
confirmed that changes to these values led to imperceptible differences to the
results obtained. The discretised variables are

$\displaystyle\eta_{i}$ $\displaystyle=\eta_{l}+i\Delta\eta,$ $\displaystyle
i$ $\displaystyle=0\ldots L,$ $\displaystyle\Delta\eta$
$\displaystyle=\frac{\eta_{r}-\eta_{l}}{L},$ (97) $\displaystyle\theta_{j}$
$\displaystyle=j\Delta\theta,$ $\displaystyle j$ $\displaystyle=0\ldots M-1,$
$\displaystyle\Delta\theta$ $\displaystyle=\frac{2\pi}{M},$ (98)
$\displaystyle\tau_{k}$ $\displaystyle=k\Delta\tau,$ $\displaystyle k$
$\displaystyle=0\ldots N,$ $\displaystyle\Delta\tau$
$\displaystyle=\frac{\tau_{1}}{N},$ (99)

and we write $\mathcal{H}_{i,j,k}$ and $\mathcal{P}_{i,j,k}$ to denote the
approximations of the dependent variables.

The numerical procedure is as follows. First, the given initial condition for
$\mathcal{H}$ and $\mathcal{P}$ is discretized. At each subsequent time step,
the interface shape is updated via an adapted midpoint (second-order Runge-
Kutta) method. The interface height at the midpoint between timesteps is given
by (using 18)

$\mathcal{H}_{i,j,k+1/2}=\mathcal{H}_{i,j,k}+\frac{\Delta\tau}{2}\Bigg{[}\frac{\eta}{2}\frac{\partial\mathcal{H}}{\partial\eta}+\nabla\cdot(\mathcal{H}\nabla\mathcal{P})\Bigg{]}_{i,j,k},$
(100)

where the term in square brackets is calculated using central differences in
the two spatial coordinates. The pressure at the midpoint,
$\mathcal{P}_{i,j,k+1/2}$ is obtained from (19) using
$\mathcal{H}=\mathcal{H}_{i,j,k+1/2}$ and applying a five-point difference
formula and the relaxation method to solve the steady problem. The
timestepping is completed with

$\mathcal{H}_{i,j,k+1}=\mathcal{H}_{i,j,k}+\Delta\tau\Bigg{[}\frac{\eta}{2}\frac{\partial\mathcal{H}}{\partial\eta}+\nabla\cdot(\mathcal{H}\nabla\mathcal{P})\Bigg{]}_{i,j,k+1/2},$
(101)

and $\mathcal{P}_{i,j,k+1}$ is then obtained in an identical fashion to
$\mathcal{P}_{i,j,k+1/2}$. We use $\Delta\tau=0.1\Delta\eta^{2}$.

Various checks were used to verify the accuracy of the numerical method.
First, we confirmed that mass conservation (29) was satisfied to within $0.05$
% at all times over a large range of the parameters; $0<M<10$ and $0<G<10$.
Second, we used the axisymmetric similarity solution as the initial condition
and found that the solution was steady. Finally, we confirmed that the results
always converged to the axisymmetric similarity solution.

###  A.1 Estimating the decay rate ($\sigma$)

In this section we describe how the decay rate, $\sigma$, of linear
perturbations to the axisymmetric self-similar flow can be estimated using the
numerical method. We use the following initial condition for the numerical
method

$\mathcal{H}(\eta,\theta,\tau=0)=\mathcal{H}_{0}\bigg{(}\frac{\eta}{1+0.1\cos(n\theta)}\bigg{)},$
(102)

for mode $n$ perturbations, where $\mathcal{H}_{0}(\eta)$ is the axisymmetric
self-similar solution. We integrate forward in time as described in Appendix
A. The error between the resultant numerical solution,
$\mathcal{H}(\eta,\theta,\tau)$ and the axisymmetric solution,
$\mathcal{H}_{0}(\eta)$ is calculated at each time via

$E(\tau)=\int_{0}^{\eta_{r}}\int_{0}^{2\pi}\sqrt{\mathcal{H}^{2}-\mathcal{H}_{0}^{2}}\,\mathrm{d}\theta\mathrm{d}\eta.$
(103)

At sufficiently late times, the linear stability analysis of §5 applies. To
predict the decay rate, $\sigma$, we best-fit a straight line to $\log E$ as a
function of $\tau\in[5,10]$. Some of the results for $G=0.01$ are shown as
crosses in figure 8.

## References

  * Al-Housseiny et al. (2012) Al-Housseiny, T.T., Tsai, P.A. & Stone, H.A. 2012 Control of interfacial instabilities using flow geometry. Nat. Phys. 8 (10), 747–750.
  * Bachu (2015) Bachu, S. 2015 Review of CO2 storage efficiency in deep saline aquifers. Int. J. Greenh. Gas Control. 40, 188–202.
  * Benham et al. (2022) Benham, G.P., Neufeld, J.A. & Woods, A.W. 2022 Axisymmetric gravity currents in anisotropic porous media. arXiv preprint arXiv:2202.01165 .
  * Bernoff & Witelski (2002) Bernoff, A.J. & Witelski, T.P. 2002 Linear stability of source-type similarity solutions of the thin film equation. Appl. Math. Lett. 15 (5), 599–606.
  * Chertock (2002) Chertock, A. 2002 On the stability of a class of self-similar solutions to the filtration-absorption equation. Eur. J. Appl. Math. 13 (2), 179–194.
  * De Wit & Homsy (1997) De Wit, A. & Homsy, G.M. 1997 Viscous fingering in periodically heterogeneous porous media. II. Numerical simulations. J. Chem. Phys. 107 (22), 9619–9628.
  * Debbabi et al. (2018) Debbabi, Y., Jackson, M.D., Hampson, G.J. & Salinas, P. 2018 Impact of the buoyancy–viscous force balance on two-phase flow in layered porous media. Transp. Porous Media 124 (1), 263–287.
  * Grundy & McLaughlin (1982) Grundy, R.E. & McLaughlin, R. 1982 Eigenvalues of the Barenblatt-Pattle similarity solution in nonlinear diffusion. Proc. R. Soc. A: Math. Phys. Eng. Sci. 383 (1784), 89–100.
  * Guo et al. (2016) Guo, B., Zheng, Z., Celia, M.A. & Stone, H.A. 2016 Axisymmetric flows from fluid injection into a confined porous medium. Phys. Fluids 28 (2), 022107.
  * Hesse et al. (2007) Hesse, M.A., Tchelepi, H.A., Cantwel, B.J. & Orr, F.M. 2007 Gravity currents in horizontal porous layers: transition from early to late self-similarity. J. Fluid Mech. 577, 363–383.
  * Hinton (2020) Hinton, E.M. 2020 Axisymmetric viscous flow between two horizontal plates. Phys. Fluids 32 (6), 063104.
  * Hinton & Woods (2018) Hinton, E.M. & Woods, A.W. 2018 Buoyancy-driven flow in a confined aquifer with a vertical gradient of permeability. J. Fluid Mech. 848, 411–429.
  * Homsy (1987) Homsy, G.M. 1987 Viscous fingering in porous media. Annu. Rev. Fluid Mech. 19 (1), 271–311.
  * Huppert et al. (2013) Huppert, H.E., Neufeld, J.A. & Strandkvist, C. 2013 The competition between gravity and flow focusing in two-layered porous media. J. Fluid Mech. 720, 5–14.
  * Huppert & Pegler (2022) Huppert, H.E. & Pegler, S.S. 2022 The fate of continuous input of relatively heavy fluid at the base of a porous medium. J. Fluid Mech. 932.
  * Huppert & Woods (1995) Huppert, H.E. & Woods, A.W. 1995 Gravity-driven flows in porous layers. J. Fluid Mech. 292, 55–69.
  * John et al. (2013) John, M.O., Oliveira, R.M., Heussler, F.H.C. & Meiburg, E. 2013 Variable density and viscosity, miscible displacements in horizontal Hele-Shaw cells. Part 2. Nonlinear simulations. J. Fluid Mech. 721, 295–323.
  * Juanes et al. (2010) Juanes, R., MacMinn, C.W. & Szulczewski, M.L. 2010 The footprint of the CO2 plume during carbon dioxide storage in saline aquifers: storage efficiency for capillary trapping at the basin scale. Transp. Porous Media 82 (1), 19–30.
  * Kowal (2021) Kowal, K.N. 2021 Viscous banding instabilities: non-porous viscous fingering. J. Fluid Mech. 926, A4.
  * Kowal & Worster (2019) Kowal, K.N. & Worster, M.G. 2019 Stability of lubricated viscous gravity currents. Part 2. Global analysis and stabilisation by buoyancy forces. J. Fluid Mech. 871, 1007–1027.
  * Lyle et al. (2005) Lyle, S., Huppert, H.E., Hallworth, M., Bickle, M. & Chadwick, A. 2005 Axisymmetric gravity currents in a porous medium. J. Fluid Mech. 543, 293–302.
  * Mathunjwa & Hogg (2006) Mathunjwa, J.S. & Hogg, A.J. 2006 Self-similar gravity currents in porous media: linear stability of the Barenblatt–Pattle solution revisited. Eur. J. Mech. B Fluids 25 (3), 360–378.
  * Newman (1984) Newman, W.I. 1984 A Lyapunov functional for the evolution of solutions to the porous medium equation to self-similarity. I. J. Math. Phys. 25 (10), 3120–3123.
  * Nijjer et al. (2018) Nijjer, J.S., Hewitt, D.R. & Neufeld, J.A. 2018 The dynamics of miscible viscous fingering from onset to shutdown. J. Fluid Mech. 837, 520–545.
  * Nijjer et al. (2022) Nijjer, J.S., Hewitt, D.R. & Neufeld, J.A. 2022 Horizontal miscible displacements through porous media: the interplay between viscous fingering and gravity segregation. J. Fluid Mech. 935, A14.
  * Nordbotten & Celia (2006) Nordbotten, J.M. & Celia, M.A. 2006 Similarity solutions for fluid injection into confined aquifers. J. Fluid Mech. 561, 307–327.
  * Paterson (1981) Paterson, L. 1981 Radial fingering in a Hele Shaw cell. J. Fluid Mech. 113, 513–529.
  * Pattle (1959) Pattle, R.E. 1959 Diffusion from an instantaneous point source with a concentration-dependent coefficient. Q. J. Mech. Appl. Math. 12 (4), 407–409.
  * Pegler et al. (2014) Pegler, S.S., Huppert, H.E. & Neufeld, J.A. 2014 Fluid injection into a confined porous layer. J. Fluid Mech. 745, 592–620.
  * Riaz & Meiburg (2003) Riaz, A. & Meiburg, E. 2003 Three-dimensional miscible displacement simulations in homogeneous porous media with gravity override. J. Fluid Mech. 494, 95–117.
  * Saffman & Taylor (1958) Saffman, P.G. & Taylor, G.I. 1958 The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid. Proc. R. Soc. A: Math. Phys. Eng. Sci. 245 (1242), 312–329.
  * Schoonman et al. (2017) Schoonman, C.M., White, N.J. & Pritchard, D. 2017 Radial viscous fingering of hot asthenosphere within the Icelandic plume beneath the North Atlantic Ocean. Earth Planet. Sci. Lett. 468, 51–61.
  * Sharma et al. (2020) Sharma, V., Nand, S., Pramanik, S., Chen, C.-Y. & Mishra, M. 2020 Control of radial miscible viscous fingering. J. Fluid Mech. 884.
  * Snyder & Tait (1998) Snyder, D. & Tait, S. 1998 A flow-front instability in viscous gravity currents. J. Fluid Mech. 369, 1–21.
  * Tchelepi (1994) Tchelepi, H.A. 1994 Viscous fingering, gravity segregation and permeability heterogeneity in two-dimensional and three-dimensional flows. PhD thesis, Stanford University.
  * Zheng et al. (2015a) Zheng, Z., Kim, H. & Stone, H.A. 2015a Controlling viscous fingering using time-dependent strategies. Phys. Rev. Lett. 115 (17), 174501.
  * Zheng et al. (2015b) Zheng, Z., Rongy, L. & Stone, H.A. 2015b Viscous fluid injection into a confined channel. Phys. Fluids 27 (6), 062105.
  * Zheng & Stone (2022) Zheng, Z. & Stone, H.A. 2022 The Influence of Boundaries on Gravity Currents and Thin Films: Drainage, Confinement, Convergence, and Deformation Effects. Annu. Rev. Fluid Mech. 54, 27–56.