Anisotropy distorts the spreading of a fixed volume porous gravity current

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


Introduction
Gravity-driven flows resulting from the release of a fluid within a porous medium are a common feature of environmental fluid dynamics.For example, such flows arise when groundwater responds to heavy rainfall [1], after the spillage of a contaminant [2], or during the geological storage of carbon dioxide in saline aquifers [3].Since all geological aquifers are heterogeneous and often this heterogeneity manifests as an anisotropic permeability field [4], it is important to quantify how this affects the migration speeds and shape of the current.This is particularly relevant to situations where it is desirable to minimise the volume of the aquifer contacted by the fluid (e.g. containing the spread of a contaminant) or to maximise the contacted volume (e.g.trapping residual saturation during CO 2 storage).
In the case of homogeneous and isotropic porous media, there have been numerous studies on the evolution of fixed volume gravity currents.Some of these studies have treated singlephase flows [5], whilst others have incorporated multiphase effects such as residual trapping [6,7] and dissolution within the ambient fluid [8].In some cases simple scaling laws were derived by exploiting the self-similar properties of gravity currents, as shown by [9] in the case of singlephase flow, and by [10] in the case of trapped residual saturation.Later work by [11][12][13] explored the transition to self-similarity in both confined and unconfined settings.However, less attention has been paid to the case of heterogeneous or anisotropic porous media, despite the relevance to real geological reservoirs.
Nevertheless, some progress has been made for specific types of buoyancy driven flows in heterogeneous media.For example, several studies have investigated how anisotropy affects convective dissolution within an ambient fluid phase [14][15][16].In the case of gravity currents resulting from constant injection, [17] explored how heterogeneities of different lateral and vertical scales affect the migration speed of a CO 2 plume.Likewise, [18] addressed the case of a gravity current resulting from point source injection in an anisotropic medium.In this study, it was shown that anisotropy can cause a build-up of pressure that stretches the flow into an ellipsoid shape during an early-time regime of the flow, before transitioning to a gravitydominated regime at much later times.However, no studies have addressed how heterogeneity affects the spreading of a released volume of fluid (i.e. in the absence of injection), despite the relevance to post-injection scenarios during CO 2 storage, and to post-leakage scenarios in the context of contaminant spills.
CO 2 storage in geological reservoirs is one of the key proposed technologies to reduce emissions and limit the effects of global warming [3].In such scenarios, buoyant CO 2 is injected into a brine-filled reservoir beneath an impermeable cap rock.Once the injection is switched off, the CO 2 rises and spreads out beneath the cap rock, with a fraction of its mass being lost to residual trapping (via the drainage/imbibition cycle) and dissolution within the surrounding brine [8,19].Hence, to quantify the trapping potential of different geological reservoirs (e.g. when choosing potential storage sites) it is important to understand how the anisotropy of the aquifer may affect the historical migration of the current across the pore space.In the context of a contaminant spillage, the objective is to contain and minimise the spread of a harmful fluid within an aquifer.Therefore, in a similar manner to the CO 2 storage problem, it is necessary to quantify how and where the contaminant fluid will spread in response to the heterogeneity of the aquifer, once the leak has been closed off.
In this study we demonstrate that anisotropy restricts the flow of the gravity current to thin finger-like regions spreading near the impermeable boundary, qualitatively similar to those predicted by other studies [17].Due to this flow distortion, the swept volume of the current is reduced for anisotropic aquifers.This indicates that isotropic aquifers may have better potential for certain forms of CO 2 trapping that depend on the contacted volume of pore space.By contrast, in the case of a contaminant spill, anisotropic aquifers may help contain the spread of the fluid by restricting the flow to a reduced fraction of the porous medium.

Finite release in two-dimensional anisotropic media (a) Release and subsequent dynamics
We consider an anisotropic porous medium in which the horizontal permeability k H is much larger than the vertical permeability k V .Such anisotropic flow properties are a common feature in geological reservoirs and may result from the deposition of successive layers of fine and coarse material or from post-depositional compaction of the formation [20].As we will show, the anisotropy of the medium restricts the vertical flow of the bulk of the fluid, which we denote region I (see figure 1), resulting in a slow migration towards the impermeable boundary.Meanwhile, gravity-driven spreading is limited to thin finger-like regions near the base, which we denote region II .
To start with we restrict our attention to two-dimensional flows (although radially symmetric flows will be addressed later in Section 5) and we consider the release of a volume (per unit width) of fluid with constant density ρ.The surrounding porous medium is initially saturated with an ambient fluid with relatively smaller density ρa < ρ.Due to the Boussinesq approximation [21], these results also apply in the case of a lighter fluid (e.g.CO 2 ) released within a porous medium saturated with a heavier fluid (e.g.brine), with the impermeable boundary located above rather than below.
For the sake of simplicity we take the viscosity of the released fluid and the ambient fluid to be the same (µ), and we consider that the initial shape of the released fluid is rectangular, with dimensions 2L × H.It should be noted, however, that these results would apply to any similar convex shape, as shown in Figure 7 and discussed in more detail in Appendix A.
The flow in the released volume of fluid is subject to the two-dimensional Darcy equations, where u is the Darcy velocity vector, p is the pressure, and k = diag(k H , k V ) is the anisotropic permeability field in the x, z directions.Combining (2.1) and (2.2), the pressure satisfies Laplace's equation with anisotropic coefficients, where the anisotropy is given by (2.4) The boundary conditions for the flow in region I are as follows.The left hand and bottom boundaries, x = 0, z = 0, are assumed to be symmetric and impermeable, respectively, so we prescribe no normal flow, u = 0 : x = 0, (2.5) Likewise, at the fluid interface, which we denote z = h(x, t), we impose the dynamic and kinematic boundary conditions: ) The former condition matches the pressure in the fluid with the ambient hydrostatic pressure (note that the reference pressure pa is the ambient value at z = 0), whilst the latter condition imposes that a particle at the interface remains at the interface [22,23].
In the limit ǫ → 0, mass conservation (2.2) indicates that if there is no vertical flow (i.e.k V = 0) then there cannot be any horizontal flow either.Essentially, the governing equation (2.3) implies that the pressure (which is continuous) is set by the ambient fluid, such that p = pa − ρagz, and the resulting velocities are u = w = 0. Hence, the interface remains at the initial position z = h 0 (x) for all time.Now let's consider the case of a strongly anisotropic porous medium, such that 0 < ǫ ≪ 1.In this case, the solution can be found by performing an asymptotic expansion of the pressure in powers of ǫ.Within this expansion, the leading order contribution to the pressure is simply the solution to the isotropic problem (i.e. with ǫ = 0 exactly).Hence, for the same reasons as described above, the leading order pressure is hydrostatic and set by the ambient fluid, p = pa − ρagz.However, by inserting this pressure into Darcy's law (2.1),we now derive a small but finite (i.e.first order) vertical velocity within region I, such that where is the buoyancy velocity and ∆ρ = ρ − ρa.Since (2.9) does not satisfy the impermeability condition (2.6), it is necessary to reevaluate the solution near z ≈ 0 using boundary layer theory.Mathematically speaking, (2.3) is a singular perturbation problem since it is a second order Partial Differential Equation (PDE) with a small parameter in front of a second derivative.This indicates that not all vertical boundary conditions can be satisfied by the leading order solution.Specifically, the dynamic boundary condition (2.7) is imposed at z = h to ensure continuity of pressure, whilst the impermeability condition (2.6) at z = 0 is left unsatisfied.To correct this requires rescaling the solution to investigate changes over a small vertical distance near z ≈ 0, also known as a boundary layer.By inspection of (2.3) it is clear that z must be rescaled by a factor of ǫ 1/2 to recover all terms in the governing equation.Hence, an appropriate choice of rescaled dimensionless variables is where ξ, ζ and P are variables which are O(1) in magnitude.Note that the pressure in (2.11) must also be rescaled by a factor ǫ 1/2 so that the boundary condition (2.6) (i.e.∂p/∂z + ρg = 0 at z = 0) balances all terms at leading order.In this way, within the boundary layer region the governing equations and boundary conditions (2.3),(2.5)and (2.6) become ) ) In addition, we require that the inner solution (within the boundary layer) matches with the outer solution (far outside the boundary layer), such that (2.15) The system is not yet complete since the governing equation (2.12) is a second order elliptic PDE which requires four boundary conditions.The fourth and final boundary condition needs more careful thought.Since there is a vertical velocity (2.9) descending through the outer region, this induces an arrival of flux ǫu b L within the inner region.However, since this flux can go neither downwards nor leftwards (due to impermeable/symmetric boundaries), it must instead exit through the right hand boundary, x = L.In other words, the right hand interface must move outwards to conserve mass, creating a new finger-like region of vertical size z ∼ O(ǫ 1/2 ), which we denote region II.Hence, the final boundary condition for region I is given by an integral constraint of the form The descent of the upper interface in region I and the resulting finger-like region II are both illustrated in figure 1b.
Before addressing these details further, we first note that (2.12)-(2.16)can be solved exactly by separation of variables.Hence, the composite solution (valid across both inner and outer regions) is given by where λn = (2n + 1)π/2.Conservation of mass within region I indicates that Hence, inserting (2.17) into (2.18)results in the governing equation for the evolution of the thickness of the current, which is This can be further simplified by ignoring exponentially small terms when h/L is larger than O(ǫ 1/2 ), and by using the fact that the infinite sum converges to −x/L.Hence, we see that the thickness within region I is given by which is valid for 0 < x < L and for h ≫ O(ǫ 1/2 ).
Next, we address the fluid flow in region II, which is the finger-like region of escaped fluid near the base of the current, which is defined for L < x < x N (t), where x N (t) is the maximum extent of the finger.This flow region is long and thin (like a classical gravity current) such that Times are given in terms of the transition time (2.37),where (a) t/t * = 0.07, (b) 0.17 and (c) 0.34.Streamlines are evaluated using velocities derived from (2.17 the horizontal velocity is much larger than the vertical velocity.Consequently, the pressure within region II is hydrostatic to good approximation, such that Similarly to (2.18), conservation of mass within this region gives which is sometimes called the Dupuit approximation.This is accompanied by boundary conditions that correspond with imposing the input flux from region I, and imposing zero thickness and zero flux at the moving front, By introducing dimensionless coordinates, we get the same system as [5] for a constant input flux (see Appendix B for further details).The solution is well known and is given in terms of the similarity variables (2.27) The self-similar shape function f (η) is defined for η ∈ [0, η N = 1.482], and is monotone decreasing from f 0 := f (0) = 1.296 to f (η N ) = 0.The solution for regions I and II is plotted in figure 2a,b,c, at several different times.Streamlines confirm that the flux into region II is fed by the shrinking of region I.Comparison is also made to a numerical solution described later in Section 3. showing both numerical and approximate solutions.Early time behaviour is given by (2.20),(2.28),whereas late time behaviour is given by (2.34),(2.35).
The maximum vertical extent of the flow is given by (2.20), whereas the maximum horizontal extent is determined through the above scalings as (2.28) These are plotted in figure 3 with dashed blue lines, thereby indicating the early-time behaviour of the released fluid.

(b) Transition to self-similarity
After a long time the flow is expected to eventually transition to the self-similar behaviour of a finite release gravity current [5,12].As discussed by [18], the late-time dynamics of a twodimensional gravity current are independent of the anisotropy of the medium since the bulk flow decouples from the ambient.Hence, the anisotropy only affects the late-time behaviour by delaying the time to transition to self-similarity.Hence, to study the late-time behaviour we first analyse the isotropic case, which determines the late-time dynamics, and then use these dynamics to derive the time t * to transition between the early and late solutions, where t * depends on the anisotropy ǫ.
As before, the thickness of the gravity current (which now occupies a single region 0 ≤ x ≤ x N (t)) satisfies the Dupuit approximation (2.22).The boundary conditions are similar to the case of constant input flux, as described above, except the left hand boundary condition is replaced with the zero flux condition ∂h ∂x = 0 : x = 0. (2.29) Consequently, mass conservation indicates that (2.30) By introducing dimensionless coordinates, and switching to similarity variables we arrive at a system of equations that can be solved analytically to give where η N = 3 2/3 , as shown by [5] (see Appendix B for further details).
The maximum vertical and horizontal extent of the flow are given by respectively.These scalings are plotted in figure 3 with dashed blue lines, thereby indicating the late-time behaviour of the gravity current.Next, we discuss the time taken to transition from the early flow regime involving two fluid regions to the late flow regime with a single region which is self-similar.In the early regime the vertical extent of the flow descends according to (2.20).By contrast, the vertical extent of region II (ǫ 1/2 LH in (2.26)) increases like ∼ f 0 (u b ǫ 2 L 2 t/φ) 1/3 .Hence, it is expected that the transition to self-similarity will occur once these two thickness scalings approach each other.Thus, the transition time t * satisfies the cubic equation (2.36) The solution has a complicated form but can be expanded in powers of ǫ ≪ 1 to give where α = H/L is the aspect ratio of the initial current shape.Hence, anisotropy delays the transition to a classical self-similar regime, which is consistent with other studies [18].At the transition time t * the maximum thickness of the current (which we denote H * ) is given by

.38)
This indicates that, by the time transition occurs for very anisotropic media, the bulk of the current has shrunk significantly.This hints towards a reduced swept volume, which we analyse further in Section 4.
The transition time t * and the transition thickness H * (given in dimensionless terms) are plotted in figure 4.These importantly depend on both the anisotropy ǫ as well as the initial aspect ratio of the flow α.It should be noted that in practice strong anisotropy ǫ ≪ 1 may cause a stretching of the flow before contact with the impermeable cap rock.Therefore, ǫ ≪ 1 is likely to be correlated with α ≪ 1.

Finite difference computations of Darcy flow
Next, in this section we compare our analytical predictions to finite difference computations of two-dimensional Darcy flow.The flow is modelled with the Darcy equations (2.1)-(2.2),accounting for different permeability values, k H and k V , in the horizontal and vertical directions.The governing equations are accompanied by boundary conditions (2.5)-(2.6),corresponding with impermeable/symmetric walls at x = 0 and z = 0.The dynamic boundary condition (2.7) is imposed on the interface z = h(x, t), which is interpolated over a gridded mesh of 150×150 points.Likewise, a spatial domain of finite size 4L × H is chosen in the x × z directions.The fluid flow is not resolved beyond the interface z > h(x, t) since pressure is assumed to be hydrostatic  in the ambient fluid.Mass conservation dictates that the interface evolves according to with suitable initial conditions h(x, 0) = h 0 (x).Boundary conditions for h are given by (2.24), (2.25) and (2.29) (see earlier discussion for further explanation).The time-dependent equation for the fluid-fluid interface (3.1) is solved using an explicit forward Euler scheme in time, and a backward eighth order scheme in space.At each time-step the Darcy equations, (2.1)-(2.2),are solved using a second order central finite difference scheme.The code used for these computations is available in the Supplementary Materials.
Results from the finite difference computations are compared with the analytical model in figures 2 and 3 for an anisotropy value ǫ = 0.01. Figure 2a,b,c shows the distorted spreading of the gravity current via a thin finger near the base.Good agreement is observed everywhere except near x ≈ L where the two regions connect.Here, the interface transitions smoothly between regions I and II, which is a second order feature that is missing from the simple analytical model.The numerical solution enables computation of the gravity current shape up to and beyond the transition time t * , as is displayed in figure 2d.This demonstrates clearly how the flow transitions from a two-region (bulk/finger) structure at early times to a slumping single-region structure at late times.
Figure 3 displays numerical computations of the vertical and horizontal extents of the gravity current across early and late time regimes.Clearly, the numerical computations reflect the transition between the early and late time analytical scalings, shown with dashed lines.Likewise, the transition between these different regimes corresponds with t/t * ≈ 1, indicating the accuracy of our prediction for t * (2.37).
As a benchmark test, we also compare these finite difference computations with a gravity current slumping in an isotropic medium, ǫ = 1.Since our analytical model only applies for ǫ ≪ 1, the initial dynamics involving regions I and II are not relevant for this case.Instead, we set the initial shape as which is simply the similarity solution (2.33).Since the initial shape satisfies the similarity conditions, this ensures that the solution remains self-similar for all time.In figure 7a,b,c, in A, we plot the self-similar evolution of the gravity current shape at various times.Excellent agreement is attained between the numerical model and the exact self-similar solution, indicating the reliability of the finite difference approach.In figure 7d,e,f, we display similar computations for the same initial shape released in an anisotropic medium, ǫ = 0.01.In this case, the flow is decomposed into regions I and II, as before, demonstrating how our simple model can be extended to account for other released shapes.Further details and discussion of this case are given in Appendix A.

Swept shape and swept volume
As described earlier, it is useful for applications (e.g.contaminant transport or CO 2 storage) to quantify the total volume contacted by the gravity current, also known as the swept volume 1 .At early times t ≪ t * region I shrinks uniformly downwards whilst region II grows upwards and outwards.Hence, the swept volume V is simply equal to the initial volume plus the instantaneous volume of region II, such that At much later times t ≫ t * , once the gravity current has transitioned to self-similarity, the thickness z = h(x, t) has some parts which are growing and other parts which are shrinking, so this requires more careful attention.To deal with this we first define the swept shape S as the maximum thickness that the current ever reached at a given value of x, such that The swept volume is then given in terms of S as 3) The swept shape (4.2) is calculated by finding the time at which the thickness is maximal, which is equivalent to η = 3 1/6 in the similarity solution f (2.33).Inserting this into (4.2) we get which is only valid for x ≫ L. Hence, the swept volume (up to a constant C) is which clearly diverges like ∼ log t as t → ∞.The constant of integration is found by equating (4.1) and (4.5) at the transition time t * .Hence, the late time behaviour of the swept volume is given by In figure 5a,b,d, the swept shape S(x) and swept volume V(t) are plotted for different values of the anisotropy ǫ.For ǫ ≪ 1, the anisotropy restricts the swept shape (for x > L) to fingerlike regions near the base of the current.The maximum thickness of the fingers is set by H * (2.38), which scales like H * ∼ ǫ 1/3 (i.e. the more anisotropic the medium, the narrower the fingers).Furthermore, anisotropy delays the transition time (2.37), since t * ∝ ǫ −1 .Meanwhile, before transition occurs, t ≪ t * , anisotropy causes the swept volume to grow slowly (since V ∝ ǫt in (4.1)), such that the gravity current contacts a smaller fraction of pore space at early times.Note that the kinks at the corners x = ±L (in S) and at the transition time t = t * (in V) would be smoothed out by higher order asymptotics or numerical computations of two-dimensional Darcy flow.However, such features do not affect the overall leading order behaviour displayed here.
For comparison, we also plot the swept shape S in the isotropic case, ǫ = 1, in figure 5c.As described at the end of Section 3, the isotropic case uses (3.2) as the initial shape and remains self-similar for all times.It should be noted that the modified initial shape results in a different initial swept volume for this case, V(0) = 4HL/3.Clearly, the classical (isotropic) self-similar solution has a larger swept shape and a faster growing swept volume than the anisotropic cases.In Section 6, we discuss how these results can be applied when modelling contaminant spills and CO 2 sequestration, noting the possible limitations of our model.

Finite release with radial symmetry
The above results can be easily extended to account for radially symmetric flows if the anisotropy remains aligned with the vertical coordinate, i.e. with permeability k H , k V in the radial/vertical directions.To extend the model from the two-dimensional case, we take the impermeable boundary to be the horizontal plane, z = 0. Likewise, we consider the released shape to be a cylinder of initial height H and radius R. Since the initial shape of the current is radially symmetric, it will spread out and remain radially symmetric for all times.
Following the derivation in Section 2(a), the thickness of the current evolves according to (2.20) at early times.This results in a radial flux of magnitude ǫu b πR 2 exiting region I (0 ≤ r ≤ R) into a growing annular region II (R ≤ r ≤ r N (t)) at the base of the current.Conservation of mass within region II gives ( Boundary conditions correspond with imposing the input flux from region I, and imposing zero thickness and zero flux at the moving front, By introducing dimensionless coordinates, we get a similar system to [24] (see Appendix B for further details).The solution is given in terms of the similarity variables The shape function f (η) is defined for 0 < η ≤ η N = 1.155, but has an unphysical singularity at the origin.This singularity, which is due to a breakdown of the hydrostatic assumption, can be addressed by introducing a non-hydrostatic (i.e.source-driven) region near the origin, as discussed in [18].However, for the sake of simplicity, we ignore such details for the present study.For our purposes, (5.6) serves as a good approximation for the finger-like growth of region II (i.e. for r ≫ R).The transition time t * is found by matching the two thicknesses of regions I and II, such that Hence, the transition time and thickness are given by where in this case the initial aspect ratio is α = H/R.Much later than the transition time t ≫ t * , the flow continues to slump as a single region, similar to Section 2(b).In this case, the entire thickness satisfies (5.1).The boundary condition (5.By introducing dimensionless coordinates, and switching to similarity variables we arrive at a system of equations that can be solved analytically to give where η N = 2 (see Appendix B for further details).
Next, let's briefly discuss the swept shape and swept volume for this case.At early times t ≪ t * , the swept volume V is equal to the initial volume plus the instantaneous volume of region II, such that (5.15) At much later times t ≫ t * , the swept shape S(r) is calculated by finding the time at which the current thickness is maximal, which is equivalent to η = 2 1/2 in the similarity solution f (5.14).Hence, the swept shape is which is only valid for r ≫ R. The swept volume, which is now defined as is calculated (following similar steps as in Section 4) for late times as . (5.18) To illustrate the radially symmetric case, the swept shape S(r) is plotted in figure 6a,b, for ǫ = 10 −2 and 10 −3 .As with the two-dimensional case, stronger anisotropy results in a reduced swept shape (i.e.smaller values of S for r > R).Likewise, although not plotted, the swept volume increases slowly like V ∝ ǫt before transition to self-similarity, indicating that anisotropy reduces the contacted volume of pore space.

Discussion and applications
The gravity-driven spreading of a finite volume of fluid in anisotropic porous media differs from the isotropic case since the vertical flow is restricted by the permeability k V ≪ k H .In the initial dynamics the bulk of the flow descends slowly with a boundary layer near the impermeable base that diverts the flow into thin finger-like regions growing slowly in the lateral direction.This partition of the flow into bulk and finger regions reduces the swept volume of the gravity current compared to the isotropic case.This indicates that released volumes in anisotropic aquifers contact a smaller fraction of the available pore space.Hence, the spread of a contaminant in an anisotropic aquifer may be easier to contain since the contacted volume is reduced.By contrast, in the case of CO 2 sequestration, where the aim is to trap as much CO 2 in the pore space as possible, isotropic aquifers may have better potential in terms of residual trapping (which is a function of the contacted volume of pore space), though this ignores other trapping mechanisms such as dissolution, structural trapping and mineralisation [19].
It is important to consider the possible limitations of this model for such realistic scenarios.First and foremost, it must be noted that the anisotropic permeability values k H , k V , are upscaled quantities which attempt to capture the macroscopic effect of small-scale heterogeneities on the flow.These upscaled quantities are only a good approximation when the vertical length scale of the flow H is much larger than the heterogeneity length scale (e.g. the width of sedimentary layers) or when the permeable interval is inherently anisotropic due to compaction effects.If this is not the case, more complex flow models are required to treat the spread of fluid beneath and through successive layers [3,25,26].
Another consideration for the case of CO 2 sequestration is the effect of trapped saturation on the dynamics and spreading of the current.As the current moves, a fraction of its mass is lost to residual trapping due to small scale capillary forces [19] and dissolution within the surrounding brine [8].According to some trapping models [7] this can arrest the spread of the current altogether.Likewise, the spreading could also be arrested by lateral heterogeneities in the capillary pressure pinning the nose of the gravity current.
Whilst the conclusion of this study is that anisotropy reduces the swept volume of the gravity current, this doesn't take into account the enhanced trapping potential due to changes in the capillary pressure across heterogeneities, also known as capillary heterogeneity trapping [27].Specifically, during the imbibition cycle small-scale capillary forces induce a build-up of saturation beneath each sedimentary layer that can account for as much as 14% of the overall trapped saturation.Hence, the optimum anisotropy will no doubt strike a balance between the trapping associated with these heterogeneities and the reduction in swept volume that they induce.
It is also worth mentioning the flow of the ambient fluid, which we have so far ignored for this study.In particular, wherever there are sharp changes in the interface shape, this may be associated with non-negligible displacement of the local ambient fluid, such that the dynamic boundary condition (2.7) cannot be applied.Likewise, in the case where there is a viscosity contrast between the released and ambient fluids, this can modify the shape of the current [28] and cause fingering instabilities [29].In both of these cases a full numerical model would be necessary to resolve such flow details.In the case of carbon sequestration, CO 2 is typically 20-30 times less viscous than brine.As shown by [28], the viscosity contrast causes an enhanced spreading of the CO 2 in the shape of a thin finger along the cap rock.Hence, this viscosity contrast distorts the spreading of CO 2 in a similar manner to anisotropy, as studied here.Therefore, the expected effect of a viscosity contrast in anisotropic media is an extremely pronounced finger-like intrusion of CO 2 .Data Accessibility.The finite difference code used in this study can be found in the Supplemental Materials or on the personal website of the author: https://people.maths.ox.ac.uk/benham/openacces.zip Competing Interests.The author reports no competing interests.
Funding.There is no funding to report.
Acknowledgements.The author wishes to thank the anonymous peer reviewers for time taken to review the manuscript.

A. Additional plots
In this section we present some additional plots comparing the numerical solution from Section 3 with various approximate and analytical solutions.Whilst in the previous sections a rectangular shape was chosen for the initial released shape of dense fluid, here we consider a curved profile with initial shape h 0 (x) given by (3.2).In the case of an isotropic medium ǫ = 1, this results in immediate self-similar behaviour, as described in Section 2(b).We use the analytical selfsimilar solution (2.33) as a means of verifying our numerical method.In figure 7a,b,c, profiles of the gravity current are shown at three different times.Overall, very good agreement is found, indicating the soundness of the numerical method.
Unfortunately, no benchmark analytical solution exists for anisotropic porous media, but we can nevertheless compare against the approximate solution derived in this study.In figure 7d,e plots are shown for the same curved initial shape (3.2) released in an anisotropic medium with ǫ = 0.01.It is straightforward to extend the approximate solution derived earlier to this initial shape.As such, the bulk fluid region I initially evolves according to Meanwhile, region II initially evolves according to the self-similar dynamics (2.27) which correspond with a finger-like region fed by a constant input flux ǫu b L from region I (see Section 2(a)).These two solutions are then simply joined together to make the plots in figure 7d,e,f.Overall, good agreement is achieved indicating that our model can be extended to this and other such similar initial shapes.

B. Further details on similarity solutions
In this section, we summarise the equations that define the different similarity solutions used in the main text.Let's start with the two-dimensional equations for region II at early times.In this case the system of equations deriving from (2.22)-(2.25)after applying coordinate transformations (2.26)-(2.27) is −f df dη = 1 : η = 0, (A 2) −f df dη = 0 : η = η N .
(A 4) These can be solved numerically for the shape function f (η) and prefactor η N = 1.482.
Next we summarise the two-dimensional equations at late times, once the gravity current has transitioned to a single slumping region.In this case the system of equations

Figure 1 .
Figure 1.Schematic diagram.(a) A finite volume 2LHφ of dense fluid is released into a porous medium with an anisotropic permeability field k H ≫ k V above an impermeable base.(b)The medium is initially saturated with a less dense fluid so that the bulk region (I) descends slowly at speed −ǫu b , where ǫ = k V /k H and u b is the buoyancy velocity (2.10), whilst a thin finger region (II) develops near the base of the current.

Figure 2 .
Figure 2. (a,b,c) Approximate solution for regions I and II compared with numerical solution in the case of ǫ = 10 −2 .

Figure 3 .
Figure 3. Evolution of the vertical (a) and horizontal (b) extents of the gravity current shown in figure 2 (i.e. for ǫ = 10 −2 ),

Figure 4 .
Figure 4. Transition time t * and the transition thickness H * , determined by solving (2.36) for different values of the anisotropy ǫ and the initial aspect ratio α.Leading order scalings for small ǫ (see (2.37),(2.38))are shown with dotted lines.

Figure 5 .
Figure 5. (a,b) Swept shape of the gravity current S(x) (4.2) in the case of ǫ = 10 −3 and 10 −2 .(c) Self-similar isotropic case ǫ = 1 (note the modified initial conditions).(d) Evolution of the swept volume V(t) for each of these cases (note that time is stretched by a factor ǫ). Transition times t * are indicated with coloured dots.The initial aspect ratio is set as α = 1 in all cases.

2 )
is replaced by a zero flux condition at the origin

Figure 6 .
Figure 6.Swept shape of the gravity current S(r) in the radially symmetric case.Initial aspect ratio is set to α = 0.2 and the anisotropy is set to (a) ǫ = 10 −2 and (b) ǫ = 10 −3 .

Figure 7 .
Figure 7. Evolution of a gravity current with initial shape given by (3.2) in the case of an isotropic medium ǫ = 1 (a,b,c) and an anisotropic medium ǫ = 0.01 (d,e,f).Dimensionless times are given by u b t/Lφ = 0.87, 1.74, 8.68 in (a,b,c) and 1.70, 5.09, 17.0 in (d,e,f).