Two-phase gravity currents in layered porous media

We examine the effects of horizontally layered heterogeneities on the spreading of two-phase gravity currents in a porous medium, with application to numerous environmental flows, most notably geological carbon sequestration. Geological heterogeneities, which are omnipresent within all reservoirs, affect the large-scale propagation of such flows through the action of small-scale capillary forces, yet the relationship between these small and large scales is poorly understood. Here, we derive a simple upscaled model for a gravity current under an impermeable cap rock, which we use to investigate the effect of a wide range of heterogeneities on the macroscopic flow. By parameterising in terms of different types of archetypal layering, we assess the sensitivity of the gravity current to the distribution and magnitude of these heterogeneities. Furthermore, since field measurements of heterogeneities are often sparse or incomplete, we quantify how uncertainty in such measurements is manifest as uncertainty in the macroscale flow predictions. Finally we apply our model to the Sleipner case study in the North Sea, positing how heterogeneities may have played a role in the observed migration of CO2.


Introduction
Injection of CO 2 into underground reservoirs to reduce greenhouse gas emissions, also known as geological carbon sequestration, is one of the major proposed technological solutions to meet future global temperature targets (Bickle 2009;Bui et al. 2018). During this process, the injected CO 2 rises as a buoyant plume within the porous aquifer, encountering impermeable cap rocks which force it to spread laterally as a gravity current. As the flow spreads out, capillary forces play a key role in determining the saturation distribution and consequent flow properties via the relative permeability and capillary pressure (Nordbotten & Celia 2011). Heterogeneities in rock properties at the 10−100 cm scales substantially amplify and complicate the effect of variations in pore-scale capillary forces, and are manifest in the large-scale saturation distributions within the CO 2 current. Hence, in order to ensure safe and efficient sequestration, it is imperative to be able to model how small-scale heterogeneities, which are ubiquitous in all subsurface reservoirs, affect spreading rates at the macroscale (Benham et al. 2021;Jackson & Krevor 2020).
The only previous attempts to model such flows in heterogeneous media involve using upscaled relative permeability data, often acquired using numerical calculations or core flooding experiments (Jackson et al. 2018), and applying these to reservoir simulators (Braun et al. 2005;Cavanagh & Haszeldine 2014;Li & Benson 2015;Trevisan et al. 2017). However, these studies are often computationally demanding and focus on a specific type of heterogeneity (e.g. over some horizontal/vertical length scale, as investigated by Jackson & Krevor (2020)). In particular, it is not currently understood how generic smallscale heterogeneities affect the propagation of such large-scale flows. Furthermore, since heterogeneities are usually measured through isolated and sparsely distributed bore hole samples, such measurements often come with a large degree of uncertainty, yet there are very few studies which discuss how such uncertainty at small-scales translates to largescale predictions. Indeed, whilst the studies of Hinton & Woods (2020a,b) investigated how variations in permeability cause shear dispersion for miscible flows, the corresponding effects due to capillary forces within immiscible flows have not yet been addressed. However, as discussed by Jackson & Krevor (2020), these capillary forces associated with the heterogeneities play a critical role during CO 2 migration, and therefore require careful modelling.
In this study, we derive a simple upscaled model for the evolution of an axisymmetric two-phase gravity current beneath an impermeable cap rock, where the structure and distribution of layered heterogeneities is a model input. This simplified approach not only greatly reduces the computational demand of modelling such small-scale details, but also allows us to study a large range of heterogeneities by parameterising them as different types of archetypal sedimentary layering. In this way, we assess how the properties of the heterogeneities affect the macroscale flow, as well as how uncertainty in the measurements is manifest in such flow predictions. Our simple model provides other key insights, such as the self-similar nature of the upscaled gravity current, scaling laws for the speed of propagation and an understanding of where and when the flow transitions between the so-called capillary and viscous limiting behaviours.
In modelling subsurface migration of CO 2 , one of the key difficulties is resolving the complex heterogeneous rock properties of the porous medium. This is a two-fold challenge, since on the one hand a resolution which captures all the details of the rock geometry is very computationally intensive, and on the other hand it is almost impossible to obtain such fine-scaled resolution over field scales from field measurements. Hence, there is a strong motivation for an upscaled modelling approach which describes the bulk, or average effect of heterogeneities on the large-scale flow features. There are many possible levels of upscaling (from the pore scale upwards, as discussed by Krevor et al. (2015)), and this depends on the desired amount of detail as well as the available data, but here we focus on length scales between the size of the heterogeneities and the size of the aquifer. Therefore, in the context of this study small-scale heterogeneities refer to variations on the scale of 10 −2 -1 m (as opposed to pore-scale heterogeneities which occur on the scale of 10 −6 -10 −3 m), and the large-scale flow refers to a gravity current which is typically around 1 -10 m thick (Cowton et al. 2016).
Heterogeneity types range from variation within the pore structure of a rock to variation in the rock type itself. Owing to the complexity of geological processes, these heterogeneities arise from many different causes, including sedimentary layering, subsequent diagenetic changes in the mineral fabrics and tectonic fracturing and faulting. Each type of heterogeneity affects the flow in a different way, via the action of smallscale capillary forces, thereby presenting a significant challenge for generic upscaling approaches. However, the low computational cost of our approach allows us to investigate a wide variety of heterogeneity types via simple parameterisations of archetypal cases. Among these, we study the effects of lithostatic compaction as well as sedimentary strata with permeability sampled from a probability distribution. The latter case is particularly useful since it captures how the uncertainty in field studies, due to a lack of measurements, is manifest in the uncertainty of modelling predictions.
When upscaling the effect of heterogeneities, a key parameter is the capillary number (Jackson et al. 2018), which is defined as the ratio between horizontal flow-driving pressure gradients (associated with Darcy flow) and vertical capillary pressure gradients (associated with the capillary forces). Hence, in the limit of small capillary number, known as the capillary limit, the capillary forces due to heterogeneities dominate the flow behaviour, while in the limit of large capillary number, known as the viscous limit, they have a negligible effect. Many previous studies focus on each of these cases separately, though recently semi-analytical approaches have been derived by Benham et al. (2021) and Boon & Benson (2021) that capture the transition between the viscous and capillary regimes, demonstrating which regions of a confined aquifer are in each of these limits, and which regions are in between the limits. However, gravity was neglected in that study, restricting the applicability to very thin aquifers. For CO 2 sequestration sites in large aquifers, gravity plays a dominant role in the rise and spreading of the buoyant plume of injected fluid (Nordbotten & Celia 2011). The role of buoyancy is characterised by the ratio between the strength of gravitational forces and capillary forces, and may be characterised by a dimensionless Bond number (Golding et al. 2013), which we define later in Section 2.3. As discussed by Benham et al. (2021), the Bond number is greater than unity for aquifers larger than around ∼ 1 m thick, in which case gravity alters the upscaled flow properties significantly. Hence, in this study we focus on the upscaled modelling of such gravity currents, so that more general injection scenarios in larger aquifers can be addressed.
There is a long history of studying gravity currents in porous media, from early work which explored the first fundamentals (Huppert & Woods 1995), to later studies which investigated the effect of confinement (Pegler et al. 2014), permeability variations (Hinton & Woods 2018), and capillary forces (Golding et al. 2013). Recently, Jackson & Krevor (2020) showed that small-scale capillary heterogeneities can significantly modify the large-scale migration of a buoyant plume within an aquifer. However, their numerical study is both computationally intensive, and does not provide general scalings for different types of heterogeneity.
The aim of the present study is to quantify the macroscopic effect of a wide range of heterogeneities on the axisymmetric injection of CO 2 beneath an impermeable cap rock. The low computational cost of our simple approach allows us to explore different parameterisations of heterogeneities, providing insights into the dominant controls on the evolution of the gravity current. Similar to Benham et al. (2021) (though focussing on gravitational effects), we investigate both the viscous limit, the capillary limit and the transition between these limits using a locally defined capillary number that determines where and when heterogeneities play an important role. We show that away from this transition zone the upscaled gravity current is self-similar, where the front moves like the square root of time (like the homogeneous case discussed by Golding et al. (2013)) and the prefactor varies significantly depending on the type and strength of the heterogeneity, as well as the Bond number. In addition, we provide a framework for managing real permeability data with uncertainty in the measurements, illustrating how this uncertainty is manifest in modelling predictions. Finally, we use our upscaled approach to investigate how heterogeneities may have affected the injection of CO 2 at the Sleipner site in the North Sea (Bickle et al. 2007).
Our paper is laid out as follows. In Section 2 we derive a simplified model for the upscaled gravity current, discussing different types of heterogeneities. Section 3 presents numerical and analytical results in the viscous and capillary limits, as well as numerical results accounting for the dynamic transition between these limits. In Section 4 we apply our results to the case study of the Sleipner project, and finally we close with some concluding remarks in Section 5.  Figure 1: Schematic diagram of an axisymmetric gravity current (with constant injection Q) in both the homogeneous case (a) and the heterogeneous case (with sedimentary strata) (c), also illustrating the corresponding vertical non-wetting saturation profiles (b,d), given by (2.10),(2.12) (note the heterogeneity wavelength is exaggerated for illustration purposes). (e) Relationship between mean non-wetting saturations (2.13) and gravity current thickness h.

Upscaled modelling of two-phase gravity currents
In this section, we outline the assumptions used to model an upscaled two-phase gravity current in a heterogeneous porous medium, making note of how the saturation of phases varies within the current. Then, we derive the upscaled governing equations and boundary conditions which describe the macroscopic dynamics. Subsequently, we discuss a variety of different types of heterogeneity and how these are manifest in the upscaled properties. Finally, despite the added complexity of the heterogeneities, we demonstrate the selfsimilar nature of the gravity current, thereby greatly reducing the complexity of the problem.

Fundamentals of two-phase flow in heterogeneous porous media
The flow scenario we consider is illustrated in figure 1 with a radial coordinate system (r, z). A buoyant non-wetting phase (e.g. CO 2 ) is injected at a point source at the origin with flow rate Q into a surrounding porous medium saturated with a denser wetting phase (e.g. water). The resulting current spreads out radially under gravity with thickness z = h(r, t) beneath a horizontal impermeable cap rock located at z = 0. Motivated by the dominant heterogeneity arising from sedimentary layering, we consider a porous medium which has vertically varying permeability k(z) and porosity φ(z).
We model this scenario using conservation of mass and Darcy's law for two-phase flow (Bear 2013). Hence, the governing equations are where subscripts i = n, w indicate non-wetting and wetting phases, and S i , p i , u i , ρ i , µ i , k ri (S i ) are the saturations, pressures, Darcy velocities, densities, viscosities, and relative permeabilities of the two phases. We assume that the pore spaces are filled, such that S n +S w = 1. Furthermore, due to capillary forces, the pressure difference between phases satisfies where p c (S n ) is the capillary pressure. As is often done, we assume that both k ri and p c depend on the saturation only for simplicity (though in general they may have more complex dependencies). These are usually approximated with empirical parameterised formulae, such as those proposed by Corey (1954); Brooks & Corey (1964); Chierici (1984). Here we use the Brooks-Corey and Corey relationships, which are given by where p e (z) is the pore entry pressure, s = S n /(1 − S w0 ) is the reduced saturation of non-wetting phase, λ represents the pore size distribution, k rn0 is the end-point relative permeability, and α is a fitting parameter. The irreducible wetting saturation S w0 is the amount of wetting phase that is permanently stored in the pores during drainage flows, and consequently the end-point relative permeability corresponds to k rn0 = k rn (1). In this new formulation, the reduced saturation s conveniently varies between 0 and 1. The pore entry pressure p e (z) is the minimum pressure difference required to allow the non-wetting phase to enter the largest pore spaces at a given position. Likewise, as the pressure difference between phases increases, the non-wetting phase is able to enter smaller and smaller pore spaces. Hence, clearly the pore entry pressure depends on the size and geometry of the pores (and hence varies vertically), and the same is true for the permeability and porosity. However, whilst this dependence has been measured for specific rock types, it is not fully understood in general. Hence, as is often done, we use power laws to relate these different quantities, such that for some constants a, b > 0, which we take to be positive since large pore size corresponds to large porosity, large permeability, and small pore entry pressure. As discussed by Cloud (1941); Nelson (1994), we do not expect these constants to be the same for different rock types. Therefore, we keep them in general form for this analysis. However, we note a commonly used scaling proposed by Leverett (1941), p e ∼ (φ/k) 1/2 which implies b = 1/2(1 − a). Motivated by field observations of gravity currents (e.g. see Cowton et al. (2016), where the aspect ratio of the gravity current at Sleipner was calculated to be less than ∼ 1/1000) and following Golding et al. (2011Golding et al. ( , 2013, we assume that the gravity current is long and thin, such that the vertical velocity is much smaller than the horizontal velocity w i u i . In this case, the pressure within each phase is approximately hydrostatic where p 0 is the pressure at the edge of the gravity current (z = h) and ∆ρ = ρ w − ρ n . The saturation is calculated by combining (2.4) and (2.9), enforcing the physical lower bound on s, such that (2.10) To determine the value of p 0 we consider that the edge of the gravity current is defined as the boundary below which no saturation of non-wetting phase exists. Hence, from (2.9),(2.10), it is sufficient to ensure that s = 0 for all z > h if we choose p 0 = min p e (z). In other words, by setting the capillary pressure at the edge of the gravity current as the smallest required pressure difference to invade any pores in the aquifer, we guarantee that anywhere below the edge of the gravity current (p c < p 0 ) no saturation will be found. Therefore, even though there may be disconnected regions of non-wetting phase within 0 z h (e.g. see figure 1c,d), there will never be such regions for z > h. The saturation distribution (2.10) represents a balance between capillary forces (due to heterogeneities) and gravitational forces. However, this is only valid for situations where capillary forces are large enough to drive the saturation into regions of larger pore space, or equivalently when the capillary number is small. Therefore, in general the saturation distribution must depend on the local capillary number N c , which is given as the ratio between the horizontal flow-driving pressure gradient and the typical vertical gradient in pore entry pressure (Benham et al. 2021). For the former, we use the pressure gradient of non-wetting phase ∂p n /∂r, and for the latter we use ∆p e /h, where ∆p e = max p e (z) − min p e (z) is the maximum difference in pore entry pressure across the aquifer (constant), and the gravity current thickness h is used as the vertical length scale. Hence, the capillary number is given by (2.11) In the limit of very small capillary number N c 1, also known as the capillary limit, the saturation distribution (2.10) remains accurate. However, when the capillary number is very large, also known as the viscous limit, capillary forces due to heterogeneities are effectively negligible (i.e. we can ignore pore entry pressure variations p e (z) = p 0 ), and the saturation distribution becomes which is identical to the homogeneous case addressed by Golding et al. (2013) †. † We note that although the upscaled description of the viscous limit is mathematically identical to the homogeneous case, the model would still have to account for flow variations due to permeability gradients through an effective permeability. The saturation would, however, be identical to the homogeneous case.
For intermediate capillary number (i.e. when the flow is neither in the viscous limit nor the capillary limit), the saturation distribution lies in between (2.10) and (2.12), and therefore the expression for the saturation must contain the capillary number itself s = s(z, h, N c ). Typically, the dependence of the saturation on the capillary number is logarithmic (Benham et al. 2021), as with the upscaled flow properties, and we will return to address this later in Section 3.3.
In figure 1c we illustrate a radially symmetric gravity current in a heterogeneous layered medium composed of sedimentary strata. We contrast this to the classic homogeneous case in figure 1a (as studied by Golding et al. (2013)), which is equivalent to the upscaled viscous limit (N c 1) for a heterogeneous medium (see earlier footnote). For each case we plot typical vertical saturation profiles in figure 1b,d. In the homogeneous case the saturation distribution satisfies a balance between capillary and buoyancy forces, so that lighter regions (of high saturation) are pushed towards the cap rock. In the heterogeneous case, the same overall balance is sustained, but within that balance capillary forces push the saturation into layers where the pores are larger. Hence, significant oscillatory behaviour is observed within the vertical saturation profile, including patches where the saturation drops to zero. This corresponds with regions where the pore spaces are too small to allow any non-wetting phase (i.e. the zero value is chosen in (2.10)). Hence, one interesting consequence of heterogeneities is that they modify the mean saturation value in the gravity current. In figure 1e we plot the mean saturation, defined as whilst varying the gravity current thickness h for both the homogeneous and heterogeneous cases. In both casess is an increasing function of h, but the heterogeneous case always has a lower mean value. This is due to the substantial fraction of the gravity current with zero saturation.

Upscaled model: Governing equation and boundary conditions
Having discussed the flow scenario and laid down the key assumptions, now we outline the upscaling procedure, deriving a single governing equation and accompanying boundary conditions for the macroscopic description of the gravity current.
To do so, (2.8) is integrated to obtain the pressure, and then the conservation of mass equation for the non-wetting phase (2.1) is integrated between z = 0 and z = h(r, t), such that where ϕ = (1 − S w0 )φ 0 is the reduced porosity scaling, φ 0 and k 0 are typical scalings for the porosity and permeability (whereφ = φ/φ 0 ), and u b = k 0 k rn0 ∆ρg/µ n is the buoyancy velocity. We note that the flow of wetting phase is ignored in this analysis under the longthin approximation. Essentially, the flow of non-wetting phase within the gravity current decouples from the flow of wetting phase, which is not present at leading order. Nevertheless, multiphase effects are still manifest at leading order via the multiphase properties, such as the relative permeability and capillary pressure. However, as we discuss later in Section 2.3, this assumption breaks down if the permeability ratio between layers becomes very small. In particular, if there are regions of very low permeability, these will act as a vertical obstruction to the flow. In such situations, as discussed by Pegler et al.
(2014), the flow must be treated as confined, where the flow of the ambient fluid plays an important role on the dynamics and therefore can no longer be ignored. We give more details of this consideration in the next section.
We also note that (2.14) is already an upscaled description of the flow, since the heterogeneities are only manifest within the integrals. Hence, (2.14) represents how the heterogeneities affect the evolution of the gravity current in a spatially averaged sense. Such an upscaling approach is desirable, since we wish to avoid having to resolve all the heterogeneities, both to reduce computation time, and also because realistically the low resolution of field measurements means that such details are uncertain anyway.
It is convenient to write (2.14) in a more standard diffusion equation form to render it amenable to conventional analysis. Therefore, by switching variables to the integrated saturation S(h, N c ), which is defined as (2.15) (2.14) is rewritten as where the flux is given by Further details of this coordinate transformation are presented in Appendix B. We note that S has dimensions of length, and F has dimensions of length squared over time. Therefore, (2.16) is just a standard diffusion equation for the total volume of non-wetting phase (per unit area), where the flux is a non-linear function that represents how capillary forces modify the flow. Hence, there is an interesting analogy between our scenario and a viscous gravity current, where the flux function represents how viscous forces modify the flow (e.g. plug flow, Poiseuille flow, etc...). As we will find out later, F is sometimes well-approximated by a power law of S, and the solutions to such equations are detailed by a large historical body of literature (see Huppert (1982) for example). In general, (2.16) must be solved in tandem with the equation for the capillary number (2.11). Therefore, writing (2.11) in terms of the integrated saturation S, we arrive at the transcendental equation for the capillary number where h is written in terms of S under the assumption that (2.15) has a uniquely-defined inverse (which we later find to be the case).
For the remainder of this study (up until Section 3.3) we restrict our attention to the two limiting cases of small and large capillary number (capillary and viscous limits), where the saturation is given by (2.10) or (2.12) and the flux is just given by F = F(S), Sedimentary strata (Hlow/Hhigh = 1)

Turbidites
Compacted Spectrum Figure 2: Illustrations of the different types of heterogeneity we consider, where the heterogeneity is characterised by variation of the permeability with depth. (a-f) represent the deposition of sediments through various geological mechanisms, whereas (g) represents compaction due to lithostatic pressure. In (a,b,c) we illustrate the case of sedimentary strata with greyscale permeability maps for three different values of the width ratio between low/high permeability regions (H low /H high ). In the spectrum case (f) we display the probability density function (PDF) of the permeability which is randomly sampled from a uniform distribution on a logarithmic scale. thereby decoupling (2.16) and (2.20). However, later in Section 3.3 we address the case of intermediate capillary number, for which the equations must be solved in tandem.
The governing equation (2.16) must be accompanied by some initial and boundary conditions to create a well-posed system. Firstly, we define the nose of the gravity current at position r = r N (t), at which the thickness is zero, such that S| r=r N (t) = 0.
(2.21) Secondly, following Golding et al. (2013), we impose global conservation of mass of the non-wetting phase, such that 2πˆr or equivalently, we impose the input flux at the origin and zero flux through the nose of the current, The finite flux value Q in (2.23) indicates that the gradient ∂S/∂r must become infinite Heterogeneity type Functional form Parameters Table 1: Definitions of the different types of heterogeneity (characterised by the permeability), as displayed in figure 2. Sedimentary strata take binary permeability values k low /k high with the width ratio of low/high regions given by H low /H high . Turbidites, the deposits of turbidity currents, consist of a periodic array of layers with linearly varying permeability, where the wavenumber n, is considered in the limit nh → ∞. In the spectrum case, permeability is a series of strata, where each layer has permeability taken from a uniform random distribution, distributed logarithmically across range [k low , k high ]. Likewise, the widths of the layers are taken from a uniform random distribution on a linear scale. Compacted rock corresponds to a permeability profile which decreases with depth under a power law β, starting with a finite value at z = 0.
as r → 0. Therefore, it is expected that the long-thin approximation made earlier may become invalid very close to injection. Furthermore, near the nose of the gravity current r = r N , where the gravity current becomes thinner than the heterogeneity length scale, we do not expect our upscaled approximation to be accurate.

Incorporating heterogeneity
To close the system, we must choose a type of vertical heterogeneity. Since we have used power laws a, b to relate the porosity and pore entry pressure to the permeability, we need only choose a functional form for k(z). Whilst in general this function may vary in three dimensions (i.e. k(x)), here we restrict our attention to pure vertical variation since, not only does this capture the leading order behaviour for sedimentary layering, but also because this is consistent with the long-thin approximation of a gravity current made earlier.
To model the permeability we have a variety of different physically-motivated choices which we list in table 1 and plot in figure 2. Firstly, sedimentary strata represent a periodic deposition of two different types of sediments, such that the permeability alternates between two values, k low and k high , in a periodic array of layers, where the width ratio of each of these is given by H low /H high (see figure 2a,b,c). Unlike the sedimentary strata, which are uniformly deposited in each layer, turbidites represent the deposition of sediments from the continuous arrival of turbidity currents, such that within each layer the permeability varies linearly. The sign of the linear slope indicates that layers become more permeable as one descends deeper, since this corresponds to the early/late arrival of large/small particles in a turbidity current. Thirdly, we consider a permeability profile which is generated by randomly sampling from a distribution, or spectrum, of permeability values, spread out logarithmically. This case is motivated by realistic measurements of sedimentary strata which are often noisy and aperiodic. Finally, we consider a compacted rock, where the permeability decreases with depth due to the buildup of lithostatic pressure over time.
Although there are many other possible choices for the permeability, we restrict our investigation to these four examples since they are canonical cases from which we may learn about the fundamental effects of heterogeneities. Each case is parameterised, either by the ratio of the permeabilities and widths of the lowest-highest permeability regions k low /k high , H low /H high , or by the compaction power law β, which represents the strength of the compaction effect.
It is important to note the possible limitations on these parameters. In particular, sufficiently low permeability layers may cause a vertical obstruction, such that an unconfined description of the gravity current is no longer applicable. To investigate the limitations on the permeability ratio k low /k high we have performed a set of numerical simulations of the two-dimensional miscible Darcy equations using the DarcyLite finite element package in Matlab, adapted to account for gravity (Liu et al. 2016;Harper et al. 2021). The miscible Darcy equations are equivalent to the immiscible equations (2.1)-(2.2) in the limit where the relative permeabilities become independent of the phases, k rn , k rw → 1, and the phase pressures equalise such that p c → 0. Studying the miscible flow problem allows us to investigate the applicability of upscaling for small values of the permeability ratio without accounting for the more complex effects of immiscible phase saturations. We do not display the numerical results here, since a rigorous analysis of this query is outside the scope of this paper, but we describe our findings here in writing instead.
For very small values of the permeability ratio (e.g. k low /k high = 0.001) there are several important observations from these numerical simulations. At early times, due to the effective obstruction from the low permeability layers, the injection is focussed within the nearest high permeability layers instead, and behaves approximately like a confined flow in which the pressure has significant streamwise gradients (i.e. deviating away from the hydrostatic condition (2.8)). As a result, the shape of the gravity current is highly distorted and loses its self-similar structure. At later times, once the gravity current has invaded a sufficient number of vertical layers, it begins to assume self-similar dynamics and the pressure becomes hydrostatic to good approximation. Therefore, there is no strict lower bound on the permeability ratio k low /k high for an upscaling procedure, but rather this becomes a question of temporal and spatial scales. In other words, for any positive permeability ratio, given enough time and spatial extent, such an injection will eventually resemble a self-similar gravity current and is therefore amenable to upscaling. However, the overall aspect ratio of the gravity current is significantly reduced (i.e. longer and thinner) due to the delay in the vertical flow caused by the obstruction of the low permeability layers. Therefore, to avoid dealing with the prolonged transient effects that precede self-similarity in the case of very low permeability ratios, for the remainder of this study we restrict our attention to k low /k high 0.01.
Continuing our upscaling analysis we note that, given a particular type of heterogeneity k(z) and power laws a, b for the porosity and pore entry pressure, the integrals (2.15),(2.18)-(2.19) must first be calculated before we can solve (2.16). For general values of a, b, these integrals must be calculated numerically, using a trapezoidal integration rule, for example. In the layered cases, we wish to remove the dependence of these integrals on the heterogeneity wavelength, since it is undesirable to have upscaled properties like F that oscillate depending on the gravity current thickness. Therefore, instead of resolving all of the details of the flow, we build a macroscopic picture of the gravity current, which is consistent with other upscaling approaches (see Rabinovich et al. (2016) for example).
In the case of sedimentary strata, since k (and therefore p e and φ) takes either one of two possible values, integrals can be simply decomposed into bulk fractionŝ thereby removing the need to resolve individual layers. A similar approach can be taken in the case of the permeability spectrum, although in that case (2.25) is replaced by a sum over the number permeability values sampled from the random distribution †. However, in the case of the turbidites, to remove the dependence on the wavenumber n (as defined in table 1), the integrals must be calculated with a fine numerical scheme for a large but finite value of nh 1. The most salient features of this analysis are the integrated saturation S(h) and the flux F(h), since F allows us to solve the diffusion equation (2.16), and S allows us to calculate the gravity current thickness by way of inversion. These both depend on a number of non-dimensional parameters. Ignoring the capillary number (since for now we restrict our attention to N c 1 or N c 1) there are a total of 8 non-dimensional parameters which govern the problem. These consist of the heterogeneity parameters k low /k high , H low /H high , β (if compaction present); the power laws relating porosity and pore entry pressure to the permeability a, b; the Brooks-Corey parameters λ, α; and finally the Bond number, which is defined as The Bond number, which can alternatively be written as Bo= ∆ρgH/p 0 , where H = Q/u b is the buoyancy length scale, is interpreted as the ratio between buoyancy forces and capillary forces. This quantity largely controls the saturation distribution (2.10),(2.12), which is evident upon dimensional analysis. For example, when Bo 1 the saturation, written in dimensionless form, approximates to regardless of which type of heterogeneity p e (z) we consider (so long as Bop 0 /p e 1). We note that some of the above parameters have already been studied by other authors. For example, the power laws a, b were already addressed by Benham et al. (2021) and the Brooks-Corey parameter λ was studied by Golding et al. (2013). Therefore, for the remainder of the current study we focus on the heterogeneity parameters k low /k high , H low /H high , β, and the Bond number as the key parameters of interest. We use the homogeneous case k low /k high = 1 as a proxy to study the viscous limit behaviour N c 1, since they are equivalent (see footnote in Section 2.1). For k low /k high < 1, however, we assume capillary limit behaviour N c 1. We fix the remaining parameters at typical values a = 1/7, b = 3/7 (using the Leverett scaling), λ = 3 and α = 4, which we have extracted from a variety of different sources (Golding et al. 2011;Berg et al. 2013;Bickle et al. 2017). For the remainder of this study (up until Section 4), we keep the same values for these parameters so that we can focus on the effect of the heterogeneities instead, but we note that our approach is by no means restricted to these values. Nevertheless, continuing with these parameter values, we illustrate how the flux F depends on the type of heterogeneity and the Bond number in figure 3a,b,c. For each of the layered cases, we use a permeability ratio value of k low /k high = 0.1, whereas in the compacted case we use a power law of β = 1. In all cases (except the spectrum case) the flux is well approximated by a power law F ∝ S ψ , for some value of ψ between 1/2 and 2. In some cases, as we will show later in Section 3.2, these power laws can be derived analytically. Figure 11 in Appendix A displays the integrated saturation S, as well as the velocity distribution u n ∝ ∆ρgk(z)k rn (s)/µ n within the gravity current. There are several interesting observations to make. Firstly, no matter which type of heterogeneity nor which Bond number we choose, the integrated saturation S(h) is always a monotone Methodology Figure 4: Schematic illustration of our methodology, with stages going from left to right (a-e). We start by parameterising the heterogeneity k(z), p e (z), φ(z); then we use (2.10) to determine the saturation distribution s(z, h); then we obtain the velocity distribution u n ∝ ∆ρgk(z)k rn (s)/µ n (velocities for high and low permeability regions are illustrated as well as the mean); then from this we calculate the integrals comprising the flux F(h(S)) (2.17); then finally we use (2.16) to calculate the gravity current thickness h (via S(h)).
increasing function, such that the inversion h = S −1 (S) is always well-posed. Secondly, we note that there is an interesting interpretation to the value of the flux power law ψ, by way of analogy to viscous gravity currents. In the governing diffusion equation for a classic viscous gravity current, the flux power law relates to the velocity distribution within that current. For example, the velocity distribution for Poiseuille flow, which varies quadratically in the vertical coordinate, when integrated gives a cubic flux power law. Likewise, a uniform plug flow, when integrated gives a linear flux power law. In general, any viscous gravity current flux power law can be achieved by considering a shear thinning/thickening power law viscosity µ ∝ (∂u/∂z) κ . Then, it is easy to show that the lubrication balance µ∂ 2 u/∂z 2 ∼ ∂p/∂x can be integrated to give a flux F = h 0 u dz which obeys the power law F ∝ h 2+1/(1+κ) . This is illustrated in figure 3d, indicating specific cases with dashed lines. For example, a shear-thinning fluid with power law κ = −5/3 will produce a flux with power law F ∝ h 1/2 . Whilst our upscaling problem is very different from a non-Newtonian viscous gravity current, the analogy is nevertheless useful in helping to relate the flux functions observed in figures 3a,b,c, to the velocity distributions within our gravity current (which are displayed in figure 11b,d,f, in Appendix A). Now that all the steps in our approach have been outlined, we summarise our methodology for analysing the gravity current in figure 4. This illustrates the steps between initially choosing a heterogeneity type and finally solving for the gravity current thickness h. We have provided some example code in the supplemental materials to demonstrate these steps in the case of sedimentary strata, including how to numerically calculate the flux functions.

Discussion of self-similarity and the numerical scheme
There is a final simplification that can be made owing to a coordinate invariance, which allows calculation of the solution using a simple numerical scheme. In particular, much like the classic single phase axisymmetric † gravity current discussed by Lyle et al. (2005), the heterogeneous case is self-similar. Upon inspection, for constant flux our governing equation (2.16) (under the assumption of viscous N c 1 or capillary N c 1 limits) admits the similarity variables where the nose of the gravity current is located at η = η N for some constant η N which we will determine shortly. To further simplify the equations, and to convert to a unit interval domain, we write our system in terms of the variables y = η/η N andf (y) = f (η). In this way, the governing equation for the integrated saturation (2.16) and the boundary conditions (2.21),(2.22) become The system can be solved numerically using a simple finite difference scheme, starting at y = 1 and marching back towards y = 0, where the second order ODE (2.30) is conveniently written as a set of two first order ODE's with boundary conditions where L is the total dimensionless flux, and 1 is a small but finite numerical value (we cannot choose = 0 since it will generate infinite gradients). To find the constant η N , we start with an initial guess η N0 , and then use Newton's method to iteratively solve the flux condition (2.31) (see our code which we have uploaded in the supplemental material).
We make the key observation that, independent of the form ofF(f ), the gravity current is self-similar, with a nose that moves like the square root of time. Hence, the heterogeneities are only capable of modifying the prefactor η N for the nose speed, not the power law (which is always r ∼ t 1/2 ). However, the heterogeneities may also change the shape of the gravity current via F and S.
As discussed later in Section 3.3, in the case where the capillary number is neither small nor large, the flux must depend on the capillary number itself. In this case, since the capillary number involves derivatives of S with respect to r, this modifies the form of the governing equation (2.16), rendering such similarity variables inadmissible. Over † Note that the two-dimensional case is not necessarily self-similar. The two-dimensional gravity current thickness scales like h ∼ t 1/3 , such that the flux function F (h) cannot be written in a general self-similar form. However, this becomes possible in certain specific cases (e.g. a linear power law F ∝ h, as discussed by Huppert & Woods (1995)). time the solution changes from a viscous limit regime (self-similar) to a capillary limit regime (self-similar), but the transition itself, therefore, cannot be self-similar.

Results: viscous limit, capillary limit and transition
Our results comprise the following three different cases. (i) The capillary number is large throughout the aquifer (viscous limit), in which case the upscaled problem is equivalent to the homogeneous case.
(ii) The capillary number is small throughout the aquifer (capillary limit), in which case the upscaled problem is dominated by the effect of the heterogeneities.
(iii) The capillary number varies across the aquifer, such that different regions simultaneously lie within the viscous limit and the capillary limit, and other regions lie between these limits.
We start by addressing the former two cases (viscous and capillary limits), for which the problem is self-similar. The homogeneous case is used to study the viscous limit (since they are equivalent in an upscaled sense) and the heterogeneous cases are used to study the capillary limit. The gravity current thickness must be calculated numerically by solving the ODE system (2.31),(2.33)-(2.36). In specific limiting cases, such as when the Bond number is very large or very small, analytical progress can be made. We first present our numerical results which we use to understand the broad effect of the heterogeneities across large parameter regimes. Then, we use asymptotic analysis to help interpret the results in certain limiting cases. Finally, we address the transition between the viscous and capillary limits, for which the full PDE system (2.16),(2.21),(2.22) must be solved in tandem with an algebraic equation for the local capillary number (2.20).

Numerical solution in the viscous and capillary limits
The capillary limit numerical solution for different types of heterogeneities is plotted in figures 5 and 6, for Bond number values between Bo = 10 −2 and Bo = 10 2 . Typical saturation profiles are also displayed as inserts in each plot. To compare the different profiles, we have normalised the thickness by twice its mid range value h half = h(r N (t)/2, t). The radial coordinate is normalised by the nose position r N (t) so that the shape remains on a fixed domain for all time. For now, we do not include plots of the viscous limit numerical solutions, since these are very similar to the study by Golding et al. (2013). However, shortly we will use these as a reference when comparing the different types of heterogeneities.
Let us first focus on the non-compacted cases in the capillary limit (figures 5 and 6a,c,e). For small Bond number Bo 1, the saturation becomes near-zero s ≈ 0, but with spikes of linearly increasing magnitude that represent the thin regions where the permeability is near its maximum value k ≈ k high . As we increase the Bond number towards unity, the saturation becomes larger, with an overall curved profile and significant oscillations. At high Bond number, the saturation tends towards s ≈ 1 everywhere except very close to z = h, where it rapidly drops to s = 0. The shape of the gravity current changes from having a sharp nose at high Bond number to having a rounded blunt nose at low Bond number. There is not a noticeable difference in the shape of the gravity current between the different types of heterogeneity.
Apart from the shape and the saturation distribution, there are two other important metrics which are useful for describing the current. Firstly, the prefactor η N relates to the speed of the advancing nose, and secondly the mid-range thickness h half = h(r N (t)/2, t) indicates the approximate size of the current. Following Golding et al. (2013), we use the classic single phase limit values as a useful reference. Using a subscript notation, these (e) (f) Figure 5: Numerical results for the capillary limit in the case of sedimentary strata (a,c,e) and turbidites (b,d,f) (with k low /k high = 1/3, H low /H high = 1). From top to bottom, capillary forces become less important with respect to gravitational forces. The radius r is given in terms of the nose position r N (t), and the thickness h is normalised by the reference value 2h half = 2h(r N (t)/2, t) for the sake of comparison. The heterogeneity wavelength is exaggerated for illustration purposes. In each plot inserts illustrate the vertical saturation profile, normalised by the uppermost value s 0 = s(0).
are given by η N0 = 1.155 and h half0 = 0.348H (Lyle et al. 2005). In figures 7a,b,c, 8a we plot these quantities for different values of the Bond number. In the limit Bo 1 all cases converge to the single phase limits, which is expected due to (2.27). Likewise, in every case the flux converges to a linear power law, corresponding to a uniform velocity profile, as can be seen in figures 3c and 11f. In the limit Bo 1 the mid-range thickness h half behaves similarly for all three layered cases, growing approximately like h half ∼Bo −1/2 , as described by Golding et al. (2013) for the homogeneous case. On the other hand, the prefactor η N behaves rather differently. In all cases, we see an increase in η N for stronger heterogeneities, indicating that capillary forces accelerate the gravity current. However, each heterogeneity affects the prefactor η N differently, as can be seen in the different shaped curves in figure 7. This reflects the complex nature of the velocity distributions and flux functions depicted in figures 3 and 11. It is interesting to note that despite having a permeability profile with the same mean value, the different variations within each layer for each heterogeneity type are sufficient to alter the flux of saturation, and hence modify the speed of propagation of the gravity current. This sheds light on both the need for detailed bore hole measurements to infer as much information about the heterogeneities as possible, as well as the usefulness of such an upscaling approach as we have taken here. For each of the different types of strata, we compare the capillary limit curves against the homogeneous case (solid black line), which is equivalent to the viscous limit. This allows us to quantify the effect of the heterogeneities on the prefactor more clearly. The strongest effect on the prefactor occurs when there are thin regions of high permeability (H low /H high = 10) in which the non-wetting saturation concentrates. This focussing of the saturation feeds into the nonlinearity of the relative permeability, thereby amplifying the effect on the flux function F. By contrast, thin regions of low permeability (H low /H high = 0.1) produce results which are very close to the homogeneous case, so we do not display them here.
In the case of the permeability spectrum, we choose a permeability distribution whose standard deviation divided by the mean is σ(k)/k 0 = 1, and whose mean permeability ratio (between lowest and highest values) is µ(k low /k high ) = 0.04. We calculate the prefactor η N for 50 different realisations of this distribution and plot the results in figure 7c. For each Bond number we display the mean value, as well as one standard deviation on either side µ(η N ) ± σ(η N ). The mean result is reminiscent of the previous sedimentary strata cases. However, it is interesting to note that the standard deviation is largest for low-medium Bond number and shrinks as the Bond number gets larger. Hence, predictions are particularly prone to uncertainty if the Bond number is less than

Spectrum Compacted
Figure 7: Nose growth prefactor η N given in terms of the single phase limit η N0 = 1.155 for all heterogeneity types, parameterised by the permeability ratio k low /k high , the width ratio H low /H high , and the compaction power law β. In the case of the permeability spectrum, we show the mean result alongside one standard deviation above and below. Limiting behaviours are illustrated with dashed lines. Solid black curves correspond to the homogeneous case, which is equivalent to the viscous limit, whereas all other curves correspond to the capillary limit.
order unity. For CO 2 sequestration applications, this indicates that particular attention towards measurements of heterogeneities should be paid for injection sites with low flow rates. Next, we move on to describe the compacted case in the capillary limit ( figure 6b,d,f). The presence of compaction is most noticeable in the small Bond number cases. By comparing figures 6b and 6a, we see that compaction significantly increases the saturation within the gravity current, which is due to the permeability gradient forcing the nonwetting phase upwards. This is accompanied by an increase in the prefactor η N and a decrease in h half , as seen in figures 7d, 8b. This is expected since, if a larger saturation is maintained, volume conservation indicates that the gravity current must be thinner. By varying the compaction power law β from 0 to 1, we illustrate a fairly uniform transition of the values of η N and h half between those of a uniform rock and those of a strongly compacted rock.
For Bo 1 the saturation becomes near uniform s ≈ 1, as before, and this is accompanied by both η N and h half converging to their single phase limits. Hence, at such large Bond numbers the effect of compaction on the saturation distribution and gravity current shape is not particularly noticeable. This is expected, since compaction forces the saturation upwards, just like gravity.
For applications to CO 2 storage, it is useful to infer from the above results how much we can expect heterogeneities to affect the speed of propagation of a gravity current. Such information allows one to make efficiency predictions for CO 2 storage that help pinpoint the best sites for injection, as well as safety predictions that ensure the CO 2 does not spread beyond the desired perimeter. Using the homogeneous case, k low /k high = 1 (which is equivalent to the viscous limit), as the base case, we define the efficiency parameter ν as the relative difference we can expect heterogeneities to make on the prefactor η N , such that ν(Bo, k(z)) = η N het /η N hom − 1. (3.1) Clearly ν depends on a number of parameters, but here we focus on the different types of heterogeneity k(z) and the Bond number. Restricting our attention to the layered cases (ignoring compaction), we plot the heterogeneity efficiency ν in figure 9 for different values of the Bond number. The largest heterogeneity efficiency is observed for the sedimentary strata with thin bands of high permeability (H low /H high = 10). As we mentioned earlier, this can be explained by a nonlinear focussing of the saturation into these high-speed bands. The most we can expect heterogeneities to accelerate the gravity current at high Bond number is around ν = 10 − 50%, whereas at low Bond number the speedup can be as much as ν = 100 − 200%. In the case of the permeability spectrum we also illustrate the standard deviation of the predictions, indicating that the results are particularly sensitive to uncertainty at low Bond number.
It is interesting to note that whilst fluid is injected at constant flow rate Q, the Bond number is held constant, but if the flow were to stop suddenly this would no longer be the case. In such situations, the buoyancy length scale H = Q/u b , which was previously used to define the Bond number, would be rendered meaningless. Instead, the appropriate length scale for the flow would be the gravity current thickness itself h which, after the cessation of Q, would gradually decrease towards zero. Hence, the Bond number P o s ti n j e c t i o n Figure 9: Heterogeneity efficiency ν (3.1), describing the relative increase in prefactor value η N due to heterogeneities, given as a ratio of the prefactor value for the homogeneous case. Here we focus on the layered cases (S.S. stands for sedimentary strata), ignoring compaction. We illustrate how in post-injection scenarios the Bond number typically decreases over time, such that heterogeneities play an increasingly important role. In the case of the permeability spectrum we plot the mean value as well as one standard deviation on either side. The permeability ratio for all cases is k low /k high = 0.04. of the flow would decrease accordingly, as illustrated in figure 9, causing the effect of the heterogeneities to be increasingly amplified with time. This is particularly relevant for CO 2 storage applications, where the injection of gas is switched off once the aquifer is deemed to have reached maximum safe capacity. Hence, in such situations, it is clear that modelling heterogeneities is essential for understanding the post-injection spread of the CO 2 . However, it is important to note that such situations involve imbibition flows, as opposed to drainage flows, as we have studied here (Woods 2015). Imbibition flows typically have different capillary pressure and relative permeability curves than drainage flows, though the approach studied here is still applicable.
In any case, we have shown here that heterogeneities have the potential to significantly alter the growth of the gravity current, and consequently careful upscaled modelling is required. In the next section, we will show that some of these limiting behaviours can be explained using asymptotic analysis.

Limiting cases and analytical solutions
Some simplifications can be made in the limits of strong and weak capillary forces (i.e. small and large Bond number). Here we address these and derive analytical solutions which we use to explain some of our earlier numerical results. We split the analysis into situations without compaction β = 0 and with compaction β > 0. As in the previous section, here we restrict our attention to viscous limit and capillary limit behaviour only.

Weak capillary forces without compaction
We already showed earlier that in the limit of large Bond number the saturation distribution (2.27) is approximately uniform s ≈ 1. Inserting this into the integrals (2.15),(2.18), we see that which allows us to calculate the dimensionless flux This linear power law matches with our numerical observations in figure 3c. Comparing (3.2),(3.4) with the study by Lyle et al. (2005), we see that the limit of large Bond number is identical to the classic single phase limit. This is of course expected, since in the limit of weak capillary forces the flow of phases decouples, such that a single phase model becomes appropriate. This explains the convergence behaviour for both η N and h half for Bo 1, as we illustrate with dashed lines in figures 7 and 8.

Strong capillary forces without compaction
To address the limit of strong capillary forces, we first consider the homogeneous case k low /k high = 1 (equivalent to the viscous limit) since this makes the analysis for the subsequent heterogeneous cases easier. Hence, in the limit of small Bo 1 in the homogeneous case, the saturation distribution (2.12) approximates to a linear function Inserting this into the integrals (2.15),(2.18) we get from which we calculate the dimensionless flux In this case, the flux has an α/2 power law, which matches with our numerical calculations in figure 3a (for which α = 4). The gravity current thickness is given by inverting (3.6), such that Clearly, the thickness grows like h ∼ Bo −1/2 as Bo→ 0, which we illustrate in figure 8 with dashed lines. Our numerical results show good agreement, indicating their robustness. We also note that the square root in (3.9) explains why we see a blunting of the gravity current nose at low Bond number in figures 5,6.
If we now consider a finite heterogeneity k low /k high < 1, then in the case of sedimentary strata, the saturation distribution takes one of two possible values (3.10) Since the low/high permeability layers are distributed according to the ratio H low /H high , the integrals (2.15),(2.18) approximate to (3.12) Hence, the dimensionless flux iŝ Like the homogeneous case, the heterogeneous flux has an α/2 power law, which also matches with our numerical calculations for sedimentary strata in figure 3a (for which α = 4). We note that in the above analysis, the saturation approximation (3.10), and consequently the flux (3.13), are only valid for heterogeneities which have a small enough permeability ratio that (1 − k low /k high )/(1 + k low /k high ) Bo. However, for this study we restrict our attention to significantly heterogeneous media (rather than weakly heterogeneous).
We also note that in both the homogeneous case and the heterogeneous case, the coefficients in (3.8),(3.13) depend on the Bond number itself. Therefore, we do not expect a constant value asymptote for the prefactor η N in the limit Bo→ 0 (except in the specific case α = 2), which is consistent with our numerical observations in figure 7.
In the case of the permeability spectrum, a similar analysis is possible since the only contribution to the integrals will come from the regions with the largest permeability value k = k high . However, since there are potentially many more than just two permeability values in the spectrum, one needs to replace the factor H high /(H high + H low ) in (3.11)-(3.13) with the fraction of the aquifer occupied by such high permeability layers, which we denote H high /H total . In the case of the turbidites, a similar analysis to the above is much more difficult, since we cannot approximate the saturation distribution as simply as (3.10).

Weak capillary forces with compaction
In the case where the rock is compacted, the saturation (2.27) must be approximately uniform s ≈ 1 in the limit of large Bond number, as before. Inserting this into the integrals (2.15),(2.18), we see that Note that if either a or aβ equals unity, we get analytical expressions with logarithms instead of (3.14),(3.15). Assuming β = 1, aβ = 1, we then calculate the dimensionless flux, which we write in terms of h for now: (3.16) Clearly the flux (3.16) is not a linear function of S (3.14) as in the homogeneous case. However, we note that for weak porosity-permeability power laws a 1 (e.g. a = 0.14 for the Salt Creek case study (Bickle et al. 2017)) the integrated saturation (3.14) approximates to a linear function S ≈ ϕh, as we observe in our numerical calculations in figure 11e. Furthermore, in situations where the compaction law is weak β 1, or alternatively where the gravity current thickness is small h/H 1, the flux (3.16) reduces toF ≈ h/H, thereby recovering the single phase limit (3.4). This is in accordance with our numerical observations in figure 3c.
However, we note that by choosing sufficiently strong compaction/porosity power laws a, β, the single phase limit is no longer recovered in the limit of large Bond number. Physically speaking, this is because even with weak capillary forces, if the medium is sufficiently compacted the velocity distribution within the gravity current becomes dominated by the permeability variation, which upon integration creates a flux power law that is not equal to one.

Strong capillary forces with compaction
To address the compacted case at small Bond number, we first consider the form of the saturation (2.10). In particular, we note that since k is monotone decreasing with depth, the saturation s is also monotone decreasing. In particular, we can calculate that s will intercept zero at some depth z * which satisfies ∆ρg(h − z * ) + p 0 = p e (z * ), and will remain zero for all larger depths than this z > z * . Therefore, without loss of generality, we take p 0 = p e (h), such that there are no regions of zero saturation within the gravity current. In this way, for small Bond number (2.10) approximates to a finite expression which is independent of Bo, ( 3.17) This explains why the properties of the gravity current (e.g. η N , h half ) asymptote to constant values for Bo→ 0 in figures 7d,8b. For general power laws a, b, λ, α, β, inserting (3.17) into the integrals (2.15),(2.18) leads to very complicated analytical expressions, which we do not display here. As we can see in figure 3a, the flux for this case does not obey a fixed power law, but the exponent varies roughly between 1 and 2.

Viscous-capillary transition
Up until now we have only discussed situations where the capillary number (2.20) is either very large (viscous limit), in which capillary forces due to heterogeneities can be ignored (i.e. effectively the homogeneous case), or very small (capillary limit), in which the heterogeneities play a dominant role on the flow behaviour (i.e. the heterogeneous cases studied above). However, in general the capillary number may vary between small and large values throughout the aquifer. In this case, neither the viscous nor the capillary limit can be applied to the flux function in (2.16), and instead the flux must depend on the local capillary number, which is effectively a measure of the local pressure gradients within the gravity current. In this section, we discuss how to model this using numerical simulations of the full PDE (2.16) coupled to the transcendental equation (2.20), thereby determining which regions of the gravity current are within the viscous and capillary limits, and which regions lie in between these limits.
Following the same approach as Benham et al. (2021), we formulate composite functions for the upscaled properties of the flow, which capture both the viscous and capillary limit regimes, as well as the transition between these limits. The two upscaled quantities of interest are the integrated saturation S and the flux F. For each of these upscaled quantities {S, F}, the transition behaviour is given in terms of the mean saturation (2.13)  Figure 10: Numerical solution of the evolution of the gravity current (2.16), accounting for transition behaviour between viscous and capillary limits using composite expressions (3.18) for the upscaled flow properties, where the capillary number is given implicitly by (3.19) (with Bo= 1). The gravity current shape, shaded to illustrate the saturation distribution (using the same colour scale as in figures 5,6), is illustrated in (a,b,c), whereas the local capillary number N c is illustrated in (d,e,f). For all plots we shade regions with capillary number larger than the transition value N c >N ct = 394 in green, and regions with N c <N ct in blue. We also illustrate one folding scale ∆ on either side of N ct with dashed lines in (d,e,f). The evolution of the gravity current nose position r N (t) is shown on a log-log plot in (g). and capillary number (2.20) by the formula where {S, F} ± = {S, F} visc ± {S, F} cap is given in terms of the viscous limit (homogeneous) and capillary limit (heterogeneous) upscaled properties derived earlier, and N ct , ∆, are two fitting parameters which represent the transition capillary number and folding scale, respectively. The precise values of these fitting parameters depend on numerous factors, including the type of heterogeneity and the specific power laws used, but here we use the same values as calculated by Benham et al. (2021), which are N ct = 394, ∆ = 5.5, and these were shown to give good comparison with numerical simulations across a large range of capillary numbers with mean relative error of ∼ 1%. By accounting for the dependence of the upscaled flux on the capillary number, the system loses its self-similarity. In mathematical terms, this can be seen by noticing that the capillary number (2.20) in the governing equation (2.16) contains derivatives of S with respect to r, compromising the self-similar structure of the equations. Therefore, since we are no longer able to convert to a single governing similarity ODE, we must instead solve the full PDE (2.16) numerically, in conjunction with the algebraic equation for the capillary number (2.20). Written in dimensionless terms, we observe that the capillary number is related to the Bond number according to For the sake of the calculations in this section, we use a mid-range value of Bo= 1 for the Bond number, since this represents an even balance between capillary and gravity forces.
In figure 10 we display the numerical solution in the case of sedimentary strata (k low /k high = 1/3, H low /H high = 1), plotted at three different times after injection, illustrating the evolution of the gravity current, as well as the spatial variation in the local capillary number (3.19). At early times the majority of the gravity current lies within the viscous limit, except the very tip of the gravity current nose. Indeed, the nose of the gravity current lies in the capillary limit for all times, since the flux vanishes there (since N c ∼ h∂p n /∂r ∼ h∂h/∂r → 0 as r → r N (t) due to (2.24)). At later times, a substantial portion of the gravity current lies within the capillary limit, whilst the viscous-capillary transition (N c =N ct ) tends to a constant position of r/H ≈ 0.2. The position of the nose of the gravity current r N (t) is plotted on a logarithmic scale in figure 10g, also indicating the viscous and capillary limit curves.
We make several observations. Firstly, we note the dynamic transition between viscouslike behaviour at early times and capillary like-behaviour at late times. This is evident from figure 10g, where we see that the nose position initially grows like ∼ t 1/2 with prefactor corresponding to the viscous limit, and later grows like ∼ t 1/2 with a capillary limit prefactor. There is a transition period at intermediate times where the behaviour does not obey a power law, as observed by Benham et al. (2021) in the absence of gravity. This can be explained by the fact that, as the gravity current spreads out radially, the capillary number reduces everywhere except near the origin (where N c → ∞). In this way, the composite expressions (3.18) in the majority of the gravity current switch between viscous limit behaviour to capillary limit behaviour over time. Secondly, we note that since the viscous-capillary transition (N c =N ct ) tends to a constant position, if we continue the simulation indefinitely, the fraction of the gravity current in the capillary limit tends towards unity, indicating that at late times a capillary limit model is a good approximation for the whole aquifer.

Discussion of applications to CO 2 sequestration
In this section we discuss the implications of the upscaled description of gravity currents in the context of CO 2 sequestration. There are many different sites that we could choose as case studies, but probably the one with the most available data is at the Sleipner project in the North Sea (Bickle et al. 2007). In particular, we focus on the top layer, which has been investigated by numerous authors (Chadwick et al. 2009;Williams & Chadwick 2017). Whilst various attempts have been made to model the CO 2 migration at Sleipner, and comparisons have been made with seismic measurements of the extent of the plume (Bickle et al. 2007;Golding et al. 2013;Cowton et al. 2016Cowton et al. , 2018, these attempts have often fitted certain parameters (e.g. the permeability) to match observations, without accounting for the possible effect of heterogeneities. Therefore, here we use our previous upscaling procedure to perform this analysis, calculating the effect that different types of heterogeneities could have had on the speed of the gravity current.
Since little is known about the geological heterogeneities at Sleipner, we investigate all the different types of layering we studied earlier, and we account for how uncertainty in such layering may propagate to uncertainty in the modelling predictions.
As described earlier in Section 2.3, there are 8 non-dimensional parameters in our model (excluding the capillary number). Of these parameters, 5 relate to type of heterogeneities, about which little is known for Sleipner. Therefore, instead we take the following values from other similar studies: the porosity-permeability power law used by Bickle et al.
(2017) a = 0.14; the pore-entry pressure power law taken from the scaling proposed by Leverett (1941) b = 0.43; the permeability ratio taken from Bickle et al. (2017) k low /k high = 0.01, no compaction β = 0, and we vary the layer ratio H low /H high between 1 and 10. Since it is unknown whether sedimentary strata, turbidites or a permeability spectrum is most appropriate, we investigate all of these different heterogeneity types.
There are 2 parameters relating to the multiphase flow properties λ and α, which define the capillary pressure and relative permeability. We note that not all relative permeability curves are parameterised as simply as (2.5). For example, other more nonlinear functional forms have been proposed by Chierici (1984). However, we note that Chadwick et al. (2009);Williams & Chadwick (2017) use the Brooks-Corey formulation to model Sleipner, which is given by (4.1) We could change our formulation to account for this modified parameterisation, but instead we notice that (4.1) can be approximated by (2.5) with mean relative error 3% using a power law value of α = 2.32. Therefore, we stick with our original formulation (2.5) for the sake of consistency, without sacrificing much accuracy. Meanwhile, the pore size distribution is estimated by Chadwick et al. (2009) as λ = 2/3. The final parameter we need to describe the problem is the Bond number, which is defined by (2.26) in terms of 6 other physical parameters (excluding gravity, g = 9.81 m/s 2 ). For the top layer at Sleipner, Williams & Chadwick (2017) estimate the temperature between 28-31 • C and pressure between 8.2-8.9 MPa, which gives a density difference of ∆ρ = 232-309 kg/m 3 . Meanwhile the viscosity of CO 2 is taken as µ n = 54.7-65.5 × 10 −6 Pa·s. The input flux is best estimated by Cowton et al. (2016), which for the first few years of injection is Q = 1.5-3 × 10 −3 m 3 /s. The mean permeability is estimated by Bickle et al. (2007) as k 0 = 1.1 − 5 × 10 −12 m 2 . Finally, the pore entry pressure and end-point relative permeability are given by Williams & Chadwick (2017) as p 0 = 1.3 kPa and k rn0 = 0.28. Putting these all together, we estimate the Bond number as Bo= 8.9 − 35.5. Now that we have determined all the parameter values (up to a given uncertainty), we follow the procedure outlined earlier to calculate the prefactor η N for the gravity current using the various different heterogeneity types. In this way, we can measure the heterogeneity efficiency ν (3.1), which tells us how much we can expect the heterogeneities to modify the speed of propagation. In the low/high Bond number estimate Bo=8.9/35.5, we find ν = 36%/31% for equally distributed sedimentary strata (H low /H high = 1), ν = 147%/108% for sedimentary strata with thin streaks of high permeability (H low /H high = 10), ν = 11%/6% for turbidites, and ν = 23 ± 11%/9 ± 8% for a permeability spectrum (mean value). In the latter case, we include the standard deviation values to indicate how these predictions vary due to the uncertainty of the heterogeneity measurements. Indeed, the large degree of uncertainty in these predictions not only illustrates the need for more detailed bore hole measurements at Sleipner, but also demonstrates the importance of accounting for such uncertainty in any modelling approach.
We note that the permeability ratio in the Sleipner field may not be as small as the value we have taken from Salt Creek, k low /k high = 0.01 (Bickle et al. 2017). Therefore we also perform the same analysis as above using a permeability ratio ten times larger. We find that the efficiencies ν calculated for k low /k high = 0.1 are between 1/5 and 3/5 of their values for k low /k high = 0.01. This indicates that, even if we have vastly overestimated the permeability ratio, the effect of heterogeneities is still likely to be significant.
From this analysis, it is clear that the possible effect that heterogeneities may have had on the injection of CO 2 at Sleipner largely depends on the type of heterogeneities present. In particular, thin sedimentary strata with very high permeability could have caused a potential speedup of more than 100%. However, for more moderate permeability profiles, such as evenly distributed strata or turbidites, these heterogeneities may have only caused a 10 − 30% speedup.

Concluding remarks
We have studied the upscaled effect of several different types of heterogeneity on the evolution of an axisymmetric gravity current under an impermeable cap rock. The four heterogeneity types considered, which were all given in terms of vertical variations in the rock properties, were each motivated by different physical mechanisms for the nonuniform deposition or compaction of sediments. We developed a general method for calculating the gravity current shape and growth rate in either the viscous or capillary limits, which involves solving a single ordinary differential equation that depends on an upscaled flux term evaluated either via numerical integration, or using analytical expressions which we derived in certain limiting cases. This simplified approach not only reduces the computation time significantly, but also provides key insights into the role of small-scale heterogeneities on the propagation of the large-scale flow.
In particular, we showed that heterogeneities have the ability to speed up the growth rate of the gravity current by means of a nonlinear focussing mechanism (via the relative permeability and capillary pressure) into high permeability layers. The degree of speedup depends on the type of heterogeneity, and most importantly the fraction of high/low permeability regions within the medium. The largest effect was seen in the case of sedimentary strata with thin regions of high permeability. Using a permeability profile composed of randomly sampled layers, we demonstrated how uncertainty in heterogeneity measurements can propagate to uncertainty in field predictions, an effect which is particularly pronounced for small values of the Bond number. We also investigated modelling the transition from the viscous limit regime to the capillary limit regime, shedding light into which regions of the gravity current are most affected by heterogeneities, and when.
The main motivation for this study was to create an upscaled description of how smallscale heterogeneities affect large-scale CO 2 migration, for safe and efficient sequestration in porous aquifers. To assess the risks associated with any CO 2 storage scheme, examining the effect of heterogeneities is essential. To illustrate this, we used the case study of CO 2 injection at the Sleipner project in the North Sea. In this case, for realistic parameter values, we demonstrated that heterogeneities may have sped up the gravity current by more than 100% during injection. However, we also illustrated that this figure depends greatly on the heterogeneity type, indicating the need for detailed core measurements from bore holes.
For future work, variations in the heterogeneities along the length of the aquifer could be studied (in addition to the vertical heterogeneities investigated here), similarly to Jackson & Krevor (2020). In such cases, the upscaled flow properties would depend on both the horizontal and vertical correlation length scales of the heterogeneities. In addition, using the upscaled results presented here, predictions of trapping efficiencies could be calculated for various sequestration sites. This could shed light onto which aquifers have heterogeneities that could potentially enhance their capability for CO 2 storage. Furthermore, as we discussed earlier, it would be interesting to investigate the evolution of the gravity current after injection has stopped (see figure 9). However, in this case both imbibition and drainage relative permeability/capillary pressure curves must be considered, depending on whether the CO 2 is invading or being withdrawn from different regions of the aquifer, similarly to Golding et al. (2017).
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 CO 2 through North Sea Geological Carbon Storage Sites" (grant no. NE/N016084/1).
Declaration of Interests. The authors report no conflict of interest.
Code written to numerically evaluate the flux integrals and calculate the similarity solution can be found on the personal website of G.P. Benham: https://yakari. polytechnique.fr/people/benham/gravity_current/upscale.m

Appendix A. Extra plots
In this Appendix we present extra plots in figure 11 for the integrated saturation S (2.15) and the Darcy velocity of the non-wetting phase u n ∝ ∆ρgk(z)k rn (s)/µ n at different Bond numbers and for different types of heterogeneity.
The relationship between the integrated saturation S and the gravity thickness h, as shown in figure 11a,c,e, is useful for understanding how to invert the solution of the governing equation (2.16). In all cases S has approximate power law dependence on h, where the power is between linear and cubic, as illustrated with dotted lines.
The velocity profiles in figure 11b,d,f, are useful for understanding the form of the flux function F (which is the integrated velocity), as displayed in figure 3a,b,c. We present the scaled velocity U = (u n (0) − u n (z))/(u n (0) − u n (h)) so that U varies between U = 0 at z = 0 and U = 1 at z = h. In this way we can see the functional form of U more clearly. For example, if U is constant (as in figure 11f), when integrated this will give a linear dependence between the flux F and the gravity current thickness h.

Appendix B. Derivation of the governing equation (2.16)
In this Appendix we provide the details for the derivation of the governing equation for the integrated saturation S, which is given by (2.16). We start by considering the governing equation for the gravity current thickness (2.14). To re-write this in terms of the the integrated saturation (2.15), we first need to transform the derivatives of h into derivatives of S. As illustrated by figure 11a,c,e, the integrated saturation is always a monotone increasing function of the gravity current thickness, so that we can define a unique inverse function h = S −1 (S).
(B 1) Then, using the chain rule, the gradient is given by In each plot we indicate power law values of 1/2, 1, 2 and 3 with dotted lines for comparison. (b,d,f) Corresponding scaled velocity profiles U = (u n (0) − u n (z))/(u n (0) − u n (h)) where u n ∝ ∆ρgk(z)k rn (s)/µ n . We plot the ensemble average of the velocity, rather than the velocity within each layer, so as not to display oscillatory behaviour between layers.
We use the inverse derivative identity to calculate derivatives of (B 1), such that