\section{Introduction}
\label{sec:intro}
Displacement flows in porous media occur in many important environmental settings including geological CO\textsubscript{2} 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 \citep{saffman1958penetration,paterson1981radial,homsy1987viscous}. 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 CO\textsubscript{2}, which in turn reduces the storage efficiency \citep{bachu2015review}. Multiple strategies have been developed to suppress viscous fingering including controlling the flow geometry and varying the volume flux of the input fluid \citep{al2012control,zheng2015controlling}. Significant miscibility between the fluids has also been shown to shut down fingering at later times \citep{nijjer2018dynamics,sharma2020control}.

In the case that the two fluids have different densities, as is relevant to CO\textsubscript{2} 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' \citep{huppert1995gravity}. There has been much research into gravity currents in porous media, particularly in confined layers in which the ambient fluid must be displaced \citep{hesse2007gravity,juanes2010footprint,zheng2022influence}. 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 \citep{nordbotten2006similarity,pegler2014fluid,guo2016axisymmetric}. It has also been shown that the combination of buoyancy and mixing can suppress the Saffman-Taylor instability at long times \citep{tchelepi1994viscous,riaz2003three,nijjer2022horizontal}. 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 \ref{fig:setupcartoon}). 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 \citep{paterson1981radial}. 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
\begin{equation} \label{eq:classicalflux}
q_r = \frac{-2 \pi r k}{\mu}  \frac{\partial p}{\partial r},
\end{equation}
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 \textit{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 \citep{kowal2019stability}. 
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 \citep{kowal2021viscous}.

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 \citep{pattle1959diffusion,grundy1982eigenvalues,newman1984lyapunov,bernoff2002linear,chertock2002stability,mathunjwa2006self}. 
 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 \citep{nordbotten2006similarity}. 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 \S \ref{sec:model}. The coupled system that governs the background pressure and the interface shape is integrated numerically in \S \ref{sec:numerics} 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 \S \ref{sec:stabilmech}, the physical mechanism that suppresses the Saffman-Taylor instability is analysed. A linear stability analysis is carried out in \S \ref{sec:linstabil} for the case of buoyancy-segregated fluids with an arbitrarily small density difference. In \S \ref{sec:vertvarperm}, 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 \S \ref{sec:conc}.


\begin{figure}
\centering
\includegraphics[width=0.95\columnwidth]{setupcartoon.eps}
\caption{(a) Schematic of the setup. (b) Cross-section of the setup.}
\label{fig:setupcartoon}
\end{figure}

\section{Model} \label{sec:model}
The setup is shown in figure \ref{fig:setupcartoon}. A similar geometry was studied by \citet{nordbotten2006similarity} and \citet{guo2016axisymmetric}. 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,
\begin{equation}
P= \bar{P} + \rho g Z \qquad \text{for} \quad 0<Z<H,
\end{equation}
\begin{equation}
P= \bar{P} -\Delta \rho g H + (\rho+\Delta \rho) g Z \qquad \text{for} \quad H<Z<H_0,
\end{equation}
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
\begin{equation}
\mathbf{U}_i = -\frac{k}{\mu_i} \nabla \bar{P},
\end{equation}
\begin{equation}
\mathbf{U}_a = -\frac{k}{\mu_a}\Big( \nabla \bar{P} - \Delta \rho g \nabla H \Big),
\end{equation}
in the input and ambient fluids, respectively.
The governing equations of the flow are
\begin{equation} \label{eq:dimgov1}
\phi \frac{\partial H}{\partial T} + \nabla \cdot \big( H \mathbf{U}_i \big) = 0
\end{equation}
and 
\begin{equation} \label{eq:dimgov2}
\nabla \cdot \big[ H \mathbf{U}_i + (H_0-H) \mathbf{U}_a  \big] = 0,
\end{equation}
which represent mass conservation integrated over the thickness of the input fluid, and mass conservation integrated over the thickness of the layer, respectively.
\par 
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
\begin{equation} \label{eq:dimglobalmass}
\int_0^{2\pi} \int_0^{R_f(\theta,T)} \phi R H \, \mathrm{d} R\, \mathrm{d} \theta = QT,
\end{equation}
where $R_f(\theta,T)$ denotes the frontal contact line where $H=0$ (see figure \ref{fig:setupcartoon}b).

Eventually, the fluid-fluid interface will touch the lower boundary \citep{nordbotten2006similarity}. 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 \citep{lyle2005axisymmetric,guo2016axisymmetric}. Recently, \citet{benham2022axisymmetric} 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 \citep[see also][]{huppert2022fate}. 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 \ref{fig:setupcartoon}b.

When the input fluid has filled the layer near the source, the boundary conditions are as follows.
At the source,
\begin{equation} \label{eq:dimBCsource}
H \to H_0, \qquad \frac{\partial \bar{P}}{\partial R}  \to \frac{-\mu_i Q}{2 \pi k H_0 R} \qquad \text{as} \quad R \to 0,
\end{equation}
where the latter arises from volume conservation. In the far field, the boundary condition is
\begin{equation} \label{eq:dimBCfar}
H \to 0, \qquad \frac{\partial \bar{P}}{\partial R}  \to \frac{-\mu_a Q}{2 \pi k H_0 R} \qquad \text{as} \quad R \to \infty.
\end{equation}
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 \eqref{eq:dimgov1} and \eqref{eq:dimgov2} combined with global mass conservation \eqref{eq:dimglobalmass} and the boundary conditions \eqref{eq:dimBCsource}, \eqref{eq:dimBCfar} 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 \S \ref{sec:numerics}), 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 \S \ref{sec:intro}). 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
\begin{equation} \label{eq:BCsRbfirst}
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,
\end{equation}
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
\begin{equation} \label{eq:BCsRbsec}
-\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^+}.
\end{equation}
At the leading contact line, $R=R_f(\theta,t)$, the corresponding boundary conditions are
\begin{equation}
H(R=R_f^+)=H(R=R_f^-)=H_0, \qquad \big[\bar{P}\big]^+_-=0,
\end{equation}
\begin{equation}
\Bigg[ \frac{\partial (\bar{P}-\Delta \rho g H)}{\partial R} - \frac{\partial R_f}{\partial \theta} \frac{1}{R} \frac{\partial  (\bar{P}-\Delta \rho g H)}{\partial \theta} \Bigg]^+_-=0,
\end{equation}
\begin{equation}
-\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^-},
\end{equation}
where $[\cdot]^+_-$ denotes the discontinuity in the quantity in square brackets at the contact line $R=R_f(\theta, t)$.


\subsection{Transformed coordinate system} \label{sec:transformed}
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
\begin{equation} \label{eq:etaandtau}
\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),
\end{equation}
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
\begin{equation}
\mathcal{H}(\eta,\theta,\tau) = \frac{H(R,\theta,T)}{H_0}, \qquad \mathcal{P}(\eta,\theta,\tau)=\frac{2 \pi k H_0 \bar{P}(R,\theta,T)}{\mu_i Q}.
\end{equation}
The two mass conservation equations (\ref{eq:dimgov1}, \ref{eq:dimgov2}) become
\begin{equation} \label{eq:localmasscons_sim}
\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,
\end{equation}
\begin{equation} \label{eq:fullmasscons_sim}
\tilde{\nabla} \cdot \Big[\big(M+(1-M)\mathcal{H}\big) \tilde{\nabla} \mathcal{P} - G M (1-\mathcal{H}) \tilde{\nabla} \mathcal{H}   \Big] = 0,
\end{equation}
where
\begin{equation}
\tilde{\nabla}=\Bigg(\frac{\partial}{\partial (\eta \cos \theta)}, \, \frac{\partial}{\partial (\eta \sin \theta)} \Bigg)
\end{equation}
is the gradient operator with respect to the transformed coordinates and we have introduced the following two dimensionless groups
\begin{equation} \label{eq:GandM}
M=\frac{\mu_i}{\mu_a}, \qquad G= \frac{2 \pi k \Delta \rho g H_0^2}{\mu_i Q},
\end{equation}
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
\begin{equation} \label{eq:bc1_sim}
\mathcal{H} \to 1, \qquad \frac{\partial \mathcal{P}}{\partial \eta}  \to \frac{-1}{\eta} \qquad \text{as} \quad \eta \to 0,
\end{equation}
and
\begin{equation} \label{eq:bc2_sim}
\mathcal{H} \to 0, \qquad \frac{\partial \mathcal{P}}{\partial \eta}  \to \frac{-M^{-1}}{\eta} \qquad \text{as} \quad \eta \to \infty.
\end{equation}
The trailing contact line in the transformed coordinates is denoted by $\eta=\eta_b(\theta,\tau)$ and the boundary conditions there (\ref{eq:BCsRbfirst}, \ref{eq:BCsRbsec}) are recast as
\begin{equation} \label{eq:transBCtrail1}
\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,
\end{equation}
\begin{equation}\label{eq:transBCtrail2}
-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^+},
\end{equation}
whilst at the leading contact line, $\eta=\eta_f(\theta,\tau)$,
\begin{equation}\label{eq:transBCfront1}
\mathcal{H}(\eta=\eta_f^+)=\mathcal{H}(\eta=\eta_f^-)=0, \qquad \big[\mathcal{P} \big]^+_-=0,
\end{equation}
\begin{equation}\label{eq:etafbc3}
\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,
\end{equation}
\begin{equation} \label{eq:transBCfront2}
-\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^-}.
\end{equation}
Mass conservation \eqref{eq:dimglobalmass} takes the form
\begin{equation} \label{eq:transformglmasscon}
\frac{1}{2 \pi} \int_0^{2 \pi} \int_0^{\eta_f} \eta \mathcal{H}_0 \, \mathrm{d} \eta \, \mathrm{d} \theta = 1.
\end{equation}



\subsection{Self-similar axisymmetric solution}
The flow has a self-similar axisymmetric solution with \citep{nordbotten2006similarity,guo2016axisymmetric}
\begin{equation}
\mathcal{H}(\eta,\theta,\tau)=\mathcal{H}_0(\eta), \qquad \mathcal{P}(\eta,\theta,\tau) = \mathcal{P}_0(\eta),
\end{equation}
and
\begin{equation}
\eta_b(\theta,\tau) = \eta_{b_0}, \qquad \eta_f(\theta,\tau) = \eta_{f_0},
\end{equation}
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 \S\S \ref{sec:numerics}, \ref{sec:stabilmech} and \ref{sec:linstabil}. 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 \citet{nordbotten2006similarity} and \citet{guo2016axisymmetric}.
For axisymmetric flow, the pressure at the top boundary, $\mathcal{P}_0(\eta)$, can be eliminated from the problem as follows. We integrate equation \eqref{eq:fullmasscons_sim} with respect to $\eta$ and apply the source flux boundary condition \eqref{eq:bc1_sim} to obtain
\begin{equation} \label{eq:P0axi}
\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}.
\end{equation}
Upon substituting \eqref{eq:P0axi} into \eqref{eq:localmasscons_sim}, we obtain the following ordinary differential equation for $\mathcal{H}_0(\eta)$,
\begin{equation} \label{eq:axisymmetricODE}
-\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{G M \mathcal{H}_0(1-\mathcal{H}_0)}{M+(1-M) \mathcal{H}_0}\frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta} \Bigg).
\end{equation}
\begin{figure}
\centering
\includegraphics[width=0.95\columnwidth]{viscshape24jun.eps}
\caption{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 (\ref{eq:h0smallG}, black dashed line).  (b) Small-$G$ approximate interface shapes \eqref{eq:h0smallG} for four values of $M$.}
\label{fig:viscshape}
\end{figure}
This axisymmetric system can be solved numerically for $\mathcal{H}_0(\eta)$ as described in \citet{guo2016axisymmetric} 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$, \eqref{eq:axisymmetricODE} has an approximate analytic solution. We neglect the right-hand side of \eqref{eq:axisymmetricODE}. Although this removes the highest derivative of $\mathcal{H}_0$, the boundary conditions \eqref{eq:transBCfront2} and \eqref{eq:transBCtrail2} 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 \eqref{eq:axisymmetricODE} furnishes the following interface shape \citep{guo2016axisymmetric}
\begin{align}
\mathcal{H}_0&=1,  &   0<&\eta< \eta_{b_0},\\
\mathcal{H}_0 &= \frac{(2M)^{1/2}}{(1-M)\eta} -\frac{M}{1-M}, &  \eta_{b_0}<&\eta< \eta_{f_0}, \label{eq:h0smallG} \\
\mathcal{H}_0&=0, &  \eta_{f_0}<&\eta,
\end{align}
where the contact lines are given by 
\begin{equation}
\eta_{b_0}=\sqrt{2M} \qquad and \qquad \eta_{f_0} =\sqrt{2/M}.
\end{equation}
Note that these provide inner bounds on the locations of the contact lines for general $G> 0$ provided that $M<1$:
\begin{equation} \label{eq:contactpointineq}
\eta_{b_0} \leq \sqrt{2M}, \qquad \eta_{f_0} \geq \sqrt{2/M},
\end{equation}
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
\begin{align}
\mathcal{P}_0&=-1-\log\Bigg(\frac{\eta}{(2M)^{1/2}}\Bigg)  &   0<&\eta< \eta_{b_0},\\
\mathcal{P}_0 &= -\frac{\eta}{(2M)^{1/2}}, &  \eta_{b_0}<&\eta< \eta_{f_0}, \\
\mathcal{P}_0&=-M^{-1} -M^{-1} \log\Bigg(\frac{\eta}{(2/M)^{1/2}}\Bigg) &  \eta_{f_0}<&\eta,
\end{align}
where $\mathcal{P}_0(\eta)$ is defined up to addition of an arbitrary constant.
\par
The numerical solution to \eqref{eq:axisymmetricODE} for $\mathcal{H}_0$ is plotted as a blue line in figure \ref{fig:viscshape}a for $G=0.1$ and $M=0.2$. The small-$G$ analytical solution \eqref{eq:h0smallG} is shown as a black dashed line and there is good agreement with the numerical result. Figure \ref{fig:viscshape}b 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 \eqref{eq:h0smallG} 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 \eqref{eq:h0smallG} becomes a more accurate approximation of the numerical solution, with the input fluid preferentially flowing near the upper boundary (figure \ref{fig:viscshape}). 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 \citep{paterson1981radial}. The derivation of the axisymmetric analytical solution \eqref{eq:h0smallG} 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.
\par
For $G$ small but non-zero, the interfaces in figure \ref{fig:viscshape}b 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
\begin{equation} \label{eq:segregatetime}
T \gg \frac{\mu_a H_0}{\Delta \rho g k},
\end{equation}
whilst the latter requires
\begin{equation} \label{eq:shallowtime}
T \gg \frac{2 \pi \phi H_0^3}{Q}.
\end{equation}
Prior to buoyancy segregation, Saffman-Taylor fingering can occur. The present paper is concerned with the post-segregation stability of the axisymmetric solutions to \eqref{eq:axisymmetricODE} to angular-dependent perturbations with $\Delta \rho>0$.
For $G\ll 1$, the first condition \eqref{eq:segregatetime} implies the second \eqref{eq:shallowtime}. 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 \eqref{eq:segregatetime} will depend only on the vertical permeability, whereas the gravity number (and hence the self-similar solutions) depends only on horizontal permeability \citep[cf.][]{benham2022axisymmetric}. 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 \S \ref{sec:vertvarperm}. The case of more complicated anisotropy such as different horizontal permeabilities in different directions is beyond the scope of the present work.
\par
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 \eqref{eq:segregatetime} and \eqref{eq:shallowtime} apply, the exact source location becomes unimportant.

\section{Numerical results} \label{sec:numerics}
\begin{figure}
\centering
\includegraphics[width=0.95\columnwidth]{fingersuppressionmode6_24jun.eps}
\caption{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$.}
\label{fig:fingersuppression}
\end{figure}


\begin{figure}
\centering
\includegraphics[width=0.95\columnwidth]{fingersuppression_mode3_27jun.eps}
\caption{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 \ref{fig:radialsupressn3}.}
\label{fig:fingersupressn3}
\end{figure}


\begin{figure}
\centering
\includegraphics[width=0.95\columnwidth]{suppression_mode3_h_p_27jun.eps}
\caption{Suppression of mode three fingers for $G=0.5$ and $M=0.5$ (same parameters and initial condition as figure \ref{fig:fingersupressn3}). 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$.}
\label{fig:radialsupressn3}
\end{figure}


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 \eqref{eq:localmasscons_sim} and \eqref{eq:fullmasscons_sim} with boundary conditions \eqref{eq:bc1_sim} and \eqref{eq:bc2_sim}, and an appropriate initial condition. We solve this system numerically using a finite-difference scheme, details of which are given in Appendix \ref{app:numerics}.

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 \citep[cf.][]{mathunjwa2006self}. Examples are shown in figures \ref{fig:fingersuppression} and \ref{fig:fingersupressn3} 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 \ref{fig:radialsupressn3}). 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 \S \ref{sec:stabilmech}). However, this discontinuity in the hydrostatic pressure is not required for the interface to be stable. Indeed, in \S \ref{sec:linstabil} we show that the interface is stable even in the limit $G \to 0$. 


\section{Stability mechanism} \label{sec:stabilmech}
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 \S \ref{sec:intro}), 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 \ref{fig:radialsupressn3}. 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 \ref{fig:viscshape}a).

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
\begin{equation} \label{eq:loc1}
 \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}}.
\end{equation}
The latter is calculated from \eqref{eq:P0axi}.
The pressure gradient driving the ambient fluid at (i) is 
\begin{equation} \label{eq:dpdetaambienti}
\frac{\mathrm{d} (\mathcal{P}_0-G \mathcal{H}_0)}{\mathrm{d} \eta}=-\frac{1}{M \eta_{f_0}}.
\end{equation}
At location (ii), just behind the leading contact line, $\eta=\eta_{f_0}^-$, the interface and pressure gradients are
\begin{equation} \label{eq:dpdetaii}
\frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta} =\frac{1}{\eta_{f_0} M G} -\frac{\eta_{f_0}}{2G}, \qquad \frac{\mathrm{d} \mathcal{P}_0}{\mathrm{d} \eta} = -\frac{\eta_{f_0}}{2},
\end{equation}
The former is calculated by taking $\mathcal{H} \to 0$ in \eqref{eq:axisymmetricODE}.
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 \ref{fig:radialsupressn3}. 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 \eqref{eq:dpdetaii}b, and its magnitude is greater than (or equal to) the magnitude of the pressure gradient in the ambient fluid \eqref{eq:dpdetaambienti} because $\eta_{f_0} \geq \sqrt{2/M}$ (see \ref{eq:contactpointineq}). 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 \S \ref{sec:linstabil}.

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
\begin{equation}
\frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta} =\frac{\eta_{b_0}}{2 M G} - \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},
\end{equation}
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
\begin{equation} \label{eq:loc4}
	\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}}.
\end{equation}
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 \eqref{eq:loc1}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 \eqref{eq:loc4}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 \ref{eq:classicalflux}). 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 \ref{fig:viscshape}a) 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 \ref{sec:linstabil}, 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.


\section{Linear stability for small density difference ($G \ll 1$)} \label{sec:linstabil}
In the present section, we demonstrate the linear stability of the $G \ll 1$ axisymmetric self-similar solutions (given by \ref{eq:h0smallG}). 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 \S \ref{sec:stabilmech}. 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
\begin{align}
\mathcal{H}&=\mathcal{H}_0(\eta) + \epsilon \mathcal{H}_1(\eta) e^{\sigma \tau + i n \theta}, \\
\eta_f &=\eta_{f_0} + \epsilon \eta_{f_1} e^{\sigma \tau + i n \theta}, \\
\eta_b &=\eta_{b_0} + \epsilon \eta_{b_1} e^{\sigma \tau + i n \theta}, \\
\mathcal{P} &= \mathcal{P}_0(\eta) +\epsilon \mathcal{P}_1(\eta) e^{\sigma \tau + i n \theta},
\end{align}
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 \ref{eq:etaandtau}).
\par
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 \citep[e.g.][]{mathunjwa2006self,kowal2019stability}. 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 \eqref{eq:transformglmasscon} 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
\begin{equation}
\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.
\end{equation}
Given that the pressure perturbation remains finite as $\eta \to \infty$ and as $\eta \to 0$, the solution in the single-phase regions is
\begin{equation} \label{eq:p1ups}
\mathcal{P}_1=c_n \eta^n \qquad \text{for} \qquad \eta< \eta_{b_0}
\end{equation}
\begin{equation} \label{eq:p1down}
\mathcal{P}_1=d_n \eta^{-n} \qquad \text{for} \qquad \eta > \eta_{f_0}
\end{equation}
where $c_n$ and $d_n$ are constants. 

After some algebra, the linearised governing equations \eqref{eq:localmasscons_sim}, \eqref{eq:fullmasscons_sim} in the interface region ($\eta_{b_0}<\eta<\eta_{f_0}$) become
\begin{equation}
(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,
\end{equation}
\begin{equation} \label{eq:h1p1rel}
(1-M)(\sigma +1/2) \mathcal{H}_1 = -\frac{M}{\eta} \frac{\mathrm{d} \mathcal{P}_1}{\mathrm{d} \eta}.
\end{equation}
We can then eliminate $\mathcal{H}_1$ and obtain a single equation for $\mathcal{P}_1$ in the interface region between the contact lines:
\begin{equation} \label{eq:p1govsimp}
\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.
\end{equation}
\begin{figure}
\centering
\includegraphics{dispersionrelationfig24jun.eps}
\caption{The dispersion function, $\mathcal{F}(\sigma,n,M)$ \eqref{eq:dispersionfundef} 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$ \eqref{eq:dispersion1} 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).}
\label{fig:dispersiondifference}
\end{figure}
To determine the pressure perturbation in the interface region, we solve \eqref{eq:p1govsimp} 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 \ref{eq:transBCtrail1}, \ref{eq:transBCtrail2}, \ref{eq:transBCfront1}, \ref{eq:transBCfront2}). At $\eta=\eta_{b_0}^+$, these boundary conditions take the form
\begin{align}
\mathcal{P}_1\rvert_{\eta_{b_0}^+} &= c_n \eta_{b_0}^n, \\
\frac{\mathrm{d} \mathcal{P}_1}{\mathrm{d} \eta} \bigg\rvert_{\eta_{b_0}^+} &= \frac{\eta_{b_1}}{\eta_{b_0}^2} + n c_n \eta_{b_0}^{n-1},\\
(\sigma+1/2) \eta_{b_1} &= - 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} + n c_n \eta_{b_0}^{n-1}\Bigg), \label{eq:kinlintrail}
\end{align}
where we have used the upstream behaviour \eqref{eq:p1ups}.
At $\eta=\eta_{f_0}^-$, the analogous boundary conditions take the form
\begin{align}
\mathcal{P}_1\rvert_{\eta_{f_0}^-} &= d_n \eta_{f_0}^{-n}, \\
\frac{\mathrm{d} \mathcal{P}_1}{\mathrm{d} \eta} \bigg\rvert_{\eta_{f_0}^-}&=\frac{\eta_{f_1}}{2} - n d_n \eta_{f_0}^{-n-1}, \\
(\sigma+1/2) \eta_{f_1} &= -\frac{\mathrm{d} \mathcal{P}_1}{\mathrm{d} \eta} \bigg\rvert_{\eta_{f_0}^-}=-\Bigg(\frac{\eta_{f_1}}{2} - n d_n \eta_{f_0}^{-n-1}\Bigg), \label{eq:kinlinlead}
\end{align}
where we have used the downstream behaviour \eqref{eq:p1down}.
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 \eqref{eq:p1govsimp}. The system is an eigenvalue problem for the growth rate $\sigma$. The general solution to equation  \eqref{eq:p1govsimp} is characterised by the value of $\sigma$ relative to the critical value
\begin{equation}
\sigma_c = \frac{-2n^2 -1}{4n^2 +1}.
\end{equation}
\begin{figure}
\centering
\includegraphics[width=0.95\columnwidth]{eigenfunctions_03March.eps}
\caption{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 \eqref{eq:h1p1rel} and \eqref{eq:eigefuns}. 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$.}
\label{fig:eigen}
\end{figure}
This critical value of $\sigma$ corresponds to a repeated root in the auxiliary equation for the ordinary differential equation \eqref{eq:p1govsimp}. 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
\begin{equation} \label{eq:eigefuns}
\mathcal{P}_1 = \eta^{1/2} \big[ a_1 \cos( \alpha \log \eta) + a_2 \sin( \alpha \log \eta)\big], 
\end{equation}
where
\begin{equation}
\alpha= \sqrt{-\frac{1}{4}-\bigg(\frac{\sigma+1/2}{\sigma +1} \bigg) n^2},
\end{equation}
and $a_1$ and $a_2$ are constants. We apply the six boundary conditions to obtain the following dispersion relation for $\sigma$
\begin{equation} \label{eq:dispersion1}
\mathcal{F}(\sigma,n,M) =0,
\end{equation}
where
\begin{equation} \label{eq:dispersionfundef}
\begin{aligned} 
\mathcal{F}&(\sigma,n,M)= \\
\Big[&\big(\sigma+ \textstyle{\frac{1}{2}}\big)n s_b -(\sigma+1)\big(\textstyle{\frac{1}{2}}s_b+\alpha c_b\big)\Big] \Big[ \big(\sigma+\textstyle{\frac{1}{2}}\big)n c_f+(\sigma+1)\big(\textstyle{\frac{1}{2}}c_f-\alpha s_f\big)\Big]-\\
\Big[&\big(\sigma+\textstyle{\frac{1}{2}} \big)n s_f +(\sigma+1)\big(\textstyle{\frac{1}{2}}s_f+\alpha c_f\big)\Big] \Big[ \big(\sigma+\textstyle{\frac{1}{2}}\big)n c_b-(\sigma+1)\big(\textstyle{\frac{1}{2}}c_b-\alpha s_b\big)\Big], 
\end{aligned}
\end{equation}
and
\begin{align}
s_b=&\sin(\alpha \log \eta_{b_0}), & s_f=&\sin(\alpha \log \eta_{f_0}), \\
c_b=&\cos(\alpha \log \eta_{b_0}), & c_f=&\cos(\alpha \log \eta_{f_0}).
\end{align}
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 \eqref{eq:dispersion1} also depends on $M$.
\par
The dispersion equation \eqref{eq:dispersion1} can be solved to obtain the growth rate, $\sigma$ for any $n \geq 1$ and $0<M<1$. The corresponding eigenfunctions are given by \eqref{eq:eigefuns}. Since $\sigma < \sigma_c<0$, the growth rate is always negative and hence the axisymmetric similarity solutions are linearly stable. 

The dispersion equation \eqref{eq:dispersion1} can have multiple solutions. The function $\mathcal{F}(\sigma,n,M)$ is plotted against $\sigma$ in figure \ref{fig:dispersiondifference} 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 \ref{fig:dispersiondifference}). 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 \eqref{eq:dispersion1} corresponds to the slowest decaying non-trivial solution and these values are indicated by circles in figure \ref{fig:dispersiondifference}. We expect this to be the decay rate that is observed and we take this solution as the correct value of $\sigma$.

The eigenfunctions \eqref{eq:eigefuns} 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 \eqref{eq:eigefuns} can be calculated from \eqref{eq:h1p1rel}. These interfacial perturbation eigenfunctions are shown in figure \ref{fig:eigen} for different values of the viscosity ratio, $M$ and the mode, $n$. The shapes are more oscillatory at higher values of $n$.


\begin{figure}
\centering
\includegraphics{decayratesM0109.eps}
\caption{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 \ref{app:numericsdecay}).}
\label{fig:decayrates}
\end{figure}

The growth rate, $\sigma$, calculated for different values of $M$ and $n$ using the dispersion relation \eqref{eq:dispersion1} is plotted using continuous lines in figure \ref{fig:decayrates}. The rate of decay is slower for larger values of $n$ and smaller viscosity ratios, $M$. 
The crosses in figure \ref{fig:decayrates} 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 \ref{app:numericsdecay}. 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 (\ref{eq:kinlintrail}, \ref{eq:kinlinlead}) furnish $\sigma=-1$. 
The thickness perturbation takes the form (see \ref{eq:h1p1rel})
\begin{equation}
\mathcal{H}_1 = \frac{2 M b}{1-M} \eta^{-1}.
\end{equation}
The numerical results for axisymmetric perturbations to the self-similar solutions corroborated that $\sigma =-1$ (shown as a cross at $n=0$ in figure \ref{fig:dispersiondifference}). 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 \citep{grundy1982eigenvalues,mathunjwa2006self}. For a different stability argument in the case of axisymmetric perturbations, see the Appendix of \citet{nordbotten2006similarity}.



\section{Vertically varying permeability} \label{sec:vertvarperm}
\begin{figure}
\centering
\includegraphics[width=0.95\columnwidth]{varpermG0shapes.eps}
\caption{Self-similar axisymmetric interface shapes in a layer with vertically varying permeability. (a) Small-$G$ analytic solutions \eqref{eq:G0intrusionvark} for different linear permeability variations \eqref{eq:linearperm}. (b) Small-$G$ analytic solutions in the case that there is a shock-like region.}
\label{fig:varpermG0shapes}
\end{figure}

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 \S \ref{subsec:shock}). 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
\begin{equation}
\mathcal{K}(Z/H_0) = \frac{k(Z)}{\frac{1}{H_0} \int_0^{H_0} k(Z)\,\mathrm{d}Z},
\end{equation}
where the denominator is the mean permeability. We also define the dimensionless depth-integrated permeability \citep{hinton2018buoyancy},
\begin{equation}
\Phi(\mathcal{H}) = \int_0^\mathcal{H} \mathcal{K}(s) \, \mathrm{d} s.
\end{equation}
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 (\ref{eq:localmasscons_sim}, \ref{eq:fullmasscons_sim}) generalise to
\begin{equation} \label{eq:localmasscons_simvar}
\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,
\end{equation}
\begin{equation} \label{eq:fullmasscons_simvar}
\tilde{\nabla} \cdot \Big[\big(M+(1-M)\Phi(\mathcal{H})\big) \tilde{\nabla} \mathcal{P} - G M (1-\Phi(\mathcal{H})) \tilde{\nabla} \mathcal{H}   \Big] = 0.
\end{equation}
This model applies provided that the fluids have become segregated by buoyancy. In a uniform layer this applies at times given by \eqref{eq:segregatetime} 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
\begin{equation} \label{eq:segtimevark}
T \gg \frac{\mu_a}{\Delta \rho g} \int_0^{H_0} \frac{1}{k_v(Z)} \mathrm{d} Z,
\end{equation}
where $k_v(Z)$ is the vertical permeability. In an isotropic, vertically heterogeneous layer, $k_v(Z)=k(Z)$. The time \eqref{eq:segtimevark} 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 \eqref{eq:shallowtime}.

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 \eqref{eq:localmasscons_simvar}, \eqref{eq:fullmasscons_simvar}, with the pressure gradient at the top boundary given by
\begin{equation}
\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},
\end{equation}
and $\mathcal{H}_0$ satisfies the following ordinary differential equation
\begin{equation} \label{eq:mch0odevark}
-\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{G M \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).
\end{equation}
The interface gradient and pressure gradient at the leading contact line are
\begin{equation} \label{eq:gradleadvark}
\frac{\mathrm{d} \mathcal{H}_0}{\partial \eta} = \frac{1}{\eta_{f_0} M G} - \frac{\eta_f}{2 G \mathcal{K}(0)}, \qquad \frac{\mathrm{d} \mathcal{P}_0}{\partial \eta} = \frac{-\eta_{f_0}}{2 \mathcal{K}(0)},
\end{equation}
respectively. Whilst at the trailing contact line
\begin{equation} \label{eq:gradtrailvark}
\frac{\mathrm{d} \mathcal{H}_0}{\mathrm{d} \eta} = \frac{\eta_{b_0}}{2 M G \mathcal{K}(1)} - \frac{1}{\eta_{b_0} G}, \qquad \frac{\mathrm{d} \mathcal{P}_0}{\mathrm{d} \eta} = \frac{-1}{\eta_{b_0}}.
\end{equation}
The small-$G$ analytic solutions are given by neglecting the right-hand side of \eqref{eq:mch0odevark} from which we obtain the following implicit solution
\begin{align}
\mathcal{H}_0&=1,  &   0<&\eta< \eta_{b_0},\\
\eta&=\frac{\sqrt{2M \mathcal{K}(\mathcal{H}_0)}}{M+(1-M)\Phi(\mathcal{H}_0)}, &  0<&\mathcal{H}_0< 1, \label{eq:G0intrusionvark}  \\
\mathcal{H}_0&=0, &  \eta_{f_0}<&\eta,
\end{align}
where the contact lines are
\begin{equation}
\eta_{f_0} = \sqrt{\frac{2  \mathcal{K}(0)}{M}}, \qquad \eta_{b_0} =\sqrt{2 M  \mathcal{K}(1)}.
\end{equation}
The associated pressure gradient in the interface region is given by
\begin{equation} \label{eq:dpdetavark1}
\frac{\mathrm{d} \mathcal{P}_0}{\mathrm{d} \eta} = -\frac{1}{\sqrt{2M \mathcal{K}(\mathcal{H}_0)}}.
\end{equation}
In a uniform layer, the interface \eqref{eq:G0intrusionvark} 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 \citep{hinton2018buoyancy}. The solution \eqref{eq:G0intrusionvark} 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$,
\begin{equation} \label{eq:linearperm}
\mathcal{K}(s) = 1+ \Delta k (s-1/2),
\end{equation}
the interface \eqref{eq:G0intrusionvark} is monotonic if and only if \citep{hinton2018buoyancy}
\begin{equation}
M < \frac{(2-\Delta k)^2}{4 - 2 \Delta k + \Delta k^2}.
\end{equation}
Provided \eqref{eq:G0intrusionvark} 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 \ref{fig:varpermG0shapes}a). 
In the case that \eqref{eq:G0intrusionvark} has a turning point, a shock must be introduced; this regime is studied in \S \ref{subsec:shock}.

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 \S \ref{sec:stabilmech} to show how the pressure gradients suppress instabilities in a layer with vertically varying permeability (but no shock-like regions; for that case, see \S \ref{subsec:shock}). 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 \S \ref{sec:stabilmech}). 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 \eqref{eq:G0intrusionvark} is monotonic, the contacts lines for $G > 0$ satisfy
\begin{equation} \label{eq:ineqconvark}
\eta_{f_0} \geq \sqrt{\frac{2  \mathcal{K}(0)}{M}}, \qquad \eta_{b_0} \leq \sqrt{2 M  \mathcal{K}(1)},
\end{equation}
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 \eqref{eq:gradleadvark} and \eqref{eq:gradtrailvark}. Thus, the pressure gradient in the ambient fluid just ahead of the leading contact line is given by
\begin{equation}
\frac{\mathrm{d}(\mathcal{P}_0-G \mathcal{H}_0)}{\mathrm{d} \eta} = - \frac{1}{M \eta_{f_0}},
\end{equation}
whilst in the input fluid just upstream of the leading contact line, the pressure gradient is
\begin{equation}
\frac{\mathrm{d} \mathcal{P}_0}{\mathrm{d} \eta} = - \frac{\eta_{f_0}}{2 \mathcal{K}(0)},
\end{equation}
and the inequality \eqref{eq:ineqconvark}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
\begin{equation}
\frac{\mathrm{d}(\mathcal{P}_0-G \mathcal{H}_0)}{\mathrm{d} \eta} = - \frac{\eta_{b_0}}{2 M \mathcal{K}(1)},
\end{equation}
whilst in the input fluid just upstream of the trailing contact line, the pressure gradient is
\begin{equation}
\frac{\mathrm{d} \mathcal{P}_0}{\mathrm{d} \eta} = - \frac{1}{\eta_{b_0}},
\end{equation}
and the inequality \eqref{eq:ineqconvark}b ensures that the magnitude of the pressure gradient is larger in the input fluid, which stabilises the interface.


\subsection{Shock-like regions} \label{subsec:shock}
In the case in which the small-$G$ interface \eqref{eq:G0intrusionvark} 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 \ref{fig:varpermG0shapes}b).  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 \ref{fig:varpermG0shapes}b). The shock at the top of the layer is at location $\eta=\eta_s$ and is of thickness $\mathcal{H}_s$, which satisfy
\begin{equation}
\eta_s=\sqrt{2F'(\mathcal{H}_s)}, \qquad \eta_s^2 \mathcal{H}_s = 2 F(\mathcal{H}_s), \quad \text{where} \quad F(\mathcal{H}) = \frac{\Phi(\mathcal{H})}{M+(1-M)\Phi(\mathcal{H})}.
\end{equation}
These arise from continuity with the solution \eqref{eq:G0intrusionvark}, 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 \ref{fig:varpermG0shapes}b). For more complicated permeability profiles, there can be multiple shock-like regions separated by rarefaction-like regions \citep{hinton2018buoyancy}.

\begin{figure}
\centering
\includegraphics[width=0.95\columnwidth]{varpermfig24jun.eps}
\caption{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 \S \ref{sec:numerics}). 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).}
\label{fig:varperm}
\end{figure}

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 \ref{fig:varperm} 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 \ref{fig:varperm}c; 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 \citep{hinton2018buoyancy}. 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
\begin{equation}
\frac{L}{H_0} \sim G \Bigg(\frac{QT}{2 \pi \phi H_0^3} \Bigg)^{1/2},
\end{equation}
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 \ref{fig:varpermG0shapes}b). 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
\begin{equation}
\eta_{f_0}\geq \eta_s \geq \sqrt{\frac{2 \mathcal{K}(0)}{M}}.
\end{equation}
The first inequality arises from buoyancy smoothing the interface beyond the shock location (see figure \ref{fig:varperm}c) whilst the second arises because the profile \eqref{eq:G0intrusionvark} 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 \eqref{eq:ineqconvark}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 \eqref{eq:ineqconvark}b remains correct. 

The gradient of the interface and the pressure at the upper boundary at the contact lines are given by \eqref{eq:gradleadvark} and \eqref{eq:gradtrailvark}, even when there are shock-like regions. Hence the stability argument of \S \ref{sec:vertvarperm} 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 \citep[see also][]{debbabi2018impact}. We note that there may be significant intermingling of the fluids in the transient evolution to the buoyancy-segregated flow \citep[e.g.][]{huppert2013competition}. There may also be a competition between viscous fingering and heterogeneity-driven fingering prior to buoyancy segregation \citep{de1997viscous}.




\section{Discussion and conclusion} \label{sec:conc}
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 \citep{schoonman2017radial}. The base self-similar flow was found for a line source and a point source by \citet{zheng2015viscous} and \citet{hinton2020axisymmetric} 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 \citep{snyder1998flow,john2013variable}. 

In the context of CO\textsubscript{2} 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. 


\section*{Acknowledgements}
E.M.H. is grateful to the University of Melbourne for the award of a Harcourt-Doig research fellowship.

\section*{Declaration of interests}
The authors report no conflict of interest.


\section{How to submit to the \textbf{\textit{Journal of Fluid Mechanics}}}
Authors must submit using the online submission and peer review system, Scholar One Manuscripts (formerly Manuscript Central) at http://mc.manuscriptcentral.com/jfm. If visiting the site for the first time, users must create a new account by clicking on `register here'. Once logged in, authors should click on the `Corresponding Author Centre', from which point a new manuscript can be submitted, with step-by-step instructions provided. Authors must at this stage specify whether the submission is a standard paper, or a {\it JFM Rapids} paper (see \S\ref{sec:filetypes} for more details). In addition, authors must specify an editor to whom the paper will be submitted from the drop-down list provided. Note that all editors exclusively deal with either standard papers or {\it JFM Rapids} (clearly indicated on the list), so please ensure that you choose an editor accordingly. Corresponding authors must provide a valid ORCID ID in order to submit a manuscript, either by linking an existing ORCID profile to your ScholarOne account or by creating a new ORCID profile. Once your submission is completed you will receive an email confirmation.  Book reviews should not be submitted via the online submission site and should instead be submitted by email to c.p.caulfield@damtp.cam.ac.uk.
 
\section{Rules of submission}\label{sec:rules_submission}
Submission of a paper implies a declaration by the author that the work has not previously been published, that it is not being considered for publication elsewhere and that it has not already been considered by a different editor of the Journal. Note that a report on a conference must be submitted within 3 months of the meeting.

\section{Types of paper}\label{sec:types_paper}
\subsection{Standard papers}
Regular submissions to JFM are termed `standard papers'. Note that such papers must contain original research, as review articles are not published in JFM. Papers should be written in a concise manner; though JFM has no page limit, each paper will be judged on its own merits, and those deemed excessive in length will be rejected or will require significant revision.  
\subsection{JFM Rapids}
{\it JFM Rapids} is devoted to the rapid publication of short, high-impact papers across the full range of fluid mechanics. Manuscripts submitted as {\it Rapids}  must be strictly 10 or fewer printed pages, and must be submitted in {\LaTeX} using the jfm.cls class file, so as to ensure that they meet the page limit and to expedite their production.  The principal and over-riding objective is to publish papers of the highest scientific quality. 

Once a paper is submitted, referees are asked to provide reports with a short turnaround.  In order to be accepted for publication in {\it JFM Rapids}, such papers must require only minor revisions to improve clarity, usually without recourse to a second round of refereeing. In this case, and at the discretion of the editor, some additional pages may be allowed to address specific points raised by the referees, such as the addition of an extra figure or some explanatory text.  

In cases where the editor, guided by the referees, judges that a paper has merit but requires substantial revision that will require significant refereeing, a decision of `revise and resubmit' will be given. On re-submission, such papers will be handled as standard JFM papers and if accepted will not subsequently appear as {\it JFM Rapids}.

{\it Rapids} will be published online within one month of final acceptance.  They will appear within a designated section on the Journal of Fluid Mechanics website.  Each {\it Rapid} will be cited and indexed as a JFM article but with a distinctive {\it Rapids} identifier, and will be assigned to a JFM volume. {\it Rapids} will not be included in the regular print versions of JFM, and accordingly they incur no colour figure charges.
 
\section{File types}\label{sec:filetypes}
Authors are strongly encouraged to compose their papers in {\LaTeX}, using the jfm.cls style file and supporting files provided at\\ https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/information/\\ instructions-contributors, with the jfm-instructions.tex file serving as a template (note that this is mandatory for {\it JFM Rapids}). A PDF of the {\LaTeX} file should then be generated and submitted via the submission site. Please note that PDFs larger than 10MB are not acceptable for review. There is no need to submit the {\LaTeX} source file alongside the PDF, but upon provisional acceptance of the paper, the {\LaTeX} source file, along with individual figure files and a PDF of the final version, will need to be submitted for typesetting purposes. 
Authors may also compose standard papers in Word, though this will lead to the paper spending a longer period in production. If using Word, please note that equations must NOT be converted to picture format and the file must be saved with the option `make equation editable'. 
All submitted video abstract files should be formatted as MP4 (H.264). MP4 has full compatibility across commonly used browsers, whereas other video formats will only work on selected browsers. This will ensure the greatest possible dissemination of this work. 
\section{Preparing your manuscript}
Authors should write their papers clearly and concisely in English, adhering to JFM's established style for notation, as provided in \S\ref{notstyle}. We encourage the submission of online supplementary material alongside the manuscript where appropriate (see section \ref{online}). Metric units should be used throughout and all abbreviations must be defined at first use, even those deemed to be well known to the readership. British spelling must be used, and should follow the \textit{Shorter Oxford English Dictionary}.


\subsection{Figures}
All authors need to acquire the correct permissions and licences to reproduce figures. Figures should be as small as possible while displaying clearly all the information required, and with all lettering readable. Every effort should be taken to avoid figures that run over more than one page. Figures submitted in colour will appear online in colour but, with the exception of {\it JFM Rapids}, all figures will be printed in black and white unless authors specify during submission that figures should be printed in colour, for which there is a charge of \pounds200 plus VAT per figure (i.e. \pounds240) with a cap of \pounds1000 plus VAT per article (i.e. \pounds1200) (colour is free for {\it JFM Rapids}).  Note that separate figures for online and print will {\bf not} be accepted and it is the author's responsibility to ensure that if a figure is to appear in colour online only, that same figure will still be meaningful when printed in black and white (for example, do not rely upon colours to distinguish lines if those colours will just appear as similar shades of grey when printed). If using colour, authors should also ensure that some consistency is applied within the manuscript. For review purposes figures should be embedded within the manuscript. Upon final acceptance, however, individual figure files will be required for production. These should be submitted in EPS or high-resolution TIFF format (1200 dpi for lines, 300 dpi for halftone and colour in CMYK format, and 600 dpi for a mixture of lines and halftone). The minimum acceptable width of any line is 0.5pt. Each figure should be accompanied by a single caption, to appear beneath, and must be cited in the text. Figures should appear in the order in which they are first mentioned in the text and figure files must be named accordingly (`Figure 1.eps', `Figure 2a.tiff', etc) to assist the production process (and numbering of figures should continue through any appendices). The word \textit {figure} is only capitalized at the start of a sentence. For example see figures \ref{fig:ka} and \ref{fig:kd}. Failure to follow figure guidelines may result in a request for resupply and a subsequent delay in the production process. Note that {\em all} figures are relabelled by the typesetter, therefore the files supplied should have editable labelling. Please ensure all figure labels are carefully checked against your originals when you receive your proofs. 

\begin{figure}
  \centerline{\includegraphics{trapped}
  \caption{Trapped-mode wavenumbers, $kd$, plotted against $a/d$ for
    three ellipses:\protect\\
    ---$\!$---,
    $b/a=1$; $\cdots$\,$\cdots$, $b/a=1.5$.}
\label{fig:ka}
\end{figure}

\begin{figure}
  \centerline{\includegraphics{modes}}
  \caption{The features of the four possible modes corresponding to
  (\textit{a}) periodic\protect\\ and (\textit{b}) half-periodic solutions.}
\label{fig:kd}
\end{figure}

\subsection{Tables}
Tables, however small, must be numbered sequentially in the order in which they are mentioned in the text. The word \textit {table} is only capitalized at the start of a sentence. See table \ref{tab:kd} for an example.

\begin{table}
  \begin{center}
\def~{\hphantom{0}}
  \begin{tabular}{lccc}
      $a/d$  & $M=4$   &   $M=8$ & Callan \etal \\[3pt]
       0.1   & 1.56905 & ~~1.56~ & 1.56904\\
       0.3   & 1.50484 & ~~1.504 & 1.50484\\
       0.55  & 1.39128 & ~~1.391 & 1.39131\\
       0.7   & 1.32281 & ~10.322 & 1.32288\\
       0.913 & 1.34479 & 100.351 & 1.35185\\
  \end{tabular}
  \caption{Values of $kd$ at which trapped modes occur when $\rho(\theta)=a$}
  \label{tab:kd}
  \end{center}
\end{table}

\subsection{Online supplementary material}\label{online}
Relevant material which is not suitable for print production, such as movies or numerical simulations/animations, can be uploaded as part of the initial submission. Movies should be designated as `Movie' and each individual file must be accompanied by a separate caption and a suitable title (eg Movie 1). Accepted formats are .mov, .mpg, .mp4, and .avi, though they should be archived as a .zip or .tar file before uploading. Each movie should be no more than 10MB. Upon publication these materials will then be hosted online alongside the final published article. Likewise, should there be detailed mathematical relations, tables or figures which are likely to be useful only to a few specialists or take up excessive space in the printed journal, these can also be published online as supplementary material [designated as `Other supplementary material']. Note that supplementary material is published `as is', with no further production performed.

\section{Editorial decisions}
\subsection{Revision}
If a revision is requested, you should upload revised files following the same procedure as for submitting a new paper. You begin by clicking on `Manuscripts with decision' in your Corresponding Author Center, and then on `Create a revision'. (Note that if you abandon the process before completing the submission, to continue the submission, you must click on `Revised manuscripts in draft'.) There is a new first page showing the decision letter and a space for your reply to the referee's/editor's comments. You also have the opportunity at this stage to upload your reply to the comments as a separate file. All the values filled in on original submission are displayed again. The ID number of the paper will be appended `.R1'. Also note that if a manuscript is submitted as a {\it JFM Rapid}, but requires substantial revision, it will be re-designated as a standard paper, and the ID and paper type will be amended to reflect this.

\subsection{Provisional acceptance}
If the paper is accepted as suitable for publication you will be sent a provisional acceptance decision. This enables you to upload the final files required for production:
(1) the final PDF or word version of the paper, designated as a `main document';
(2) any source files (see section \ref{sec:filetypes}) which must be designated as `production (not for review)' and uploaded as a single .zip or .tar file.

\subsection{Acceptance}
On receipt of the production files you will be sent an email indicating completion of the acceptance process.

\section{Publication process}
Once a paper has been accepted for publication and the source files have been uploaded, the manuscript will be sent to Cambridge University Press for copyediting and typesetting, and will be assigned a digital object identifier (doi). When the proof is ready, authors will receive an email alert containing a link to the PDF of the proof, and instructions for its correction and return. It is imperative that authors check their proofs closely, particularly the equations and figures, which should be checked against the accepted file, as the production schedule does not allow for corrections at a later stage. Once ready for printing, papers will be published online on Cambridge Journals Online in the current `open' volume. {\it JFM Rapids} will also be assigned to the current open volume, but will receive an article number in place of traditional pagination, as they will not appear in the print version of that volume. Each volume will be kept open for approximately two weeks, at which point it will be considered `closed', and then printed and distributed. The following volume will be immediately opened. Note that the PDF published online is the Version of Record, matching the print version, and no further alterations/corrections to this document will be allowed. The corresponding author is emailed a link to the published article when it is first published online.

\section{Corrigenda}
The Journal will publish corrigenda that alter significant conclusions made in a paper.  Such corrigenda should be submitted to an associate editor, who will consider the submission similarly to a new paper and may consult referees if appropriate.  When published, corrigenda are clearly linked with the original articles to which they refer, and the articles to them.

The Journal does not normally publish corrigenda to amend typographical errors, so it is extremely important that authors make very careful checks of their manuscript at every stage, including the reading of proofs, prior to publication.



\section{Obtaining help}
Technical support for the online submission system is available by clicking on the `Get Help Now' link at the top-right corner of each page of the submission site. Any other questions relating to the submission or publication process should be directed to the JFM Editorial Assistant, Mrs Amanda Johns, at ajohns@cambridge.org.

\section{Cambridge Language Editing Service} 
We suggest that authors whose first language is not English have their manuscripts checked by a native English speaker before submission. This is optional but will help to ensure that any submissions that reach peer review can be judged exclusively on academic merit. We offer a Cambridge service which you can find out more about at https://www.cambridge.org/core/services/authors/language-services, and suggest that authors make contact as appropriate. Please note that use of language editing services is voluntary and at the author’s own expense. Use of these services does not guarantee that the manuscript will be accepted for publication nor does it restrict the author to submitting to a Cambridge-published journal


\section{Notation and style}\label{notstyle}
Generally any queries concerning notation and journal style can be answered by viewing recent pages in the Journal. However, the following guide provides the key points to note. It is expected that Journal style will be followed, and authors should take care to define all variables or entities upon first use. Also note that footnotes are not normally accepted.

\subsection{Mathematical notation}
\subsubsection{Setting variables, functions, vectors, matrices etc}
{\bf Italic font} should be used for denoting variables, with multiple-letter symbols avoided except in the case of dimensionless numbers such as $\Rey$, $\Pran$ and $\Pen$ (Reynolds, Prandtl, and P\'eclet numbers respectively, which are defined as \verb}\Rey}, \verb}\Pran} and {\verb}\Pen} in the template).\\
{\bf Upright Roman font} (or upright Greek where appropriate) should be used for:\\
Operators: sin, log, d, $\Delta$, e etc.\\
Constants: i ($\sqrt{-1}$), $\upi$ (defined as \verb}\upi}), etc.\\
Functions: $\Ai$, $\Bi$ (Airy functions, defined as \verb}\Ai} and \verb}\Bi}), $\Real$ (real part, defined as \verb}\Real}), $\Imag$ (imaginary part, defined as \verb}\Imag}), etc.\\
Physical units: cm, s, etc\\
Abbreviations: c.c. (complex conjugate), h.o.t. (higher-order terms), DNS, etc.\\
{\bf Bold italic font} (or bold sloping Greek) should be used for:\\
Vectors (with the centred dot for a scalar product also in bold): $\boldsymbol{i \cdot j}$\\
{\bf Bold sloping sans serif font}, defined by the \verb}\mathsfbi} macro, should be used for:\\
Tensors and matrices: $\mathsfbi{D}$ \\
{\bf Script font} (for example $\mathcal{G}$, $\mathcal{R}$) can be used as an alternative to italic when the same letter denotes a different quantity (use \verb}\mathcal} in \LaTeX)\\
The product symbol ($\times$) should only be used to denote multiplication where an equation is broken over more than one line, to denote a cross product, or between numbers (the $\cdot$ symbol should not be used, except to denote a scalar product specifically).\\ 

\subsubsection{Other symbols}
A centred point should be used only for the scalar product of vectors.
Large numbers that are not scientific powers should not include commas, but have the
form 1600 or 16 000 or 160 000.
Use \textit{O} to denote `of the order of', not the \LaTeX\ $\mathcal{O}$.

\section{Citations and references}
All papers included in the References section must be cited in the article, and vice versa. Citations should be included as, for example ``It has been shown \citep{Rogallo81} that...'' (using the {\verb}\citep} command, part of the natbib package) ``recent work by \citet{Dennis85}...'' (using {\verb}\citet}).
The natbib package can be used to generate citation variations, as shown below.\\
\verb#\citet[pp. 2-4]{Hwang70}#:\\
\citet[pp. 2-4]{Hwang70} \\
\verb#\citep[p. 6]{Worster92}#:\\
\citep[p. 6]{Worster92}\\
\verb#\citep[see][]{Koch83, Lee71, Linton92}#:\\
\citep[see][]{Koch83, Lee71, Linton92}\\
\verb#\citep[see][p. 18]{Martin80}#:\\
\citep[see][p. 18]{Martin80}\\
\verb#\citep{Brownell04,Brownell07,Ursell50,Wijngaarden68,Miller91}#:\\
\citep{Brownell04,Brownell07,Ursell50,Wijngaarden68,Miller91}\\
The References section can either be built from individual \verb#\bibitem# commands, or can be built using BibTex. The BibTex files used to generate the references in this document can be found in the zip file at http://journals.cambridge.org/\linebreak[3]data/\linebreak[3]relatedlink/\linebreak[3]jfm-ifc.zip.\\
Where there are up to ten authors, all authors' names should be given in the reference list. Where there are more than ten authors, only the first name should appear, followed by et al.\\

Acknowledgements should be included at the end of the paper, before the References section or any appendicies, and should be a separate paragraph without a heading. Several anonymous individuals are thanked for contributions to these instructions.

