Two-phase gravity currents in layered porous media
TThis draft was prepared using the LaTeX style file belonging to the Journal of Fluid Mechanics Two-phase gravity currents in layeredporous media
Graham P. Benham † , Mike J. Bickle , Jerome A. Neufeld Department of Earth Sciences, University of Cambridge, Cambridge CB3 0EZ, UK BP Institute, University of Cambridge, Cambridge CB3 0EZ, UK Department of Applied Mathematics and Theoretical Physics, University of Cambridge,Cambridge CB3 0WA, UK(Received xx; revised xx; accepted xx)
We examine the effects of horizontally layered heterogeneities on the spreading of two-phase gravity currents in a porous medium, with application to numerous environmentalflows, most notably geological carbon sequestration. Geological heterogeneities, whichare omnipresent within all reservoirs, affect the large-scale propagation of such flowsthrough the action of small-scale capillary forces, yet the relationship between thesesmall and large scales is poorly understood. Here, we derive a simple upscaled modelfor a gravity current under an impermeable cap rock, which we use to investigate theeffect of a wide range of heterogeneities on the macroscopic flow. By parameterising interms of different types of archetypal layering, we assess the sensitivity of the gravitycurrent to the distribution and magnitude of these heterogeneities. Furthermore, sincefield measurements of heterogeneities are often sparse or incomplete, we quantify howuncertainty in such measurements is manifest as uncertainty in the macroscale flowpredictions. 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 CO .
1. Introduction
Injection of CO into underground reservoirs to reduce greenhouse gas emissions, alsoknown as geological carbon sequestration, is one of the major proposed technologicalsolutions to meet future global temperature targets (Bickle 2009; Bui et al. 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 saturationdistribution and consequent flow properties via the relative permeability and capillarypressure (Nordbotten & Celia 2011). Heterogeneities in rock properties at the 10 −
100 cmscales substantially amplify and complicate the effect of variations in pore-scale capillaryforces, and are manifest in the large-scale saturation distributions within the CO current.Hence, in order to ensure safe and efficient sequestration, it is imperative to be able tomodel how small-scale heterogeneities, which are ubiquitous in all subsurface reservoirs,affect spreading rates at the macroscale (Benham et al. et al. et al. et al. † Email address for correspondence: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] F e b type of heterogeneity (e.g. over some horizontal/vertical length scale, as investigated byJackson & Krevor (2020)). In particular, it is not currently understood how generic small-scale heterogeneities affect the propagation of such large-scale flows. Furthermore, sinceheterogeneities are usually measured through isolated and sparsely distributed bore holesamples, such measurements often come with a large degree of uncertainty, yet there arevery few studies which discuss how such uncertainty at small-scales translates to large-scale predictions. Indeed, whilst the studies of Hinton & Woods (2020 a , b ) investigatedhow variations in permeability cause shear dispersion for miscible flows, the correspondingeffects due to capillary forces within immiscible flows have not yet been addressed.However, as discussed by Jackson & Krevor (2020), these capillary forces associatedwith the heterogeneities play a critical role during CO migration, and therefore requirecareful modelling.In this study, we derive a simple upscaled model for the evolution of an axisymmetrictwo-phase gravity current beneath an impermeable cap rock, where the structure anddistribution of layered heterogeneities is a model input. This simplified approach notonly greatly reduces the computational demand of modelling such small-scale details,but also allows us to study a large range of heterogeneities by parameterising themas different types of archetypal sedimentary layering. In this way, we assess how theproperties of the heterogeneities affect the macroscale flow, as well as how uncertaintyin the measurements is manifest in such flow predictions. Our simple model providesother key insights, such as the self-similar nature of the upscaled gravity current, scalinglaws for the speed of propagation and an understanding of where and when the flowtransitions between the so-called capillary and viscous limiting behaviours.In modelling subsurface migration of CO , one of the key difficulties is resolvingthe complex heterogeneous rock properties of the porous medium. This is a two-foldchallenge, since on the one hand a resolution which captures all the details of the rockgeometry is very computationally intensive, and on the other hand it is almost impossibleto obtain such fine-scaled resolution over field scales from field measurements. Hence,there is a strong motivation for an upscaled modelling approach which describes thebulk, or average effect of heterogeneities on the large-scale flow features. There are manypossible 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 sizeof the aquifer. Therefore, in the context of this study small-scale heterogeneities referto variations on the scale of 10 − - 1 m (as opposed to pore-scale heterogeneities whichoccur on the scale of 10 − - 10 − m), and the large-scale flow refers to a gravity currentwhich is typically around 1 - 10 m thick (Cowton et al. et al. capillary limit , the capillary forces due to heterogeneities dominate the flowbehaviour, while in the limit of large capillary number, known as the viscous limit , theyhave 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 capillaryregimes, demonstrating which regions of a confined aquifer are in each of these limits, andwhich regions are in between the limits. However, gravity was neglected in that study,restricting the applicability to very thin aquifers.For CO sequestration sites in large aquifers, gravity plays a dominant role in therise 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 gravitationalforces and capillary forces, and may be characterised by a dimensionless Bond number(Golding et al. et al. (2021), the Bond number is greater than unity for aquifers larger than around ∼ et al. et al. beneath an impermeable caprock. The low computational cost of our simple approach allows us to explore differentparameterisations of heterogeneities, providing insights into the dominant controls on theevolution of the gravity current. Similar to Benham et al. (2021) (though focussing ongravitational effects), we investigate both the viscous limit, the capillary limit and thetransition between these limits using a locally defined capillary number that determineswhere and when heterogeneities play an important role. We show that away from thistransition zone the upscaled gravity current is self-similar, where the front moves like thesquare root of time (like the homogeneous case discussed by Golding et al. (2013)) and theprefactor 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 realpermeability data with uncertainty in the measurements, illustrating how this uncertaintyis manifest in modelling predictions. Finally, we use our upscaled approach to investigatehow heterogeneities may have affected the injection of CO at the Sleipner site in theNorth Sea (Bickle et al. HeterogeneousHomogeneous rrzz zzhh ssρ n , µ n ρ w , µ w φ ( z ) , k ( z ) , p e ( z ) Q s = 0 h ( r, t )(a) (b)(c) (d)(e) Figure 1: Schematic diagram of an axisymmetric gravity current (with constant injection Q ) in both the homogeneous case (a) and the heterogeneous case (with sedimentarystrata) (c), also illustrating the corresponding vertical non-wetting saturation profiles(b,d), given by (2.10),(2.12) (note the heterogeneity wavelength is exaggerated forillustration purposes). (e) Relationship between mean non-wetting saturation ¯ s (2.13)and gravity current thickness h .
2. Upscaled modelling of two-phase gravity currents
In this section, we outline the assumptions used to model an upscaled two-phase gravitycurrent in a heterogeneous porous medium, making note of how the saturation of phasesvaries within the current. Then, we derive the upscaled governing equations and boundaryconditions which describe the macroscopic dynamics. Subsequently, we discuss a varietyof different types of heterogeneity and how these are manifest in the upscaled properties.Finally, despite the added complexity of the heterogeneities, we demonstrate the self-similar nature of the gravity current, thereby greatly reducing the complexity of theproblem. 2.1.
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 ) is injected at a point source at the originwith flow rate Q into a surrounding porous medium saturated with a denser wettingphase (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 thedominant heterogeneity arising from sedimentary layering, we consider a porous mediumwhich 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 φ ( z ) ∂S i ∂t + ∇ · u i = 0 , i = n, w, (2.1) u i = − k ( z ) k ri ( S i ) µ i ∇ ( p i − ρ i g z ) , i = n, w, (2.2)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 relativepermeabilities 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 phasessatisfies p n − p w = p c ( S n ) , (2.3)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 morecomplex dependencies). These are usually approximated with empirical parameterisedformulae, 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 p c = p e ( z )(1 − s ) − /λ , (2.4) k rn = k rn s α , (2.5)where p e ( z ) is the pore entry pressure, s = S n / (1 − S w ) is the reduced saturation ofnon-wetting phase, λ represents the pore size distribution, k rn is the end-point relativepermeability, and α is a fitting parameter. The irreducible wetting saturation S w is theamount of wetting phase that is permanently stored in the pores during drainage flows,and consequently the end-point relative permeability corresponds to k rn = k rn (1). Inthis 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 allowthe non-wetting phase to enter the largest pore spaces at a given position. Likewise, asthe pressure difference between phases increases, the non-wetting phase is able to entersmaller and smaller pore spaces. Hence, clearly the pore entry pressure depends on thesize and geometry of the pores (and hence varies vertically), and the same is true forthe permeability and porosity. However, whilst this dependence has been measured forspecific rock types, it is not fully understood in general. Hence, as is often done, we usepower laws to relate these different quantities, such that φ ∝ k a , (2.6) p e ∝ k − b , (2.7)for some constants a, b >
0, which we take to be positive since large pore size correspondsto 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 differentrock types. Therefore, we keep them in general form for this analysis. However, we notea commonly used scaling proposed by Leverett (1941), p e ∼ ( φ/k ) / which implies b = 1 / − a ).Motivated by field observations of gravity currents (e.g. see Cowton et al. (2016), wherethe aspect ratio of the gravity current at Sleipner was calculated to be less than ∼ / et al. (2011, 2013), we assume that the gravity current is long andthin, such that the vertical velocity is much smaller than the horizontal velocity w i (cid:28) u i .In this case, the pressure within each phase is approximately hydrostatic ∂p i ∂z = ρ i g, i = n, w, (2.8)and consequently, (2.8) is integrated to match the capillary pressure (2.3), giving p c = − ∆ρg ( z − h ) + p , (2.9)where p 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 lowerbound on s , such that s = max (cid:40) − (cid:20) p p e ( z ) + ∆ρg ( h − z ) p e ( z ) (cid:21) − λ , (cid:41) . (2.10)To determine the value of p we consider that the edge of the gravity current is definedas 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 = min p e ( z ).In other words, by setting the capillary pressure at the edge of the gravity current asthe smallest required pressure difference to invade any pores in the aquifer, we guaranteethat anywhere below the edge of the gravity current ( p c < p ) no saturation will befound. Therefore, even though there may be disconnected regions of non-wetting phasewithin 0 (cid:54) z (cid:54) 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 (dueto heterogeneities) and gravitational forces. However, this is only valid for situationswhere capillary forces are large enough to drive the saturation into regions of largerpore space, or equivalently when the capillary number is small. Therefore, in generalthe saturation distribution must depend on the local capillary number N c , which isgiven as the ratio between the horizontal flow-driving pressure gradient and the typicalvertical gradient in pore entry pressure (Benham et al. ∂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 acrossthe aquifer (constant), and the gravity current thickness h is used as the vertical lengthscale. Hence, the capillary number is given byN c = (cid:12)(cid:12)(cid:12)(cid:12) h∆p e ∂p n ∂r (cid:12)(cid:12)(cid:12)(cid:12) . (2.11)In the limit of very small capillary number N c (cid:28)
1, also known as the capillary limit ,the saturation distribution (2.10) remains accurate. However, when the capillary numberis very large, also known as the viscous limit , capillary forces due to heterogeneities areeffectively negligible (i.e. we can ignore pore entry pressure variations p e ( z ) = p ), andthe saturation distribution becomes s = 1 − (cid:20) ∆ρg ( h − z ) p (cid:21) − λ , (2.12)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 mathematicallyidentical to the homogeneous case, the model would still have to account for flow variations dueto permeability gradients through an effective permeability. The saturation would, however, beidentical to the homogeneous case.
For intermediate capillary number (i.e. when the flow is neither in the viscous limitnor 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 islogarithmic (Benham et al. et al. (2013)), which is equivalent to the upscaledviscous limit (N c (cid:29)
1) for a heterogeneous medium (see earlier footnote). For each casewe plot typical vertical saturation profiles in figure 1b,d. In the homogeneous case thesaturation distribution satisfies a balance between capillary and buoyancy forces, so thatlighter regions (of high saturation) are pushed towards the cap rock. In the heterogeneouscase, the same overall balance is sustained, but within that balance capillary forcespush the saturation into layers where the pores are larger. Hence, significant oscillatorybehaviour is observed within the vertical saturation profile, including patches where thesaturation drops to zero. This corresponds with regions where the pore spaces are toosmall to allow any non-wetting phase (i.e. the zero value is chosen in (2.10)). Hence, oneinteresting consequence of heterogeneities is that they modify the mean saturation valuein the gravity current. In figure 1e we plot the mean saturation, defined as¯ s ( h ) = 1 h ˆ h s d z, (2.13)whilst varying the gravity current thickness h for both the homogeneous and heteroge-neous cases. In both cases ¯ s is an increasing function of h , but the heterogeneous casealways has a lower mean value. This is due to the substantial fraction of the gravitycurrent with zero saturation.2.2. Upscaled model: Governing equation and boundary conditions
Having discussed the flow scenario and laid down the key assumptions, now weoutline the upscaling procedure, deriving a single governing equation and accompanyingboundary 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 massequation for the non-wetting phase (2.1) is integrated between z = 0 and z = h ( r, t ),such that ϕ ∂∂t ˆ h ˆ φ ( z ) s d z − u b r ∂∂r (cid:32) rk k rn ∂h∂r ˆ h k ( z ) k rn ( s ) d z (cid:33) = 0 , (2.14)where ϕ = (1 − S w ) φ is the reduced porosity scaling, φ and k are typical scalingsfor the porosity and permeability (where ˆ φ = φ/φ ), and u b = k k rn ∆ρg/µ n is thebuoyancy velocity.We note that the flow of wetting phase is ignored in this analysis under the long-thin approximation. Essentially, the flow of non-wetting phase within the gravity currentdecouples from the flow of wetting phase, which is not present at leading order. Neverthe-less, multiphase effects are still manifest at leading order via the multiphase properties,such as the relative permeability and capillary pressure. However, as we discuss laterin Section 2.3, this assumption breaks down if the permeability ratio between layersbecomes very small. In particular, if there are regions of very low permeability, these willact 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 playsan important role on the dynamics and therefore can no longer be ignored. We give moredetails of this consideration in the next section.We also note that (2.14) is already an upscaled description of the flow, since theheterogeneities are only manifest within the integrals. Hence, (2.14) represents how theheterogeneities 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 theheterogeneities, both to reduce computation time, and also because realistically the lowresolution 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 renderit amenable to conventional analysis. Therefore, by switching variables to the integratedsaturation S ( h, N c ), which is defined as S = ϕ ˆ h ˆ φ ( z ) s d z, (2.15)(2.14) is rewritten as ∂ S ∂t = 1 r ∂∂r (cid:20) r F ( S , N c ) ∂ S ∂r (cid:21) , (2.16)where the flux is given by F = u b k k rn (cid:20) K ( h, N c ) S h ( h, N c ) (cid:21) , (2.17)and the two functions K ( h, N c ) and S h ( h, N c ) are defined as K = ˆ h k ( z ) k rn ( s ) d z, (2.18) S h = ϕ ˆ h ˆ φ ( z ) ∂s/∂h d z = ∂ S ∂h . (2.19)Further details of this coordinate transformation are presented in Appendix B. We notethat 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-wettingphase (per unit area), where the flux is a non-linear function that represents how capillaryforces modify the flow. Hence, there is an interesting analogy between our scenario anda viscous gravity current, where the flux function represents how viscous forces modifythe flow (e.g. plug flow, Poiseuille flow, etc...). As we will find out later, F is sometimeswell-approximated by a power law of S , and the solutions to such equations are detailedby 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 thetranscendental equation for the capillary numberN c = ∆ρg∆p e (cid:12)(cid:12)(cid:12)(cid:12) h ( S , N c ) S h ( h ( S , N c ) , N c ) ∂ S ∂r (cid:12)(cid:12)(cid:12)(cid:12) , (2.20)where h is written in terms of S under the assumption that (2.15) has a uniquely-definedinverse (which we later find to be the case).For the remainder of this study (up until Section 3.3) we restrict our attention to thetwo 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 ), H low /H high = 0 . H low /H high = 1 H low /H high = 10 (a) (b) (c) rz rz rz HighLow z k (d) z k (e)
PDFlog k z k (f) z k (g) Sedimentary strata ( H low /H high = 1 ) Turbidites CompactedSpectrum
Figure 2: Illustrations of the different types of heterogeneity we consider, wherethe 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 caseof sedimentary strata with greyscale permeability maps for three different values of thewidth 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 randomlysampled from a uniform distribution on a logarithmic scale.thereby decoupling (2.16) and (2.20). However, later in Section 3.3 we address the caseof intermediate capillary number, for which the equations must be solved in tandem.The governing equation (2.16) must be accompanied by some initial and boundaryconditions to create a well-posed system. Firstly, we define the nose of the gravity currentat 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 thenon-wetting phase, such that 2 π ˆ r N ( t )0 r S d r = Qt, (2.22)or equivalently, we impose the input flux at the origin and zero flux through the nose ofthe current, − π (cid:20) r F ∂ S ∂r (cid:21) r =0 = Q, (2.23)2 π (cid:20) r F ∂ S ∂r (cid:21) r = r N ( t ) = 0 . (2.24)The finite flux value Q in (2.23) indicates that the gradient ∂ S /∂r must become infinite0 Heterogeneity type Functional form ParametersSedimentary strata k = (cid:40) k low k high with ratio H low /H high k low /k high , H low /H high Turbidites k = k (cid:2) − (2 A/π ) tan − cot (2 nπz ) (cid:3) A = − k low /k high k low /k high Spectrum k = exp ( U [log k low , log k high ]) k low /k high Compacted k = k (1 + z ) − β β (cid:62) Table 1: Definitions of the different types of heterogeneity (characterised by thepermeability), 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 linearlyvarying permeability, where the wavenumber n , is considered in the limit nh → ∞ . In thespectrum case, permeability is a series of strata, where each layer has permeability takenfrom 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 alinear scale. Compacted rock corresponds to a permeability profile which decreases withdepth 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 maybecome 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.2.3. Incorporating heterogeneity
To close the system, we must choose a type of vertical heterogeneity. Since we haveused 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 varyin three dimensions (i.e. k ( x )), here we restrict our attention to pure vertical variationsince, 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 currentmade earlier.To model the permeability we have a variety of different physically-motivated choiceswhich we list in table 1 and plot in figure 2. Firstly, sedimentary strata representa periodic deposition of two different types of sediments, such that the permeabilityalternates between two values, k low and k high , in a periodic array of layers, where thewidth ratio of each of these is given by H low /H high (see figure 2a,b,c). Unlike thesedimentary strata, which are uniformly deposited in each layer, turbidites representthe deposition of sediments from the continuous arrival of turbidity currents, such thatwithin each layer the permeability varies linearly. The sign of the linear slope indicatesthat layers become more permeable as one descends deeper, since this corresponds tothe early/late arrival of large/small particles in a turbidity current. Thirdly, we considera permeability profile which is generated by randomly sampling from a distribution, or spectrum , of permeability values, spread out logarithmically. This case is motivated by1realistic measurements of sedimentary strata which are often noisy and aperiodic. Finally,we consider a compacted rock, where the permeability decreases with depth due to thebuildup of lithostatic pressure over time.Although there are many other possible choices for the permeability, we restrict ourinvestigation to these four examples since they are canonical cases from which we maylearn about the fundamental effects of heterogeneities. Each case is parameterised, eitherby 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 strengthof the compaction effect.It is important to note the possible limitations on these parameters. In particular, suffi-ciently low permeability layers may cause a vertical obstruction, such that an unconfineddescription of the gravity current is no longer applicable. To investigate the limitations onthe permeability ratio k low /k high we have performed a set of numerical simulations of thetwo-dimensional miscible Darcy equations using the DarcyLite finite element package in
Matlab , adapted to account for gravity (Liu et al. et al. k rn , k rw →
1, and thephase pressures equalise such that p c →
0. Studying the miscible flow problem allowsus to investigate the applicability of upscaling for small values of the permeability ratiowithout accounting for the more complex effects of immiscible phase saturations. We donot display the numerical results here, since a rigorous analysis of this query is outsidethe 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 . k low /k high for an upscaling procedure,but rather this becomes a question of temporal and spatial scales. In other words, for anypositive permeability ratio, given enough time and spatial extent, such an injection willeventually 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. longerand thinner) due to the delay in the vertical flow caused by the obstruction of the lowpermeability layers. Therefore, to avoid dealing with the prolonged transient effects thatprecede self-similarity in the case of very low permeability ratios, for the remainder ofthis study we restrict our attention to k low /k high (cid:62) . 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 valuesof 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 onthe heterogeneity wavelength, since it is undesirable to have upscaled properties like F that oscillate depending on the gravity current thickness. Therefore, instead of resolvingall of the details of the flow, we build a macroscopic picture of the gravity current, whichis 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 of2two possible values, integrals can be simply decomposed into bulk fractions ˆ · dz = H low ´ k low · dz + H high ´ k high · dz H low + H high , (2.25)thereby removing the need to resolve individual layers. A similar approach can be takenin the case of the permeability spectrum, although in that case (2.25) is replaced by a sumover 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 definedin table 1), the integrals must be calculated with a fine numerical scheme for a large butfinite value of nh (cid:29) S ( h ) and theflux F ( h ), since F allows us to solve the diffusion equation (2.16), and S allows us tocalculate the gravity current thickness by way of inversion. These both depend on anumber of non-dimensional parameters. Ignoring the capillary number (since for nowwe restrict our attention to N c (cid:28) c (cid:29)
1) there are a total of 8 non-dimensionalparameters 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 andpore entry pressure to the permeability a , b ; the Brooks-Corey parameters λ , α ; andfinally the Bond number, which is defined asBo = (cid:18) Q∆ρgµ n k k rn p (cid:19) / . (2.26)The Bond number, which can alternatively be written as Bo= ∆ρgH/p , where H = (cid:112) Q/u b is the buoyancy length scale, is interpreted as the ratio between buoyancyforces and capillary forces. This quantity largely controls the saturation distribution(2.10),(2.12), which is evident upon dimensional analysis. For example, when Bo (cid:29) s = 1 − (cid:20) p e ( z ) /p + Bo( h/H − z/H ) p e ( z ) /p (cid:21) − λ ≈ , (2.27)regardless of which type of heterogeneity p e ( z ) we consider (so long as Bo p /p e (cid:29) a, b were already addressed by Benham et al. (2021) andthe Brooks-Corey parameter λ was studied by Golding et al. (2013). Therefore, for theremainder 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 thehomogeneous case k low /k high = 1 as a proxy to study the viscous limit behaviour N c (cid:29) k low /k high <
1, however, weassume capillary limit behaviour N c (cid:28) a = 1 / b = 3 / λ = 3 and α = 4, which we have extracted from a variety of different sources(Golding et al. et al. et al. † Note that in the case of the permeability spectrum, we sample N pairs of values { k i , H i } ( i = 1 , . . . , N ) from a random distribution of permeabilities and layer widths. Once sampled, itdoes not matter how these values are arranged. Therefore, a bulk decomposition like (2.25) isstill possible. (a) Bo= 10 − (b) Bo= 1 (c) (d)
Bo= 10 Figure 3: (a,b,c) Variation of the flux F (2.17) of the integrated saturation S in (2.16)for different values of the Bond number Bo. Both F and S are normalised by referencevalues (measured at twice the mid range value of the gravity current thickness, h half = h ( r N ( t ) / , t )) for illustration purposes. In each plot we indicate power law values of 1 /
2, 1and 2 with dotted lines for comparison. (d) Analogy between a two-phase gravity currentin a heterogeneous porous medium and a non-Newtonian viscous gravity current withviscosity power law µ ∝ ( ∂u/∂z ) κ . The resultant flux power law is given by ´ h u d z ∝ h / (1+ κ ) , as indicated with the blue curve. Red dashed lines illustrate particular powerlaw values of interest.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 ofthe layered cases, we use a permeability ratio value of k low /k high = 0 .
1, whereas in thecompacted case we use a power law of β = 1. In all cases (except the spectrum case) theflux is well approximated by a power law F ∝ S ψ , for some value of ψ between 1 / S , as well as the velocitydistribution u n ∝ ∆ρgk ( z ) k rn ( s ) /µ n within the gravity current. There are severalinteresting observations to make. Firstly, no matter which type of heterogeneity norwhich Bond number we choose, the integrated saturation S ( h ) is always a monotone4 hz k hz s hz u n F h h r (a) (b) (c) (d) (e) Heterogeneity Saturation Velocity Flux Thickness
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 illustratedas 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 − ( 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 fora classic viscous gravity current, the flux power law relates to the velocity distributionwithin that current. For example, the velocity distribution for Poiseuille flow, which variesquadratically 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 consideringa shear thinning/thickening power law viscosity µ ∝ ( ∂u/∂z ) κ . Then, it is easy to showthat the lubrication balance µ∂ u/∂z ∼ ∂p/∂x can be integrated to give a flux F = ´ h u d z which obeys the power law F ∝ h / (1+ κ ) . This is illustrated in figure 3d,indicating specific cases with dashed lines. For example, a shear-thinning fluid with powerlaw κ = − / F ∝ h / . Whilst our upscalingproblem is very different from a non-Newtonian viscous gravity current, the analogy isnevertheless useful in helping to relate the flux functions observed in figures 3a,b,c, to thevelocity 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 methodol-ogy for analysing the gravity current in figure 4. This illustrates the steps between initiallychoosing a heterogeneity type and finally solving for the gravity current thickness h . Wehave provided some example code in the supplemental materials to demonstrate thesesteps in the case of sedimentary strata, including how to numerically calculate the fluxfunctions. 2.4. Discussion of self-similarity and the numerical scheme
There is a final simplification that can be made owing to a coordinate invariance, whichallows calculation of the solution using a simple numerical scheme. In particular, much5like 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 governingequation (2.16) (under the assumption of viscous N c (cid:29) c (cid:28) η = ( ϕ /Qu b ) / rt − / , (2.28) S = Hϕf ( η ) , (2.29)where the nose of the gravity current is located at η = η N for some constant η N which wewill determine shortly. To further simplify the equations, and to convert to a unit intervaldomain, we write our system in terms of the variables y = η/η N and ˆ f ( y ) = f ( η ). Inthis way, the governing equation for the integrated saturation (2.16) and the boundaryconditions (2.21),(2.22) become (cid:104) y ˆ F ( ˆ f ) ˆ f (cid:48) (cid:105) (cid:48) + 12 η N y ˆ f (cid:48) = 0 , (2.30) η N (cid:20) π ˆ y ˆ f ( y ) d y (cid:21) = 1 , (2.31)ˆ f (1) = 0 , (2.32)where ˆ F = F ϕ/u b H . The system can be solved numerically using a simple finitedifference scheme, starting at y = 1 and marching back towards y = 0, where the secondorder ODE (2.30) is conveniently written as a set of two first order ODE’s with boundaryconditions L (cid:48) ( y ) = − η N y L F ( ˆ f ) , (2.33)ˆ f (cid:48) ( y ) = L y ˆ F ( ˆ f ) , (2.34) L (1) = (cid:15), (2.35)ˆ f (1) = (cid:15), (2.36)where L is the total dimensionless flux, and (cid:15) (cid:28) (cid:15) = 0 since it will generate infinite gradients). To find the constant η N , we start with an initial guess η N , and then use Newton’s method to iterativelysolve the flux condition (2.31) (see our code which we have uploaded in the supplementalmaterial).We make the key observation that, independent of the form of ˆ F ( ˆ f ), the gravitycurrent is self-similar, with a nose that moves like the square root of time. Hence, theheterogeneities are only capable of modifying the prefactor η N for the nose speed, notthe power law (which is always r ∼ t / ). However, the heterogeneities may also changethe shape of the gravity current via F and S .As discussed later in Section 3.3, in the case where the capillary number is neithersmall nor large, the flux must depend on the capillary number itself. In this case, sincethe capillary number involves derivatives of S with respect to r , this modifies the formof the governing equation (2.16), rendering such similarity variables inadmissible. Over † Note that the two-dimensional case is not necessarily self-similar. The two-dimensionalgravity current thickness scales like h ∼ t / , such that the flux function F ( h ) cannot be writtenin a general self-similar form. However, this becomes possible in certain specific cases (e.g. alinear power law F ∝ h , as discussed by Huppert & Woods (1995)).
3. 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 casethe upscaled problem is equivalent to the homogeneous case.(ii) The capillary number is small throughout the aquifer (capillary limit), in whichcase the upscaled problem is dominated by the effect of the heterogeneities.(iii) The capillary number varies across the aquifer, such that different regions simulta-neously lie within the viscous limit and the capillary limit, and other regions lie betweenthese limits.We start by addressing the former two cases (viscous and capillary limits), for whichthe 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 tostudy the capillary limit. The gravity current thickness must be calculated numerically bysolving the ODE system (2.31),(2.33)-(2.36). In specific limiting cases, such as when theBond number is very large or very small, analytical progress can be made. We first presentour numerical results which we use to understand the broad effect of the heterogeneitiesacross large parameter regimes. Then, we use asymptotic analysis to help interpret theresults in certain limiting cases. Finally, we address the transition between the viscousand capillary limits, for which the full PDE system (2.16),(2.21),(2.22) must be solvedin tandem with an algebraic equation for the local capillary number (2.20).3.1.
Numerical solution in the viscous and capillary limits
The capillary limit numerical solution for different types of heterogeneities is plottedin figures 5 and 6, for Bond number values between Bo = 10 − and Bo = 10 . Typicalsaturation profiles are also displayed as inserts in each plot. To compare the different pro-files, we have normalised the thickness by twice its mid range value h half = h ( r N ( t ) / , t ).The radial coordinate is normalised by the nose position r N ( t ) so that the shape remainson a fixed domain for all time. For now, we do not include plots of the viscous limitnumerical 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 ofheterogeneities.Let us first focus on the non-compacted cases in the capillary limit (figures 5 and6a,c,e). For small Bond number Bo (cid:28)
1, the saturation becomes near-zero s ≈
0, butwith spikes of linearly increasing magnitude that represent the thin regions where thepermeability is near its maximum value k ≈ k high . As we increase the Bond numbertowards unity, the saturation becomes larger, with an overall curved profile and significantoscillations. At high Bond number, the saturation tends towards s ≈ z = h , where it rapidly drops to s = 0. The shape of the gravity currentchanges from having a sharp nose at high Bond number to having a rounded blunt noseat low Bond number. There is not a noticeable difference in the shape of the gravitycurrent between the different types of heterogeneity.Apart from the shape and the saturation distribution, there are two other importantmetrics which are useful for describing the current. Firstly, the prefactor η N relates to thespeed of the advancing nose, and secondly the mid-range thickness h half = h ( r N ( t ) / , t )indicates the approximate size of the current. Following Golding et al. (2013), we use theclassic single phase limit values as a useful reference. Using a subscript notation, these7 s/s z s/s zs/s z s/s zs/s z s/s z Sedimentary strata ( H low /H high = 1) Turbidites B o = − B o = B o = (a) (b)(c) (d)(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 / , 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 thereference value 2 h half = 2 h ( r N ( t ) / , t ) for the sake of comparison. The heterogeneitywavelength is exaggerated for illustration purposes. In each plot inserts illustrate thevertical saturation profile, normalised by the uppermost value s = s (0).are given by η N = 1 .
155 and h half = 0 . H (Lyle et al. (cid:29) (cid:28) h half behaves similarly for all three layeredcases, growing approximately like h half ∼ Bo − / , as described by Golding et al. (2013) forthe homogeneous case. On the other hand, the prefactor η N behaves rather differently. Inall cases, we see an increase in η N for stronger heterogeneities, indicating that capillaryforces 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 thecomplex nature of the velocity distributions and flux functions depicted in figures 3 and11. It is interesting to note that despite having a permeability profile with the same meanvalue, the different variations within each layer for each heterogeneity type are sufficientto alter the flux of saturation, and hence modify the speed of propagation of the gravitycurrent. This sheds light on both the need for detailed bore hole measurements to infer8 s/s z s/s zs/s z s/s zs/s z s/s z CompactedSpectrum B o = − B o = B o = (a) (b)(c) (d)(e) (f) Figure 6: Numerical results for the capillary limit in the case of spectrum permeability(a,c,e) (with mean permeability ratio k low /k high = 0 .
04) and compacted rock (b,d,f)(with compaction power law β = 1). In each plot inserts illustrate the vertical saturationprofile, normalised by the uppermost value s = s (0). The heterogeneity wavelength isexaggerated for illustration purposes.as much information about the heterogeneities as possible, as well as the usefulness ofsuch an upscaling approach as we have taken here.For each of the different types of strata, we compare the capillary limit curves againstthe homogeneous case (solid black line), which is equivalent to the viscous limit. Thisallows us to quantify the effect of the heterogeneities on the prefactor more clearly. Thestrongest 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 ofthe saturation feeds into the nonlinearity of the relative permeability, thereby ampli-fying 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 wedo not display them here.In the case of the permeability spectrum, we choose a permeability distribution whosestandard deviation divided by the mean is σ ( k ) /k = 1, and whose mean permeabilityratio (between lowest and highest values) is µ ( k low /k high ) = 0 .
04. We calculate theprefactor η N for 50 different realisations of this distribution and plot the results infigure 7c. For each Bond number we display the mean value, as well as one standarddeviation on either side µ ( η N ) ± σ ( η N ). The mean result is reminiscent of the previoussedimentary strata cases. However, it is interesting to note that the standard deviationis 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 than9 (a) (b) Sedimentary strata Turbidites (c) (d)
Spectrum Compacted
Figure 7: Nose growth prefactor η N given in terms of the single phase limit η N = 1 . k low /k high , the widthratio H low /H high , and the compaction power law β . In the case of the permeabilityspectrum, we show the mean result alongside one standard deviation above and below.Limiting behaviours are illustrated with dashed lines. Solid black curves correspond tothe homogeneous case, which is equivalent to the viscous limit, whereas all other curvescorrespond to the capillary limit.order unity. For CO sequestration applications, this indicates that particular attentiontowards measurements of heterogeneities should be paid for injection sites with low flowrates.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. Bycomparing figures 6b and 6a, we see that compaction significantly increases the saturationwithin the gravity current, which is due to the permeability gradient forcing the non-wetting phase upwards. This is accompanied by an increase in the prefactor η N and adecrease in h half , as seen in figures 7d, 8b. This is expected since, if a larger saturation ismaintained, volume conservation indicates that the gravity current must be thinner. Byvarying the compaction power law β from 0 to 1, we illustrate a fairly uniform transitionof the values of η N and h half between those of a uniform rock and those of a stronglycompacted rock.For Bo (cid:29) s ≈
1, as before, and this is0 (a) (b)
CompactedNon-compacted
Figure 8: Mid-range thickness of the gravity current h half = h ( η N / h half = 0 . H for the non-compacted cases (a) and the compactedcase (b). Limiting behaviours are illustrated with dashed lines.accompanied by both η N and h half converging to their single phase limits. Hence, atsuch large Bond numbers the effect of compaction on the saturation distribution andgravity current shape is not particularly noticeable. This is expected, since compactionforces the saturation upwards, just like gravity.For applications to CO storage, it is useful to infer from the above results how muchwe can expect heterogeneities to affect the speed of propagation of a gravity current. Suchinformation allows one to make efficiency predictions for CO storage that help pinpointthe best sites for injection, as well as safety predictions that ensure the CO does notspread beyond the desired perimeter. Using the homogeneous case, k low /k high = 1 (whichis 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 − . (3.1)Clearly ν depends on a number of parameters, but here we focus on the different types ofheterogeneity 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 valuesof the Bond number. The largest heterogeneity efficiency is observed for the sedimentarystrata 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-speedbands. The most we can expect heterogeneities to accelerate the gravity current at highBond number is around ν = 10 − ν = 100 − Q , the Bondnumber is held constant, but if the flow were to stop suddenly this would no longerbe the case. In such situations, the buoyancy length scale H = (cid:112) Q/u b , which waspreviously used to define the Bond number, would be rendered meaningless. Instead, theappropriate 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 number1 P o s t - i n j e c t i o n Figure 9: Heterogeneity efficiency ν (3.1), describing the relative increase in prefactorvalue η N due to heterogeneities, given as a ratio of the prefactor value for thehomogeneous case. Here we focus on the layered cases (S.S. stands for sedimentarystrata), ignoring compaction. We illustrate how in post-injection scenarios the Bondnumber typically decreases over time, such that heterogeneities play an increasinglyimportant role. In the case of the permeability spectrum we plot the mean value aswell as one standard deviation on either side. The permeability ratio for all cases is k low /k high = 0 . storage applications, where the injection of gasis 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 forunderstanding the post-injection spread of the CO . However, it is important to notethat such situations involve imbibition flows, as opposed to drainage flows, as we havestudied here (Woods 2015). Imbibition flows typically have different capillary pressureand relative permeability curves than drainage flows, though the approach studied hereis still applicable.In any case, we have shown here that heterogeneities have the potential to significantlyalter the growth of the gravity current, and consequently careful upscaled modelling isrequired. In the next section, we will show that some of these limiting behaviours can beexplained using asymptotic analysis.3.2. 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 solutionswhich we use to explain some of our earlier numerical results. We split the analysis intosituations without compaction β = 0 and with compaction β >
0. As in the previoussection, here we restrict our attention to viscous limit and capillary limit behaviour only.23.2.1.
Weak capillary forces without compaction
We already showed earlier that in the limit of large Bond number the saturationdistribution (2.27) is approximately uniform s ≈
1. Inserting this into the integrals(2.15),(2.18), we see that
S ≈ ϕh, (3.2)
K ≈ k k rn h, (3.3)which allows us to calculate the dimensionless fluxˆ F = ˆ f . (3.4)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 Bondnumber is identical to the classic single phase limit. This is of course expected, since inthe limit of weak capillary forces the flow of phases decouples, such that a single phasemodel becomes appropriate. This explains the convergence behaviour for both η N and h half for Bo (cid:29)
1, as we illustrate with dashed lines in figures 7 and 8.3.2.2.
Strong capillary forces without compaction
To address the limit of strong capillary forces, we first consider the homogeneouscase k low /k high = 1 (equivalent to the viscous limit) since this makes the analysis forthe subsequent heterogeneous cases easier. Hence, in the limit of small Bo (cid:28) s ≈ λ Bo ( h − z ) H . (3.5)Inserting this into the integrals (2.15),(2.18) we get
S ≈ ϕλ Bo2
H h , (3.6) K ≈ k k rn Hλ α Bo α α + 1 (cid:18) hH (cid:19) α +1 , (3.7)from which we calculate the dimensionless fluxˆ F = (cid:20) α/ ( λ Bo) α/ − α + 1 (cid:21) ˆ f α/ . (3.8)In this case, the flux has an α/ α = 4). The gravity current thickness is given by inverting (3.6),such that h = H (cid:16) f /λ Bo (cid:17) / . (3.9)Clearly, the thickness grows like h ∼ Bo − / as Bo →
0, which we illustrate in figure 8 withdashed 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 gravitycurrent 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 sedimentarystrata, the saturation distribution takes one of two possible values s ≈ (cid:40) k = k low λ Bo( h − z ) /H : k = k high . (3.10)3Since the low/high permeability layers are distributed according to the ratio H low /H high ,the integrals (2.15),(2.18) approximate to S ≈ (cid:18) H high H high + H low (cid:19) (cid:18) k high k (cid:19) a ϕλ Bo2
H h , (3.11) K ≈ (cid:18) H high H high + H low (cid:19) k high k rn Hλ α Bo α α + 1 (cid:18) hH (cid:19) α +1 . (3.12)Hence, the dimensionless flux isˆ F = (cid:34)(cid:18) H high H high + H low (cid:19) α/ (cid:18) k high k (cid:19) − a (1+ α/ α/ ( λ Bo) α/ − α + 1 (cid:35) ˆ f α/ . (3.13)Like the homogeneous case, the heterogeneous flux has an α/ α = 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 smallenough permeability ratio that (1 − k low /k high ) / (1 + k low /k high ) (cid:29) Bo. However, for thisstudy we restrict our attention to significantly heterogeneous media (rather than weaklyheterogeneous).We also note that in both the homogeneous case and the heterogeneous case, thecoefficients in (3.8),(3.13) depend on the Bond number itself. Therefore, we do not expecta constant value asymptote for the prefactor η N in the limit Bo → α = 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 onlycontribution to the integrals will come from the regions with the largest permeabilityvalue k = k high . However, since there are potentially many more than just two perme-ability 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 theabove is much more difficult, since we cannot approximate the saturation distribution assimply as (3.10).3.2.3. Weak capillary forces with compaction
In the case where the rock is compacted, the saturation (2.27) must be approximatelyuniform s ≈ S ≈ ϕH (1 + h/H ) − aβ − − aβ , (3.14) K ≈ k k rn H (1 + h/H ) − β − − β . (3.15)Note that if either a or aβ equals unity, we get analytical expressions with logarithmsinstead of (3.14),(3.15). Assuming β (cid:54) = 1, aβ (cid:54) = 1, we then calculate the dimensionlessflux, which we write in terms of h for now:ˆ F = (1 + h/H ) β ( a − − (1 + h/H ) aβ − β . (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 (cid:28) a = 0 . et al. S ≈ ϕh , as we observe in our numerical calculationsin figure 11e. Furthermore, in situations where the compaction law is weak β (cid:28)
1, oralternatively where the gravity current thickness is small h/H (cid:28)
1, the flux (3.16) reducesto ˆ
F ≈ h/H , thereby recovering the single phase limit (3.4). This is in accordance withour 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 mediumis sufficiently compacted the velocity distribution within the gravity current becomesdominated by the permeability variation, which upon integration creates a flux powerlaw that is not equal to one.3.2.4.
Strong capillary forces with compaction
To address the compacted case at small Bond number, we first consider the form ofthe saturation (2.10). In particular, we note that since k is monotone decreasing withdepth, 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 = p e ( z ∗ ), and willremain zero for all larger depths than this z > z ∗ . Therefore, without loss of generality,we take p = p e ( h ), such that there are no regions of zero saturation within the gravitycurrent. In this way, for small Bond number (2.10) approximates to a finite expressionwhich is independent of Bo, s ≈ − (cid:18) z h (cid:19) λβb . (3.17)This explains why the properties of the gravity current (e.g. η N , h half ) asymptote toconstant values for Bo → 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 notobey a fixed power law, but the exponent varies roughly between 1 and 2.3.3. Viscous-capillary transition
Up until now we have only discussed situations where the capillary number (2.20) iseither very large (viscous limit), in which capillary forces due to heterogeneities can beignored (i.e. effectively the homogeneous case), or very small (capillary limit), in whichthe heterogeneities play a dominant role on the flow behaviour (i.e. the heterogeneouscases studied above). However, in general the capillary number may vary between smalland large values throughout the aquifer. In this case, neither the viscous nor the capillarylimit can be applied to the flux function in (2.16), and instead the flux must depend onthe local capillary number, which is effectively a measure of the local pressure gradientswithin the gravity current. In this section, we discuss how to model this using numericalsimulations of the full PDE (2.16) coupled to the transcendental equation (2.20), therebydetermining which regions of the gravity current are within the viscous and capillarylimits, and which regions lie in between these limits.Following the same approach as Benham et al. (2021), we formulate composite func-tions for the upscaled properties of the flow, which capture both the viscous and capillarylimit regimes, as well as the transition between these limits. The two upscaled quantitiesof interest are the integrated saturation S and the flux F . For each of these upscaledquantities {S , F} , the transition behaviour is given in terms of the mean saturation (2.13)5 (a) (b) (c)(d) (e) (f) tu b /H = 0 . tu b /H = 0 . tu b /H = 0 . V i s c o u s C a p ill a r y V i s c o u s C a p ill a r y V i s c o u s C a p ill a r y N ct × ∆ N ct N ct × ∆ − (g) Figure 10: Numerical solution of the evolution of the gravity current (2.16), accountingfor transition behaviour between viscous and capillary limits using composite expressions(3.18) for the upscaled flow properties, where the capillary number is given implicitlyby (3.19) (with Bo= 1). The gravity current shape, shaded to illustrate the saturationdistribution (using the same colour scale as in figures 5,6), is illustrated in (a,b,c), whereasthe local capillary number N c is illustrated in (d,e,f). For all plots we shade regions withcapillary number larger than the transition value N c > N c t = 394 in green, and regionswith N c < N c t in blue. We also illustrate one folding scale ∆ on either side of N c t withdashed lines in (d,e,f). The evolution of the gravity current nose position r N ( t ) is shownon a log-log plot in (g).and capillary number (2.20) by the formula {S , F} trans = 12 (cid:20) {S (¯ s ) , F (¯ s ) } − tanh (cid:18) log N c − log N c t log ∆ (cid:19) + {S (¯ s ) , F (¯ s ) } + (cid:21) , (3.18)where {S , F} ± = {S , F} visc ± {S , F} cap is given in terms of the viscous limit (homoge-neous) and capillary limit (heterogeneous) upscaled properties derived earlier, and N c t , ∆ , are two fitting parameters which represent the transition capillary number and foldingscale, respectively. The precise values of these fitting parameters depend on numerousfactors, including the type of heterogeneity and the specific power laws used, but here weuse the same values as calculated by Benham et al. (2021), which are N c t = 394, ∆ = 5 . ∼ S 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 mustinstead solve the full PDE (2.16) numerically, in conjunction with the algebraic equationfor the capillary number (2.20). Written in dimensionless terms, we observe that thecapillary number is related to the Bond number according toN c = Bo p ∆p e (cid:12)(cid:12)(cid:12)(cid:12) h ( S , N c ) /H S h ( h ( S , N c ) , N c ) ∂ S ∂r (cid:12)(cid:12)(cid:12)(cid:12) . (3.19)For the sake of the calculations in this section, we use a mid-range value of Bo= 1 forthe Bond number, since this represents an even balance between capillary and gravityforces.In figure 10 we display the numerical solution in the case of sedimentary strata( k low /k high = 1 / , 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 localcapillary number (3.19). At early times the majority of the gravity current lies withinthe viscous limit, except the very tip of the gravity current nose. Indeed, the nose of thegravity current lies in the capillary limit for all times, since the flux vanishes there (sinceN c ∼ h∂p n /∂r ∼ h∂h/∂r → r → r N ( t ) due to (2.24)). At later times, a substantialportion of the gravity current lies within the capillary limit, whilst the viscous-capillarytransition (N c =N c t ) tends to a constant position of r/H ≈ .
2. The position of the noseof the gravity current r N ( t ) is plotted on a logarithmic scale in figure 10g, also indicatingthe viscous and capillary limit curves.We make several observations. Firstly, we note the dynamic transition between viscous-like behaviour at early times and capillary like-behaviour at late times. This is evidentfrom figure 10g, where we see that the nose position initially grows like ∼ t / withprefactor corresponding to the viscous limit, and later grows like ∼ t / with a capillarylimit prefactor. There is a transition period at intermediate times where the behaviourdoes 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 → ∞ ). Inthis way, the composite expressions (3.18) in the majority of the gravity current switchbetween viscous limit behaviour to capillary limit behaviour over time. Secondly, we notethat since the viscous-capillary transition (N c =N c t ) tends to a constant position, if wecontinue the simulation indefinitely, the fraction of the gravity current in the capillarylimit tends towards unity, indicating that at late times a capillary limit model is a goodapproximation for the whole aquifer.
4. Discussion of applications to CO sequestration In this section we discuss the implications of the upscaled description of gravitycurrents in the context of CO sequestration. There are many different sites that wecould choose as case studies, but probably the one with the most available data is atthe Sleipner project in the North Sea (Bickle et al. et al. migration at Sleipner, and comparisons have been made with seismic measurements of theextent of the plume (Bickle et al. et al. et al. et al. (2017) a = 0 .
14; the pore-entry pressure power law taken from the scaling proposedby 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 between1 and 10. Since it is unknown whether sedimentary strata, turbidites or a permeabilityspectrum is most appropriate, we investigate all of these different heterogeneity types.There are 2 parameters relating to the multiphase flow properties λ and α , which definethe capillary pressure and relative permeability. We note that not all relative permeabilitycurves are parameterised as simply as (2.5). For example, other more nonlinear functionalforms 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 k rn = k rn s (cid:2) − (1 − s ) (cid:3) . (4.1)We could change our formulation to account for this modified parameterisation, butinstead 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 poresize distribution is estimated by Chadwick et al. (2009) as λ = 2 / g =9 .
81 m/s ). For the top layer at Sleipner, Williams & Chadwick (2017) estimate thetemperature between 28-31 ◦ C and pressure between 8 . . ∆ρ = 232-309 kg/m . Meanwhile the viscosity of CO is taken as µ n = 54 . . × − Pa · s. The input flux is best estimated by Cowton et al. (2016), which forthe first few years of injection is Q = 1 . × − m /s. The mean permeability isestimated by Bickle et al. (2007) as k = 1 . − × − m . Finally, the pore entrypressure and end-point relative permeability are given by Williams & Chadwick (2017) as p = 1 . k rn = 0 .
28. Putting these all together, we estimate the Bond numberas Bo= 8 . − . η N for the gravitycurrent using the various different heterogeneity types. In this way, we can measure theheterogeneity efficiency ν (3.1), which tells us how much we can expect the heterogeneitiesto modify the speed of propagation. In the low/high Bond number estimate Bo=8 . / . ν = 36% /
31% for equally distributed sedimentary strata ( H low /H high = 1), ν =147% / H low /H high =10), ν = 11% /
6% for turbidites, and ν = 23 ± / ±
8% for a permeability spectrum(mean value). In the latter case, we include the standard deviation values to indicate howthese 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 moredetailed bore hole measurements at Sleipner, but also demonstrates the importance ofaccounting for such uncertainty in any modelling approach.8We note that the permeability ratio in the Sleipner field may not be as small as thevalue we have taken from Salt Creek, k low /k high = 0 .
01 (Bickle et al. ν calculated for k low /k high = 0 . / / k low /k high = 0 .
01. This indicates that, even if we have vastly overestimatedthe 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 havehad on the injection of CO at Sleipner largely depends on the type of heterogeneitiespresent. In particular, thin sedimentary strata with very high permeability could havecaused a potential speedup of more than 100%. However, for more moderate permeabilityprofiles, such as evenly distributed strata or turbidites, these heterogeneities may haveonly caused a 10 −
30% speedup.
5. Concluding remarks
We have studied the upscaled effect of several different types of heterogeneity on theevolution of an axisymmetric gravity current under an impermeable cap rock. The fourheterogeneity types considered, which were all given in terms of vertical variations inthe rock properties, were each motivated by different physical mechanisms for the non-uniform deposition or compaction of sediments. We developed a general method forcalculating the gravity current shape and growth rate in either the viscous or capillarylimits, which involves solving a single ordinary differential equation that depends onan upscaled flux term evaluated either via numerical integration, or using analyticalexpressions which we derived in certain limiting cases. This simplified approach not onlyreduces the computation time significantly, but also provides key insights into the roleof small-scale heterogeneities on the propagation of the large-scale flow.In particular, we showed that heterogeneities have the ability to speed up the growthrate of the gravity current by means of a nonlinear focussing mechanism (via the relativepermeability and capillary pressure) into high permeability layers. The degree of speedupdepends on the type of heterogeneity, and most importantly the fraction of high/lowpermeability regions within the medium. The largest effect was seen in the case ofsedimentary strata with thin regions of high permeability. Using a permeability profilecomposed of randomly sampled layers, we demonstrated how uncertainty in heterogeneitymeasurements can propagate to uncertainty in field predictions, an effect which is partic-ularly pronounced for small values of the Bond number. We also investigated modellingthe transition from the viscous limit regime to the capillary limit regime, shedding lightinto 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 small-scale heterogeneities affect large-scale CO migration, for safe and efficient sequestrationin porous aquifers. To assess the risks associated with any CO storage scheme, examiningthe effect of heterogeneities is essential. To illustrate this, we used the case study of CO injection at the Sleipner project in the North Sea. In this case, for realistic parametervalues, we demonstrated that heterogeneities may have sped up the gravity current bymore than 100% during injection. However, we also illustrated that this figure dependsgreatly on the heterogeneity type, indicating the need for detailed core measurementsfrom bore holes.For future work, variations in the heterogeneities along the length of the aquifer couldbe studied (in addition to the vertical heterogeneities investigated here), similarly toJackson & Krevor (2020). In such cases, the upscaled flow properties would dependon both the horizontal and vertical correlation length scales of the heterogeneities. In9addition, using the upscaled results presented here, predictions of trapping efficienciescould be calculated for various sequestration sites. This could shed light onto whichaquifers have heterogeneities that could potentially enhance their capability for CO storage. Furthermore, as we discussed earlier, it would be interesting to investigate theevolution of the gravity current after injection has stopped (see figure 9). However, inthis case both imbibition and drainage relative permeability/capillary pressure curvesmust be considered, depending on whether the CO is invading or being withdrawn fromdifferent regions of the aquifer, similarly to Golding et al. (2017).This research is funded in part by the GeoCquest consortium, a BHP-supportedcollaborative project between Cambridge, Stanford and Melbourne Universities, and bya NERC consortium grant “Migration of CO through North Sea Geological CarbonStorage 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 similaritysolution 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 atdifferent Bond numbers and for different types of heterogeneity.The relationship between the integrated saturation S and the gravity thickness h , asshown in figure 11a,c,e, is useful for understanding how to invert the solution of thegoverning 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 fluxfunction F (which is the integrated velocity), as displayed in figure 3a,b,c. We presentthe scaled velocity U = ( u n (0) − u n ( z )) / ( u n (0) − u n ( h )) so that U varies between U = 0at z = 0 and U = 1 at z = h . In this way we can see the functional form of U moreclearly. For example, if U is constant (as in figure 11f), when integrated this will give alinear 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 equationfor the integrated saturation S , which is given by (2.16). We start by considering thegoverning equation for the gravity current thickness (2.14). To re-write this in terms ofthe the integrated saturation (2.15), we first need to transform the derivatives of h intoderivatives of S . As illustrated by figure 11a,c,e, the integrated saturation is always amonotone increasing function of the gravity current thickness, so that we can define aunique inverse function h = S − ( S ) . (B 1)Then, using the chain rule, the gradient is given by ∂h∂r = ∂h∂ S ∂ S ∂r . (B 2)0 (a) (b)(c) (d)(e) (f) Bo= 10 − Bo= 1Bo= 10 Figure 11: (a,c,e) Variation of the integrated saturation S (2.15) for different values of theBond number Bo and different types of heterogeneity. Both S and h are normalised byreference values (at h = 2 h half ) for illustration purposes. In each plot we indicate powerlaw values of 1 /
2, 1, 2 and 3 with dotted lines for comparison. (b,d,f) Corresponding scaledvelocity 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.1We use the inverse derivative identity to calculate derivatives of (B 1), such that ∂h∂ S = (cid:18) ∂ S ∂h (cid:19) − . (B 3)Finally, by defining the flux function K as (2.18), we arrive at the governing equation for S , which is ∂ S ∂t = u b k k rn r ∂∂r (cid:18) r (cid:20) K ∂ S /∂h (cid:21) ∂ S ∂r (cid:19) . (B 4) REFERENCESBear, J
Dynamics of fluids in porous media . Courier Corporation.
Benham, GP, Bickle, MJ & Neufeld, JA
J. Fluid Mech. (in press) . Berg, S, Oedai, S & Ott, H –brine systems in sandstone. Int. J. Greenh. Gas Con. , 478–492. Bickle, MJ
Nat. Geosci. (12), 815–818. Bickle, M, Chadwick, A, Huppert, HE, Hallworth, M & Lyle, S
EarthPlanet. Sci. Lett. (1-2), 164–176.
Bickle, M, Kampman, N, Chapman, H, Ballentine, C, Dubacq, B, Galy, A,Sirikitputtisak, T, Warr, O, Wigley, M & Zhou, Z , brine and silicate minerals during geological carbon storage: Modelling based on afield CO injection experiment. Chem. Geol. , 17–31.
Boon, M & Benson, SM plume migration. Sub Judice . Braun, C, Helmig, R & Manthey, S
J. Contam. Hydrol. (1-2), 47–85. Brooks, R & Corey, T
Hydrology Papers, ColoradoState University , 37. Bui, M, Adjiman, CS, Bardow, A, Anthony, EJ, Boston, A, Brown, S, Fennell, PS,Fuss, S, Galindo, A & Hackett, LA
Energy Environ. Sci. (5), 1062–1176. Cavanagh, AJ & Haszeldine, RS plume requires fractured shale barriers within the utsira formation. Int. J.Greenh. Gas Con. , 101–112. Chadwick, RA, Noy, DJ & Holloway, S as a greenhouse gas mitigation measure. Pet. Geosci. (1), 59–73. Chierici, GL
Soc.Petrol. Eng. J. (03), 275–276. Cloud, WF
OilWeekly (8), 26.
Corey, AT
Prod. Monthly (1), 38–41. Cowton, LR, Neufeld, JA, White, NJ, Bickle, MJ, White, JC & Chadwick, RA -filledlayer at the Sleipner field, North Sea. J. Geophys. Res. Solid Earth (7), 5068–5085.
Cowton, LR, Neufeld, JA, White, NJ, Bickle, MJ, Williams, GA, White, JC &Chadwick, RA flow simulations at theSleipner field, North Sea. Earth Planet. Sci. Lett. , 121–133.
Golding, MJ, Huppert, HE & Neufeld, JA
Phys. Fluids (3), 036602. Golding, MJ, Huppert, HE & Neufeld, JA
J. Fluid Mech. , 550–577.
Golding, MJ, Neufeld, JA, Hesse, MA & Huppert, HE
J. Fluid Mech. , 248–270.
Harper, G, Liu, J, Tavener, S & Wildey, T
Comput. Methods inAppl. Mech. Eng. , 113469.
Hinton, EM & Woods, AW
J. Fluid Mech. , 411–429.
Hinton, EM & Woods, AW a Shear dispersion in a porous medium. Part 1. an intrusionwith a steady shape.
J. Fluid Mech. . Hinton, EM & Woods, AW b Shear dispersion in a porous medium. Part 2. an intrusionwith a growing shape.
J. Fluid Mech. . Huppert, HE
J. Fluid Mech. , 43–58.
Huppert, HE & Woods, AW
J. Fluid Mech. ,55–69.
Jackson, SJ, Agada, S, Reynolds, CA & Krevor, S
Water Resour. Res. (4), 3139–3161. Jackson, SJ & Krevor, S storage. Geophys. Res. Lett. p. e2020GL088616.
Krevor, S, Blunt, MJ, Benson, SM, Pentland, CH, Reynolds, C, Al-Menhali, A &Niu, B
Int. J. Greenh. Gas Con. , 221–237. Leverett, MC
Trans. AIME (01), 152–169.
Li, B & Benson, SM plumemigration in storage aquifers. Adv. Water Resour. , 389–404. Liu, J, Sadre-Marandi, F & Wang, Z
Procedia Comput. Sci. , 1301–1312. Lyle, S, Huppert, HE, Hallworth, M, Bickle, M & Chadwick, A
J. Fluid Mech. , 293–302.
Nelson, PH
The Log Analyst (03). Nordbotten, JM & Celia, MA
Geological storage of CO : modeling approaches forlarge-scale simulation . John Wiley & Sons. Pegler, SS, Huppert, HE & Neufeld, JA
J. Fluid Mech. , 592–620.
Rabinovich, A, Li, B & Durlofsky, LJ
Water Resour. Res. (10), 7645–7667. Trevisan, L, Krishnamurthy, PG & Meckel, Tip A saturation forbuoyant flow in clastic aquifers. Int. J. Greenh. Gas Con. , 237–249. Williams, GA & Chadwick, RA
Energy Procedia , 2856–2870.
Woods, AW