Axisymmetric gravity currents in anisotropic porous media

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 CO$_2$ migration through permeable aquifers, as occurs in the context of carbon capture and storage.


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/τ ) 1/3 f (x/L)(t/τ ) −2/3 , (1.1) where L = Q x µ/k∆ρg and τ = Q x (µ/k∆ρ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, ∆ρ, and the viscosity of the injected fluid, µ (Huppert & Woods 1995). If one extends the analysis to account † Email address for correspondence: graham.benham@ladhyx.polytechnique.fr arXiv:2202.01165v2 [physics.flu-dyn] 3 Feb 2022 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 For the solution to be valid we require that (1.2) is small compared to the corresponding hydrostatic pressure gradient This requires that 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 CO 2 sequestration and illustrate how anisotropy of the porous medium leads to the ratio k V /k H influencing the dispersal of the CO 2 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.

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, k = 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 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 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 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 = O(10 4 ).

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 earlytime 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 φ, 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 ρ 2 < ρ 1 , whereas the viscosities of the two fluids are assumed to be equal µ 1 = µ 2 = µ. 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 0, z 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.
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 nonuniform coefficients, such that 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 Likewise, conservation of mass (within a stretched hemisphere) dictates that 2πφHR 2 /3 = Qt. (2.5) This indicates that the early time dynamics are described by 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 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 ∂p/∂r ≈ ∆ρg∂h/∂r, so that the integrated gravity-driven flux scales like (2.9) By inserting approximate scalings, h ∼ H, r ∼ R, ∂h/∂r ∼ −H/R, this provides the balance where the buoyancy velocity u b = k H ∆ρg/µ. This immediately provides an expression for the thickness of the current at the transition between pressure-to buoyancy-driven flow, (2.11) and by equating this to the early time behaviour (2.6) we find that the transition time is 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 = O(10 4 ) (Martinius et al. 1999;Bergmo et al. 2017). The critical time t * distinguishes two distinct regimes. At early times t t * , or when H H * , the effects of gravity are negligible and injection pressures dominate the flow, whereas at late times t t * , or when H 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.

Late-time dynamics: classical analysis
At much later times t t * , once the flow has grown to a thickness H 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 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 (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.
The system of equations (2.13)-(2.16) admits the similarity variables where the dimensionless shape f and prefactor η N satisfy The solution can be calculated numerically, giving a prefactor value η 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.
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).

Back to the future: late-time dynamics revisited
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 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 where G is the dimensionless shape function, ξ 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 ∼ t 1/2 power law described earlier. In particular, conservation of mass (in conjunction with (2.22)) dictates that (2.24) If a > 0, it follows that R grows as ∼ t b , where b < 1/2. Hence, we introduce a corresponding set of rescaled outer variables where F is the dimensionless shape profile, the nose position ζ 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 (hydro-static) region becomes ( 2.27) The second two terms on the left hand side of (2.27) are of the order O(t −2a ) and so become vanishingly small at large times. Assuming the time derivative of the scaled shape decreases as ∂F/∂t ∼ 1/t for large times, then the first term on the left hand side is also of the order O(t −2a ). Hence, after sufficiently late times (2.27) yields the result where 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 We note that the integrated gravity-driven flux, which is given by grows according to 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 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πφˆR 0 rh dr = Qt. (2.32) Hence, by inserting (2.29) into (2.32), we arrive at a relationship involving ζ N and ξ 0 , which is ). Noting that becomes vanishingly small at late times (since b > a), and noting that´1 0 y(− log y) 1/2 dy = (π/32) 1/2 , we can rearrange (2.33) to give A further final condition is required to fully determine the coefficients ζ N , ξ 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 In addition, the thickness and radial scalings are given by 42) 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 gravitycurrent, 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/π) log R p /R) 1/2 ∼ (2/(7π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πK/2 ≈ (5.96 × 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.

Comparison with numerical simulations
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 where ρ = ρ 1 inside the injected fluid domain (which we denote Ω 1 ) and ρ = ρ 2 outside the injected fluid domain (which we denote Ω 2 ). The anisotropic permeability matrix  The numerical solution at r = 0 is not available due to the singular nature of the divergence in the continuity equation ∇ · 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 Ω 1 and Ω 2 is given by a set of points ∂Ω = r(t), which is treated as a passive scalar and is therefore advected at the fluid velocity, such that dr dt = 1 φ u(r), (2.46) with initial conditions given by a small surface r(0) = 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 Ω 1 and Ω 2 are updated via interpolation over a gridded mesh of 150×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 Ω 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 t * , the flow leaves the origin in straight lines, and the shape of the current is given by (2.8). Later, when t ≈ 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 ∼ t −2/7 . Much later, when t 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 hyd | ∞ , normalised by the weight of the injected fluid ∆ρgH. This deviation does not fall to zero over time, but instead appears to tend to a constant value of around |p − p hyd | ∞ /∆ρgH ≈ 3. This indicates that the flux is not well approximated by the gravitydriven 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 hyd | ∞ /∆ρ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∂h/∂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 ∼ t 1/2−3/7 ∼ t 1/14 (in the case of radius) and like ∼ t 0−1/7 ∼ 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 ≈ 2.68 and 10 −6/7 ≈ 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).

Comparison with experimental data
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 cm × 61 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 = ∆ρg/ρ w (where ρ w = 1 g/cm 3 is the density of water) is given in Table 1 walls). Hence, the porosity and permeability are estimated using the values calculated for randomly packed beads, φ = 0.37 and k = 6.8 × 10 −9 m 2 (Lyle et al. 2005). However, the experiments conducted at low flow rates Q = 0.020 − 0.209 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 where we have used an average flow depth value z = 2d, and the permeability is estimated using the Kozeny-Carman relationship (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 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).
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 ≈ 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 ± 0.0176 which is close to 1/7 ≈ 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 cm 3 /s would result in current thicknesses much smaller than a single bead size (H * 0.3 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 ≈ 7.2 times larger (i.e. with a 439 cm × 439 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 ∼ (Dt) 1/2 , ∼ (dU t) 1/2 and ∼ (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 cm is the diameter of the Ballotini beads. As an estimate for the diffusivity we take D = 10 −10 m 2 /s, and for the velocity we take the value at the top of the current U ∼ dH/dt. Inputting these, we calculate diffusion/dispersion scales of ∼ 0.06 cm, ∼ 0.13 cm and ∼ 0.04 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.

Relevance to carbon storage sites
In this section we situate our results within the context of carbon storage applications, in which CO 2 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 CO 2 plume is in a pressuredriven 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 CO 2 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 Regime Variable 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 = ζ N (Q 2 u b /φ 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. 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 ∆ρg/µ, porosity φ and anisotropy K = k H /k V .
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 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 CO 2 .
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 In Salah Sleipner 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 (pointsource) injection. Corresponding data for the Sleipner and In Salah carbon sequestration sites (assuming radial injection) are illustrated, together with uncertainty estimates for parameter values.
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.
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 φ. For comparison, we also display data points for two CO 2 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 CO 2 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 CO 2 migration rather than gravitydriven spreading.

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 ∼ (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 CO 2 sequestration sites injection occurs over a vertical interval of around 1 − 10 m, such that the flow is delivered uniformly across the current depth at source. In such cases the vertical velocities associated with the inner pressuredriven 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 ∼ 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 CO 2 in brine. For example, we have neglected the viscosity contrast between these phases for simplicity (CO 2 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 = µ 1 /µ 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µ/∆ρg) 1/4 and t * = (2π/3)(µ 3 /(∆ρ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.