c 2008 Cambridge University Press J. Fluid Mech. (2008), vol. 611, pp. 35–60.  doi:10.1017/S002211200800219X Printed in the United Kingdom 35 Gravity currents with residual trapping M. A. H E S S E, F. M. O R R JR A N D H. A. T C H E L E P I† Department of Energy Resources Engineering, Stanford University, Stanford CA 94305, USA (Received 14 December 2007 and in revised form 28 April 2008) Motivated by geological carbon dioxide (CO2 ) storage, we present a verticalequilibrium sharp-interface model for the migration of immiscible gravity currents with constant residual trapping in a two-dimensional confined aquifer. The residual acts as a loss term that reduces the current volume continuously. In the limit of a horizontal aquifer, the interface shape is self-similar at early and at late times. The spreading of the current and the decay of its volume are governed by power-laws. At early times the exponent of the scaling law is independent of the residual, but at late times it decreases with increasing loss. Owing to the self-similar nature of the current the volume does not become zero, and the current continues to spread. In the hyperbolic limit, the leading edge of the current is given by a rarefaction and the trailing edge by a shock. In the presence of residual trapping, the current volume is reduced to zero in finite time. Expressions for the up-dip migration distance and the final migration time are obtained. Comparison with numerical results shows that the hyperbolic limit is a good approximation for currents with large mobility ratios even far from the hyperbolic limit. In gently sloping aquifers, the current evolution is divided into an initial near-parabolic stage, with power-law decrease of volume, and a later near-hyperbolic stage, characterized by a rapid decay of the plume volume. Our results suggest that the efficient residual trapping in dipping aquifers may allow CO2 storage in aquifers lacking structural closure, if CO2 is injected far enough from the outcrop of the aquifer. 1. Introduction 1.1. Carbon dioxide storage in saline aquifers This study is motivated by carbon capture and storage (CCS), also known as carbon dioxide (CO2 ) sequestration, in geological formations (for a summary see Metz et al. 2006). CCS involves capturing the CO2 arising from industrial processes and transporting it to a suitable storage site, where most of it is expected to be sequestered from the atmosphere permanently. Geological storage of CO2 has been proposed in saline aquifers, deep-sea sediments, depleted oil and gas reservoirs, and in combination with enhanced oil recovery or enhanced coal bed methane recovery. To make a significant contribution to the mitigation of climate change, gigatonnes of CO2 must be stored in the subsurface every year. CO2 storage on this scale is currently viable only in saline aquifers. The term saline aquifer is used for deep sedimentary rocks saturated with formation waters or brines, which contain high † Author to whom correspondence should be addressed. 36 M. A. Hesse, F. M. Orr Jr and H. A. Tchelepi (a) –15 (e) 10 –5 Trapped CO2: 4% 0 5 10 15 km –5 (b) 0 10 20 Trapped CO2: 6%  = 0.3 t = 112 years 30 40 50 km  = 0.3 t = 112 years (f ) 20%  = 10 t = 3730 years 19%  = 1.5 t = 560 years 44% =4 t = 1500 years 69% =8 t = 3000 years (g) (c) 32%  = 100 t = 37300 years 41%  = 1000 t = 373000 years (d) (h) Figure 1. Migration of the mobile CO2 current (dark grey) due to buoyancy and the formation of a zone containing residual saturation in its wake. The dimensional times correspond to the physical parameters discussed in § 6.2. (a–d ) Horizontal aquifer; (e–h) aquifer with 5◦ slope. concentrations of dissolved salts. Owing to their high salinity, these brines are not used as drinking water or for irrigation. Saline aquifers occur in all sedimentary basins and are not restricted to coal, oil, or gas provinces. CO2 is injected in the supercritical state to decrease the storage volume (Holloway & Savage 1993). For a typical geothermal gradient of 25 ◦ C km−1 , the depth required for CO2 to be in the supercritical state is approximately 800 m. Supercritical CO2 is less dense than the brine under all continental and shallow marine storage conditions. An impermeable seal, i.e. a shale or clay layer, overlying the storage formation is therefore always necessary to prevent direct upward migration of the buoyant supercritical CO2 . In addition, a concave downward structure, a structural trap or closure, is thought to be necessary to prevent upward migration of the CO2 along an inclined seal towards the surface. Much of the saline aquifer storage volume is in gently dipping regional aquifers that do not have a structural closure. Along-dip migration and eventual leakage is therefore a major concern for CO2 storage in such aquifers. However, trapping mechanisms continually reduce the volume of the CO2 plume. Trapping therefore limits both the time scale over which leakage is possible and the maximum up-dip migration distance (Ennis-King & Paterson 2002; Hesse, Tchelepi & Orr 2006). The two primary trapping mechanisms are the dissolution of CO2 into the brine (Lindeberg & Wessel-Berg 1997; Ennis-King, Preston & Paterson 2005) and the formation of residual saturation in the wake of the CO2 plume (Kumar et al. 2005; Mo et al. 2005); the latter is illustrated in figure 1. Dissolved CO2 is trapped, because it is negatively buoyant and will migrate downward. The large regional aquifers considered here are generally deep and have relatively low permeability. The onset time for convective dissolution of CO2 in low-permeability aquifers is large, and the convective fluxes after onset will be small (Riaz et al. 2006; Hassanzadeh, PooladiDarvish & Keith 2007). We therefore expect dissolution to be slow, and we assume that the formation of residual CO2 is the main trapping mechanism in deep aquifers. Gravity currents with residual trapping (a) 1 mm (b) Sc 1–Sbr Scr 0 37 (c) ∆x Q c(x, t) h x, t) Q c (x+∆ H h x,t) Q b (x+∆ hres Q b(x, t) z x θ Figure 2. (a) Microtomogram of non-wetting CO2 bubbles (black) and water (white) in the pore space between grains (grey), from Benson et al. (2006). (b) Comparison between non-wetting phase saturation profile after imbibition (– – –), inferred from a wire-line log (Koperna & Kuuskraa 2006), and the simplified step function saturation profile corresponding to the sharp-interface approximation (——). (c) The geometry of the porous layer and the variables used in the derivation are shown. 1.2. Immiscible displacements and residual saturations In many sedimentary rocks, supercritical CO2 is typically the non-wetting phase relative to the ambient brine. At the front of the CO2 plume, the CO2 saturation increases, and the brine is drained from the pore space. The capillary entry pressure prevents the drainage of the brine from the smallest pores, resulting in an incomplete displacement. We refer to the brine left behind the advancing CO2 front as the residual brine Sbr . Bachu & Bennion (2007) showed that Sbr can range from 0.2 to 0.68 at storage conditions in saline aquifers. The high end of these values is surprising and may in part be due to heterogeneity and gravity segregation in the experiments. They also show that the presence of residual water reduces the apparent permeability of the CO2 to approximately 1/5 of the single-phase permeability. We refer to this value as the relative permeability of CO2 , denoted krc . If the CO2 plume is migrating laterally as a gravity current, the CO2 saturation decreases at the trailing edge of the plume (figure 1), and the ambient brine imbibes into the pore space previously occupied by CO2 . Preferential imbibition of the brine into the smaller pores and interfacial instabilities leave CO2 behind as disconnected bubbles and ganglia of CO2 (figure 2a), which are effectively immobile. We refer to this immobile CO2 saturation as the residual CO2 saturation, Scr , and to the process as residual trapping. Bachu & Bennion (2007) report values of Scr from 0.1 to 0.35 for saline aquifers in the Alberta basin, indicating that they will trap CO2 efficiently. Most work on residual trapping during CO2 storage has focused on the effect of hysteresis on the magnitude of Scr , and the design of injection strategies that maximize residual trapping during, or shortly after, the injection period (Mo et al. 2005; Juanes et al. 2006; Ide, Jessen & Orr 2007). We focus on the effect of the vertical sweep efficiency, S, which may be defined as the fraction of the aquifer that is contacted by the CO2 plume during its buoyancy-driven migration in the post-injection period. A volume balance argument shows that the up-dip migration distance, x↑ , up to a constant is given by Vc x↑ ≈ , (1.1) φH SScr where Vc is the volume of CO2 injected, H the thickness of the aquifer, and φ 38 M. A. Hesse, F. M. Orr Jr and H. A. Tchelepi the porosity. Ennis-King & Paterson (2002) used this relationship to highlight the sensitivity of the migration distance to the magnitude of Scr . Equation (1.1) shows that the effect of S on the migration distance of the CO2 plume is comparable to that of Scr . The sweep is expected to be less than unity, because gravity segregation and viscous instabilities lead to the formation of a gravity tongue along the top of the aquifer. 1.3. Previous work on gravity currents and CO2 storage While the residual saturation is an upscaled parameter determined by viscous and capillary forces on the pore scale, the vertical sweep is a macroscopic quantity that can be studied using fluid mechanical models of gravity currents in porous media. Nordbotten, Celia & Bachu (2005), Lyle et al. (2005), Vella & Huppert (2006), Verdon & Woods (2007) and Bickle et al. (2007) used sharp-interface models to study the development of the CO2 plume during the injection period in confined and unconfined aquifers. During injection, the CO2 saturations increase and no residual CO2 forms. After several decades, the injection of CO2 ends, and the CO2 plume migrates as a gravity current for several hundreds or thousands of years. Thus, Hesse et al. (2006, 2007b) argued that the long-term evolution of the CO2 plume can be modelled as an instantaneous release of finite volume. In the post-injection period, the CO2 saturation decreases in the wake of the gravity current and residual trapping reduces the current volume (figure 1). Kochina, Mikhailov & Filinov (1983) extended the sharp-interface model to allow for a constant residual saturation, and studied the decay of the plume volume after a finite release in a horizontal aquifer. Hesse et al. (2006) used a similar model to study the effect of slope on the migration distance in unconfined aquifers. We extend their work to confined sloping aquifers and consider the limiting cases. 2. Governing equations 2.1. Derivation We consider the flow of supercritical CO2 with density ρc = ρ and viscosity μc and of brine with density ρb = ρ + ρ and viscosity μb in a sloping porous layer of constant thickness H , with dip angle θ > 0, and infinite lateral extent (figure 2c). Both fluids are considered incompressible, and we assume that the porous medium is homogeneous and isotropic with permeability k and porosity φ, and that the top and bottom boundaries are impermeable. We consider an aquifer and a gravity current with a large aspect ratio (length : height ≫ 1) so that the velocities in the z-direction are negligible over most of the current, and the pressure in each phase is hydrostatic. We also assume that the displacement is slow enough that gravity–capillary equilibrium is maintained in any vertical cross-section of the current (figure 2b). This assumption of vertical-equilibrium has been used extensively in petroleum engineering and in hydrology (for examples see Yortsos 1995; Bear & Ryzhik 1998). If the thickness of the released fluid is large compared to the capillary transition zone, the saturation in each fluid is approximately constant (figure 2b). As a first approximation we replace the transition zone by a sharp interface and assume that the saturations above and below the interface are constant. We refer to the resulting models as a vertical-equilibrium sharp-interface model, to emphasize the two main assumptions underlying the model. We denote the thickness of the CO2 plume by hc = h so that the depth of the brine is given by hb = H − h. In this case the pressure Gravity currents with residual trapping 39 in the aquifer is given by p=  pI − gρ(z − hb ) + Pc for z > hb , pI − g(ρ + ρ)(z − hb ) for z 6 hb , (2.1) where pI is the unknown pressure at the interface, Pc is the constant capillary pressure, and g is the acceleration due to gravity. The volume flux per unit width qp of phase p ∈ {c, b} is given by the multiphase extension of Darcy’s law qp = −kλp ∂φp /∂x, where φp = p −gρp (x sin θ +z cos θ) is the potential of phase p, and λp = krp /μp is the mobility of phase p, assumed to be constant. The flow rate per unit width of phase p is given by Qp = hp qp . In the absence of a source term, and with the assumption of incompressibility, the global conservation of volume is given by Qc + Qb = 0. Using this constraint we can eliminate ∂pI /∂x from the expressions for the flow rates, and we obtain   hλc (H − h)λb ∂h sin θ − cos θ Qc = −Qb = kgρ . (2.2) hλc + (H − h)λb ∂x To obtain an equation for the evolution of the interface, we consider the conservation of the CO2 volume Vc over region x and time t as shown in figure 2(c). Vc = φ(1 − Sbr )hx = (Qc |x − Qc |x+x )t + Rc . (2.3) The source term Rc accounts for the volume of CO2 that is lost as residual saturation Scr in the wake of the plume. Following Kochina et al. (1983) we assume a constant residual saturation, Scr , so that  −hxScr for ∂h/∂t < 0, Rc = (2.4) 0 for ∂h/∂t > 0. Taking limits for small x and t the equation for the evolution of the interface is given by    ∂ h(H − h) ∂h ∂h =κ cos θ , (2.5) − sin θ + ∂t ∂x h(M − 1) + H ∂x where we have introduced two parameters: the conductivity κ of the CO2 is given by ⎧ kλc ρg ⎪ ⎪ ⎨κ1 = φ(1 − S − S ) for ∂h/∂t < 0, br cr (2.6) κ= ⎪ kλc ρg ⎪ ⎩κ0 = for ∂h/∂t > 0, φ(1 − Sbr ) where κ1 > κ0 , and the mobility ratio is given by M = λc /λb . The mobility ratio M measures the change in the mobilities across the advancing part of the interface, because we neglect the effect of the residual CO2 on the mobility of the brine. For geological CO2 storage the ambient brine is less mobile and M > 1, and we restrict the discussion to this range. Hesse et al. (2007b) showed that this choice of parameters allows a simple reduction of (2.5) in the limits max(h) ≪ H or M ≪ 1. For Scr = 0, the coefficient (2.6) becomes continuous and (2.5) reduces to the form given by Bear (1972). The thickness of the zone saturated with residual saturation is given by hr [x, t1 ] = hmax [x, t1 ] − h[x, t1 ], where the vertical sweep hmax is the largest thickness the CO2 40 M. A. Hesse, F. M. Orr Jr and H. A. Tchelepi plume has achieved at a particular location up to the time of consideration hmax [x, t1 ] = max(h[x, t]), t6t1 ∀x. We also refer to hr as the residual surface that identifies the fraction of the aquifer swept by the CO2 plume. We study the evolution of a finite release of CO2 into an aquifer saturated with brine. Owing to the diffusive nature of the governing equations, the details of the particular shape of the initial condition have a small effect on the intermediate and long-term evolution. We consider initial conditions with compact support  h0 [x], |x| 6 L/2, (2.7) h [x, t = 0] = 0, |x| > L/2. In all cases discussed in this paper we simply choose h0 = H , so that the initial aspect ratio of the current is given by A0 = L/H . The initial saturations in the CO2 and in the brine are Sc = 1 − Sbr and Sb = 1. This corresponds to a porous medium that is invaded by the released fluid for the first time, which is appropriate for CO2 storage in saline aquifers. 2.2. Dimensionless form We choose the following dimensionless variables η = h/H, ξ = x/L, σ = κ/κ1 6 1, t/t ∗ , where L is the width of the initial condition. We may choose either a diffusive or an advective characteristic time scale, t ∗ , given by t a = L(κ1 sin θ)−1 a or t d = L2 (κ1 H cos θ)−1 , (2.8a, b) d so that τ = t/t and ϑ = t/t . Substituting these definitions into (2.5) we obtain two dimensionless equations σ −1 ητ + fξ = Pe −1 (f ηξ )ξ or σ −1 ηϑ + Pefξ = (f ηξ )ξ , (2.9a, b) depending on the choice of the characteristic time scale. The flux function is defined as η(1 − η) , (2.10) f = η(M − 1) + 1 and the dimensionless discontinuous coefficient is given by σ = 1, ητ < 0 or ηϑ < 0, 1 − ǫ, ητ > 0 or ηϑ > 0, (2.11) depending on the choice of the characteristic time. Equations (2.9a) and (2.9b) have the same dimensionless governing parameters M= krc μb , μc krb Pe = A0 tan θ, ǫ= Scr , 1 − Sbr the mobility ratio, M > 0, the Péclet number, Pe > 0, and the residual, 0 6 ǫ < 1. The advective and diffusive time scales are related by τ = Peϑ. The dimensionless initial condition (2.7) is given by  1, |ξ | < 12 , (2.12) η(ξ, 0) = 0, |ξ | > 12 . Gravity currents with residual trapping 41 3. Parabolic limit: horizontal aquifers Hesse et al. (2007b) discussed the evolution of a finite release in a confined horizontal aquifer, Pe = 0, without residual saturations, ǫ = 0. They showed that the evolution of a rectangular initial condition, given by (2.12), can be understood in terms of an early self-similar solution, corresponding to a tilting interface in a confined aquifer with the similarity variable ξ/ϑ 1/2 , followed by a transition to a late self-similar solution, corresponding to a finite release in an unconfined aquifer with the similarity variable ξ/ϑ 1/3 . For M > 1, the late similarity solution, which neglects the ambient fluid, is a good approximation only when the current is very thin relative to the thickness of the aquifer. In the presence of residual trapping an analogous transition from an early to a late scaling law is observed. However, we show below that the exponents of the early scaling laws are independent of the residual, whereas Bear & Ryzhik (1998) showed that the exponents of the late scaling laws are a function of the residual. In the terminology of Barenblatt (1996), the evolution of the current transitions from a similarity solution of the first kind to a similarity solution of the second kind. 3.1. Early similarity solutions with loss Initially, the two interfaces of the current evolve independently and symmetrically (figure 1a). Until they interact, they evolve as tilting interfaces described by a selfsimilar solution given by a step function initial condition  a, ξ < ξ̃ , η= (3.1) b, ξ > ξ̃ , where ξ̃ is the initial position of the interface. Here we consider only the cases a = 1 and b = 0, but self-similar solutions are not limited to this choice. The outward and inward propagating tips of the interface are denoted by ξo and ξi respectively. For ξ̃ < ξ 6 ξo the interface is advancing, and the thickness of the current is increasing, ηϑ > 0, so that σ = 1 − ǫ. For ξi 6 ξ < ξ̃ the interface is receding, and the thickness of the current is decreasing, ηϑ < 0, so that σ = 1. The location ξ̃ of the discontinuity in the coefficient is constant at early times. Equation (2.9b) and the initial condition (3.1) are self-similar in the variable ζ = ξ/ϑ 1/2 . Therefore, the propagation of the plume tip at early time is given by the scaling law ξo = a[ǫ, M]ϑ 1/2 , where a is a constant that decreases with increasing ǫ. Inserting Pe = 0 and the similarity variable into (2.5b), we obtain the ordinary differential equation   d η(1 − η) dη ζ dη = σ̃ , (3.2) − 2 dζ dζ η(M − 1) + 1 dζ where σ [ηζ ζϑ ] = σ̃ [ηζ ]. The mobility ratio, M, and the residual, ǫ, determine the shape of the similarity solution at early times. The inner and outer boundaries of integration ζi and ζo are unknown, and generally are functions of M and ǫ. These eigenvalues must be determined as part of the solution, requiring two additional constraints. The boundary conditions on the inward and outward propagating tips are given by η [ζi ] = 1, dη dζ = ζi ζi M , 2 (3.3a, b) 42 M. A. Hesse, F. M. Orr Jr and H. A. Tchelepi (a) 1.0 (b) 1.0 ζi ζo 0 0 –0.5 1 ǫ= 0 –1.0 = = 1 ǫ= 0 –1.5 ǫ 0.5 ǫ η 0.5 ζi ζo 0.5 0 –0.50 –0.25 1.0 0 ζ 0.25 0.50 0.75 ζ (c) 20 (d ) 20 0.2 0.1 –0.3 0.3 15 15 0.4 M 10 10 0.5 –0.4 0.6 5 1 0 0.2 –0.5 –0.6–0.7 5 0.7 0.8 0.4 0.6 0.8 1 1.0 0 0.2 ǫ 0.4 0.6 –0.8 0.8 1.0 ǫ Figure 3. (a) Interface shapes for M = 1 and ǫ increasing from 0 to 1 in 0.2 increments. (b) Interface shapes for M = 10 and ǫ increasing from 0 to 1 in 0.2 increments. (c) Contours of the outward propagating tip position, ζo . (d ) Contours of the inward propagating tip position, ζi . η [ζo ] = 0, dη dζ = ζo −ζo , 2(1 − ǫ) (3.3c, d ) respectively. The conditions on η specify the vertical extent of the interface. The conditions on ηζ have no direct physical interpretation, but they are required to satisfy (3.2) at the boundaries ζi and ζ0 , and they provide the two additional constraints necessary to determine the eigenvalues. Together (3.2) and (3.3 a–d ) form a nonlinear eigenvalue problem for the shape of the interface and the two unknown tip positions, ζi and ζo . This eigenvalue problem has been solved numerically as a function of the governing parameters M and ǫ using a shooting method. Figures 3(a–b) show that increasing M leads to the formation of a gravity tongue, while increasing ǫ tends to shorten the gravity tongue and may lead to an inflection in the interface. For M > 1, the outward propagating tip position, ζo , is mostly a function of ǫ, while the inward propagating tip position, ζi , is more strongly affected by M (figure 3c, d ). In the limit where all fluid is trapped, ǫ = 1, the interface does not propagate, ξo = ξ̃ , but the interface does recede, ξi = ξ̃ . However, the distinction between the fluids across the receding interface is lost, and the interface loses it physical significance. The evolution of this imaginary interface is still given by (3.2), but it does not correspond to any fluid movement. The early similarity solution shows that residual trapping has only a weak influence on the evolution of the receding part of the interface. Instead residual trapping slows down the advancing interface, because less fluid is supplied from the receding part. A similar behaviour is observed in the hyperbolic model discussed in § 4. Gravity currents with residual trapping (a) 43 (b) 102 1/2 early late 101 1/3 β ξˆo 1/4 100 ǫ=0 ǫ = 0.25 ǫ = 0.50 ǫ = 0.75 1/8 0 0.25 0.50 ǫ 0.75 1.0 10–1 –2 10 100 102 104  Figure 4. Change of the solution from the early to the late scaling law (M = 1). (a) The scaling laws for ξo ∝ ϑ β at early and late times as a function of the residual, ǫ. The symbols at ǫ of 0, 1/4, 1/2 and 3/4 show the exponents obtained by fitting the numerical data shown in (b). For the early scaling law the data was fit for ϑ between 0.1 and 10 and for the late scaling law the data was fit for ϑ > 1000. (b) The shifted tip position, ξ̂o = ξo − ξ̃ , is shown as a function of ϑ, for increasing residuals, where ξ̃ = 1. 3.2. Regime transition and late solution Hesse et al. (2007b) showed that f → η at late times when the current is much thinner than the thickness of the aquifer, so that (2.9b) reduces to ηϑ = σ (ηηξ )ξ , (3.4) for all M. The only parameter in this equation is the residual, ǫ, in the discontinuous coefficient, σ . Semi-analytic solutions for the interface shape have been obtained in the radial (Kochina et al. 1983; Barenblatt 1996), and the one-dimensional case (Bear & Ryzhik 1998). For a compact initial condition, the solution of this parabolic equation with a moving discontinuous coefficient is a similarity solution of the second kind (Barenblatt 1996). These solutions are characterized by an anomalous exponent in the scaling laws that govern the solution. In this case we obtain ξo = c1 [ǫ, M, η0 ] ϑ β(ǫ) , ηmax = c2 [ǫ, M, η0 ] ϑ 2β(ǫ)−1 , V = c3 [ǫ, M, η0 ] ϑ 3β(ǫ)−1 , (3.5a) (3.5b) (3.5c) for the tip of the plume, the height of the plume, and the volume of the plume. Figure 4(a) shows the decrease of the anomalous exponent β with increasing ǫ. For (3.4), the exponent must be determined numerically from a nonlinear eigenvalue problem (see Bear & Ryzhik 1998). As expected, increasing residual trapping leads to a faster decrease of the plume volume and slower tip propagation. Another characteristic of this type of solution is that the constants of proportionality (c1 –c3 ) in these scaling laws depend on the initial condition. Figure 4(b) shows the transition from an exponent of 1/2, to the late anomalous exponent β 6 1/3. Initially, only the constant in the scaling law is affected by ǫ, but at late times, both the coefficient and the exponent depend on ǫ. The sharp-interface model presented here predicts that the volume of the current always decays as a power law in horizontal aquifers. In this case, the current is never exhausted and continues to spread. Capillary forces, which are neglected in this model, will arrest the movement of the current once it becomes thin. 44 M. A. Hesse, F. M. Orr Jr and H. A. Tchelepi 4. Hyperbolic limit: sloping aquifers For confined sloping aquifers the current evolution is not self-similar, and no general solutions to (2.5) are known. For large Pe numbers, (2.9a) reduces to a simpler equation that allows an analysis of the current evolution. For Pe → ∞, we obtain the following quasi-linear hyperbolic equation ητ + σ ληξ = 0, (4.1) where the derivative of the flux function is denoted by λ ≡ fη . 4.1. Riemann problems First we consider two Riemann problems to understand the effect of trapping, i.e. the discontinuous coefficient σ , on rarefactions and shocks. Then we consider rectangular initial condition (2.12), representing a finite release of fluid, and the resulting wave interaction. For non-zero values of loss, the wave interaction will lead to an extinction of the current in finite time. 4.1.1. Rarefaction Consider the piecewise constant initial data  1, ξ 6 ξ̃ , η[ξ, τ = 0] = 0, ξ > ξ̃ . (4.2) The weak solution, a spreading wave or rarefaction, is obtained by the method of characteristics. This rarefaction is a tilting interface analogous to the self-similar solutions presented in § 3.1. Following the notation introduced there, we refer to the tip propagating to the left as ξi and to the tip propagating to the right as ξo . The characteristics of the solution originating at the discontinuity are given by ξ = σ [ητ ] λ[η]τ + ξ̃ , (4.3) and therefore all solutions are self-similar in the coordinate ξ/τ . The discontinuity in σ leads to a discontinuity of the slope of the interface at ξ̃ , but the interface itself remains continuous, because λ = 0 at ξ̃ . For M = 1 and ǫ < 1 the equation for the evolution of η[ξ, τ ], given by inverting (4.3) is ⎧ 1, ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ −τ + M̂(ξ0 − ξ ) + Mτ (τ + M̂(ξ − ξ0 )) ⎪ ⎪ ⎪ , ⎪ ⎪ ⎪ )) M̂(τ + M̂(ξ − ξ 0 ⎪ ⎨ 1 η[ξ, τ ] = √ , ⎪ ⎪ (1 + M) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ −τ ǫ̂ + M̂(ξ0 − ξ ) + Mτ ǫ̂(τ ǫ̂ + M̂(ξ − ξ0 )) ⎪ ⎪ ⎪ , ⎪ ⎪ ⎪ M̂(τ ǫ̂ + M̂(ξ − ξ0 )) ⎩ 0, ξ 6 ξi , ξi < ξ < ξ̃ , ξ = ξ̃ , ξ̃ < ξ 6 ξo , ξo < ξ, (4.4a) Gravity currents with residual trapping (a) (b) 0.06 ǫ=0 0.04 σf 45 1.00 0.25 0.75 0.50 η 0.50 0.02 0.75 ǫ=1 0 0.25 1.00 0 0.25 0.50 η 0.75 1.0 0 –0.5 ǫ=0 0 0.5 1.0 ξ/τ Figure 5. (a) The flux function f is shown for M = 10 and increasing trapping ǫ. (b) The corresponding solution profiles given by (4.4) in self-similar coordinates. where M̂ = M − 1, ǫ̂ = 1 − ǫ. For M = 1 and ǫ < 1, the solution becomes piecewise linear and is given by ⎧ 1, ξ 6 ξi , ⎪ ⎪ ⎪   ⎪ ⎪ ξ0 − ξ ⎪ 1 ⎪ +1 , ξi < ξ < ξ̃ , ⎪ ⎪ 2 ⎪ 2τ ⎨ (4.4b) η [ξ, τ ] = 12 , ξ = ξ̃ , ⎪ ⎪   ⎪ ⎪1 ξ0 − ξ ⎪ ⎪ ⎪ ⎪ 2 2τ (1 − ǫ) + 1 , ξ̃ < ξ 6 ξo , ⎪ ⎪ ⎩ 0, ξo < ξ, where the positions of the outward and inward propagating tips of the rarefaction are given by ξo = (1 − ǫ)τ + ξ̃ , ξi = − τ + ξ̃ . M (4.5) For tilting interfaces in the parabolic limit, we have shown that ξi is a stronger function of M than ǫ. In the hyperbolic limit the entire receding interface has lost its dependence on ǫ. Where the interface is advancing, f is scaled by σ = 1−ǫ (figure 5a), and the residual reduces the propagation speed of the advancing interface (figure 5b). This behaviour is analogous to the self-similar solutions in the parabolic limit. The solution for ǫ = 1 is still given by (4.4), but the advancing interface is replaced by a discontinuity at ξ̃ with a jump from 0 to ηmax . A comparison of figures 3(b) and 5(b) shows that M and ǫ have a similar effect on the shape of the interface in both the parabolic and hyperbolic limits. For M = 0, (4.1) reduces to the linear advection equation ητ + σ ηξ = 0, (4.6) with a moving discontinuity in the velocity, σ . In this limit an ambiguity arises in σ (ητ ), because ητ is singular at the discontinuity. The solution to (4.6) can be obtained taking the limit of (4.4) as M → 0 and requiring that the solution changes continuously with changes in the parameter M. In this limit, the rarefaction degenerates to a contact 46 M. A. Hesse, F. M. Orr Jr and H. A. Tchelepi discontinuity η [ξ, τ ; M = 0] =  1, ξ 6 Λc τ + ξ̃ , 0, Λc τ + ξ̃ < ξ, (4.7) where Λc = 1 − ǫ is the speed of the discontinuity. 4.1.2. Shock Consider the piecewise constant initial data  ηL , ξ 6 ξ̃ , η[ξ, τ = 0] = ηR , ξ > ξ̃ , (4.8) where ηL < ηR . The weak solution of (4.1) with (4.8) is a shock of strength, ηs = ηR −ηL , given by  ηR , ξ 6 Λτ + ξ̃ , η [ξ, τ ; M] = (4.9) ηL , Λτ + ξ̃ < ξ, where Λ is the shock speed. Similarly to (4.6) an ambiguity arises, because ητ is singular at the shock. Therefore we define the shock speed considering (4.1) regularized by a diffusion term ητ + σ [ητ ] fη ηξ = νηξ ξ , (4.10) where ν > 0. For initial data (4.8) the entropy-satisfying solution to (4.1) is defined by the solution of (4.10) in the limit of small ν, the so-called vanishing viscosity solution. We assume the solution of (4.10) with (4.8) takes the form of a travelling wave with constant speed Λν , so that η [ξ, τ ] = η̂ [ζ ], where the travelling-wave coordinate is ζ = ξ − Λν τ . For a steady travelling wave, the time derivative that determines the value of σ is given by ητ = ζτ η̂ζ = −Λν η̂ζ . (4.11) For initial data (4.8) we expect the solution to be monotonically increasing, η̂ζ > 0, so that the sign of ητ is determined by the direction of wave propagation, given by the sign of Λν . The constant wave speed implies that ητ does not change sign, and hence σ is continuous, but its magnitude depends on the sign of the wave speed, σ = σ̂ [Λν ]. In the moving coordinate system, (4.10) reduces to the ordinary differential equation with the boundary conditions η̂ → ηR for ζ → +∞, and η̂ → ηL for ζ → −∞. After integrating once and eliminating the constant of integration with either boundary condition, the following expression is obtained η̂ζ = σ M(η̂ − ηL )(ηR − η̂) ν(M̂ η̂ + 1)(M̂ηL + 1)(M̂ηR + 1) . (4.12) This confirms our assumption that η̂ζ > 0 for ηL < η̂ < ηR . Equation (4.12) does not depend on the wave speed, because consistency with both boundary conditions requires that the wave speed be given by Λν = σ 1 − M̂ηR ηL − ηL − ηR (M̂ηL + 1)(M̂ηR + 1) , (4.13) which has already been substituted into (4.12). For M = 1 the solution to (4.12) is    1 (ζ − Cν)σ (ηR − ηL ) η̂ = ηL + ηR + (ηR − ηL ) tanh , (4.14) 2 2ν Gravity currents with residual trapping (a) 6 (b) 6 (×10–2) 4 (×10–2) 4 σf σf 2 (c)1.0 (×10–2) 0.5 Λ 2 0.5 η 1.0 rec. 0 adv. –0.5 adv. rec. 0 47 0 0.5 η 1.0 –1.0 0 0.5 ηL 1.0 Figure 6. Shock construction for M = 10. (a) In the absence of trapping, ǫ = 0, the construction for advancing (adv.) and receding (rec.) shocks is the same. (b) In the presence of trapping, ǫ = 1/2, the shock construction for advancing shocks is modified. (c) Shock speed as a function of ηL , for ηL > ηc the shock speed depends on the amount of trapping. where C is a constant of integration that can be determined by the initial position of the discontinuity, ξ̃ , and σ is given by (2.11). Equations (4.11) and (4.12) allow us to write an equivalent definition of σ in terms of the sign of the wave speed σ [ητ ] = σ̂ [Λν ] = 1, Λν > 0, 1 − ǫ, Λν < 0. (4.15) This definition is not ambiguous at discontinuities in the limit ν → 0. With this definition of σ , the Rankine–Hugoniot jump condition for the shock is Λ = σ̂ f [ηL ] − f [ηR ] = Λν , ηL − ηR (4.16) so that the shock speed of the limiting solution is identical to the wave speed of the viscous solution. The geometric interpretation of the jump condition (4.16), is a chord joining f [ηL ] and f [ηR ], as shown in figure 6(a, b). For a given ηR and M the direction of wave propagation as a function of ηL is given by ⎧ ⎪ ⎨< 0, ηc < ηL < ηR , (4.17) Λν = 0, ηL = ηc , ⎪ ⎩> 0, 0 < η < η , L c where 1 − ηR , (4.18) (M − 1)ηR + 1 which corresponds to a horizontal tangent and a stationary shock. The speed of shocks moving to the right is not affected by trapping, while the speed of waves travelling to the left is reduced owing to trapping (figure 6c). This behaviour of the shock is analogous to the behaviour of the rarefaction. In both cases, the speed of the interface is reduced, by a factor 1 − ǫ due to trapping, when it is advancing, but it is not affected if the interface is receding. However, the parameter that identifies whether the interface recedes or not is different for the rarefaction and the shock. For an interface with a finite slope, i.e. a rarefaction, the change of the height of the interface with time, ητ , determines which parts of the interface recede. For a vertical interface, i.e. a shock, the sign of shock speed, Λ, determines whether the interface recedes across the shock or not. ηc = 48 M. A. Hesse, F. M. Orr Jr and H. A. Tchelepi 50 16 (a) (b) D D 40 12 30 τ 8 τ 20 4 C C 10 B 0 A –1.5 0 F 2 4 ξ 6 8 10 0 –5 B A 0F 10 20 30 ξ (c) 1 16 η 8 τ 0 –1 0 2 4 ξ 6 8 10 0 Figure 7. Characteristic portraits of the evolution of initial data (2.12), in the (ξ, τ )-plane. The current is bounded by the shock, A to D, in the back and the leading characteristic of the rarefaction, F to D, in the front. (a) Analytic solution for M = 1 and ǫ = 0.4. (b) Semi-analytic solution for M = 5 and ǫ = 0.4. (c) Evolution of the current interface (black) and the residual surface (grey). 4.2. Wave interaction The evolution of a finite release of fluid is governed by the interaction of the shock at the back of the current with the rarefaction at the front of the current. In the presence of trapping, this interaction leads to a rapid reduction of current volume, and gives rise to a finite migration distance and time. Consider the piecewise constant initial data given by (2.12). Initially, the solution of the front is a rarefaction given by (4.4) with ξ̃ = 1/2, for all M > 0, and the evolution of the back is a stationary shock (A to B in figure 7) with constant shock strength (figure 8b), because Λ = 0 for ηL = 0 and ηR = 1 for all M > 0. When the back end of the rarefaction reaches the shock, the shock strength decreases and it begins to move. The solution, η, will then be given by (4.4), but η [ξ < ξs ] = 0, where the ξs is the shock position. 4.2.1. Shock evolution The waves begin to interact at the transition time, τt = M (B in figure 7), and the shock speed begins to increase as the shock strength decreases. Considering the shock strength, ηs , as a parameter along the shock path (B to D), so that τ = τ [ηs ], and ξ = ξ [ηs ], we differentiate (4.3) and (4.16) to obtain two differential equations ξηs = ληs [ηs ] τ + λ [ηs ] τηs , ηs ξηs = f [ηs ] τηs . (4.19) Gravity currents with residual trapping 16 (b) 49 (a) 12 τ 8 4 0. C 0 B A 1 0.5 ηs ηr 0 ǫ= ǫ=1 τt 0 ǫ= τp 4 ǫ= 0.7 ǫ= 0 .6 0.5 ǫ= F 1 (c) 0.5 0 –1 0 2 4 6 8 10 ξ Figure 8. (a) Portrait of the shock path and the leading characteristic of the rarefaction, for   increasing residual. The locus of the termination points, τf = τf ξ↑ , is shown as a dotted line. (b) The corresponding temporal of the shock strength, ηs . (c) The corresponding  evolution  final residual surfaces, ηr = ηr τf . Eliminating ξηs , we obtain the differential equation f [ηs ] − ηs λ [ηs ] dηs = , dτ ηs λη [ηs ] τ (4.20) for the evolution of the shock strength. Together with the initial condition ηs [τt ] = 1, it determines the evolution of the shock strength (figure 8b). Once ηs is known, the shock position is given by ξs = λ [ηs ] τ + 1/2. We solve for ηs instead of ξs , because (4.20) is always separable. If trapping occurs, ǫ = 0, the evolution of ηs is divided into an early period, τt 6 τ 6 τp , during which the shock interacts with the receding portion of the rarefaction (B to C in figure 7), and a later period τp < τ 6 τf , where the shock interacts with the advancing portion of the rarefaction (C to D). Hence, the problem has three time scales: the transition time, τt , when interaction between waves begins (B ), the passing time, τp , when the shock passes the stationary point of the rarefaction and starts to interact with the advancing section of the rarefaction (C ), and the final migration time, τf , when the volume of the current goes to zero (D). In the early period, the shock strength and position are independent of ǫ and are given by √ √ 1 √ ηs = , ξs = − 21 + ( M − τ )2 . (4.21a, b) 1 − M + Mτ 50 M. A. Hesse, F. M. Orr Jr and H. A. Tchelepi We have shown in § 4.1 that the receding section of the interface is not affected by trapping. Therefore, the shock strength and position are independent of the amount of trapping. The early period extends until the shock migrates past the initial extent of the current, ξs < 1/2. The passing time and the shock strength at the passing time are given by √ 1 √ . (4.22a, b) τp = (1 + M)2 , ηs [τp ] = 1+ M In the absence of trapping, (4.21a, b) remain valid for all τ > τt (dashed line in figure 8). After τp , an explicit expression for ηs and ξs is obtained only for M = 1. We give the full details for the case M = 1 to illustrate the structure of the solution, and for M > 1 we integrate (4.20) numerically, using (4.22) as an initial condition. For τp 6 τ 6 τf , the shock strength for M = 1 is given by ⎧ 1+1/(ǫ−1) −1−1/(2(ǫ−1)) t ⎨ǫ − 2 , ǫ = 21 , 2ǫ − 1 (4.23) ηs = ⎩  √ 1 1 ǫ = 2. log 2/ t + 2 , For ǫ = 1/2, the right-hand side of (4.20) loses the explicit dependence on ηs , leading to a simplified expression. The expression for the shock path is given by ⎧ ⎨ (ǫ − 1) (t − 22+1/(ǫ−1) t 1/(2(1−ǫ)) ) + 1 , ǫ = 1 , 2 2 ξs = 2ǫ − 1 (4.24) ⎩1 1 (log(t/4)t + 1), ǫ = . 2 2 In the second phase, the shock speed and position depend strongly on the amount of trapping. Figure 8(a) illustrates the effect of increasing trapping on the current evolution in the (ξ, τ )-plane for M = 1, and similar behaviour is observed numerically for other M (figure 7b). From (4.24), (4.23) and (4.4), the aspect ratio of the current is given by ⎧ ⎪ ⎨1 + (1 − ǫ)τ, τ < 1, ξo − ξs 1 < τ 6 4, = 2τ − ǫτ 3/2 , (4.25) A= ⎪ ηs ⎩2(1 − ǫ)τ, 4 < τ < τf , for M = 1. Even in the presence of trapping, A is always increasing for M > 1. The final residual surface, ηr , is the interface that divides the part of the aquifer that has been invaded by the current, and now contains trapped fluid, from the region of the aquifer that has not been invaded. The final residual surface for ξ < 1/2 is the initial condition. For ξ > 1/2, the final residual surface is defined as ηr = ηs [ξs ] and can be calculated by eliminating τ between (4.23) and (4.24) to obtain ξ = ξ [ηr ]. However, it is not possible to invert this relationship to obtain an explicit expression for ηr ⎧ ⎨ 1 + (ǫ − 1) (ϕ(η ) − 22+1/(ǫ−1) ϕ(η )1/(2−2ǫ) ), ǫ = 1 . r r 2 2ǫ − 1 ξ= 2 (4.26) ⎩ 1−2ηr 1 1 ǫ = 2. (2 − 4ηr ) + 2 , e where ϕ[ηr ] = (2−ǫ/(ǫ−1) (−2ǫηr + ηr + ǫ))1/(2ǫ−1)−1 . The shape of the final residual surface is shown in figure 8(c) for increasing residual saturations. 4.2.2. Final migration time and distance The finite migration distance is an important characteristic of the hyperbolic limit, in contrast to the self-similar solutions for the horizontal case. The maximum migration Gravity currents with residual trapping (a) 20 51 (b) 20 15 15 8 8 6 M 10 6 10 4 4 2 5 5 2 0 1 0 0.2 0.4 0.6 0.8 ǫ 1.0 1 0 0.2 0.4 0.6 0.8 1.0 ǫ Figure 9. Contour plots of (a) the up-dip migration distance, log(ξ↑ ), and (b) of the total migration time, log(τf ). The data for several potential storage aquifers in Alberta, Canada, are shown as circles to indicate typical parameter values (Bachu & Bennion 2007). distance and the corresponding total migration time are amongst the most important time and length scales for CO2 storage, and the analytic and semi-analytic results (figure 9) presented here show how these scales depend on the governing parameters. The total migration distance, ξ↑ , is given by 1 − 21+1/(2ǫ−1) (ǫ − 1)ǫ 1/(2ǫ−1)−1 , ǫ = 12 , 2 ξ↑ = 1 (4.27) + 2e, ǫ = 12 . 2 Increasing mobility increases ξ↑ , because a gravity tongue forms at the top of the aquifer and only a small fraction of the aquifer is swept and contributes to residual trapping. The final migration time for M = 1 is given by  −1−1/(ǫ−1) −2(ǫ−1)/(2(ǫ−1)+1) (2 ǫ) , ǫ = 21 , τf = (4.28) 4e, ǫ = 21 . The dimensionless final migration time increases strongly with increasing M. The dimensionless time is based only on the properties of the released fluid. As M increases, the ambient fluid becomes increasingly less mobile, which retards the movement of the released fluid. 4.2.3. Volume evolution The final migration time, τf , bounds the time available for leakage from the storage reservoir. For smaller times, the volume of mobile CO2 is a basic quantity necessary to estimate potential leakage over time. Analytic expressions for the evolution of the current volume are given by integrating (4.4) over the extent of the plume. The normalized current volume V = V (τ )/V (τ0 ) is given by  ξo √ η dξ = 1 − τ ǫ(1 + M)−2 V= (4.29a) ξs for τ < τp . In this early period, the evolution of both the shock and the receding part of the rarefaction is independent of ǫ. Therefore the swept area containing residual saturation is also independent of ǫ, and the current volume decreases in proportion to ǫ. The decay of the current volume slows down with increasing M, 52 M. A. Hesse, F. M. Orr Jr and H. A. Tchelepi (a) 1.0 (b) 1.0 0.8 0.8 0.6 0.6  0.4 0.2 0 0.4 0.1 0.3 0.4 1 0.2 0.2 25 τ 50 0 5 10 15 20 25 40 100 200 300 τ Figure 10. The effect of the governing parameters on the evolution of the current volume. The passing and final migration times, are shown as circles and triangles. (a) Increasing the residual, ǫ, from 0.1 to 0.9 at constant M = 1. (b) Increasing the mobility ratio, M, from 1 to 40 at constant ǫ = 0.4. because the receding interface slows down. This behaviour is also observed in semianalytic solutions for τ > τp . At later times the lower limit of the integral is known only for M = 1, and in this case the current volume is given by ⎧ 2 1 − ǫ  −ǫ/(1−ǫ) 1/(2(1−ǫ)) ⎪ ⎪ t − ǫt , ǫ = 21 , 2 ⎨ 2 t (1 − 2ǫ) (4.29b) V=    2 ⎪ ⎪ 1 ⎩ t log t − 1 , ǫ = 2, 8 4 for τ < τf = 4e. Figure 10(a) shows how strongly the current volume is affected by the magnitude of ǫ. This evolution is distinctly different from the power-law decay obtained from the self-similar solutions for horizontal aquifers. It suggests that sloping aquifers reduce the current volume particularly effectively. An evolution similar to (4.29) is observed for a wide range of M (figure 10b). The residual volume is given by Vr = 1 − V, and from (4.29) we see that the volume of trapped CO2 increases linearly initially, but more slowly for τ > τp . The efficiency of residual trapping decreases with time, because the size of the receding interface, ηs , decreases rapidly, once the shock interacts with the advancing section of the rarefaction. 5. Numerical results In this section, we present numerical solutions to (2.9) and show that the current evolution for intermediate values of Pe can be understood in terms of an early essentially parabolic period followed by a late near hyperbolic regime. We also discuss important length and time scales of the solution as a function of increasing Pe. The numerical solutions were computed using standard finite-volume techniques and explicit time integration (Leveque 2002). 5.1. Near hyperbolic behaviour for large M In CO2 storage, the Pe number is not necessarily large, but the mobility ratio is generally larger than unity, M ≈ 5. Here we show that the particular shape of f leads to near hyperbolic behaviour of (2.9a) for small Pe and large M. Consider the Gravity currents with residual trapping 53 (a) 0.2 0 1 η 0.1 M = 20 0 15 10 10 5 20 30 40 50 hyp. num. sim. 60 (b) 0.1 hyp. num. η 0 10 ǫ=0 0.1 0.2 0.3 20 30 40 50 60 ξ Figure 11. Several comparisons between the numerical solution (ξ = 0.01) to (2.9a) and the hyperbolic solution for Pe = 2 at τ = 50. (a) Increasing M at constant ǫ = 0. (b) Increasing ǫ at constant M = 10. differentiated form of (2.9a) given by σ −1 ητ + ληξ = Pe −1 (λ(ηξ )2 + f ηξ ξ ). (5.1) In the absence of source terms, the thickness η of the current decreases monotonically with time, and for M ≫ 1, the current takes the shape of a triangle composed of a very elongate, flat, advancing tongue in the front and a short, steep, receding tongue in the back (figure 1g–h). This asymmetry increases with Pe, M and τ . The hyperbolic approximation introduced in § 4 replaces the short, steep, receding interface of the current by a shock. In the advancing tongue ηξ ξ ≪ ηξ ≪ 1, and hence both the diffusive, f ηξ ξ , and the second-order term, ληξ2 are small. The diffusive term is reduced further over the √ entire current, because max(f ) = 1/(1 + M)2 decreases with increasing M. The second-order term ηξ2 is not small in the steep receding back of the current, and figure 11 shows that the error in the hyperbolic approximation is localized there. Figure 11(a) illustrates how the error introduced by the hyperbolic approximation decreases with increasing M. For M > 5 and τ > 2τp , the only significant error in the hyperbolic approximation is due to a small off-set in the shock position. In the limit M = 0, a similarity solution to the full equation is available (Huppert & Woods 1995), and it is shown to provide an indication of the accuracy of the numerical method used. Once loss is introduced, ǫ > 0, the small differences in the sweep of the full solution and the hyperbolic solution lead to differences in plume volume at any given time. This volume error accumulates over time and leads to differences between the shock location in the hyperbolic model and the back of the plume in the numerical solution (figure 11b). 5.2. Transition from parabolic to hyperbolic behaviour In gently sloping media, the current transitions from initial parabolic behaviour to hyperbolic behaviour at late times (figure 12). Owing to the distinctly different behaviour of these limiting cases for ǫ > 0, the transition is apparent in the evolution of the current volume. An initial power-law decay of the current volume is followed by a much more rapid decay characteristic of the hyperbolic limit. 54 M. A. Hesse, F. M. Orr Jr and H. A. Tchelepi 37 373 t (years) 3730 37 300 373 000 100 100 10–1 10–1  10–2 10–3 10–1 Full equation: Pe = 0 (θ = 0°) Pe = 0.44 (θ = 12 °) Pe = 0.87 (θ = 1°) Pe = 4.37 (θ = 5°) Pe = 8.82 (θ = 10°) Pe = 13.4 (θ = 15°) Limiting cases: similarity soln hyperbolic soln 100 10–2 101 102 10–3 103  = t/t d = tκ1H cos θ/L2 Figure 12. Evolution of the current volume, V, as a function of diffusive dimensionless time, ϑ, for M = 5 and increasing Pe. The time axis at the top gives dimensional values corresponding to the parameter choice discussed in § 6, assuming cos θ ≈ 1. In the limit of large Pe, the plume volume is given by (4.29) and its semi-analytic extension to M = 1. This limiting solution is plotted in figure 12 as a function of the diffusive dimensionless time ϑ = τ/Pe, so that changes in Pe correspond to a stretching of the time axis. We observe that the limiting solution describes the rapid decay of the current volume at late times, even for Pe < 1. 5.3. Length and time scales for gravity currents with residual The length and time scales that are important for gravity currents in sloping layers are the propagation distance against gravity, ξ↓ , and the time of the associated reversal in the direction of the interface movement, τt . For a buoyant current the down-dip propagation of the interface is due to the across-slope component of gravity described by the parabolic part of the equation. The reversal of the interface movement indicates the transition from early near-parabolic to late near-hyperbolic behaviour. The continued volume loss and the resulting extinction of the current give rise to additional length and time scales, the up-dip migration distance, ξ↑ , and the final migration time, τf . Numerical results for all four scales as a function of increasing Pe are shown in figure 13. To determine the scales associated with the end of the current, the numerical simulation was terminated at V = 10−3 . This value was chosen to allow a reasonable detection of the gravity current tips, which become strongly affected by numerical diffusion for smaller cutoff volumes. The rapid decay of the plume volume at late times (figure 12) suggests that the results are not very sensitive to the cutoff, and the reasonable agreement between numerical and limiting analytic results in figure 9 confirms this. The down-dip migration √ distance, ξ↓ , is not very sensitive to M or ǫ. For small Pe, it follows the relation τt = 1/ 6Pe obtained for Gravity currents with residual trapping (b) 20 (a) 2 Numerical: M=0 M=1 M=5 M = 10 M = 15 M = 20 –ξ↓ 1 τt 10 similarity soln hyperbolic soln 0 5 10 15 0 20 (c) 100 (d) 160 75 120 ξ↑ 50 τf 0 5 10 15 20 5 10 Pe 15 20 80 40 25 0 55 5 10 Pe 15 20 0 Figure 13. Time and length scales obtained from the numerical solution of (2.9) are shown as symbols (ξ = 0.01). All results are for a residual of ǫ = 0.4. The limiting hyperbolic and self-similar solutions are indicated by full and dashed lines, respectively. (a) Down-dip migration distance, ξ↓ . (b) Transition time, τt . (c) Up-dip migration time, ξ↑ . (d ) Final migration time, τf . the similarity solution for M = ǫ = 0. As τt increases with decreasing Pe, the current thickness decreases as well, and the similarity solution for the unconfined aquifer becomes a better approximation. As Pe increases, ξ↓ approaches the position of the initial condition, which is also the down-dip migration distance in the hyperbolic limit. The transition time, τt , appears to approach the limiting hyperbolic solution much slower than the other parameters. The up-dip migration distance, ξ↑ , and the final migration time, τf , show a more complicated behaviour as Pe increases. As Pe goes to zero both become very large, which is consistent with the late similarity solution for the limiting case Pe = 0. As Pe increases, the hyperbolic limit is approached quickly, except for M = 0. However, for Pe ≈ 0.2 − 0.5, a minimum in both quantities is observed for M > 1. Hesse et al. (2006) showed that both ξ↑ and τf decrease with increasing Pe for M = 0. Figure 9(c, d ) shows a similar behaviour of ξ↑ and τf for M < 1, but we observe an increase of ξ↑ and τf for M > 1 for Pe > 1. 6. Discussion 6.1. Assumptions The most limiting assumptions in the derivation of the governing equations are the homogeneity of the porous medium, the uniform saturation in the CO2 plume, and the incompressibility of the CO2 . Although simple forms of heterogeneity in 56 M. A. Hesse, F. M. Orr Jr and H. A. Tchelepi the across-slope direction can be included (Huppert & Woods 1995), general forms of heterogeneity require a full numerical solution. Ide et al. (2007) investigated gravitational spreading of a CO2 plume in a horizontal two-dimensional heterogeneous aquifer with a reservoir simulator. They showed that the effect of capillary forces on residual trapping is larger than the effect of moderate permeability heterogeneity. In the absence of capillary forces, their numerical simulations resemble figure 1(a–d ). Simulations including capillary forces lead to a thicker current with a non-uniform vertical saturation profile and an increase in capillary trapping (Mo et al. 2005). It may therefore be interesting to extend the model presented here to non-uniform saturation distributions. The infinitely spreading solution obtained for horizontal aquifers is a direct result of neglecting capillary forces. The current will stop propagating once the hydrostatic potential at the tip of the plume decreases below the entry pressure necessary for the CO2 to invade the largest pores. Below we show that even in very gently dipping aquifers, the CO2 may rise vertically by several hundred metres, if the up-dip migration distance is long. In these cases, the expansion of the CO2 will lead to an increase in plume volume that will be significant, in particular if the CO2 becomes gaseous. In the hyperbolic model for sloping aquifers, we neglect the second-order term, λ(ηξ )2 . This simplification leads to the unphysical vertical interface at the trailing end of the current. Despite this simplification, we find that the resulting quasi-linear hyperbolic model describes many characteristics of the current evolution very well. The minimum in the migration distance and time observed in figure 13 is not explained by the limiting cases. The most likely explanation is the down-dip migration of the current against gravity. This down-dip migration is represented by the second-order term, λ(ηξ )2 . Inclusion of this term leads to a fully nonlinear hyperbolic model σ −1 ητ + ληξ = Pe −1 λ(ηξ )2 . (6.1) This model could still be solved by the method of characteristics allowing fast and accurate solutions to be obtained by numerical integration. In terms of geological CO2 storage, we have neglected the dissolution of CO2 into the brine. Two types of dissolution must be distinguished, direct dissolution of CO2 into the residual brine and the enhanced dissolution of CO2 at the interface between CO2 and brine owing to convective currents in the underlying brine (Ennis-King et al. 2005; Riaz et al. 2006). The direct dissolution into the residual brine can be modelled in this framework, but the effect is likely to be small compared to the residual trapping. The enhanced dissolution due to convective transport may become important at later times, because the gravity tongue provides a large interface. This process could be modelled conceptually through an effective loss term, but only preliminary data are available to estimate the magnitude of this loss term (Hesse, Riaz & Tchelepi 2007a). An extension of the model presented here would allow a comparison of the relative contribution of dissolution and residual trapping. The leading edge of the CO2 plume is both gravitationally and viscously unstable and this leads to the formation of a gravity tongue in our solutions. Riaz & Tchelepi (2006) have shown that this gravity tongue is also the dominant feature in high-resolution numerical simulations of two-dimensional vertical displacements. The displacement may also be viscously unstable in the transverse direction, similar to the Saffman–Taylor instability in a Hele-Shaw cell (Saffman & Taylor 1958), or the instability of viscous currents on inclined planes discussed by Lister (1992). The linear instability of an immiscible displacement in a porous medium has been studied by Yortsos & Hickernell (1989) and Riaz & Meiburg (2004). However, these Gravity currents with residual trapping 57 studies consider advection-dominated one-dimensional base flows and not the twodimensional base flow in gravity–capillary equilibrium that is of interest here. The study of the stability of the two-dimensional flows presented here is an interesting direction for future work, but is beyond the scope of this paper. 6.2. Implications for CO2 storage in saline aquifers Nicot (2007) considered the storage of roughly a fifth of the CO2 emissions from coal power stations in the state of Texas over the next 50 years. This would require the injection of approximately 370 × 106 m3 per year of supercritical CO2 , equivalent to 50 MtCO2 per year, into the central section of the Carrizo–Wilcox aquifer. He envisaged a line of 50 wells aligned perpendicularly to the dip of the aquifer with a 1.6 km spacing. These wells act as a 80 km line source, and our one-dimensional solution is a reasonable approximation of the displacement away from the edges. To obtain first-order estimates of the migration distance of the injected CO2 , we assume the following constant aquifer properties estimated from the data given in Nicot (2007): H ≈ 200 m, k ≈ 500 mD, φ ≈ 0.15, the average depth of injection 2.7 km, the average distance to the outcrop 100 km, and therefore an average dip angle, θ ≈ 1.5◦ . We assume the following fluid properties: ρ ≈ 300 kg m−3 , Scr =Sbr = 0.2, μb /μc ≈ 10, krc ≈ 0.2, krb = 1 (at the advancing interface). Given these data, the appropriate length scale at the end of 50 years of injection is L ≈ 10 km, so that the initial aspect ratio is A = L/H ≈ 50. The three governing parameters are therefore M ≈ 5, ǫ ≈ 0.25 and Pe = A tan θ ≈ 1.4. Figure 13 shows that the migration distance and time for M = 5 are already close to the hyperbolic limit for Pe = 1.4, so that we obtain ξ↑ ≈ 80 and τf ≈ 110 from figure 9. The along-slope migration distance would be approximately ξ↑ L = 800 km, so that the migration distance exceeds the up-dip extent of the aquifer significantly. Under the simplifying assumptions made here, this particular injection scenario does not ensure that all CO2 would be trapped by residual saturation alone. The initial length scale L would have to be reduced to 100 km/ξ↑ = 1.25 km to achieve residual trapping of all injected CO2 . This could be achieved by increasing the width of the injection zone from 80 to approximately 620 km. Another option would be to reduce the volume injected into this aquifer from 50 to roughly 6.5 MtCO2 per year, assuming the original injection width of 80 km. In this case, the CO2 would remain mobile for 1550 years after the end of injection. The Carrizo–Wilcox formation has a relatively high permeability for deep saline aquifers and is relatively thick. Therefore, the time estimated for trapping all injected CO2 is shorter and the storage volume is larger than expected for other deep aquifers. In many cases, the properties of deep saline aquifers are closer to: k ≈ 30 mD, φ ≈ 0.1, H ≈ 20 m, and the dip angle may be as small as 1/2◦ (Michael et al. 2006). In this case, the storage volume is less than 10 % of that discussed above for the same width of injection. To assist comparison, we assume all other properties are the same, so that P e ≈ 1/2. The volume evolution for these cases and various larger angles is shown in figure 12, and the evolution of the interface shape for θ = 0◦ and 5◦ is also shown in figure 1. We see that the migration time of the current is reduced from approximately 40 000 years for a slope of 1/2◦ to approximately 2000 years for a slope of 15◦ . The up-dip migration distance decreases from 85.6 km at 1/2◦ to a minimum of 82.8 km at 1◦ and slowly increases again to 84.5 km at 15◦ . As the slope of the aquifer increases, the distance the plume can migrate before it reaches the surface decreases. Assuming an injection depth of 2000 m, only angles less than 2◦ prevent leakage in the example considered above, and even for θ = 1/2◦ the 58 M. A. Hesse, F. M. Orr Jr and H. A. Tchelepi (b) 20 (a) 20 0.1° 15 0.5° 0.01% 15 1° 0.1% 2° M 10 4° 10 6° 1% 10° 1 10% 20° 35° 90° 5 0 0.2 0.4 0.6 ǫ 0.8 5 1.0 1 50% 90% 0 0.2 0.4 0.6 0.8 1.0 ǫ Figure 14. Contour plots of (a) the critical angle, θc , assuming d0 = 2 and (b) the sweep, S, as functions of M and ǫ. The data for several potential storage aquifers in Alberta, Canada, is shown as circles to indicate typical parameter values (Bachu & Bennion 2007). CO2 rises 750 m vertically. Therefore, the depth of the current at the end of migration is the important criterion for storage security. Given the dimensionless initial depth d0 = D0 /L, the dimensionless final depth is d = −d0 + ξ↑ sin θ. We can define the critical angle at which the injected CO2 reaches the surface as   d0 . (6.2) θc = arcsin ξ↑ Figure 14(a) shows that the critical angle, θc , increases with ǫ and decreases with M. For a given M and ǫ the CO2 will not reach the surface for θ < θc . Despite rapid trapping compared to the horizontal case, CO2 storage in saline aquifers is limited by the poor vertical sweep, S, and the resulting long migration distances. The vertical sweep can be defined as S = hr  /H = ηr , where hr  is the average residual surface. The initial volume of mobile CO2 is given by φ(1 − Sbr ) HL, and the residual volume of CO2 at the end of the migration is given by φScr hr  (x↑ − x↓ ). An expression for the vertical sweep, S, can be obtained from the conservation of volume and is given by ǫ , (6.3) S= ξ↑ − ξ↓ where ξ↓ ≈ 0 in the hyperbolic limit. Figure 14(b) shows that the sweep is a strong function of ǫ, increasing monotonically from zero at ǫ = 0 to unity for ǫ = 1. For the data from Bachu & Bennion (2007), the sweep can vary from less than 10−3 to more than 10−1 , roughly three orders of magnitude, while the residual, ǫ, itself only varies between 0.1 and 0.7. Therefore the residual CO2 saturation, Scr , influences the up-dip migration distance given by (1.1) more strongly through its influence on the sweep than through its direct volumetric effect. 7. Conclusion We have derived a vertical-equilibrium sharp-interface model governing the migration of a gravity current with residual saturations in a two-dimensional confined aquifer. This model allows an analysis of the effect of the imperfect vertical sweep on the migration distance and time of the current. Gravity currents with residual trapping 59 The parabolic limit of the governing equations, corresponding to a horizontal aquifer, allows self-similar solutions at early and late times. The spreading of the current and the decay of the volume are given by power laws. At late times, the exponent of the power law depends on the residual, ǫ. Despite continued volume loss, the current volume is never reduced to zero, and the current continues to spread. In the hyperbolic limit of the governing equations, the current is described by a travelling wave with hysteresis. The leading edge of the interface is a rarefaction, and the trailing edge is a shock. The current volume is reduced continuously by the residual, and the volume is reduced to zero in finite time. We have obtained expressions for the up-dip migration distance and the final migration time of the current. In both the parabolic self-similar solutions and the hyperbolic rarefaction and shock, we observe that the residual trapping affects the advancing section of the interface more strongly than the receding section of the interface. Although trapping occurs at the receding interface, the advancing interface is affected, because less fluid is supplied. In gently sloping aquifers, an initial near-parabolic regime with power-law decay of the volume is followed by a near-hyperbolic regime with very rapid volume decay. Increasing the slope of the aquifer or the initial aspect ratio of the current reduces the duration of initial parabolic period. Our results suggest that lateral migration of the injected CO2 along the seal will trap the CO2 relatively quickly as residual saturation. Residual trapping is quite effective in sloping aquifers with small mobility ratios and high residual CO2 saturations. However, the long migration distances of CO2 , owing to the formation of a gravity tongue, may limit the volume of CO2 that can be stored in sloping regional aquifers. This research was supported by the Global Climate and Energy Project (GCEP) and the SUPRI-B reservoir simulation affiliates program at Stanford University. REFERENCES Bachu, S. & Bennion, B. 2007 Effects of in-situ conditions on relative permeability characteristics of CO2 –brine systems. Environ. Geol. doi 10.1007/s00254-007-0946-9. Barenblatt, G. I. 1996 Scaling, Self-Similarity, and Intermediate Asymptotics. Cambridge University Press. Bear, J. 1972 Dynamics of Fluids in Porous Media. Elsevier. Bear, J. & Ryzhik, V. 1998 On displacement of NAPL lenses and plumes in a phreatic aquifer. Transport Porous Med. 33, 227–255. Benson, S. M., Tomutsa, L., Silin, D. & Kneafsey, T. 2006 Core scale and pore scale studies of carbon dioxide migration in saline formations. In Proc. GHGT-8, Trondheim, Norway. Bickle, M., Chadwick, A., Huppert, H. E., Hallworth, M. & Lyle, S. 2007 Modelling carbon dioxide accumulation at Sleipner: implications for underground carbon storage. Earth Planet. Sci. Lett. 255, 164–176. Ennis-King, J. & Paterson, L. 2002 Engineering aspects of geological sequestration of carbon dioxide. In Asia Pacific Oil and Gas Conf. and Exhibition, Melbourne, Australia. Ennis-King, J., Preston, I. & Paterson, L. 2005 Onset of convection in anisotropic porous media subject to a rapid change in boundary conditions. Phys. Fluids 17(8), 084107. Hassanzadeh, H., Pooladi-Darvish, M. & Keith, D. W. 2007 Scaling behavior of convective mixing, with application to geological storage of CO2 . AIChE J. 53(5), 1121–1131. Hesse, M. A., Tchelepi, H. A. & Orr, Jr, F. M. 2006 Scaling analysis of the migration of CO2 in saline aquifers. In SPE Annu. Tech. Conf. and Exhibition, San Antonio, TX. Hesse, M. A., Riaz, A. & Tchelepi, H. A. 2007a Convective dissolution of CO2 in saline aquifers. Geochim. Cosmochim. Acta 71 (15), A401. 60 M. A. Hesse, F. M. Orr Jr and H. A. Tchelepi Hesse, M. A., Tchelepi, H. A., Cantwell, B. J. & Orr, Jr, F. M. 2007b Gravity currents in horizontal porous layers: transition from early to late self-similarity. J. Fluid Mech. 577, 363–383. Holloway, S & Savage, D. 1993 The potential for aquifer disposal of carbon-dioxide in the UK. Energy Convers. Manage. 34(9–11), 925–932. Huppert, H. E. & Woods, A. W. 1995 Gravity-driven flows in porous media. J. Fluid Mech. 292, 55–69. Ide, S. T., Jessen, K. & Orr, Jr, F. M. 2007 Storage of CO2 in saline aquifers: effects of gravity, viscous, and capillary forces on amount and timing of trapping. Intl J. Greenh. Gas Control 1(4), 481–491. Juanes, R., Spiteri, E. J., Orr, Jr, F. M. & Blunt, M. J. 2006 Impact of relative permeability hysteresis on geological CO2 storage. Water Resour. Res. 42 (W12418), 1–13. Kochina, I. N., Mikhailov, N. N. & Filinov, M. V. 1983 Groundwater mound damping. Intl J. Engng Sci. 21, 413–421. Koperna, G. J. & Kuuskraa, V. A. 2006 Assessing technical and economic recovery of oil resources in residual oil zones. Tech. Rep. Advanced Resources International, 4501 Fairfax Drive, Suite 910, Arlington, VA 22203 USA. Kumar, A., Ozah, R., Noh, M., Pope, G. A., Bryant, S., Sepehrnoori, K. & Lake, L. W. 2005 Reservoir simulation of CO2 storage in deep saline aquifers. Soc. Petrol. Engrs J. September pp. 336–348. Leveque, R. J. 2002 Finite Volume Methods for Hyperbolic Problems. Cambridge University Press. Lindeberg, E. & Wessel-Berg, D. 1997 Vertical convection in an aquifer column under a gas cap of CO2 . Energy Conserv. Manage. 38, SS229–SS234. Lister, J. R. 1992 Viscous flows down an inclined plane from point and line sources. J. Fluid Mech. 242, 631–653. Lyle, S., Huppert, H. E., Hallworth, M., Bickle, M. & Chadwick, A. 2005 Axisymmetric gravity currents in a porous medium. J. Fluid Mech. 543, 293–302. Metz, B., Davidson, O., de Coninck, H., Loos, M. & Meyer, L. (ed.) 2006 Special Report on Carbon Dioxide Capture and Storage. Cambridge University Press. Michael, K., Bachu, S., Buschkuehle, B. E., Haug, K., Grobe, M. & Lytviak, A. T. 2006 Comprehensive characterization of a potential site for CO2 geological storage in central Alberta, Canada. In Proc. CO2SC Symp. 2006, Berkeley, CA. Mo, S., Zweigel, P., Lindeberg, E. & Akervoll, I. 2005 Effect of geologic parameters on CO2 storage in deep saline aquifers. In SPE Eurospec/EAGE Annu. Conf. Madrid, Spain. Nicot, J.-P. 2007 Evaluation of large-scale carbon storage on fresh-water sections of aquifers: a Texas study. Intl J. Greenh. Gas Control (in press). Nordbotten, J. M., Celia, M. A. & Bachu, S. 2005 Injection and storage of CO2 in deep saline aquifers: analytical solution for the CO2 plume evolution during plume injection. Transport Porous Med. 58, 339–360. Riaz, A. & Meiburg, E. 2004 Linear stability of radial displacements in porous media: influence of velocity-induced dispersion and concentration-dependent diffusion. Phys. Fluids 16 (10), 3592. Riaz, A. & Tchelepi, H. A. 2006 Numerical simulation of immiscible two-phase flow in porous media. Phys. Fluids 18 (014104). Riaz, A., Hesse, M., Tchelepi, H. & Orr, Jr, F. M. 2006 Onset of convection in a gravitationally unstable, diffusive boundary layer in porous media. J. Fluid Mech. 548, 87–111. Saffman, P. G. & Taylor, G. I. 1958 The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more visous fluid. Proc. R. Soc. Lond. A 245, 312–329. Vella, D. & Huppert, H. E. 2006 Gravity currents in a porous medium at an inclined plane. J. Fluid Mech. 555, 353–362. Verdon, J. & Woods, A. W. 2007 Gravity driven reacting flows in a confined porous aquifer. J. Fluid Mech. 588, 29–41. Yortsos, Y. C. 1995 A theoretical analysis of vertical flow equilibrium. Transport Porous Med. 18, 107–129. Yortsos, Y. C. & Hickernell, F. J. 1989 Linear stability of immiscible displacement in porous media. SIAM J. Appl. Maths 49, 730–748.