# Axisymmetric gravity currents in anisotropic porous media

Graham P. Benham1,2<EMAIL_ADDRESS>Jerome A.
Neufeld1,2,3  Andrew W. Woods1,2  1 Department of Earth Sciences, University
of Cambridge, Cambridge CB3 0EZ, UK 2 Centre for Environmental and Industrial
Flows, University of Cambridge, Cambridge CB3 0EZ, UK 3 Department of Applied
Mathematics and Theoretical Physics, University of Cambridge, Cambridge CB3
0WA, UK

###### Abstract

We explore the motion of an axisymmetric gravity current in an anisotropic
porous medium in which the horizontal permeability is larger than the vertical
permeability. It is well known that the classical axisymmetric gravity current
supplied by a constant point source of fluid has an unphysical singularity
near the origin. We address this by considering a pressure-dominated region
near the origin which allows for vertical flow from the source, such that the
current remains of finite depth, whilst beyond this region the flow is
gravity-dominated. At early times the inner pressure-driven region controls
the spreading of the current, but at late times the inner region occupies a
progressively smaller fraction of the current such that the radius increases
as $\sim t^{3/7}$, while the depth near the origin increases approximately as
$\sim t^{1/7}$. The presence of anisotropy highlights this phenomenon, since
the vertical permeability maintains an effect on the flow at late times
through the pressure-driven flow near the origin. Using these results we
provide some quantitative insights into the dominant dynamics which control
CO2 migration through permeable aquifers, as occurs in the context of carbon
capture and storage.

##  1 Introduction

Buoyancy-driven flows in porous media resulting from the injection of fluid of
different density to the original reservoir fluid have been studied in some
detail owing to their importance for geothermal power production, carbon
capture and storage and enhanced oil recovery, for example. The initial models
of gravity-driven flow in a porous medium by Barenblatt (1996) and Huppert &
Woods (1995) established some simple similarity solutions for two-dimensional
gravity currents resulting from steady injection, and tested these models with
analogue laboratory experiments in isotropic porous media. These models were
based on the assumption that the vertical pressure gradient is everywhere
hydrostatic and that the flow becomes long and thin, so that the predominantly
horizontal flow is driven by the horizontal pressure gradient associated with
variations in the thickness of the current. The solution takes the form

$h/L=({t}/{\tau})^{1/3}f\left[(x/L)(t/\tau)^{-2/3}\right],$ (1.1)

where $L=Q_{x}\mu/k\Delta\rho g$ and $\tau=Q_{x}(\mu/k\Delta\rho g)^{2}$ are
dimensional length and time scalings, given in terms of the input flow rate
per unit length, $Q_{x}$, the permeability, $k$ (assumed isotropic), the
density difference between the two fluids, $\Delta\rho$, and the viscosity of
the injected fluid, $\mu$ (Huppert & Woods, 1995). If one extends the analysis
to account for different permeabilities in the horizontal and vertical
directions, $k_{H}$ and $k_{V}$, the assumption of hydrostatic pressure in the
flow means that the vertical permeability does not feature in the solution,
which is the same as for the isotropic case. This is curious since the
vertical permeability of the porous medium may be much smaller than the
horizontal permeability owing to the original geological processes of
formation and compaction of the medium (Corbett & Jensen, 1992; Martinius et
al., 1999; Woods, 2015). The key to this apparent paradox is to assess the
vertical pressure gradient required to drive the vertical flow which is
implicit in the gravity current solution. With a very small vertical
permeability, this pressure gradient scales as

$\frac{\partial p}{\partial z}\sim\frac{\mu}{k_{V}}\frac{\partial h}{\partial
t}\sim\frac{\mu L}{\tau^{1/3}t^{2/3}k_{V}}.$ (1.2)

For the solution to be valid we require that (1.2) is small compared to the
corresponding hydrostatic pressure gradient

$\frac{\partial p}{\partial z}\sim{\Delta\rho g}.$ (1.3)

This requires that

$t/\tau\gg\left({k_{H}}/{k_{V}}\right)^{3/2},$ (1.4)

and we see the role of the vertical permeability in determining the time at
which the flow adjusts to the similarity solution. It is worth noting that
this result provides an extension of the analysis presented by Huppert &
Pegler (2022) who considered the early time pressure solution in isotropic
media.

For radial injection from a point source, the situation is more complex. In
the classical similarity solution (Lyle et al., 2005) there is a singularity
in the calculated depth of the current at the origin. However, this result is
unphysical, and in order that the flow remains of finite depth at the origin
there is a near-origin adjustment zone for all time. Inclusion of this near-
origin adjustment zone leads to a different set of scaling laws and a new
solution for the current that is influenced by $k_{V}$ for all time. In this
study we address these details using a theoretical analysis in conjunction
with numerical simulations, and we show the analysis is consistent with a
series of new laboratory experiments. Our study goes beyond that of Huppert &
Pegler (2022), where it was assumed that the pressure-driven flow simply
adjusts to the classical similarity solution.

We consider the application of our results to the case of CO2 sequestration
and illustrate how anisotropy of the porous medium leads to the ratio
$k_{V}/k_{H}$ influencing the dispersal of the CO2 at all times, in contrast
to the classical gravity current solution. Furthermore, we show that for
typical injection rates the transition in flow regime from a pressure-
dominated flow to a gravity-dominated flow may occur once the flow has reached
depths of $10-1000$ m at time scales of $1-10$ years, depending on the value
of $k_{V}/k_{H}$. This suggests that the initial pressure-driven flow persists
for a significant period of the injection, and in some cases the flow may
never become gravity-driven over the time and length scales relevant to those
sites.

The structure of the paper is laid out as follows. Section 2 addresses the
flow scenario, treated analytically at early and late times, and comparisons
are made with both numerical simulations and porous bead experiments. Section
3 applies the results to the context of carbon sequestration, calculating the
criteria for whether several sites are in a pressure- or gravity-driven
regime, and Section 4 closes with a discussion of these results.

##  2 Axisymmetric injection into anisotropic porous media

###  2.1 A note on anisotropic permeability

Before discussing the flow scenario, we first briefly discuss the geological
context of anisotropy in porous media. In subsurface geological reservoirs it
is common for the permeability of rocks to differ significantly depending on
the direction of the flow (Corbett & Jensen, 1992; Martinius et al., 1999).
This may result from post-depositional compaction of the formation or from the
deposition of successive layers of fine and coarse material. Hence, the
permeability field is sometimes written as a three-dimensional diagonal
matrix, $\mathbf{k}=\mathrm{Diag}(k_{x},k_{y},k_{z})$, or equivalently in a
different basis depending on the direction of compaction (e.g. see studies on
cross-bedding (Woods, 2015)).

We account for such situations, but we restrict our attention to the case
where the compaction is aligned with the vertical coordinate, and hence the
permeability is reduced to horizontal and vertical variation only,
$(k_{x},k_{y},k_{z})=(k_{H},k_{H},k_{V})$, where $k_{H}$ and $k_{V}$ are
constants. Hence, the anisotropy is characterised by the single dimensionless
parameter

$K=k_{H}/k_{V}.$ (2.1)

This formulation is relevant either for flow through a single thick
sedimentary layer that has been compacted vertically, resulting in a binary
permeability field, or for flow through a system of many horizontal
sedimentary layers, in which the values $k_{H}$, $k_{V}$, can be interpreted
as effective permeabilities (Cardwell & Parsons, 1945; Kumar et al., 1997;
Woods, 2015). For example, in the latter scenario, if the layers vary between
permeabilities $k_{1}$ and $k_{2}$ over alternating layer widths $d_{1}$ and
$d_{2}$, then the effective permeability values are given by the arithmetic
and harmonic mean values

$k_{H}\approx\frac{d_{1}k_{1}+d_{2}k_{2}}{d_{1}+d_{2}},\quad
k_{V}\approx\frac{d_{1}+d_{2}}{d_{1}/k_{1}+d_{2}/k_{2}}.$ (2.2)

This tends to be a good approximation as long as the depth of the flow is much
larger than the widths of the layers $d_{1},d_{2}$. Hence, in this system of
horizontal sedimentary layers the anisotropy is always such that $K\geq 1$. In
fact, for many geological systems the horizontal and vertical permeability
have been observed to differ by several orders of magnitude (Martinius et al.,
1999), resulting in anisotropy values as large as $K=\mathcal{O}(10^{4})$.

###  2.2 The flow scenario at early times and subsequent transition behaviour

Next we move on to describe the flow scenario we consider, and describe how
this evolves at early times, before the effects of gravity become significant.
For the isotropic case, as discussed by Huppert & Pegler (2022), the constant
input of fluid from a point source results in the current expanding in the
shape of a self-similar hemisphere at early times. The radial and vertical
extent of the current are equal and grow like the cube root of time. The flow
transitions to a gravity-driven regime once the weight of the fluid dominates
over the injection pressure. In the following analysis, we extend this early-
time analysis to the case of injection into an anisotropic medium,
demonstrating how the anisotropy $K$ affects both the early-time dynamics and
the transition behaviour.

We consider the constant axisymmetric injection of fluid (at flow rate $Q$)
into a porous medium with anisotropy $K$ and porosity $\phi$, bounded from
below by an impermeable substrate, as illustrated in figure 1. The porous
medium is initially saturated with an ambient phase which has a lighter
density than the injected fluid $\rho_{2}<\rho_{1}$, whereas the viscosities
of the two fluids are assumed to be equal $\mu_{1}=\mu_{2}=\mu$. For
simplicity, we assume that there is no mixing between the these two fluids,
such that the interface between them remains sharp. We choose a cylindrical
polar coordinate system $(r\geq 0,z\geq 0)$ and we denote the radial and
vertical extent of the current, $R(t)$ and $H(t)$, which increase over time
due to the constant injection rate $Q$ from a point at the origin.

$\boldsymbol{z}$$\boldsymbol{H}$$\boldsymbol{R}$

$\boldsymbol{Q}$

$\boldsymbol{z}$$\boldsymbol{r}$$\boldsymbol{r}$$\boldsymbol{R_{p}}$$\boldsymbol{\rho=\rho_{1}}$$\boldsymbol{\rho=\rho_{2}<\rho_{1}}$$\boldsymbol{\mu,\phi}$$\boldsymbol{k_{V}}$$\boldsymbol{k_{H}}$Anisotropy(b)
$\boldsymbol{t\gg t^{*}}$(a) $\boldsymbol{t\ll t^{*}}$ Figure 1: Schematic
diagram illustrating the transition from a pressure-driven flow at early times
(a) to a gravity-driven flow with a pressure-driven boundary layer near the
origin at late times (b) in an anisotropic porous medium. Contour lines are
indicative of the typical pressure field and approximate streamlines are
sketched.

At early times the pressure is dominated by the viscous resistance due to
injection and gravity has a negligible effect. Hence, due to conservation of
mass the pressure must satisfy Laplace’s equation. However, since the medium
is anisotropic this results in non-uniform coefficients, such that

$k_{H}\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial p}{\partial
r}\right)+k_{V}\frac{\partial^{2}p}{\partial z^{2}}=0.$ (2.3)

Since $k_{H}$ and $k_{V}$ are constants, the dynamics can be easily mapped
from an isotropic medium, and hence the lateral and vertical extents must
satisfy

$H^{2}/R^{2}=1/K.$ (2.4)

Likewise, conservation of mass (within a stretched hemisphere) dictates that

$2\pi\phi HR^{2}/3=Qt.$ (2.5)

This indicates that the early time dynamics are described by

$\displaystyle H$ $\displaystyle=K^{-1/3}(3Qt/2\pi\phi)^{1/3},$ (2.6)
$\displaystyle R$ $\displaystyle=K^{1/6}(3Qt/2\pi\phi)^{1/3}.$ (2.7)

Hence, at early times the effect of anisotropy is to stretch/squash the flow
in the more/less permeable direction. It follows that the thickness of the
current, $z=h(r,t)$, is given by

$h/H=\left[1-(r/R)^{2}\right]^{1/2},$ (2.8)

at early times (see figure 2a). It is surprising to note that (within this
early time regime) a strongly anisotropic medium results in a current which is
very long and thin, but in which the pressure is not hydrostatic. This is
contrary to the common assumption that long and thin currents automatically
imply a hydrostatic pressure.

After a significant amount of time the current grows to a thickness where the
weight of fluid begins to dominate over the injection pressures. This is
equivalent to the moment when the input flow rate $Q$ is matched by the flow
rate due to hydrostatic pressure gradients. For a hydrostatic current, we have
that $\partial p/\partial r\approx\Delta\rho g\partial h/\partial r$, so that
the integrated gravity-driven flux scales like

$Q_{g}\sim-\frac{k_{H}\Delta\rho g}{\mu}rh\frac{\partial h}{\partial r}\sim
Q.$ (2.9)

By inserting approximate scalings, $h\sim H$, $r\sim R$, $\partial h/\partial
r\sim-H/R$, this provides the balance

$Q_{g}\sim u_{b}H^{2}\sim Q,$ (2.10)

where the buoyancy velocity $u_{b}=k_{H}\Delta\rho g/\mu$. This immediately
provides an expression for the thickness of the current at the transition
between pressure- to buoyancy-driven flow,

$H^{*}=(Q/u_{b})^{1/2},$ (2.11)

and by equating this to the early time behaviour (2.6) we find that the
transition time is

$t^{*}=(2K\phi\pi/3)(Q/u_{b}^{3})^{1/2}.$ (2.12)

Therefore, anisotropy causes the transition to occur at later times, as
expected. Indeed, the significance of this delayed transition is immediately
appreciable when one considers that some geological formations have anisotropy
values as large as $K=\mathcal{O}(10^{4})$ (Martinius et al., 1999; Bergmo et
al., 2017).

The critical time $t^{*}$ distinguishes two distinct regimes. At early times
$t\ll t^{*}$, or when $H\ll H^{*}$, the effects of gravity are negligible and
injection pressures dominate the flow, whereas at late times $t\gg t^{*}$, or
when $H\gg H^{*}$, gravity plays a significant role. However, clearly the flow
must remain pressure-driven very close to the source, even at late times. This
pressure-driven region is responsible for distributing the flow across the
vertical extent of the current, whereupon it diverts to the remaining gravity-
driven region far away from the origin (see figure 1b). Whilst various authors
have alluded to the existence of this pressure-driven boundary layer, it is
not known how its lateral and vertical extent evolve over time, nor how it
couples with the remaining gravity current. Furthermore, it is unknown how
anisotropy affects the flow after transition occurs, which is the subject of
the following sections.

###  2.3 Late-time dynamics: classical analysis

At much later times $t\gg t^{*}$, once the flow has grown to a thickness $H\gg
H^{*}$, it retains much of the character of a classical porous gravity
current. Lyle et al. (2005) described the late-time dynamics in an isotropic
porous medium by assuming a hydrostatic pressure profile everywhere within the
current. This results in the corresponding radially symmetric thin-film
equation

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

for the evolution of the current thickness $z=h(r,t)$. The above equation is
accompanied by boundary conditions imposing zero thickness at the nose,
constant injection flux at the origin, and zero flux through the nose, such
that

$\displaystyle h$ $\displaystyle=0,\quad\mathrm{at}\quad r=R(t),$ (2.14)
$\displaystyle-2\pi rh\frac{\partial h}{\partial r}$ $\displaystyle\rightarrow
Q,\quad\mathrm{as}\quad r\rightarrow 0,$ (2.15) $\displaystyle-2\pi
rh\frac{\partial h}{\partial r}$ $\displaystyle\rightarrow
0,\quad\mathrm{as}\quad r\rightarrow R(t).$ (2.16)

It should be noted that, whilst these equations were formulated for the case
of an isotropic porous medium, the analysis can be easily transposed to the
case of anisotropic permeability, $k_{H}$, $k_{V}$. However, the vertical
permeability $k_{V}$ vanishes from the analysis since the injected flow and
the ambient flow are decoupled (an assumption which we revisit later),
resulting in an identical system of equations.

$\boldsymbol{R_{p}(t)}$(a)(b)(c) Figure 2: Numerical solution of the current
shape $z=h(r,t)$, illustrating analytical solutions with red lines. Early-time
profile in (a) ($t/t^{*}=10^{-3}$) is normalised by early scalings
(2.6)-(2.7). Transition regime profiles (b) ($t/t^{*}=1-10^{2}$) have no
normalisation. Late-time profile in (c) ($t/t^{*}=10^{3}$) is normalised by
late scalings (2.41)-(2.42). The position of the late-time boundary layer near
the origin $r=R_{p}(t)$ (relatively shrinking) is indicated in blue.

The system of equations (2.13)-(2.16) admits the similarity variables

$h=(Q/u_{b})^{1/2}f(\eta/\eta_{N}),\quad\eta=r(Qu_{b})^{-1/4}(t/\phi)^{-1/2},$
(2.17)

where the dimensionless shape $f$ and prefactor $\eta_{N}$ satisfy

$\displaystyle-\frac{1}{2}\eta f^{\prime}$
$\displaystyle=\frac{1}{\eta}\frac{\mathrm{d}}{\mathrm{d}\eta}\left[\eta
f\frac{\mathrm{d}f}{\mathrm{d}\eta}\right],$ (2.18) $\displaystyle f$
$\displaystyle=0,\quad\mathrm{at}\quad\eta=\eta_{N},$ (2.19)
$\displaystyle-2\pi\eta f\frac{\mathrm{d}f}{\mathrm{d}\eta}$
$\displaystyle\rightarrow 1,\quad\mathrm{as}\quad\eta\rightarrow 0,$ (2.20)
$\displaystyle-2\pi\eta f\frac{\mathrm{d}f}{\mathrm{d}\eta}$
$\displaystyle\rightarrow 0,\quad\mathrm{as}\quad\eta\rightarrow\eta_{N}.$
(2.21)

The solution can be calculated numerically, giving a prefactor value
$\eta_{N}=1.155$, and the resultant shape is plotted in figure 2c as a red
dash-dot line. The influx condition (2.20) results in singular behaviour of
the thickness $h$ near the origin, which is a consequence of the hydrostatic
pressure assumption. This is unphysical and can be addressed by considering
the pressure-driven flow near the origin. In particular, the flux within this
region is not buoyancy-driven, and so the inflow boundary condition (2.20)
(resulting in singular behaviour of the thickness) no longer applies. Instead,
the inflow condition requires a non-hydrostatic flux component within the
pressure-driven zone, and as we will see, this maintains a dominant
contribution to the flow (even at late times) such that the thin-film
approximation (2.13) cannot be used near the origin.

###  2.4 Back to the future: late-time dynamics revisited

$\boldsymbol{z}$$\boldsymbol{r}$$\boldsymbol{R_{p}\sim\xi_{0}\,t^{a}}$$\boldsymbol{R\sim\zeta_{N}\,t^{b}}$

$\boldsymbol{H\sim\xi_{0}\,t^{a}}$

Inner(pressure-driven)Outer(gravity-driven) Figure 3: Schematic diagram
showing the different variables of the late-time regime. The analysis in
Section 2.4 shows that the power law values are $a=1/7,b=3/7$. Note the aspect
ratio is exaggerated for illustration purposes, but in fact $H$ is smaller
than $R_{p}$ due to anisotropy (2.22).

Within the inner pressure-driven region strong vertical velocities deliver the
injected flow all the way to the uppermost extent of the current. The inner
pressure-driven radius and the thickness scaling are related (for the same
reasons as (2.4)) according to

$H^{2}/R_{p}^{2}=1/K.$ (2.22)

Our isotropic porous bead experiments, discussed later in Section 2.6, show
that the thickness of the current near the origin increases slowly over time
(e.g. see figure 5). Therefore, we attempt to describe the solution within
this inner region using a set of rescaled variables where we allow for the
possibility of an evolving thickness at the origin. Using a weak power law
value $0<a<1$, the variables are rescaled according to

$\xi=r/R_{p}(t),\quad h=K^{-1/2}\mathcal{G}(\xi,t)R_{p}(t),\quad
R_{p}(t)=\xi_{0}u_{b}^{\frac{3a-1}{2}}Q^{\frac{1-a}{2}}(t/\phi)^{a},$ (2.23)

where $\mathcal{G}$ is the dimensionless shape function, $\xi_{0}$ is a
dimensionless constant that we determine later, and the powers of $u_{b}$ and
$Q$ in (2.23) are chosen to be dimensionally consistent.

The consequence of having a growing inner region is that the radial extent of
the corresponding outer region does not grow according to the classical $R\sim
t^{1/2}$ power law described earlier. In particular, conservation of mass (in
conjunction with (2.22)) dictates that

$\pi\phi HR^{2}\sim\pi\phi K^{-1/2}R_{p}R^{2}\sim Qt.$ (2.24)

If $a>0$, it follows that $R$ grows as $\sim t^{b}$, where $b<1/2$. Hence, we
introduce a corresponding set of rescaled outer variables

$\zeta=r/R(t),\quad h=K^{-1/2}\mathcal{F}(\zeta,t)R_{p}(t),\quad
R(t)=\zeta_{N}u_{b}^{\frac{3b-1}{2}}Q^{\frac{1-b}{2}}(t/\phi)^{b},$ (2.25)

where $\mathcal{F}$ is the dimensionless shape profile, the nose position
$\zeta_{N}$ is a dimensionless constant to be determined, and $h$ is rescaled
to match with the inner region. Mass conservation (2.24) indicates that

$a+2b=1.$ (2.26)

Given these rescaled variables (2.25), the thin film equation (2.13) in the
outer (hydrostatic) region becomes

$\frac{R^{2}}{R_{p}}\frac{\partial\mathcal{F}}{\partial
t}-\frac{RR^{\prime}}{R_{p}}\zeta\frac{\partial\mathcal{F}}{\partial\zeta}+\frac{R^{2}R_{p}^{\prime}}{R_{p}^{2}}\mathcal{F}=\frac{u_{b}}{K^{1/2}\phi\zeta}\frac{\partial}{\partial\zeta}\left[\zeta\mathcal{F}\frac{\partial\mathcal{F}}{\partial\zeta}\right].$
(2.27)

The second two terms on the left hand side of (2.27) are of the order
$\mathcal{O}(t^{-2a})$ and so become vanishingly small at large times.
Assuming the time derivative of the scaled shape decreases as
${\partial\mathcal{F}}/{\partial t}\sim 1/t$ for large times, then the first
term on the left hand side is also of the order $\mathcal{O}(t^{-2a})$. Hence,
after sufficiently late times (2.27) yields the result

$-2\pi\zeta\mathcal{F}\frac{\mathrm{d}\mathcal{F}}{\mathrm{d}\zeta}\approx\mathcal{Q},$
(2.28)

where $\mathcal{Q}$ is the scaled flux of injected fluid, which is uniform in
magnitude across the extent of the gravity current, and which is yet unknown.
Upon further integration and applying the boundary condition (2.14), we arrive
at an expression for the outer solution, which is

$\mathcal{F}\approx\left(-\frac{\mathcal{Q}}{\pi}\log\zeta\right)^{1/2}.$
(2.29)

We note that the integrated gravity-driven flux, which is given by

$Q_{g}=-2\pi u_{b}{r}{h}\frac{\partial{h}}{\partial{r}},$ (2.30)

grows according to

$Q_{g}=K^{-1}\mathcal{Q}u_{b}R_{p}(t)^{2}.$ (2.31)

At the transition thickness, $H=K^{-1/2}R_{p}=H^{*}=(Q/u_{b})^{1/2}$, the
gravity-driven flux (2.31) equals the input flux $Q$, indicating that
$\mathcal{Q}=1$. We also note that, after the transition time $t>t^{*}$, the
gravity-driven flux (2.31) continues to grow unbounded over time. Hence, this
can satisfy neither the boundary condition at the origin (2.15), nor that at
the nose (2.16). Losing the ability to impose the boundary condition at the
nose indicates that a shock profile will form, over which the flux
discontinuously drops to zero (this is a feature of such equations, e.g. see
(Whitham, 2011)) Hence, in terms of the classification of partial differential
equations, the problem is hyperbolic, unlike the case of a two-dimensional
gravity current due to a line source (Huppert & Woods, 1995) which is
parabolic.

Next we consider conservation of mass to determine the model parameters. In
particular, (2.13), (2.15) and (2.16) can be combined to give a single volume
integral

$2\pi\phi\int_{0}^{R}rh\,\mathrm{d}r=Qt.$ (2.32)

Hence, by inserting (2.29) into (2.32), we arrive at a relationship involving
$\zeta_{N}$ and $\xi_{0}$, which is

$2\pi\zeta_{N}^{2}\xi_{0}(\pi K)^{-1/2}\int_{\epsilon(t)}^{1}y(-\log
y)^{1/2}\,\mathrm{d}y=1,$ (2.33)

where $\epsilon=R_{p}(t)/R(t)=\mathcal{O}(t^{-(b-a)})$. Noting that $\epsilon$
becomes vanishingly small at late times (since $b>a$), and noting that
$\int_{0}^{1}y(-\log y)^{1/2}\,\mathrm{d}y=(\pi/32)^{1/2}$, we can rearrange
(2.33) to give

$\zeta_{N}^{2}\xi_{0}\approx{{(8K)^{1/2}}}/{\pi}.$ (2.34)

A further final condition is required to fully determine the coefficients
$\zeta_{N}$, $\xi_{0}$, and this is provided by matching the flux of the
current between the inner and outer regions. Moving from the outer region to
the origin the thickness tends to a finite but growing value $h=H(t)$ (with
zero slope) and a balance is sustained at the matching radius $r=R_{p}$, so
the gravity-driven flux becomes

$-2\pi u_{b}{r}{h}\frac{\partial{h}}{\partial{r}}=2\pi
K^{-1}u_{b}{R_{p}^{3}}{R^{-1}}=Q,\quad\mathrm{at}\quad r=R_{p}(t).$ (2.35)

Hence, (2.35) suggests that

$b=3a,$ (2.36)

which together with (2.26), gives the power law values

$a=1/7,\quad b=3/7.$ (2.37)

Furthermore, (2.35) implies that the coefficients satisfy

$\xi_{0}^{3}/\zeta_{N}=K/2\pi.$ (2.38)

Hence, the two equations (2.34) and (2.38), can be solved simultaneously to
give

$\displaystyle\xi_{0}=(2^{1/2}\pi^{3})^{-1/7}K^{5/14}$ $\displaystyle\approx
0.58K^{5/14},$ (2.39)
$\displaystyle\zeta_{N}=(2^{11/2}\pi^{-2})^{1/7}K^{1/14}$
$\displaystyle\approx 1.24K^{1/14}.$ (2.40)

In addition, the thickness and radial scalings are given by

$\displaystyle H$
$\displaystyle=K^{-1/7}(Q^{3}u_{b}^{-2}\pi^{-3}2^{-1/2})^{1/7}(t/\phi)^{1/7},$
(2.41) $\displaystyle R$
$\displaystyle=K^{1/14}(2^{11/2}Q^{2}u_{b}\pi^{-2})^{1/7}(t/\phi)^{3/7},$
(2.42) $\displaystyle R_{p}$
$\displaystyle=K^{5/14}(Q^{3}u_{b}^{-2}\pi^{-3}2^{-1/2})^{1/7}(t/\phi)^{1/7},$
(2.43)

which now all reflect a dependence on the anisotropy $K$. We note one further
detail of this analysis. Whilst $H(t)$ (2.41) captures the overall thickness
scaling of the gravity-current, it does not reflect the precise value at the
origin. In particular, if we insert the matching radius $r=R_{p}$ into the
thickness profile in the outer region (2.29), we find that the thickness near
the origin actually grows like

$h(R_{p}(t),t)=K^{-1/2}R_{p}(-(1/\pi)\log R_{p}/R)^{1/2}\sim(2/(7\pi K)\log
t)^{1/2}t^{1/7}.$ (2.44)

However, since the relative growth of the logarithmic component is very slow,
this only becomes appreciable at time scales $t/t^{*}$ larger than around
$e^{7\pi K/2}\approx(5.96\times 10^{4})^{K}$. To summarise, the scaling for
$H(t)$ (2.41) provides the general behaviour of the thickness everywhere in
the current except near the origin, and (2.44) must be used to capture the
near-origin behaviour at very late times.

###  2.5 Comparison with numerical simulations

$\boldsymbol{R\sim t^{3/7}}$$\boldsymbol{R\sim t^{1/3}}$(a)

$\boldsymbol{H\sim t^{1/7}(\log t)^{1/2}}$$\boldsymbol{H\sim
t^{1/7}}$$\boldsymbol{H\sim t^{1/3}}$(b)

Figure 4: Analytical, numerical and experimental data for the radial (a) and
vertical (b) extents of the current, $R$, $H$. Numerical results are shown for
both an isotropic medium $K=1$ and an anisotropic medium $K=10^{4}$, whereas
experimental results are only for the isotropic case. Analytical early and
late time scalings (2.6)-(2.7),(2.41)-(2.43), are shown with solid black lines
and alternate scalings derived by Lyle et al. (2005) ($R\sim t^{1/2},H\sim
H^{*}$) are shown with dotted green lines. Best fit power law
$H/H^{*}\propto(t/t^{*})^{0.1537\pm 0.0176}$ for the late time thickness data
is shown as an insert in (b).

We compare our analytical predictions for the early-time, late-time, and
transition behaviour, to finite difference numerical simulations of
axisymmetric Darcy flow in an anisotropic porous medium. Before discussing
these results, we first briefly outline the numerical method. We model the
flow using the Darcy equations (in cylindrical polar coordinates), which are

$\displaystyle\mathbf{u}_{i}=-\frac{1}{\mu}\mathbf{k}\nabla\left(p+\rho_{i}gz\right),\quad\nabla\cdot\mathbf{u}_{i}=0,\quad\mathrm{in}\quad\Omega_{i},\quad
i=1,2,$ (2.45)

where $\rho=\rho_{1}$ inside the injected fluid domain (which we denote
$\Omega_{1}$) and $\rho=\rho_{2}$ outside the injected fluid domain (which we
denote $\Omega_{2}$). The anisotropic permeability matrix $\mathbf{k}$ is
accounted for by imposing different values, $k_{H}$ and $k_{V}$, in the radial
and vertical directions.

The numerical solution at $r=0$ is not available due to the singular nature of
the divergence in the continuity equation $\nabla\cdot\mathbf{u}=0$. Instead,
the domain begins at a small but finite value of $r_{0}=0.05H^{*}$ close to
the origin (we note that the solution is insensitive to the precise value
$0.05$). The governing equations (2.45) are accompanied by boundary conditions
corresponding to impermeable/symmetric walls at $r=r_{0}$ and $z=0$, except
for a constant input flux $Q$ injected over a small region near the origin,
and far-field (Dirichlet) pressure conditions imposed at the boundaries of a
large but finite domain.

The interface between the domains $\Omega_{1}$ and $\Omega_{2}$ is given by a
set of points $\partial\Omega=\mathbf{r}(t)$, which is treated as a passive
scalar and is therefore advected at the fluid velocity, such that

$\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t}=\frac{1}{\phi}\mathbf{u}(\mathbf{r}),$
(2.46)

with initial conditions given by a small surface
$\mathbf{r}(0)=\mathbf{r}_{0}$ close to the origin. Following the Lagrangian
approach, the interface is discretised and grid points are advected along
streamlines. Therefore, at every time-step the domains $\Omega_{1}$ and
$\Omega_{2}$ are updated via interpolation over a gridded mesh of
150$\times$150 points. Likewise the grid is stretched over time to account for
the fact that the fluid region grows across several orders of magnitude. The
dynamic equation for the fluid-fluid interface (2.46) is solved using an
explicit fourth order Runge-Kutta scheme with an adaptive time-step, where at
each time-step the Darcy equations, (2.45), (with updated domains
$\Omega_{i}$) are solved using a second order finite difference scheme. A
small amount of Gaussian smoothing is applied to the interface at every time-
step to aid stability, but this is done under the constraint of mass
conservation.

In figure 2 the gravity current shape $h(r,t)$ is shown together with
streamlines at different times. To facilitate displaying the data on similar
scales, we stretch the horizontal and vertical coordinates using the early and
late scalings. At early times $t\ll t^{*}$, the flow leaves the origin in
straight lines, and the shape of the current is given by (2.8). Later, when
$t\approx t^{*}$, the flow enters a transition regime in which gravity becomes
appreciable in the majority of the current, except in the near-origin
pressure-driven region, whose relative radial extent $R_{p}(t)/R(t)$
diminishes in size like $\sim t^{-2/7}$. Much later, when $t\gg t^{*}$, the
shape of the gravity current converges to the approximate solution (2.29), and
the streamlines away from the origin become horizontal. We also plot the
horizontal and vertical extents $R$, $H$, in figure 4. At early
($t/t^{*}<10^{-2}$) and late ($t/t^{*}>10^{2}$) times, the numerical solution
matches the analytical scalings (2.6)-(2.7),(2.41)-(2.42) closely. At very
late times, the slow logarithmic growth (2.44) can be observed as a slight
deviation in the vertical extent $H$ away from the late time scaling (2.41),
which is consistent with our predictions.

To understand the role of the inner pressure-driven region, we have also
calculated the maximum absolute value of the pressure deviation from
hydrostatic $|p-p_{\mathrm{hyd}}|_{\infty}$, normalised by the weight of the
injected fluid $\Delta\rho gH$. This deviation does not fall to zero over
time, but instead appears to tend to a constant value of around
$|p-p_{\mathrm{hyd}}|_{\infty}/\Delta\rho gH\approx 3$. This indicates that
the flux is not well approximated by the gravity-driven component (2.30)
everywhere within the current. As a result, the classical thin-film equation
(2.13) and the similarity variables (2.17) cannot be accurate in the late-time
regime. In particular, the boundary condition at the origin feeding the
gravity current with flux $Q$ cannot be supplied with a gravity-driven term of
the form (2.30) as this would cause an infinite thickness. Instead, the only
way for the thickness to remain finite is with a pressure-driven input flux at
the origin, which explains why the pressure deviation from hydrostatic may
never drop to zero over time.

We have also performed the same numerical calculation in the case of a two-
dimensional injection (from a line-source) and in this case the pressure
deviation $|p-p_{\mathrm{hyd}}|_{\infty}/\Delta\rho gH$ does drop to zero over
time, indicating that (by contrast) the two-dimensional flux is well
approximated by the corresponding gravity-driven component ($-u_{b}h\partial
h/\partial x$) everywhere within the current at late times.

Since the power laws $t^{3/7}$ and $t^{1/2}$ are not hugely different, we have
also tried comparing the numerical solution with the classical self-similar
scalings (2.17) proposed by Lyle et al. (2005), as shown with dotted green
lines in figure 4. Since the disparity between the different proposed scalings
grows larger over time like $\sim t^{1/2-3/7}\sim t^{1/14}$ (in the case of
radius) and like $\sim t^{0-1/7}\sim t^{-1/7}$ (in the case of thickness),
after a time scale $t/t^{*}=10^{6}$ this corresponds to disparaging factors of
$10^{6/14}\approx 2.68$ and $10^{-6/7}\approx 0.14$, respectively. Since these
discrepancies lie outside the 5$\%$ error observed in figure 4, this indicates
that the scalings (2.41)-(2.42) are more accurate than the self-similar
scalings (2.17).

###  2.6 Comparison with experimental data

$Q$ (cm3/s) | 0.020 | 0.021 | 0.043 | 0.125 | 0.209 | 17.36 | 17.67 | 31.62
---|---|---|---|---|---|---|---|---
$g^{\prime}$ (cm/s2) | 172.66 | 178.54 | 159.90 | 159.90 | 151.07 | 22.56 | 22.56 | 22.56
Table 1: List of experiments and corresponding parameter values.

Next we compare our results with experiments conducted in a porous media tank.
We use the same apparatus for our experiments as Lyle et al. (2005) except
that we use a peristaltic pump instead of a fluid reservoir to control the
flow rate (for improved accuracy). A square perspex tank with a
$61\,\mathrm{cm}\times 61\,\mathrm{cm}$ base is filled with 0.3 cm glass
Ballotini beads saturated with fresh water. Due to the difficulties associated
with creating an anisotropic permeability field from these beads, our
experiments are restricted to the isotropic case. Hence, the dependence on the
anisotropy parameter $K$ in the analytical scalings (2.6)-(2.7),(2.41)-(2.43)
is checked against our numerical simulations, whilst the dependence on the
remaining parameters (e.g. dimensional scalings and power laws) is checked
against both numerical simulations and experiments (in the isotropic case).

Two cameras are positioned above and alongside the tank, pointing vertically
downwards (plan view) and horizontally (side view), to capture both the radius
and the thickness of the current. Salty water dyed with red food colouring is
injected into the corner of the tank at different salt concentrations leading
to a variety of density contrasts. By injecting into one corner of the tank
the experiment represents one quarter of the full axisymmetric scenario whilst
allowing for a larger range of radial data given the tank size.

A full list of the different input flow rates and reduced gravity values
$g^{\prime}=\Delta\rho g/\rho_{w}$ (where $\rho_{w}=1\,\mathrm{g/cm^{3}}$ is
the density of water) is given in Table 1. The experiments at high flow rates
$Q=17.36-31.62\,\mathrm{cm^{3}/s}$ result in a current with an order
$\mathcal{O}(1)$ aspect ratio, with most of the injected fluid located in the
interior region of the tank (away from the walls). Hence, the porosity and
permeability are estimated using the values calculated for randomly packed
beads, $\phi=0.37$ and $k=6.8\times 10^{-9}\,\mathrm{m}^{2}$ (Lyle et al.,
2005). However, the experiments conducted at low flow rates
$Q=0.020-0.209\,\mathrm{cm^{3}/s}$ result in a gravity current no thicker than
1-3 bead diameters ($z=d-3d$) and hence the flow is subject to wall effects.
In these cases the porosity is estimated using the empirical formula derived
by Ribeiro et al. (2010) accounting for wall effects, which is

$\phi=0.373+0.917e^{-0.824\,z/d}\approx 0.55,$ (2.47)

where we have used an average flow depth value $z=2d$, and the permeability is
estimated using the Kozeny-Carman relationship

$k=\frac{d^{2}\phi^{3}}{180(1-\phi)^{2}}\approx 4.36\times
10^{-8}\,\mathrm{m}^{2}.$ (2.48)

It is worth noting how the shear in permeability due to wall effects (2.47)
may modify the flow. As discussed by Hinton & Woods (2018), whilst such shear
effects may alter the coefficients of the spreading rates, they do not modify
the power laws, which motivates the use of different constant values here.

By varying the flow rate and reduced gravity, a variety of different
transition time scalings are achieved, producing data in the range
$t/t^{*}=10^{-3}-10^{5}$. We also complement our data with the experiments of
Lyle et al. (2005) which are in the intermediate range $t/t^{*}=0.08-22$. To
reach such late time values we inject using a very slow flow rate of
$Q=0.02\,\mathrm{cm^{3}/s}$ over one hour with photos taken at logarithmically
spaced time intervals starting at 30 s (see figure 5). After post-processing
the images, the radial and vertical extents of the current are extracted by
dividing the pixels into binary values according to a threshold value (enabled
by the colour contrast due to the dye).

$\boldsymbol{H=H^{*}}$$\boldsymbol{t/t^{*}=119}$$\boldsymbol{t/t^{*}=1908}$$\boldsymbol{t/t^{*}=30532}$$\boldsymbol{t/t^{*}=1540}$$\boldsymbol{t/t^{*}=12316}$$\boldsymbol{t/t^{*}=197060}$$\boldsymbol{Q=0.209$cm$^{3}$/s,$H^{*}=0.36\,\mathrm{cm}}$$\boldsymbol{Q=0.020$cm$^{3}$/s,$H^{*}=0.11\,\mathrm{cm}}$(a)(b)
Figure 5: (a,b) Closeup snapshots of the gravity current thickness near the
origin during two injection experiments with different flow rates. Photos are
taken over 30 minutes (a) and 1 hour (b), and spaced logarithmically in time.

Data for the thickness and radius of the current are plotted together with
numerical and analytical results in figure 4. The radius data $R(t)$
approximately follows the $3/7$ power law (2.42), however one could argue that
the $1/2$ power law (2.17) is just as accurate within the error margins of the
experimental data. More clarity is achieved by observing the variation of the
thickness near the origin $H(t)$, as shown by the photos in figure 5. In
particular, the thickness is clearly observed to increase beyond the
transition value by as much as $H\approx 5.5H^{*}$ over the course of the
experiments, albeit with a very small growth rate. Indeed, without running the
experiments for such a long time it would appear that the thickness has ceased
growing entirely. By contrast, when observed at logarithmically spaced time
intervals the thickness shows no signs of abating. Post-processing of the
images to extract the interface position reveals that the thickness follows
the 1/7 power law (2.41) to good approximation, as shown in figure 4. We have
also tried using a least-squares minimisation that fits an arbitrary power law
$H=At^{B}$ to the combined experimental dataset (over four different flow
rates) at late times $t>100t^{*}$, as shown in the figure insert. This
produces a power law of $B=0.1537\pm 0.0176$ which is close to $1/7\approx
0.1429$. Moreover, this slightly larger numerical value for the power law
might hint towards the faster logarithmic growth predicted earlier.

It is not possible to discern the logarithmic correction to the thickness
(2.44) quantitatively since this would require much longer time scales than we
could achieve experimentally. Nor is it possible to accurately discern the
convergence of the shape of the current to the profile (2.29) over long times,
since the aspect ratio of the current is incredibly slender (around 200:1). To
extend the experiments to larger values of $t/t^{*}$ is challenging in
practice, since smaller flow rates than $Q=0.02\,\mathrm{cm^{3}/s}$ would
result in current thicknesses much smaller than a single bead size ($H^{*}\ll
0.3\,\mathrm{cm}$) that would be difficult to observe. Alternately one could
run the experiments for a longer time, but this requires a much larger tank
than we have available. For example to extend $t/t^{*}$ by a factor 100 would
require a tank of dimensions $100^{3/7}\approx 7.2$ times larger (i.e. with a
$439\,\mathrm{cm}\times 439\,\mathrm{cm}$ base). Nevertheless, our experiments
clearly show an increasing thickness near the origin with the approximate 1/7
power law behaviour (2.41). Hence, by conservation of mass (2.24) it follows
that the radius must grow according to a 3/7 power law, and the shape profile
(2.29) must become a good approximation over long times (e.g. see analysis in
Section 2.4).

It is important to note the possible effects of diffusion/dispersion over such
long time scales. The length scales for molecular diffusion, fluid dispersion,
and Taylor dispersion are approximately $\ell\sim(Dt)^{1/2}$,
$\ell\sim(dUt)^{1/2}$ and $\ell\sim(d^{2}U^{2}t/48D)^{1/2}$, respectively,
where $D$ is the molecular diffusivity of salt or dye (both similar values),
$U$ is a characteristic velocity scale, and $d=0.3\,\mathrm{cm}$ is the
diameter of the Ballotini beads. As an estimate for the diffusivity we take
$D=10^{-10}\,\mathrm{m^{2}/s}$, and for the velocity we take the value at the
top of the current $U\sim\mathrm{d}H/\mathrm{d}t$. Inputting these, we
calculate diffusion/dispersion scales of $\ell\sim 0.06\,\mathrm{cm}$,
$\ell\sim 0.13\,\mathrm{cm}$ and $\ell\sim 0.04\,\mathrm{cm}$, which are all
considered small within the experimental context. It is clear from the images
in figure 5 that the interface near the origin remains relatively sharp,
indicating that the effects of dispersion are negligible there. Hence, this
validates our approach of post-processing the images to extract the current
thickness at late times.

##  3 Relevance to carbon storage sites

In this section we situate our results within the context of carbon storage
applications, in which CO2 is injected into a porous reservoir saturated with
brine and bounded above by an impermeable cap rock (note that under the
Boussinesq approximation our analysis equally applies to a lighter fluid
injected into a heavier fluid which is bounded above). In particular, the
ability to determine whether an injected CO2 plume is in a pressure-driven
regime or a gravity-driven regime is useful for determining the expected
shape, and the expected plume migration speeds. This is important for ensuring
that the CO2 can be stored as safely and efficiently as possible.

For this analysis we compare our estimates for the transition time and
thickness scalings $t^{*}$, $H^{*}$, to typical parameter values from field
sites. This enables us to estimate whether the required injection time, or
whether the confines of the aquifer (e.g. the available space between
impermeable cap rocks) are sufficient for a gravity current to form. To make
this comparison we require approximate parameter ranges for the injection flux
$Q$, as well as the buoyancy velocity $u_{b}=k_{H}\Delta\rho g/\mu$, porosity
$\phi$ and anisotropy $K=k_{H}/k_{V}$.

| Regime | Variable | Prefactor | Scaling | $t^{n}$ | $K^{n}$
---|---|---|---|---|---|---
| $t\ll t_{x}^{*}$ | $R_{x}$ | $2/\pi^{1/2}$ | $(Q_{x}/\phi)^{1/2}$ | 1/2 | 1/4
| $t\ll t_{x}^{*}$ | $H_{x}$ | $2/\pi^{1/2}$ | $(Q_{x}/\phi)^{1/2}$ | 1/2 | $-1/4$
Planar | $t\gg t_{x}^{*}$ | $R_{x}$ | $\eta_{N}\approx 1.482$ | $(Q_{x}u_{b}/\phi^{2})^{1/3}$ | 2/3 | $\sim$
| $t\gg t_{x}^{*}$ | $H_{x}$ | $f_{0}\approx 1.296$ | $(Q_{x}^{2}/u_{b}\phi)^{1/3}$ | 1/3 | $\sim$
|  | $t_{x}^{*}$ | $\pi/4$ | $\phi Q_{x}/u_{b}^{2}$ | $\sim$ | $3/2$
|  | $H_{x}^{*}$ | $1$ | $Q_{x}/u_{b}$ | $\sim$ | $1/2$
| $t\ll t^{*}$ | $R$ | $(3/2\pi)^{1/3}$ | $(Q/\phi)^{1/3}$ | 1/3 | 1/6
| $t\ll t^{*}$ | $H$ | $(3/2\pi)^{1/3}$ | $(Q/\phi)^{1/3}$ | 1/3 | -1/3
Radial | $t\gg t^{*}$ | $R$ | $\zeta_{N}=(2^{11/2}/\pi^{2})^{1/7}$ | $(Q^{2}u_{b}/\phi^{3})^{1/7}$ | 3/7 | 1/14
| $t\gg t^{*}$ | $H$ | $\xi_{0}=(2^{1/2}\pi^{3})^{-1/7}$ | $(Q^{3}/u_{b}^{2}\phi)^{1/7}$ | 1/7 | $-1/7$
|  | $t^{*}$ | $2\pi/3$ | $\phi(Q/u_{b}^{3})^{1/2}$ | $\sim$ | $1$
|  | $H^{*}$ | $1$ | $(Q/u_{b})^{1/2}$ | $\sim$ | $0$
Table 2: List of asymptotic limiting behaviours for the horizontal and
vertical extent of the flow, as well as the transition scalings, for the case
of planar injection (from a line source) and radial injection (from a point
source). For example, the late time extent in the radial case is given by
$R=\zeta_{N}(Q^{2}u_{b}/\phi^{3})^{1/7}t^{3/7}K^{1/14}$. Note, the vertical
extent in the case of radial injection $H$ has a slow logarithmic dependence
(2.44) which we emit here for simplicity.

As described by Huppert & Pegler (2022) in the isotropic case, the time to
transition to a gravity current may be significantly prolonged by small
permeability values. Here we show that the transition time (2.12) may be
further prolonged by the presence of anisotropy, such that even sedimentary
systems with large horizontal permeability values $k_{H}$ may still remain in
a pressure-driven regime for a long time (e.g. hundreds or thousands of years)
if the vertical permeability $k_{V}$ is significantly smaller ($K\gg 1$). This
has significant consequences for modelling approaches, which often directly
assume a gravity-driven plume within the timeframe of the injection site, and
indicates the need for detailed measurements of heterogeneities to ensure
accurate predictions for the migration of CO2.

In addition to the case of axisymmetric radial injection, which is the focus
of the current study, we also briefly describe the analogous case of two-
dimensional planar injection (from a line-source) into anisotropic media. This
is a simple extension of the study by Huppert & Pegler (2022) for isotropic
media, incorporating different permeabilities $k_{H}$, $k_{V}$, in the
horizontal and vertical directions. As such, for the present analysis we skip
the derivation of these scalings and simply present them in Table 2 for
reference. To distinguish the planar case from the radial (axisymmetric) case,
we introduce subscript notation for the transition time $t^{*}_{x}$, the
transition thickness $H^{*}_{x}$, and the injection flux per unit width
$Q_{x}$. Likewise, the variables $R_{x}$ and $H_{x}$ in the planar case are
equivalent to the lateral and vertical extents of the current. We note that
the dynamics of the current for planar injection only depend on the anisotropy
parameter $K$ at early times, not at late times, as discussed in the
Introduction.

The references for the parameter values used comprise a variety of sources
describing different carbon storage sites around the world (Bickle et al.,
2007; Chadwick et al., 2009; Oldenburg et al., 2011; Golding et al., 2013;
Cowton et al., 2016; Williams & Chadwick, 2017; Bickle et al., 2017; Cowton et
al., 2018). The buoyancy velocity is estimated using parameter values
$\mu=\left[5.5-6.6\right]\times 10^{-5}\,\mathrm{Pa}\cdot\,\mathrm{s}$,
$\Delta\rho=232-309\,\mathrm{kg/m}^{3}$, and
$k_{H}=10^{-14}-10^{-12}\,\mathrm{m}^{2}$, resulting in
$u_{b}=\left[3.4-551\right]\times 10^{-7}\,\mathrm{m/s}$. For the injection
flux, we take estimates for a two-dimensional line source as
$Q_{x}=\left[0.4-3\right]\times 10^{-4}\,\mathrm{m}^{2}/\,\mathrm{s}$, and for
a radial point source as $Q=\left[0.1-4\right]\times
10^{-2}\,\mathrm{m}^{3}/\,\mathrm{s}$. Porosity is taken as $\phi=0.2-0.25$.
The degree of anisotropy varies considerably between CO2 sequestration sites.
For example, observations of permeability variations were around two orders of
magnitude for Salt Creek, USA (Bickle et al., 2017), and four orders of
magnitude for the Tilje formation in Norway (a potential sequestration site)
(Martinius et al., 1999; Bergmo et al., 2017). Hence, for this study we
consider anisotropy in the range $K=1-10^{4}$.

(a)(b)In SalahSleipner Figure 6: Transition time $t^{*}$, $t_{x}^{*}$ (a) and
transition thickness $H^{*}$, $H_{x}^{*}$ (b) for different anisotropy values
$K$, in the case of planar (line-source) injection and radial (point-source)
injection. Corresponding data for the Sleipner and In Salah carbon
sequestration sites (assuming radial injection) are illustrated, together with
uncertainty estimates for parameter values.

Collating all of this data, we plot the transition time $t^{*}$ and transition
thickness $H^{*}$ in figure 6 for a range of anisotropy values, in both the
planar and radial cases. Upper and lower bounds are displayed which correspond
to the range of possible values for the injection flow rates $Q_{x}$, $Q$,
buoyancy velocity $u_{b}$, and porosity $\phi$. For comparison, we also
display data points for two CO2 sequestration sites: Sleipner, Norwegian North
Sea (in operation 1996 - present), and In Salah, Algeria (in operation 2004 -
2012). Taking both of these as point source (radial) injection sites, and
using approximate anisotropy values $K=10-100$, we calculate the transition
time scale as less than a day for Sleipner and 280-700 days for In Salah.
Likewise the transition thickness for Sleipner is around 2-9 m, as opposed to
170-440 m for In Salah. This indicates that the modern day plume at Sleipner
is very likely to be in a gravity-driven regime, whereas the CO2 flow at the
less permeable site at In Salah is likely to have been in a pressure-driven
regime for a significant period of operation. In particular, the large
transition thickness at In Salah indicates that interaction with other
confining stratigraphy (i.e. within 440 m of the injection point) was a
likelier dominant control on CO2 migration rather than gravity-driven
spreading.

##  4 Discussion

The transition from a pressure- to a gravity-driven flow has been addressed
for constant axisymmetric injection into an anisotropic porous medium. We have
derived scalings for the early and late time growth of the current, as well as
the timescale at which the transition occurs, showing close comparison with
analogue experiments and numerical simulations. Early time scalings revealed
that anisotropy causes a long and thin current in which the pressure is not
hydrostatic, with a delayed transition to the gravity-driven regime. Hence
this study presents a paradigm shift for such flows, which are typically
assumed to be gravity-driven if long and thin. Analysis of the late time
regime showed that a region near the origin must remain pressure-driven.
Within this pressure-driven region, strong vertical velocities are required to
deliver a uniform horizontal flow to the rest of the gravity-driven current.
The pressure within this region cannot become hydrostatic, since this causes
an unphysical singular current thickness near the origin. Therefore, the inner
pressure-driven region serves to counteract this singularity, thereby
providing a significant non-hydrostatic contribution to the flux even at late
times.

One consequence of this pressure component is a finite but ever-growing
thickness near the origin (growing like $\sim(\log t)^{1/2}t^{1/7}$), which
contrasts previous theories that assumed a constant thickness scaling. Another
consequence is that the late time growth of the gravity current remains
affected by the vertical permeability value $k_{V}$ (and hence the anisotropy
ratio $K$) due to vertical velocities near the origin. This is in contrast to
the two-dimensional (planar) case, in which anisotropy only affects the time
to transition, not the late time behaviour. Our theoretical predictions were
confirmed by comparison with numerical simulations of anisotropic Darcy flow
run across 12 orders of magnitude in time, and porous media tank experiments
in the isotropic case. Transition timescales were calculated for anisotropy
values of carbon sequestration applications, indicating that some field sites
may never reach the gravity-driven regime over the course of operation.

It should be noted that for some CO2 sequestration sites injection occurs over
a vertical interval of around $1-10\,\mathrm{m}$, such that the flow is
delivered uniformly across the current depth at source. In such cases the
vertical velocities associated with the inner pressure-driven region described
in this study will not emerge unless the thickness scaling
$H^{*}=(Q/u_{b})^{1/2}$ is larger than this interval. Consequently, in such
cases the classical dynamics of self-similarity with a radial extent growing
like $R\sim t^{1/2}$ described by Lyle et al. (2005) remain accurate. Hence,
the results presented here are only relevant for point source injection sites,
or for sites with a vertical injection interval smaller than the transition
thickness $H^{*}$ (as is the case for anisotropic aquifers).

Whilst we have mentioned carbon sequestration applications throughout this
study, we have neglected several physical effects which are relevant to the
flow of CO2 in brine. For example, we have neglected the viscosity contrast
between these phases for simplicity (CO2 is typically around 20-30 times less
viscous than brine). By extending this study to account for such effects, it
is expected that the viscosity ratio $M=\mu_{1}/\mu_{2}$ will appear in the
early and late dynamics for $H$ and $R$ wherever the anisotropy parameter $K$
does. This is because the viscosity ratio can only affect the dynamics
whenever the flow of injected and ambient fluids are coupled together (just
like $K$). It would also be interesting to extend these results to account for
multiphase effects, similarly to Golding et al. (2013) in the case of
homogeneous media, and similarly to Benham et al. (2021) in the case of
heterogeneous media.

There are several interesting perspectives from this work that are worth
discussing. Firstly, the different dynamics between early and late times
present an opportunity for inverse modelling to help interpret petrophysical
information about the reservoir. For example, if it can be shown that the
radial extent of the flow is growing according to the early time power law
behaviour, then this can place bounds on the permeability and anisotropy of
the flow due to heterogeneities. This is useful for field sites since
information about the heterogeneities is often restricted to sparse core
measurements or coarse seismic surveys.

It is also worth noting the analogy between the current study and the
canonical case of a classical viscous gravity current. It is easy to follow
the above analysis for the case of constant injection of a viscous fluid onto
an impermeable substrate in either the planar or radial case. It transpires
that similar scalings can be derived for the time and thickness required to
transition from a pressure- to a gravity-driven flow. For example, in the
radial case these are given by $H^{*}=(Q\mu/\Delta\rho g)^{1/4}$ and
$t^{*}=(2\pi/3)(\mu^{3}/(\Delta\rho g)^{3}Q)^{1/4}$. However, it is not clear
how the inner pressure-driven region affects the viscous gravity current at
late times. In particular, the no-slip condition at the base of the current
produces a different flow pattern to that considered here, making it difficult
to translate our other results. Nevertheless, one could investigate such flow
patterns by introducing tracers at the source of the gravity current, shedding
light on the size and role of the pressure-driven region at late times.

###### Acknowledgements.

This research is funded in part by the GeoCquest consortium, a BHP-supported
collaborative project between Cambridge, Stanford and Melbourne Universities,
and by a NERC consortium grant “Migration of CO2 through North Sea Geological
Carbon Storage Sites” (grant no. NE/N016084/1).
Declaration of Interests. The authors report no conflict of interest.

## References

  * Barenblatt (1996) Barenblatt, GI 1996 Scaling, self-similarity, and intermediate asymptotics. Cambridge University Press.
  * Benham et al. (2021) Benham, GP, Bickle, MJ & Neufeld, JA 2021 Two-phase gravity currents in layered porous media. J. Fluid Mech. 922, A7.
  * Bergmo et al. (2017) Bergmo, PES, Emmel, BU, Anthonsen, KL, Aagaard, P, Mortensen, GM & Sundal, A 2017 Quality ranking of the best CO2 storage aquifers in the Nordic countries. Energy Procedia 114, 4374–4381.
  * Bickle et al. (2007) Bickle, M, Chadwick, A, Huppert, HE, Hallworth, M & Lyle, S 2007 Modelling carbon dioxide accumulation at Sleipner: Implications for underground carbon storage. Earth Planet. Sci. Lett. 255 (1-2), 164–176.
  * Bickle et al. (2017) Bickle, M, Kampman, N, Chapman, H, Ballentine, C, Dubacq, B, Galy, A, Sirikitputtisak, T, Warr, O, Wigley, M & Zhou, Z 2017 Rapid reactions between CO2, brine and silicate minerals during geological carbon storage: Modelling based on a field CO2 injection experiment. Chem. Geol. 468, 17–31.
  * Cardwell & Parsons (1945) Cardwell, WT & Parsons, RL 1945 Average permeabilities of heterogeneous oil sands. Trans. AIME 160 (01), 34–42.
  * Chadwick et al. (2009) Chadwick, RA, Noy, DJ & Holloway, S 2009 Flow processes and pressure evolution in aquifers during the injection of supercritical CO2 as a greenhouse gas mitigation measure. Pet. Geosci. 15 (1), 59–73.
  * Corbett & Jensen (1992) Corbett, PWM & Jensen, JL 1992 Variation of reservoir statistics according to sample spacing and measurement type for some intervals in the lower brent group. The Log Analyst 33 (01).
  * Cowton et al. (2016) Cowton, LR, Neufeld, JA, White, NJ, Bickle, MJ, White, JC & Chadwick, RA 2016 An inverse method for estimating thickness and volume with time of a thin CO2-filled layer at the Sleipner field, North Sea. J. Geophys. Res. Solid Earth 121 (7), 5068–5085.
  * Cowton et al. (2018) Cowton, LR, Neufeld, JA, White, NJ, Bickle, MJ, Williams, GA, White, JC & Chadwick, RA 2018 Benchmarking of vertically-integrated CO2 flow simulations at the Sleipner field, North Sea. Earth Planet. Sci. Lett. 491, 121–133.
  * Golding et al. (2013) Golding, MJ, Huppert, HE & Neufeld, JA 2013 The effects of capillary forces on the axisymmetric propagation of two-phase, constant-flux gravity currents in porous media. Phys. Fluids 25 (3), 036602.
  * Hinton & Woods (2018) Hinton, EM & Woods, AW 2018 Buoyancy-driven flow in a confined aquifer with a vertical gradient of permeability. J. Fluid Mech. 848, 411–429.
  * Huppert & Pegler (2022) Huppert, HE & Pegler, SS 2022 The fate of continuous input of relatively heavy fluid at the base of a porous medium. J. Fluid Mech. 932.
  * Huppert & Woods (1995) Huppert, HE & Woods, AW 1995 Gravity-driven flows in porous layers. J. Fluid Mech. 292, 55–69.
  * Kumar et al. (1997) Kumar, A, Farmer, CL, Jerauld, GR & Li, D 1997 Efficient upscaling from cores to simulation models. In SPE Annual Technical Conference and Exhibition. OnePetro.
  * Lyle et al. (2005) Lyle, S, Huppert, HE, Hallworth, M, Bickle, M & Chadwick, A 2005 Axisymmetric gravity currents in a porous medium. J. Fluid Mech. 543, 293–302.
  * Martinius et al. (1999) Martinius, AW, Ringrose, PS, Næss, A, Wen, R & Hentz, TF 1999 Multi-scale characterisation and modelling of heterolithic tidal systems, offshore mid-Norway. In Advanced Reservoir Characterization for the 21st Century. Proceedings of GCSSEPM Research Conference, pp. 193–204.
  * Oldenburg et al. (2011) Oldenburg, CM, Jordan, PD, Nicot, JP, Mazzoldi, A, Gupta, AK & Bryant, SL 2011 Leakage risk assessment of the In Salah CO2 storage project: Applying the certification framework in a dynamic context. Energy Procedia 4, 4154–4161.
  * Ribeiro et al. (2010) Ribeiro, AM, Neto, P & Pinho, C 2010 Mean porosity and pressure drop measurements in packed beds of monosized spheres: side wall effects. Int. Rev. Chem. Eng. 2 (1), 40–46.
  * Whitham (2011) Whitham, GB 2011 Linear and nonlinear waves, , vol. 42. John Wiley & Sons.
  * Williams & Chadwick (2017) Williams, GA & Chadwick, RA 2017 An improved history-match for layer spreading within the Sleipner plume including thermal propagation effects. Energy Procedia 114, 2856–2870.
  * Woods (2015) Woods, AW 2015 Flow in porous rocks. Cambridge University Press.