Research

applied mathematics, fluid mechanics, geoscience, hydrology

Graham P. Benham

# Anisotropy distorts the spreading of a fixed volume porous gravity current

Graham P. Benham1 1Mathematical Institute, University of Oxford, Oxford OX2
6GG, UK<EMAIL_ADDRESS>
###### Abstract

We consider the release and subsequent gravity-driven spreading of a finite
volume of fluid in an anisotropic porous medium bounded by an impermeable
substrate. When the permeability in the vertical direction is much smaller
than the horizontal direction, as is the case in many real geological
reservoirs, this restricts the spread of the current to a very thin layer near
the impermeable base. Using a combination of asymptotic analysis and finite
difference computations of Darcy flow, we show that there exist two distinct
flow regimes. At early times the bulk of the current descends slowly and
uniformly, injecting fluid into thin finger-like regions near the base. At
much later times the current transitions to the classical gravity-driven
solution and continues to spread with a self-similar shape. One interesting
consequence is that the swept volume of the current grows differently
depending on the anisotropy of the medium. This has important consequences for
managing contaminant spills, where it is important to minimise the contacted
volume of the aquifer, or during geological CO2 sequestration where a larger
contacted volume results in more CO2 being stored.

###### keywords:

gravity currents, porous media, anisotropy

##  1 Introduction

Gravity-driven flows resulting from the release of a fluid within a porous
medium are a common feature of environmental fluid dynamics. For example, such
flows arise when groundwater responds to heavy rainfall [1], after the
spillage of a contaminant [2], or during the geological storage of carbon
dioxide in saline aquifers [3]. Since all geological aquifers are
heterogeneous and often this heterogeneity manifests as an anisotropic
permeability field [4], it is important to quantify how this affects the
migration speeds and shape of the current. This is particularly relevant to
situations where it is desirable to minimise the volume of the aquifer
contacted by the fluid (e.g. containing the spread of a contaminant) or to
maximise the contacted volume (e.g. trapping residual saturation during CO2
storage).

In the case of homogeneous and isotropic porous media, there have been
numerous studies on the evolution of fixed volume gravity currents. Some of
these studies have treated single-phase flows [5], whilst others have
incorporated multiphase effects such as residual trapping [6, 7] and
dissolution within the ambient fluid [8]. In some cases simple scaling laws
were derived by exploiting the self-similar properties of gravity currents, as
shown by [9] in the case of single-phase flow, and by [10] in the case of
trapped residual saturation. Later work by [11, 12, 13] explored the
transition to self-similarity in both confined and unconfined settings.
However, less attention has been paid to the case of heterogeneous or
anisotropic porous media, despite the relevance to real geological reservoirs.

Nevertheless, some progress has been made for specific types of buoyancy
driven flows in heterogeneous media. For example, several studies have
investigated how anisotropy affects convective dissolution within an ambient
fluid phase [14, 15, 16]. In the case of gravity currents resulting from
constant injection, [17] explored how heterogeneities of different lateral and
vertical scales affect the migration speed of a CO2 plume. Likewise, [18]
addressed the case of a gravity current resulting from point source injection
in an anisotropic medium. In this study, it was shown that anisotropy can
cause a build-up of pressure that stretches the flow into an ellipsoid shape
during an early-time regime of the flow, before transitioning to a gravity-
dominated regime at much later times. However, no studies have addressed how
heterogeneity affects the spreading of a released volume of fluid (i.e. in the
absence of injection), despite the relevance to post-injection scenarios
during CO2 storage, and to post-leakage scenarios in the context of
contaminant spills.

CO2 storage in geological reservoirs is one of the key proposed technologies
to reduce emissions and limit the effects of global warming [3]. In such
scenarios, buoyant CO2 is injected into a brine-filled reservoir beneath an
impermeable cap rock. Once the injection is switched off, the CO2 rises and
spreads out beneath the cap rock, with a fraction of its mass being lost to
residual trapping (via the drainage/imbibition cycle) and dissolution within
the surrounding brine [8, 19]. Hence, to quantify the trapping potential of
different geological reservoirs (e.g. when choosing potential storage sites)
it is important to understand how the anisotropy of the aquifer may affect the
historical migration of the current across the pore space. In the context of a
contaminant spillage, the objective is to contain and minimise the spread of a
harmful fluid within an aquifer. Therefore, in a similar manner to the CO2
storage problem, it is necessary to quantify how and where the contaminant
fluid will spread in response to the heterogeneity of the aquifer, once the
leak has been closed off.

In this study we demonstrate that anisotropy restricts the flow of the gravity
current to thin finger-like regions spreading near the impermeable boundary,
qualitatively similar to those predicted by other studies [17]. Due to this
flow distortion, the swept volume of the current is reduced for anisotropic
aquifers. This indicates that isotropic aquifers may have better potential for
certain forms of CO2 trapping that depend on the contacted volume of pore
space. By contrast, in the case of a contaminant spill, anisotropic aquifers
may help contain the spread of the fluid by restricting the flow to a reduced
fraction of the porous medium.

##  2 Finite release in two-dimensional anisotropic media

###  2.1 Release and subsequent dynamics

$x$$z$$\bm{H}$

Symmetry

Impermeable$\bm{L}$$x$$z$$\bm{\rho,\mu}$$w=0$

$u=0$

$\bm{w=-\epsilon
u_{b}}$$\bm{\rho_{a}<\rho}$$\bm{k_{H}}$$\bm{k_{V}}$(a)(b)$\bm{z=h(x,t)}$III$\bm{\mathcal{O}(\epsilon^{1/2})}$
Figure 1: Schematic diagram. (a) A finite volume $2LH\phi$ of dense fluid is
released into a porous medium with an anisotropic permeability field $k_{H}\gg
k_{V}$ above an impermeable base. (b) The medium is initially saturated with a
less dense fluid so that the bulk region (I) descends slowly at speed
$-\epsilon u_{b}$, where $\epsilon=k_{V}/k_{H}$ and $u_{b}$ is the buoyancy
velocity (10), whilst a thin finger region (II) develops near the base of the
current.

We consider an anisotropic porous medium in which the horizontal permeability
$k_{H}$ is much larger than the vertical permeability $k_{V}$. Such
anisotropic flow properties are a common feature in geological reservoirs and
may result from the deposition of successive layers of fine and coarse
material or from post-depositional compaction of the formation [20]. As we
will show, the anisotropy of the medium restricts the vertical flow of the
bulk of the fluid, which we denote region I (see figure 1), resulting in a
slow migration towards the impermeable boundary. Meanwhile, gravity-driven
spreading is limited to thin finger-like regions near the base, which we
denote region II .

To start with we restrict our attention to two-dimensional flows (although
radially symmetric flows will be addressed later in Section 5) and we consider
the release of a volume (per unit width) of fluid with constant density
$\rho$. The surrounding porous medium is initially saturated with an ambient
fluid with relatively smaller density $\rho_{a}<\rho$. Due to the Boussinesq
approximation [21], these results also apply in the case of a lighter fluid
(e.g. CO2) released within a porous medium saturated with a heavier fluid
(e.g. brine), with the impermeable boundary located above rather than below.

For the sake of simplicity we take the viscosity of the released fluid and the
ambient fluid to be the same ($\mu$), and we consider that the initial shape
of the released fluid is rectangular, with dimensions $2L\times H$. It should
be noted, however, that these results would apply to any similar convex shape,
as shown in Figure 7 and discussed in more detail in Appendix A.

The flow in the released volume of fluid is subject to the two-dimensional
Darcy equations,

$\displaystyle\mathbf{u}$
$\displaystyle=-\frac{1}{\mu}\underline{\underline{\mathbf{k}}}\cdot\nabla\left[p+\rho
gz\right],$ (1) $\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=0,$ (2)

where $\mathbf{u}$ is the Darcy velocity vector, $p$ is the pressure, and
$\underline{\underline{\mathbf{k}}}=\mathrm{diag}(k_{H},k_{V})$ is the
anisotropic permeability field in the $x,z$ directions. Combining (1) and (2),
the pressure satisfies Laplace’s equation with anisotropic coefficients,

$\frac{\partial^{2}p}{\partial x^{2}}+\epsilon\frac{\partial^{2}p}{\partial
z^{2}}=0,$ (3)

where the anisotropy is given by

$\epsilon=\frac{k_{V}}{k_{H}}\ll 1.$ (4)

The boundary conditions for the flow in region I are as follows. The left hand
and bottom boundaries, $x=0$, $z=0$, are assumed to be symmetric and
impermeable, respectively, so we prescribe no normal flow,

$\displaystyle u=0:$ $\displaystyle\quad x=0,$ (5) $\displaystyle w=0:$
$\displaystyle\quad z=0.$ (6)

Likewise, at the fluid interface, which we denote $z=h(x,t)$, we impose the
dynamic and kinematic boundary conditions:

$\displaystyle p=p_{a}-\rho_{a}gh:$ $\displaystyle\quad z=h(x,t),$ (7)
$\displaystyle w=\phi\frac{\partial h}{\partial t}+u\frac{\partial h}{\partial
x}:$ $\displaystyle\quad z=h(x,t).$ (8)

The former condition matches the pressure in the fluid with the ambient
hydrostatic pressure (note that the reference pressure $p_{a}$ is the ambient
value at $z=0$), whilst the latter condition imposes that a particle at the
interface remains at the interface [22, 23].

In the limit $\epsilon\rightarrow 0$, mass conservation (2) indicates that if
there is no vertical flow (i.e. $k_{V}=0$) then there cannot be any horizontal
flow either. Essentially, the governing equation (3) implies that the pressure
(which is continuous) is set by the ambient fluid, such that
$p=p_{a}-\rho_{a}gz$, and the resulting velocities are $u=w=0$. Hence, the
interface remains at the initial position $z=h_{0}(x)$ for all time.

Now let’s consider the case of a strongly anisotropic porous medium, such that
$0<\epsilon\ll 1$. In this case, the solution can be found by performing an
asymptotic expansion of the pressure in powers of $\epsilon$. Within this
expansion, the leading order contribution to the pressure is simply the
solution to the isotropic problem (i.e. with $\epsilon=0$ exactly). Hence, for
the same reasons as described above, the leading order pressure is hydrostatic
and set by the ambient fluid, $p=p_{a}-\rho_{a}gz$. However, by inserting this
pressure into Darcy’s law (1), we now derive a small but finite (i.e. first
order) vertical velocity within region I, such that

$w=-\epsilon u_{b},$ (9)

where

$u_{b}=\frac{k_{H}\Delta\rho g}{\mu},$ (10)

is the buoyancy velocity and $\Delta\rho=\rho-\rho_{a}$. Since (9) does not
satisfy the impermeability condition (6), it is necessary to reevaluate the
solution near $z\approx 0$ using boundary layer theory.

Mathematically speaking, (3) is a singular perturbation problem since it is a
second order Partial Differential Equation (PDE) with a small parameter in
front of a second derivative. This indicates that not all vertical boundary
conditions can be satisfied by the leading order solution. Specifically, the
dynamic boundary condition (7) is imposed at $z=h$ to ensure continuity of
pressure, whilst the impermeability condition (6) at $z=0$ is left
unsatisfied. To correct this requires rescaling the solution to investigate
changes over a small vertical distance near $z\approx 0$, also known as a
boundary layer. By inspection of (3) it is clear that $z$ must be rescaled by
a factor of $\epsilon^{1/2}$ to recover all terms in the governing equation.
Hence, an appropriate choice of rescaled dimensionless variables is

$x=L\xi,\quad z=\epsilon^{1/2}L\zeta,\quad
p=p_{a}-\rho_{a}gz+\epsilon^{1/2}\Delta\rho gL\,P(\xi,\zeta),$ (11)

where $\xi$, $\zeta$ and $P$ are variables which are $\mathcal{O}(1)$ in
magnitude. Note that the pressure in (11) must also be rescaled by a factor
$\epsilon^{1/2}$ so that the boundary condition (6) (i.e. $\partial p/\partial
z+\rho g=0$ at $z=0$) balances all terms at leading order. In this way, within
the boundary layer region the governing equations and boundary conditions
(3),(5) and (6) become

$\displaystyle\frac{\partial^{2}P}{\partial\xi^{2}}+\frac{\partial^{2}P}{\partial\zeta^{2}}$
$\displaystyle=0,$ (12) $\displaystyle\frac{\partial P}{\partial\xi}$
$\displaystyle=0:\quad\xi=0,$ (13) $\displaystyle\frac{\partial
P}{\partial\zeta}$ $\displaystyle=-1:\quad\zeta=0.$ (14)

In addition, we require that the inner solution (within the boundary layer)
matches with the outer solution (far outside the boundary layer), such that

$P\rightarrow 0:\quad\zeta\rightarrow\infty.$ (15)

The system is not yet complete since the governing equation (12) is a second
order elliptic PDE which requires four boundary conditions. The fourth and
final boundary condition needs more careful thought. Since there is a vertical
velocity (9) descending through the outer region, this induces an arrival of
flux $\epsilon u_{b}L$ within the inner region. However, since this flux can
go neither downwards nor leftwards (due to impermeable/symmetric boundaries),
it must instead exit through the right hand boundary, $x=L$. In other words,
the right hand interface must move outwards to conserve mass, creating a new
finger-like region of vertical size $z\sim\mathcal{O}(\epsilon^{1/2})$, which
we denote region II. Hence, the final boundary condition for region I is given
by an integral constraint of the form

$\int_{0}^{\infty}-\left.\frac{\partial
P}{\partial\xi}\right|_{\xi=1}\,\mathrm{d}\zeta=1.$ (16)

The descent of the upper interface in region I and the resulting finger-like
region II are both illustrated in figure 1b.

Before addressing these details further, we first note that (12)-(16) can be
solved exactly by separation of variables. Hence, the composite solution
(valid across both inner and outer regions) is given by

$p=p_{a}-\rho_{a}gz+2\epsilon^{1/2}\Delta\rho
gL\sum_{n=0}^{\infty}\frac{(-1)^{n}}{\lambda_{n}^{2}}\cos\left[\frac{\lambda_{n}x}{L}\right]\exp\left[-\frac{\lambda_{n}z}{\epsilon^{1/2}L}\right],$
(17)

where $\lambda_{n}=(2n+1)\pi/2$. Conservation of mass within region I
indicates that

$\phi\frac{\partial h}{\partial t}+\frac{k_{H}}{\mu}\frac{\partial}{\partial
x}\left[\int_{0}^{h(x,t)}-\frac{\partial p}{\partial
x}\,\mathrm{d}z\right]=0.$ (18)

Hence, inserting (17) into (18) results in the governing equation for the
evolution of the thickness of the current, which is

$\phi\frac{\partial h}{\partial t}-u_{b}\epsilon L\frac{\partial}{\partial
x}\sum_{n=0}^{\infty}\frac{2(-1)^{n+1}\sin\lambda_{n}x/L}{\lambda_{n}^{2}}\left(1-\exp\left[-\frac{\lambda_{n}h(x,t)}{\epsilon^{1/2}L}\right]\right)=0.$
(19)

This can be further simplified by ignoring exponentially small terms when
$h/L$ is larger than $\mathcal{O}(\epsilon^{1/2})$, and by using the fact that
the infinite sum converges to $-x/L$. Hence, we see that the thickness within
region I is given by

$h=H-\frac{\epsilon u_{b}t}{\phi},$ (20)

which is valid for $0<x<L$ and for $h\gg\mathcal{O}(\epsilon^{1/2})$.

(a)(b)(c)(d)III Figure 2: (a,b,c) Approximate solution for regions I and II
compared with numerical solution in the case of $\epsilon=10^{-2}$. Times are
given in terms of the transition time (37), where (a) $t/t^{*}=0.07$, (b)
$0.17$ and (c) $0.34$. Streamlines are evaluated using velocities derived from
(17) and (21). (d) Further numerical solutions are shown at times
$t/t^{*}=0.42,0.69,1,1.82,3.22$.

Next, we address the fluid flow in region II, which is the finger-like region
of escaped fluid near the base of the current, which is defined for
$L<x<x_{N}(t)$, where $x_{N}(t)$ is the maximum extent of the finger. This
flow region is long and thin (like a classical gravity current) such that the
horizontal velocity is much larger than the vertical velocity. Consequently,
the pressure within region II is hydrostatic to good approximation, such that

$p=p_{a}-\rho_{a}gh-\rho g\left(z-h\right).$ (21)

Similarly to (18), conservation of mass within this region gives

$\phi\frac{\partial h}{\partial t}=u_{b}\frac{\partial}{\partial
x}\left[h\frac{\partial h}{\partial x}\right],$ (22)

which is sometimes called the Dupuit approximation. This is accompanied by
boundary conditions that correspond with imposing the input flux from region
I,

$-u_{b}h\frac{\partial h}{\partial x}=\epsilon u_{b}L:\quad x=L,$ (23)

and imposing zero thickness and zero flux at the moving front,

$\displaystyle h\rightarrow 0:\quad x\rightarrow x_{N}(t),$ (24)
$\displaystyle-u_{b}h\frac{\partial h}{\partial x}\rightarrow 0:\quad
x\rightarrow x_{N}(t).$ (25)

By introducing dimensionless coordinates,

$h=\epsilon^{1/2}L\mathcal{H},\quad x=L(1+\xi),\quad t=\frac{\phi
L}{\epsilon^{1/2}u_{b}}\tau,$ (26)

we get the same system as [5] for a constant input flux (see Appendix B for
further details). The solution is well known and is given in terms of the
similarity variables

$\eta=\xi/\tau^{2/3},\quad\mathcal{H}=\tau^{1/3}f(\eta).$ (27)

The self-similar shape function $f(\eta)$ is defined for
$\eta\in[0,\eta_{N}=1.482]$, and is monotone decreasing from
$f_{0}:=f(0)=1.296$ to $f(\eta_{N})=0$. The solution for regions I and II is
plotted in figure 2a,b,c, at several different times. Streamlines confirm that
the flux into region II is fed by the shrinking of region I. Comparison is
also made to a numerical solution described later in Section 3.

The maximum vertical extent of the flow is given by (20), whereas the maximum
horizontal extent is determined through the above scalings as

$x_{N}(t)=L(1+\eta_{N}\tau^{2/3}).$ (28)

These are plotted in figure 3 with dashed blue lines, thereby indicating the
early-time behaviour of the released fluid.

(a)(b) Figure 3:  Evolution of the vertical (a) and horizontal (b) extents of
the gravity current shown in figure 2 (i.e. for $\epsilon=10^{-2}$), showing
both numerical and approximate solutions. Early time behaviour is given by
(20),(28), whereas late time behaviour is given by (34),(35).

###  2.2 Transition to self-similarity

After a long time the flow is expected to eventually transition to the self-
similar behaviour of a finite release gravity current [5, 12]. As discussed by
[18], the late-time dynamics of a two-dimensional gravity current are
independent of the anisotropy of the medium since the bulk flow decouples from
the ambient. Hence, the anisotropy only affects the late-time behaviour by
delaying the time to transition to self-similarity. Hence, to study the late-
time behaviour we first analyse the isotropic case, which determines the late-
time dynamics, and then use these dynamics to derive the time $t^{*}$ to
transition between the early and late solutions, where $t^{*}$ depends on the
anisotropy $\epsilon$.

As before, the thickness of the gravity current (which now occupies a single
region $0\leq x\leq x_{N}(t)$) satisfies the Dupuit approximation (22). The
boundary conditions are similar to the case of constant input flux, as
described above, except the left hand boundary condition is replaced with the
zero flux condition

$\frac{\partial h}{\partial x}=0:\quad x=0.$ (29)

Consequently, mass conservation indicates that

$\int_{0}^{x_{N}(t)}h\,\mathrm{d}x=HL.$ (30)

By introducing dimensionless coordinates,

$x=(HL)^{1/2}\xi,\quad h=(HL)^{1/2}\mathcal{H},\quad
t=\frac{\phi(HL)^{1/2}}{u_{b}}\tau,$ (31)

and switching to similarity variables

$\eta=\xi/\tau^{1/3},\quad\mathcal{H}=f(\eta)/\tau^{1/3},$ (32)

we arrive at a system of equations that can be solved analytically to give

$f(\eta)=\frac{1}{6}(\eta_{N}^{2}-\eta^{2}),$ (33)

where $\eta_{N}=3^{2/3}$, as shown by [5] (see Appendix B for further
details).

The maximum vertical and horizontal extent of the flow are given by

$\displaystyle h(0,t)$
$\displaystyle=\frac{\eta_{N}^{2}(HL)^{1/2}}{6\tau^{1/3}},$ (34)
$\displaystyle x_{N}(t)$ $\displaystyle=\eta_{N}(HL)^{1/2}\tau^{1/3},$ (35)

respectively. These scalings are plotted in figure 3 with dashed blue lines,
thereby indicating the late-time behaviour of the gravity current.

Next, we discuss the time taken to transition from the early flow regime
involving two fluid regions to the late flow regime with a single region which
is self-similar. In the early regime the vertical extent of the flow descends
according to (20). By contrast, the vertical extent of region II
($\epsilon^{1/2}L\mathcal{H}$ in (26)) increases like $\sim
f_{0}(u_{b}\epsilon^{2}L^{2}t/\phi)^{1/3}$. Hence, it is expected that the
transition to self-similarity will occur once these two thickness scalings
approach each other. Thus, the transition time $t^{*}$ satisfies the cubic
equation

$\left(H-\frac{\epsilon
u_{b}t^{*}}{\phi}\right)^{3}=\frac{f_{0}^{3}u_{b}\epsilon^{2}L^{2}t^{*}}{\phi}.$
(36)

The solution has a complicated form but can be expanded in powers of
$\epsilon\ll 1$ to give

$t^{*}=\frac{L\phi}{u_{b}\epsilon}\left[\alpha-
f_{0}\alpha^{1/3}\epsilon^{1/3}+\frac{f_{0}^{2}}{3\alpha^{1/3}}\epsilon^{2/3}+\ldots\right],$
(37)

where $\alpha=H/L$ is the aspect ratio of the initial current shape. Hence,
anisotropy delays the transition to a classical self-similar regime, which is
consistent with other studies [18]. At the transition time $t^{*}$ the maximum
thickness of the current (which we denote $H^{*}$) is given by

$H^{*}=L\epsilon^{1/3}\left[f_{0}\alpha^{1/3}-\frac{f_{0}^{2}}{3\alpha^{1/3}}\epsilon^{1/3}+\ldots\right].$
(38)

This indicates that, by the time transition occurs for very anisotropic media,
the bulk of the current has shrunk significantly. This hints towards a reduced
swept volume, which we analyse further in Section 4.

The transition time $t^{*}$ and the transition thickness $H^{*}$ (given in
dimensionless terms) are plotted in figure 4. These importantly depend on both
the anisotropy $\epsilon$ as well as the initial aspect ratio of the flow
$\alpha$. It should be noted that in practice strong anisotropy $\epsilon\ll
1$ may cause a stretching of the flow before contact with the impermeable cap
rock. Therefore, $\epsilon\ll 1$ is likely to be correlated with $\alpha\ll
1$.

(a)(b) Figure 4: Transition time $t^{*}$ and the transition thickness $H^{*}$,
determined by solving (36) for different values of the anisotropy $\epsilon$
and the initial aspect ratio $\alpha$. Leading order scalings for small
$\epsilon$ (see (37),(38)) are shown with dotted lines.

##  3 Finite difference computations of Darcy flow

Next, in this section we compare our analytical predictions to finite
difference computations of two-dimensional Darcy flow. The flow is modelled
with the Darcy equations (1)-(2), accounting for different permeability
values, $k_{H}$ and $k_{V}$, in the horizontal and vertical directions. The
governing equations are accompanied by boundary conditions (5)-(6),
corresponding with impermeable/symmetric walls at $x=0$ and $z=0$. The dynamic
boundary condition (7) is imposed on the interface $z=h(x,t)$, which is
interpolated over a gridded mesh of 150$\times$150 points. Likewise, a spatial
domain of finite size $4L\times H$ is chosen in the $x\times z$ directions.
The fluid flow is not resolved beyond the interface $z>h(x,t)$ since pressure
is assumed to be hydrostatic in the ambient fluid. Mass conservation dictates
that the interface evolves according to

$\phi\frac{\partial h}{\partial t}+\frac{k_{H}}{\mu}\frac{\partial}{\partial
x}\left[\int_{0}^{h(x,t)}-\frac{\partial p}{\partial
x}\,\mathrm{d}z\right]=0,$ (39)

with suitable initial conditions $h(x,0)=h_{0}(x)$. Boundary conditions for
$h$ are given by (24), (25) and (29) (see earlier discussion for further
explanation). The time-dependent equation for the fluid-fluid interface (39)
is solved using an explicit forward Euler scheme in time, and a backward
eighth order scheme in space. At each time-step the Darcy equations, (1)-(2),
are solved using a second order central finite difference scheme. The code
used for these computations is available in the Supplementary Materials.

Results from the finite difference computations are compared with the
analytical model in figures 2 and 3 for an anisotropy value $\epsilon=0.01$.
Figure 2a,b,c shows the distorted spreading of the gravity current via a thin
finger near the base. Good agreement is observed everywhere except near
$x\approx L$ where the two regions connect. Here, the interface transitions
smoothly between regions I and II, which is a second order feature that is
missing from the simple analytical model. The numerical solution enables
computation of the gravity current shape up to and beyond the transition time
$t^{*}$, as is displayed in figure 2d. This demonstrates clearly how the flow
transitions from a two-region (bulk/finger) structure at early times to a
slumping single-region structure at late times.

Figure 3 displays numerical computations of the vertical and horizontal
extents of the gravity current across early and late time regimes. Clearly,
the numerical computations reflect the transition between the early and late
time analytical scalings, shown with dashed lines. Likewise, the transition
between these different regimes corresponds with $t/t^{*}\approx 1$,
indicating the accuracy of our prediction for $t^{*}$ (37).

As a benchmark test, we also compare these finite difference computations with
a gravity current slumping in an isotropic medium, $\epsilon=1$. Since our
analytical model only applies for $\epsilon\ll 1$, the initial dynamics
involving regions I and II are not relevant for this case. Instead, we set the
initial shape as

$h_{0}(x)=H\left[1-\left(\frac{x}{L}\right)^{2}\right],$ (40)

which is simply the similarity solution (33). Since the initial shape
satisfies the similarity conditions, this ensures that the solution remains
self-similar for all time. In figure 7a,b,c, in Appendix A, we plot the self-
similar evolution of the gravity current shape at various times. Excellent
agreement is attained between the numerical model and the exact self-similar
solution, indicating the reliability of the finite difference approach. In
figure 7d,e,f, we display similar computations for the same initial shape
released in an anisotropic medium, $\epsilon=0.01$. In this case, the flow is
decomposed into regions I and II, as before, demonstrating how our simple
model can be extended to account for other released shapes. Further details
and discussion of this case are given in Appendix A.

##  4 Swept shape and swept volume

(a)(b)(c)(d)$x/L$ Figure 5: (a,b) Swept shape of the gravity current
$\mathcal{S}(x)$ (42) in the case of $\epsilon=10^{-3}$ and $10^{-2}$. (c)
Self-similar isotropic case $\epsilon=1$ (note the modified initial
conditions). (d) Evolution of the swept volume $\mathcal{V}(t)$ for each of
these cases (note that time is stretched by a factor $\epsilon$). Transition
times $t^{*}$ are indicated with coloured dots. The initial aspect ratio is
set as $\alpha=1$ in all cases.

As described earlier, it is useful for applications (e.g. contaminant
transport or CO2 storage) to quantify the total volume contacted by the
gravity current, also known as the swept volume111Note that since our model is
in two dimensions, the swept volume is taken as per unit width. Also note that
the total volume of fluid is constant and given by
$2HL\phi=\phi\mathcal{V}(0)$.. At early times $t\ll t^{*}$ region I shrinks
uniformly downwards whilst region II grows upwards and outwards. Hence, the
swept volume $\mathcal{V}$ is simply equal to the initial volume plus the
instantaneous volume of region II, such that

$\mathcal{V}(t)=2HL+2\epsilon u_{b}Lt.$ (41)

At much later times $t\gg t^{*}$, once the gravity current has transitioned to
self-similarity, the thickness $z=h(x,t)$ has some parts which are growing and
other parts which are shrinking, so this requires more careful attention. To
deal with this we first define the swept shape $\mathcal{S}$ as the maximum
thickness that the current ever reached at a given value of $x$, such that

$\mathcal{S}(x)=\max_{t}\left\\{h(x,t)\right\\}.$ (42)

The swept volume is then given in terms of $\mathcal{S}$ as

$\mathcal{V}(t)=2\int_{0}^{x_{N}(t)}\mathcal{S}\,\mathrm{d}x.$ (43)

The swept shape (42) is calculated by finding the time at which the thickness
is maximal, which is equivalent to $\eta=3^{1/6}$ in the similarity solution
$f$ (33). Inserting this into (42) we get

$\mathcal{S}(x)=\frac{HL}{3^{1/2}x},$ (44)

which is only valid for $x\gg L$. Hence, the swept volume (up to a constant
$C$) is

$\mathcal{V}(t)=\frac{2HL}{3^{1/2}}\log x_{N}(t)+C,$ (45)

which clearly diverges like $\sim\log t$ as $t\rightarrow\infty$. The constant
of integration is found by equating (41) and (45) at the transition time
$t^{*}$. Hence, the late time behaviour of the swept volume is given by

$\mathcal{V}(t)=2HL+2\epsilon
u_{b}Lt^{*}+\frac{2HL}{3^{1/2}}\log\left(\frac{t}{t^{*}}\right)^{1/3}.$ (46)

In figure 5a,b,d, the swept shape $\mathcal{S}(x)$ and swept volume
$\mathcal{V}(t)$ are plotted for different values of the anisotropy
$\epsilon$. For $\epsilon\ll 1$, the anisotropy restricts the swept shape (for
$x>L)$ to finger-like regions near the base of the current. The maximum
thickness of the fingers is set by $H^{*}$ (38), which scales like
$H^{*}\sim\epsilon^{1/3}$ (i.e. the more anisotropic the medium, the narrower
the fingers). Furthermore, anisotropy delays the transition time (37), since
$t^{*}\propto\epsilon^{-1}$. Meanwhile, before transition occurs, $t\ll
t^{*}$, anisotropy causes the swept volume to grow slowly (since
$\mathcal{V}\propto\epsilon t$ in (41)), such that the gravity current
contacts a smaller fraction of pore space at early times. Note that the kinks
at the corners $x=\pm L$ (in $\mathcal{S}$) and at the transition time
$t=t^{*}$ (in $\mathcal{V}$) would be smoothed out by higher order asymptotics
or numerical computations of two-dimensional Darcy flow. However, such
features do not affect the overall leading order behaviour displayed here.

For comparison, we also plot the swept shape $\mathcal{S}$ in the isotropic
case, $\epsilon=1$, in figure 5c. As described at the end of Section 3, the
isotropic case uses (40) as the initial shape and remains self-similar for all
times. It should be noted that the modified initial shape results in a
different initial swept volume for this case, $\mathcal{V}(0)=4HL/3$. Clearly,
the classical (isotropic) self-similar solution has a larger swept shape and a
faster growing swept volume than the anisotropic cases. In Section 6, we
discuss how these results can be applied when modelling contaminant spills and
CO2 sequestration, noting the possible limitations of our model.

##  5 Finite release with radial symmetry

The above results can be easily extended to account for radially symmetric
flows if the anisotropy remains aligned with the vertical coordinate, i.e.
with permeability $k_{H},k_{V}$ in the radial/vertical directions. To extend
the model from the two-dimensional case, we take the impermeable boundary to
be the horizontal plane, $z=0$. Likewise, we consider the released shape to be
a cylinder of initial height $H$ and radius $R$. Since the initial shape of
the current is radially symmetric, it will spread out and remain radially
symmetric for all times.

Following the derivation in Section 2.1, the thickness of the current evolves
according to (20) at early times. This results in a radial flux of magnitude
$\epsilon u_{b}\pi R^{2}$ exiting region I ($0\leq r\leq R$) into a growing
annular region II ($R\leq r\leq r_{N}(t)$) at the base of the current.
Conservation of mass within region II gives

$\phi\frac{\partial h}{\partial t}=u_{b}\frac{1}{r}\frac{\partial}{\partial
r}\left[rh\frac{\partial h}{\partial r}\right].$ (47)

Boundary conditions correspond with imposing the input flux from region I,

$-2\pi u_{b}rh\frac{\partial h}{\partial r}=\epsilon u_{b}\pi R^{2}:\quad
r=R,$ (48)

and imposing zero thickness and zero flux at the moving front,

$\displaystyle h\rightarrow 0:\quad r\rightarrow r_{N}(t),$ (49)
$\displaystyle-2\pi ru_{b}h\frac{\partial h}{\partial r}\rightarrow 0:\quad
r\rightarrow r_{N}(t).$ (50)

By introducing dimensionless coordinates,

$h=(\pi\epsilon)^{1/2}R\mathcal{H},\quad r=R(1+\xi),\quad t=\frac{\phi
R}{(\pi\epsilon)^{1/2}u_{b}}\tau,$ (51)

we get a similar system to [24] (see Appendix B for further details). The
solution is given in terms of the similarity variables

$\eta=\xi/\tau^{1/2},\quad\mathcal{H}=f(\eta).$ (52)

The shape function $f(\eta)$ is defined for $0<\eta\leq\eta_{N}=1.155$, but
has an unphysical singularity at the origin. This singularity, which is due to
a breakdown of the hydrostatic assumption, can be addressed by introducing a
non-hydrostatic (i.e. source-driven) region near the origin, as discussed in
[18]. However, for the sake of simplicity, we ignore such details for the
present study. For our purposes, (52) serves as a good approximation for the
finger-like growth of region II (i.e. for $r\gg R$). The transition time
$t^{*}$ is found by matching the two thicknesses of regions I and II, such
that

$H-\frac{\epsilon u_{b}t^{*}}{\phi}=(\pi\epsilon)^{1/2}R.$ (53)

Hence, the transition time and thickness are given by

$\displaystyle t^{*}$ $\displaystyle=\frac{\phi R}{\epsilon
u_{b}}\left(\alpha-(\pi\epsilon)^{1/2}\right),$ (54) $\displaystyle H^{*}$
$\displaystyle=(\pi\epsilon)^{1/2}R,$ (55)

where in this case the initial aspect ratio is $\alpha=H/R$. Much later than
the transition time $t\gg t^{*}$, the flow continues to slump as a single
region, similar to Section 2.2. In this case, the entire thickness satisfies
(47). The boundary condition (48) is replaced by a zero flux condition at the
origin, which is

$\frac{\partial h}{\partial r}=0:\quad r=0.$ (56)

Hence, mass conservation indicates that

$2\pi\int_{0}^{r_{N}(t)}rh\,\mathrm{d}r=\pi HR^{2}.$ (57)

By introducing dimensionless coordinates,

$r=(HR^{2})^{1/3}\xi,\quad h=(HR^{2})^{1/3}\mathcal{H},\quad
t=\frac{\phi(HR^{2})^{1/3}}{u_{b}}\tau,$ (58)

and switching to similarity variables

$\eta=\xi/\tau^{1/4},\quad\mathcal{H}=f(\eta)/\tau^{1/2},$ (59)

we arrive at a system of equations that can be solved analytically to give

$f(\eta)=\frac{1}{8}(\eta_{N}^{2}-\eta^{2}),$ (60)

where $\eta_{N}=2$ (see Appendix B for further details).

Next, let’s briefly discuss the swept shape and swept volume for this case. At
early times $t\ll t^{*}$, the swept volume $\mathcal{V}$ is equal to the
initial volume plus the instantaneous volume of region II, such that

$\mathcal{V}(t)=\pi HR^{2}+\epsilon u_{b}\pi R^{2}t.$ (61)

At much later times $t\gg t^{*}$, the swept shape $\mathcal{S}(r)$ is
calculated by finding the time at which the current thickness is maximal,
which is equivalent to $\eta=2^{1/2}$ in the similarity solution $f$ (60).
Hence, the swept shape is

$\mathcal{S}(r)=\frac{HR^{2}}{2r^{2}},$ (62)

which is only valid for $r\gg R$. The swept volume, which is now defined as

$\mathcal{V}(t)=2\pi\int_{0}^{r_{N}(t)}r\mathcal{S}\,\mathrm{d}r.$ (63)

is calculated (following similar steps as in Section 4) for late times as

$\mathcal{V}(t)=\pi HR^{2}+\epsilon u_{b}\pi R^{2}t^{*}+\pi
HR^{2}\log\left(\frac{t}{t^{*}}\right)^{1/4}.$ (64)

To illustrate the radially symmetric case, the swept shape $\mathcal{S}(r)$ is
plotted in figure 6a,b, for $\epsilon=10^{-2}$ and $10^{-3}$. As with the two-
dimensional case, stronger anisotropy results in a reduced swept shape (i.e.
smaller values of $\mathcal{S}$ for $r>R$). Likewise, although not plotted,
the swept volume increases slowly like $\mathcal{V}\propto\epsilon t$ before
transition to self-similarity, indicating that anisotropy reduces the
contacted volume of pore space.

(a)(b) Figure 6: Swept shape of the gravity current $\mathcal{S}(r)$ in the
radially symmetric case. Initial aspect ratio is set to $\alpha=0.2$ and the
anisotropy is set to (a) $\epsilon=10^{-2}$ and (b) $\epsilon=10^{-3}$.

##  6 Discussion and applications

The gravity-driven spreading of a finite volume of fluid in anisotropic porous
media differs from the isotropic case since the vertical flow is restricted by
the permeability $k_{V}\ll k_{H}$. In the initial dynamics the bulk of the
flow descends slowly with a boundary layer near the impermeable base that
diverts the flow into thin finger-like regions growing slowly in the lateral
direction. This partition of the flow into bulk and finger regions reduces the
swept volume of the gravity current compared to the isotropic case. This
indicates that released volumes in anisotropic aquifers contact a smaller
fraction of the available pore space. Hence, the spread of a contaminant in an
anisotropic aquifer may be easier to contain since the contacted volume is
reduced. By contrast, in the case of CO2 sequestration, where the aim is to
trap as much CO2 in the pore space as possible, isotropic aquifers may have
better potential in terms of residual trapping (which is a function of the
contacted volume of pore space), though this ignores other trapping mechanisms
such as dissolution, structural trapping and mineralisation [19].

It is important to consider the possible limitations of this model for such
realistic scenarios. First and foremost, it must be noted that the anisotropic
permeability values $k_{H}$, $k_{V}$, are upscaled quantities which attempt to
capture the macroscopic effect of small-scale heterogeneities on the flow.
These upscaled quantities are only a good approximation when the vertical
length scale of the flow $H$ is much larger than the heterogeneity length
scale (e.g. the width of sedimentary layers) or when the permeable interval is
inherently anisotropic due to compaction effects. If this is not the case,
more complex flow models are required to treat the spread of fluid beneath and
through successive layers [25, 3, 26].

Another consideration for the case of CO2 sequestration is the effect of
trapped saturation on the dynamics and spreading of the current. As the
current moves, a fraction of its mass is lost to residual trapping due to
small scale capillary forces [19] and dissolution within the surrounding brine
[8]. According to some trapping models [7] this can arrest the spread of the
current altogether. Likewise, the spreading could also be arrested by lateral
heterogeneities in the capillary pressure pinning the nose of the gravity
current.

Whilst the conclusion of this study is that anisotropy reduces the swept
volume of the gravity current, this doesn’t take into account the enhanced
trapping potential due to changes in the capillary pressure across
heterogeneities, also known as capillary heterogeneity trapping [27].
Specifically, during the imbibition cycle small-scale capillary forces induce
a build-up of saturation beneath each sedimentary layer that can account for
as much as 14$\%$ of the overall trapped saturation. Hence, the optimum
anisotropy will no doubt strike a balance between the trapping associated with
these heterogeneities and the reduction in swept volume that they induce.

It is also worth mentioning the flow of the ambient fluid, which we have so
far ignored for this study. In particular, wherever there are sharp changes in
the interface shape, this may be associated with non-negligible displacement
of the local ambient fluid, such that the dynamic boundary condition (7)
cannot be applied. Likewise, in the case where there is a viscosity contrast
between the released and ambient fluids, this can modify the shape of the
current [28] and cause fingering instabilities [29]. In both of these cases a
full numerical model would be necessary to resolve such flow details. In the
case of carbon sequestration, CO2 is typically 20-30 times less viscous than
brine. As shown by [28], the viscosity contrast causes an enhanced spreading
of the CO2 in the shape of a thin finger along the cap rock. Hence, this
viscosity contrast distorts the spreading of CO2 in a similar manner to
anisotropy, as studied here. Therefore, the expected effect of a viscosity
contrast in anisotropic media is an extremely pronounced finger-like intrusion
of CO2.

The finite difference code used in this study can be found in the Supplemental
Materials or on the personal website of the author:
https://people.maths.ox.ac.uk/benham/openacces.zip The author reports no
competing interests. There is no funding to report. The author wishes to thank
the anonymous peer reviewers for time taken to review the manuscript.

##  Appendix A Additional plots

(a)(b)(c)(d)(e)(f) Figure 7: Evolution of a gravity current with initial shape
given by (40) in the case of an isotropic medium $\epsilon=1$ (a,b,c) and an
anisotropic medium $\epsilon=0.01$ (d,e,f). Dimensionless times are given by
$u_{b}t/L\phi=$ 0.87, 1.74, 8.68 in (a,b,c) and 1.70, 5.09, 17.0 in (d,e,f).

In this section we present some additional plots comparing the numerical
solution from Section 3 with various approximate and analytical solutions.
Whilst in the previous sections a rectangular shape was chosen for the initial
released shape of dense fluid, here we consider a curved profile with initial
shape $h_{0}(x)$ given by (40). In the case of an isotropic medium
$\epsilon=1$, this results in immediate self-similar behaviour, as described
in Section 2.2. We use the analytical self-similar solution (33) as a means of
verifying our numerical method. In figure 7a,b,c, profiles of the gravity
current are shown at three different times. Overall, very good agreement is
found, indicating the soundness of the numerical method.

Unfortunately, no benchmark analytical solution exists for anisotropic porous
media, but we can nevertheless compare against the approximate solution
derived in this study. In figure 7d,e,f plots are shown for the same curved
initial shape (40) released in an anisotropic medium with $\epsilon=0.01$. It
is straightforward to extend the approximate solution derived earlier to this
initial shape. As such, the bulk fluid region I initially evolves according to

$h=h_{0}(x)-\frac{\epsilon u_{b}t}{\phi}.$ (65)

Meanwhile, region II initially evolves according to the self-similar dynamics
(27) which correspond with a finger-like region fed by a constant input flux
$\epsilon u_{b}L$ from region I (see Section 2.1). These two solutions are
then simply joined together to make the plots in figure 7d,e,f. Overall, good
agreement is achieved indicating that our model can be extended to this and
other such similar initial shapes.

##  Appendix B Further details on similarity solutions

In this section, we summarise the equations that define the different
similarity solutions used in the main text. Let’s start with the two-
dimensional equations for region II at early times. In this case the system of
equations deriving from (22)-(25) after applying coordinate transformations
(26)-(27) is

$\displaystyle\frac{1}{3}f-\frac{2}{3}\eta\frac{\mathrm{d}f}{\mathrm{d}\eta}$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}\eta}\left[f\frac{\mathrm{d}f}{\mathrm{d}\eta}\right],$
(66) $\displaystyle-f\frac{\mathrm{d}f}{\mathrm{d}\eta}$
$\displaystyle=1:\quad\eta=0,$ (67) $\displaystyle f$
$\displaystyle=0:\quad\eta=\eta_{N},$ (68)
$\displaystyle-f\frac{\mathrm{d}f}{\mathrm{d}\eta}$
$\displaystyle=0:\quad\eta=\eta_{N}.$ (69)

These can be solved numerically for the shape function $f(\eta)$ and prefactor
$\eta_{N}=1.482$.

Next we summarise the two-dimensional equations at late times, once the
gravity current has transitioned to a single slumping region. In this case the
system of equations (22),(24),(29),(30), after applying coordinate
transformations (31)-(32), becomes

$\displaystyle-\frac{1}{3}f-\frac{1}{3}\eta\frac{\mathrm{d}f}{\mathrm{d}\eta}$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}\eta}\left[f\frac{\mathrm{d}f}{\mathrm{d}\eta}\right],$
(70) $\displaystyle f$ $\displaystyle=0:\quad\eta=\eta_{N},$ (71)
$\displaystyle\frac{\mathrm{d}f}{\mathrm{d}\eta}$
$\displaystyle=0:\quad\eta=0,$ (72)
$\displaystyle\int_{0}^{\eta_{N}}f\,\mathrm{d}\eta$ $\displaystyle=1.$ (73)

These can be solved analytically to give (33) and $\eta_{N}=3^{2/3}$.

In the radially symmetric case, the governing equations for region II at early
times are (47)-(50). Hence, after applying the coordinate transformations
(51)-(52), we get the system of equations

$\displaystyle-\frac{1}{2}\eta\frac{\mathrm{d}f}{\mathrm{d}\eta}$
$\displaystyle=\left(\frac{1}{1+\eta\tau^{1/2}}\right)\frac{\mathrm{d}}{\mathrm{d}\eta}\left[(1+\eta\tau^{1/2})f\frac{\mathrm{d}f}{\mathrm{d}\eta}\right],$
(74)
$\displaystyle-2\pi\left(\frac{1+\eta\tau^{1/2}}{\tau^{1/2}}\right)f\frac{\mathrm{d}f}{\mathrm{d}\eta}$
$\displaystyle=1:\quad\eta=0,$ (75) $\displaystyle f$
$\displaystyle=0:\quad\eta=\eta_{N},$ (76)
$\displaystyle-2\pi\left(\frac{1+\eta\tau^{1/2}}{\tau^{1/2}}\right)f\frac{\mathrm{d}f}{\mathrm{d}\eta}$
$\displaystyle=0:\quad\eta=\eta_{N}.$ (77)

After sufficiently long times $1\ll\tau\ll\tau^{*}$ (where $\tau^{*}$ is the
dimensionless transition time), the time-dependence is removed from the above
system. In other words, we consider when enough time has passed that the
finger has grown far from the initial shape, but not so long for transition to
occur. In this case, the system of equations simplifies to

$\displaystyle-\frac{1}{2}\eta\frac{\mathrm{d}f}{\mathrm{d}\eta}$
$\displaystyle=\frac{1}{\eta}\frac{\mathrm{d}}{\mathrm{d}\eta}\left[\eta
f\frac{\mathrm{d}f}{\mathrm{d}\eta}\right],$ (78) $\displaystyle-2\pi\eta
f\frac{\mathrm{d}f}{\mathrm{d}\eta}$ $\displaystyle=1:\quad\eta=0,$ (79)
$\displaystyle f$ $\displaystyle=0:\quad\eta=\eta_{N},$ (80)
$\displaystyle-2\pi\eta f\frac{\mathrm{d}f}{\mathrm{d}\eta}$
$\displaystyle=0:\quad\eta=\eta_{N}.$ (81)

These can be solved numerically for the shape function $f(\eta)$ and prefactor
$\eta_{N}=1.155$.

Finally we summarise the radially symmetric equations at late times, once the
gravity current has transitioned to a single slumping region. In this case the
system of equations (47),(49),(56),(57), after applying coordinate
transformations (58)-(59), becomes

$\displaystyle-\frac{1}{2}f-\frac{1}{4}\eta\frac{\mathrm{d}f}{\mathrm{d}\eta}$
$\displaystyle=\frac{1}{\eta}\frac{\mathrm{d}}{\mathrm{d}\eta}\left[\eta
f\frac{\mathrm{d}f}{\mathrm{d}\eta}\right],$ (82) $\displaystyle f$
$\displaystyle=0:\quad\eta=\eta_{N},$ (83)
$\displaystyle\frac{\mathrm{d}f}{\mathrm{d}\eta}$
$\displaystyle=0:\quad\eta=0,$ (84) $\displaystyle 2\int_{0}^{\eta_{N}}\eta
f\,\mathrm{d}\eta$ $\displaystyle=1.$ (85)

These can be solved analytically to give (60) and $\eta_{N}=2$.

## References

  * [1] Guérin A, Devauchelle O, Lajeunesse E. 2014 Response of a laboratory aquifer to rainfall. J. Fluid Mech. 759.
  * [2] Bear J, Cheng A. 2010 Modeling groundwater flow and contaminant transport. Springer.
  * [3] Huppert H, Neufeld J. 2014 The fluid mechanics of carbon dioxide sequestration. Ann. Rev. Fluid Mech. 46, 255–272.
  * [4] Woods A. 2015 Flow in porous rocks. Cambridge University Press.
  * [5] Huppert H, Woods A. 1995 Gravity-driven flows in porous layers. J. Fluid Mech. 292, 55–69.
  * [6] Hesse M, Orr F, Tchelepi H. 2008 Gravity currents with residual trapping. J. Fluid Mech. 611, 35–60.
  * [7] Golding M, Huppert H, Neufeld J. 2017 Two-phase gravity currents resulting from the release of a fixed volume of fluid in a porous medium. J. Fluid Mech. 832, 550–577.
  * [8] MacMinn C, Szulczewski M, Juanes R. 2011 CO2 migration in saline aquifers. Part 2. Capillary and solubility trapping. J. Fluid Mech. 688, 321–351.
  * [9] Barenblatt G. 1952 On some unsteady motions of a fluid and a gas in a porous medium. Prikl. Mat. Makh. 16, 67–78.
  * [10] Kochina I, Mikhailov N, Filinov M. 1983 Groundwater mound damping. Int. J. Eng. Sci. 21, 413–421.
  * [11] Hesse M, Tchelepi H, Cantwel B, Orr F. 2007 Gravity currents in horizontal porous layers: transition from early to late self-similarity. J. Fluid Mech. 577, 363–383.
  * [12] Ball T, Huppert H, Lister J, Neufeld J. 2017 The relaxation time for viscous and porous gravity currents following a change in flux. J. Fluid Mech. 821, 330–342.
  * [13] Zheng Z. 2023 The radial slump of a gravity current in a confined porous layer. Proc. R. Soc. A. 479, 20220696.
  * [14] Cheng P, Bestehorn M, Firoozabadi A. 2012 Effect of permeability anisotropy on buoyancy-driven flow for CO2 sequestration in saline aquifers. Water Resour. Res. 48.
  * [15] Green C, Ennis-King J. 2014 Steady dissolution rate due to convective mixing in anisotropic porous media. Adv. Water Resour. 73, 65–73.
  * [16] De Paoli M, Zonta F, Soldati A. 2017 Dissolution in anisotropic porous media: Modelling convection regimes from onset to shutdown. Phys. Fluids 29, 026601.
  * [17] Jackson S, Krevor S. 2020 Small-scale capillary heterogeneity linked to rapid plume migration during CO2 storage. Geophys. Res. Lett. p. e2020GL088616.
  * [18] Benham G, Neufeld J, Woods A. 2022 Axisymmetric gravity currents in anisotropic porous media. J. Fluid Mech. 952, A23.
  * [19] Krevor S, Blunt M, Benson S, Pentland C, Reynolds C, Al-Menhali A, Niu B. 2015 Capillary trapping for geologic carbon dioxide storage–From pore scale physics to field scale implications. Int. J. Greenh. Gas Con. 40, 221–237.
  * [20] Corbett P, Jensen J. 1992 Variation of reservoir statistics according to sample spacing and measurement type for some intervals in the Lower Brent Group. The Log Analyst 33.
  * [21] Soltanian M, Amooie M, Dai Z, Cole D, Moortgat J. 2016 Critical dynamics of gravito-convective mixing in geological carbon sequestration. Sci. Rep. 6, 1–13.
  * [22] Gilmore K, Sahu C, Benham G, Neufeld J, Bickle M. 2022 Leakage dynamics of fault zones: Experimental and analytical study with application to CO2 storage. J. Fluid Mech. 931.
  * [23] Benham G. 2023 The near-field shape and stability of a porous plume. J. Fluid Mech. 955, A13.
  * [24] Lyle S, Huppert H, Hallworth M, Bickle M, Chadwick A. 2005 Axisymmetric gravity currents in a porous medium. J. Fluid Mech. 543, 293–302.
  * [25] Neufeld J, Huppert H. 2009 Modelling carbon dioxide sequestration in layered strata. J. Fluid Mech. 625, 353–370.
  * [26] Hewitt D. 2022 Evolution of convection in a layered porous medium. J. Fluid Mech. 941.
  * [27] Harris C, Jackson S, Benham G, Krevor S, Muggeridge A. 2021 The impact of heterogeneity on the capillary trapping of CO2 in the Captain Sandstone.. Int. J. Greenh. Gas Con. 112, 103511.
  * [28] Pegler S, Huppert H, Neufeld J. 2014 Fluid injection into a confined porous layer. J. Fluid Mech. 745, 592–620.
  * [29] Saffman P, Taylor G. 1958 The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid. Proc. R. Soc. A. 245, 312–329.