Influence of capillary-pressure models on CO2 solubility trapping

Advances in Water Resources 62 (2013) 488–498

Contents lists available at ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Influence of capillary-pressure models on CO2 solubility trapping Boxiao Li ⇑, Hamdi A. Tchelepi, Sally M. Benson Department of Energy Resources Engineering, Stanford University, 367 Panama St, Rm 65, Stanford, CA 94305, USA

a r t i c l e

i n f o

Article history: Available online 20 August 2013 Keywords: CO2 sequestration Numerical simulation Brooks–Corey Van Genuchten Capillary entry pressure CO2 solubility trapping

a b s t r a c t The typical shape of a capillary-pressure curve is either convex (e.g., Brooks–Corey model) or S-shaped (e.g., van Genuchten model). It is not universally agreed which model reflects natural rocks better. The difference between the two models lies in the representation of the capillary entry pressure. This difference does not lead to significantly different simulation results for modeling CO2 sequestration in aquifers without considering CO2 dissolution. However, we observe that the van-Genuchten-type capillary-pressure model accelerates CO2 solubility trapping significantly compared with the Brooks–Corey-type model. We also show that the simulation results are very sensitive to the slope of the van-Genuchtentype curve around the entry-pressure region. For the representative examples we study, the differences can be so large as to have complete dissolution of the CO2 plume versus persistence of over 50% of the plume over a 5000-year period. The cause of such sensitivity to the capillary-pressure model is studied. Particularly, we focus on how the entry pressure is represented in each model. We examine the mass-transfer processes under gravitycapillary equilibrium, molecular diffusion, convective mixing, and in the presence of small-scale heterogeneities. Laboratory measurement of capillary-pressure curves and some important implementation issues of capillary-pressure models in numerical simulators are also discussed. Most CO2 sequestration simulations in the literature employ one of the two capillary-pressure models. It is important to recognize that these two representations lead to very different predictions of long-term CO2 sequestration. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Carbon dioxide capture and sequestration (CCS) is one important approach to reduce anthropogenic CO2 emissions [1]. After being captured from sources like power plants, the CO2 is transported to, and injected into underground geological formations, such as deep saline aquifers. While trapping beneath a seal with a high capillary entry pressure is the primary means of retaining CO2 underground, secondary trapping mechanisms can further increase long term storage security [1]. The dissolution of CO2 in the resident water, also known as CO2 solubility trapping, is an important secondary CO2 trapping mechanism in the postinjection period. Over periods of hundreds to thousands of years, depending on the properties of the storage reservoir, a significant fraction of the CO2 can dissolve [2,3]. Modeling CO2 solubility trapping requires accurate description of the CO2 plume migration and the saturation distribution. Capillary pressure (P c ) plays an important role, especially during the post-injection period [4]. Therefore, an accurate representation of the capillary pressure is

crucial for modeling CO2 solubility trapping, and hence long-term CO2 sequestration. Here, we briefly introduce the model for multiphase flow of CO2 and water in porous media. The mass-conservation equations for CO2 and water are:

@ ð/qw Sw Þ þ r  ðqw uw Þ þ qw ¼ 0; @t @ ð/qn Sn Þ þ r  ðqn un Þ þ qn ¼ 0; @t

E-mail address: [email protected] (B. Li). 0309-1708/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.advwatres.2013.08.005

ð1bÞ

where t is the time, / is the porosity, Si ; qi ; ui , and qi (i ¼ w; n) are the saturation, density, flow rate, and source/sink of each phase. The subscript w and n denote the wetting (water) and nonwetting (CO2) phases, respectively. The flow rate of each phase is expressed using the multiphase extension of Darcy’s Law [5]:

uw ¼ k un ¼ k

⇑ Corresponding author. Tel.: +1 650 725 0830; fax: +1 650 725 2099.

ð1aÞ

krw

ðrpw þ qw g rhÞ;

ð2aÞ

ðrpn þ qn g rhÞ;

ð2bÞ

lw krn

ln

where k is the permeability, kri ; li ; pi (i ¼ w; n) are the relative permeability, viscosity, and pressure of each phase. g is the

489

B. Li et al. / Advances in Water Resources 62 (2013) 488–498

gravitational acceleration and h is the altitude. To close the massconservation equations, two more equations are needed:

Sw þ Sn ¼ 1;

ð3aÞ

Pc ¼ pn  pw ;

ð3bÞ

where Pc is the capillary pressure, which relates the pressures of CO2 and water. It is a function of saturation, often expressed as Pc ðSw Þ. Such capillary-pressure–saturation relationship, also called the capillary-pressure curve, can be measured in laboratories. Typically, the capillary-pressure curve is either S-shaped (e.g., van Genuchten model [6]) or convex (e.g., Brooks–Corey model [7]), as illustrated in Fig. 1(a) and (b). The two models represent the entry of the nonwetting phase into the porous medium differently. The Brooks–Corey model uses a plateau that ends with a nonzero capillary entry pressure, while the van Genuchten model uses a steep slope that connects the end-point (usually zero) to the plateau region. The end-point is defined as P c ðSw ¼ 1Þ. We refer to the entry-pressure representation in Fig. 1(a) as van-Genuchtentype (VG-type), and that in Fig. 1(b) as Brooks–Corey-type (BCtype). We name the region around the end-point in VG-type models the ‘entry-slope’ region. The height of this region is the pressure difference between the plateau and the end-point, and its width is given by a ‘threshold’ nonwetting-phase saturation, Snt , as illustrated in Fig. 1(a). Note that in the VG-type model, the Pc ðSw ¼ 1Þ can be zero as in Fig. 1(a), or nonzero. The capillary-pressure curves measured in the laboratory using mercury porosimetry usually have an entry-slope region similar to Fig. 1(a), though Pc ðSw ¼ 1Þ is usually nonzero. The physical meaning and importance of the entry-slope region has been the subject of several investigations [8–11]. Consensus, however, has not been reached about whether this region should be ignored or honored when interpreting the measurements. Ignoring the entry-slope region amounts to applying a BC-type representation, while preserving it amounts to a VG-type representation. Here, we show that the representation of the capillary entry pressure can affect the rate of CO2 solubility trapping significantly, even when the width of the entry-slope region (Snt ) is made very small. Large-scale simulations of CO2 sequestration in saline aquifers have been widely studied [12–16]. The sensitivity of CO2 solubility trapping to the representation of the capillary entry pressure (i.e., VG-type versus BC-type curves) has not been investigated. We note that the van Genuchten capillary-pressure model is almost exclusively used in the simulations performed by TOUGH2, a simulator that is widely used in modeling CO2 sequestration [17]. Besides large-scale behavior, it is also important to understand the physics that happen at smaller scales during CO2 solubility trapping, for example, the density-driven convective mixing. Analytical methods and small-scale high-resolution simulations have been widely used to model the onset of convection and the development of the convective fingers [3,18,20–22]. In most of these and related works, capillary pressure is neglected. Although recent work has investigated the role of the capillary transition zone (also

known as the capillary fringe) [23], the impact of using different capillary-pressure models is not studied. In this paper, we compare the BC-type and VG-type capillarypressure models and study their influences on CO2 solubility trapping. In particular, we focus on their representations of the capillary entry pressure. The problem is analyzed at both large and small scales, and a variety of factors are studied, including gravity-capillary equilibrium, molecular diffusion, convective mixing, and small-scale heterogeneities. Laboratory measurements of capillary-pressure curves and implementation issues of capillarypressure models in simulators are also discussed. 2. Simulation model The aquifer model used in this study is a 2D vertical cross-section (in the x and z directions) of the model proposed in [24], which is used as a benchmark for long-term CO2 sequestration simulations. The 2D model is illustrated in Fig. 2, where the size of the aquifer and the location of the injection well are also marked. The aquifer has a dip of 1%, and it is deep enough such that the injected CO2 is in supercritical phase, which is immiscible with the resident water. The top and bottom are no-flow boundaries. The initial in-situ pressures at the left and right boundaries are held constant. The grid information is listed in Table 1. The drainage relative-permeability (kri ; i ¼ w; n) and capillary-pressure curves provided in [24] are used here, namely, 4

krw ¼ ðSw Þ ;

2 2 krn ¼ 0:4 ½ 1  ðS^w Þ ð1  S^w Þ ;

Sw ¼ ðSw  Swi Þ=ð1  Swi Þ; S^w ¼ ðSw  Swi Þ=ð1  Swi  Snc Þ;

Pc ¼ P c;e 



Sw  Swi 1  Swi

ð4cÞ

0:5 ð5Þ

;

where the irreducible water saturation (Swi ) is 0.2, the critical CO2 saturation (Snc ) is zero, and the capillary entry pressure (Pc;e ) is 0.2 bar. The subscripts w and n denote the wetting (water) and nonwetting (CO2) phases, respectively. Eq. (5) is a BC-type capillary-pressure curve. The relative-permeability and capillary-pressure curves are also plotted in Fig. 3. The rock properties, salinity,

2-D (x-z direction) 1 % dip x z 50 m Closed boundary Constant pressure boundary Injection well Fig. 2. Sketch of the aquifer model.

Pc

Pc Table 1 Properties of the aquifer model. plateau

Snt Pc,e

Entry-slope region

0

(a)

Sw = 1

0

(b)

Sw = 1

Fig. 1. Capillary entry-pressure representations: (a) van-Genuchten-type (VGtype). The dashed circle highlights the entry-slope region; (b) Brooks–Corey-type (BC-type).

ð4aÞ ð4bÞ

Name

Value

Grid size Gridblock size Permeability Porosity Depth at the well Temperature Salinity

N x ¼ 150; N y ¼ 1; N z ¼ 40 Dx ¼ Dy ¼ 100 m, Dz ¼ 1:25 m kx ¼ ky ¼ kz ¼ 100 md / ¼ 0:15 3025 m 84.4 °C 0

490

B. Li et al. / Advances in Water Resources 62 (2013) 488–498

we do not use the original van Genuchten model as in [6], due to its difficulty in adjusting the steepness of the entry-slope region. Fig. 4 illustrates the distributions of the CO2 saturation (Sn ) and the mole fraction of CO2 in the aqueous phase (xCO2 ) over time. Because it is lighter than water, the injected supercritical CO2 rises to the top of the aquifer and migrates upward and to the right (the dip is not shown in the figure). As CO2 dissolves into the water, the plume becomes thinner and it eventually disappears completely. The density of the water increases as more CO2 dissolves in it. This ultimately causes the heavier water to sink downward into the lighter CO2–free water in the form of convective fingers that transport CO2–saturated water deep into the aquifer.

1.4

H2O

1.2

0.8

1 Pc (bar)

Relative permeability

1

0.6 CO2 0.4

0.8 0.6 0.4

0.2

0.2

0

0 0

0.2 0.4 0.6 0.8 Sw

0

1

0.2 0.4 0.6 0.8 Sw

(a)

1

(b)

Fig. 3. (a) Relative-permeability curves; (b) Capillary-pressure curve.

and the depth of the injection well are the same as those used in the benchmark model, and they are summarized in Table 1. The aquifer is initially fully saturated with water under hydrostatic equilibrium. Pure CO2 is injected at a rate of 9000 metric tons per year. Injection lasts for 20 years. The simulated time period is 5000 years. All the simulations in this study were performed using an Equation of State (EOS) based compositional simulator, namely, the Stanford General Purpose Research Simulator (GPRS) [25,26]. The flash calculation assumes instantaneous equilibrium between CO2 and water, which is a widely-used assumption in reservoir simulators. Detailed EOS tuning parameters for CO2–water system can be found in [27]. The Fully Implicit Method (FIM) [28] was used in all the simulations. The simulations shown here do not account for hysteresis in either the capillary-pressure, or the relative-permeability curves. We have performed simulations accounting for hysteresis, and have confirmed that the findings reported here will not be altered. 3. Simulation results 3.1. Base case The base case simulation uses the VG-type capillary entry-pressure representation (as in Fig. 1(a)). Specifically, the entry-pressure representation of Eq. (5) is transformed from BC-type into VGtype:

8  0:5 w Swi > < Pc;e S1S if Swi 6 Sw 6 1  Snt ; wi Pc ¼  0:5 > : Pc;e 1Snt Swi ð1  Sw Þ if 1  Snt < Sw 6 1: Snt 1S

ð6Þ

wi

The Snt is set as 0.005. P c;e and Swi are 0.2 bar and 0.2, respectively. Pc ðSw ¼ 1Þ is zero. For most part of the curve (Swi 6 Sw 6 1  Snt ), it is exactly equal to Eq. (5). At 1  Snt < Sw 6 1, however, it uses a steep slope to connect the plateau to the zero endpoint. Note that

3.2. Capillary entry-pressure representation Keeping all other parameters unchanged, simulations were performed using different representations of the capillary entry pressure. The capillary-pressure curves are shown in Fig. 5. The curve labeled ‘BC-type’ is the original Eq. (5). The curves ‘Snt ¼ 0:005’ (Base Case) and ‘Snt ¼ 0:0005’ employ a VG-type representation based on Eq. (6) with their corresponding Snt values. The ‘Silin’ curve is a VG-type capillary-pressure model proposed in [29]:

   a a 1=a2 Pc ¼ A ðSw Þ 1  1 þ B 1  ðSw Þ 2 ;

where is defined in Eq. (4b), and the fitting parameters are A ¼ B ¼ 0:2 bar, a1 ¼ 0:5, and a2 ¼ 9. All four curves in Fig. 5 are very similar, except for the small difference near Sw ¼ 1. We have also confirmed that the ‘capillary diffusion’ terms for the four cases,   kn  dP c  DðSw Þ ¼ kkwwþk  , are bounded within Sw 2 ½0; 1, where ki is the n dSw phase mobility, ki ¼ kri =li ; i ¼ w; n, and li is the viscosity. Severe numerical difficulties would occur if D were not bounded [30]. These apparently small differences in the capillary-pressure curves lead to very large differences in the long-term predictions. Fig. 6(a) and (b) plot the travel distance of the plume tip (defined by the white arrows in Fig. 4) and the plume volume versus time. While the ‘Snt ¼ 0:005’ (Base Case) indicates that complete dissolution of the CO2 plume occurs after 4000 years, less than half of the injected CO2 has dissolved in the ‘BC-type’ case. In addition, the simulation results are sensitive to the steepness of the entry-slope region. The plume disappearance is fastest for the ‘Silin’ case, and slowest for the ‘BC-type’ case. We also observed such strong sensitivity of the long-term simulation results to the capillary entrypressure representation using industrial reservoir simulators. We note that the sensitivity to the capillary entry-pressure representation is not obvious during the injection period, but is quite important in the long post-injection period. Note that although differences in the fluid pressure due to using different capillary entry-pressure representations are present, such

Mole fraction of CO2 inaqueous phase ( x CO 2 )

CO2 Saturation ( S n ) 0m

20 years 100 years

3000 years 4000 years

Sn 0.8 0.6 0.4 0.2 0.0

50 m 0m 50 m 0m

1000 years

ð7Þ

Sw

x CO 2

50 m 0m

0.032 0.024 0.016 0.008 0.000

50 m 0m 50 m 0

5000

10000

15000 m

5000

10000

15000 m

Fig. 4. Base case simulation results: CO2 saturation (Sn ) and the mole fraction of CO2 in aqueous phase (xCO2 ) distributions over time. The up-dip direction is from left to right (the dip is not shown here). The white arrows indicate the locations of the plume tip. Hysteresis is not simulated for simplicity.

491

B. Li et al. / Advances in Water Resources 62 (2013) 488–498

250000

0.3

1 Snt = 0.005 Snt=0.005

0.6

Snt = 0.0005 Snt=0.0005 Silin Silin

Plume Volume (m3)

0.8

Pc (bar)

Pc (bar)

BC-type Entry-plateau 0.2 0.1

0.4 0 0.99

0.2

0.995 Sw

1

Coarse, Virtual PcBC-type

200000

Coarse, Snt = 0.005 Snt = 0.005

150000

Middle, Virtual PcBC-type Middle, Snt = 0.005 Snt = 0.005

100000

Fine, BC-type Virtual Pc Fine, Snt = 0.005 Snt = 0.005

50000

0

0

0.2

0.4

0.6 Sw

0.8

1

0

1000

2000 3000 Time (year)

4000

5000

Fig. 8. Simulation results under different grid resolutions. Coarse: Dx ¼ 100 m; Middle: Dx ¼ 50 m; Fine: Dx ¼ 25 m.

differences (0.2 bar at most) are negligible compared with the overall fluid pressure level (about 300 bar), and they are very unlikely to affect the thermodynamic properties of the CO2 and water, including the CO2 solubility in water. It is important to note that such sensitivity to the entry-pressure representation is not observed when dissolution is not modeled. Fig. 7(a) and (b) show the travel distance and the extent (measured from tail to tip) of the CO2 plume, if CO2 is not allowed to dissolve in water. The simulation is terminated when the CO2 plume reaches the right boundary. As illustrated by the figures, the capillary entry-pressure representation has almost no impact on the results when dissolution is not modeled. So far, all the simulations are performed using a relatively coarse grid (see Table 1), where the gridblock size in the x direction, Dx, is 100 m. The simulation results where Dx is refined by two and four times are shown in Fig. 8. As can be seen, simulation

results are sensitive to grid resolution. The finer the grid, the faster is the CO2 dissolution. However, in all three resolutions, the VG-type entry-pressure representation consistently yields much faster CO2 dissolution than the BC-type. Studies with even finer resolutions are presented in the subsequent sections, where we show that such accelerated dissolution persists with grid refinement.

4. Cause of the sensitivity to entry-pressure representation In this section, we investigate why different entry-pressure representations affect the rate of CO2 solubility trapping. We first investigate CO2 dissolution and diffusion in a 1D vertical column subject to gravity and capillarity. Following this, 2D, two-phase simulations of density-driven convective mixing are studied.

8000

250000

7000 6000 5000

Virtual Pc BC-type

4000

Snt 0.005 Snt = 0.005

3000

Snt = 0.0005 Snt 0.0005

2000

Silin Silin

1000

Plume Volume (m3)

Plume Tip Travel Distance (m)

Fig. 5. Capillary-pressure curves with different entry-pressure representations. All the curves are very similar, except for the region around the entry-pressure.

200000 150000 100000 50000

0

0 0

1000

2000 3000 Time (year)

4000

5000

0

1000

(a)

2000 3000 Time (year)

4000

5000

(b)

10000

12000

8000

10000

6000

Virtual Pc BC-type Snt = 0.005 Snt 0.005

4000

Snt = 0.0005 Snt 0.0005

2000

Silin Silin

Plume Extent (m)

Plume Tip Travel Distance (m)

Fig. 6. The sensitivity of simulation results to different capillary entry-pressure representations: (a) Distance that plume tip traveled vs. time; (b) Plume volume vs. time.

8000 Virtual Pc BC-type

6000

Snt = 0.005 Snt 0.005

4000

Snt = 0.0005 Snt 0.0005

2000

0

Silin Silin

0 0

1000

2000 3000 Time (year)

(a)

4000

0

1000

2000 3000 Time (year)

4000

(b)

Fig. 7. Simulation results for different capillary entry-pressure representations when dissolution is not modeled: (a) Distance that plume tip traveled vs. time; (b) Plume extent vs. time.

B. Li et al. / Advances in Water Resources 62 (2013) 488–498

4.1. Gravity-capillary equilibrium 4.1.1. CO2 column height In the storage formation, the injected CO2 flows upward due to buoyancy, accumulates beneath the impermeable caprock, and forms a capillary transition zone (also known as the capillary fringe). The thickness of the transition zone is called the column height (hc ), measured by the vertical distance between the caprock and the bottom of the transition zone. The column height and the saturation distribution depend on the interaction between buoyancy and capillary forces. Fig. 9(a) illustrates the pressure fields of water and CO2 under gravity-capillary equilibrium, when a VG-type representation is used. Here, P c ðSw ¼ 1Þ is zero, as in Fig. 1(a). The pressures of the CO2 and water, and the capillary pressure can be written as:

pw ðzÞ ¼ pw jz¼hc  qw gðhc  zÞ;

ð8aÞ

pn ðzÞ ¼ pn jz¼hc  qn gðhc  zÞ;

ð8bÞ

Pc ðzÞ ¼ pn ðzÞ  pw ðzÞ;

ð8cÞ

where z is the depth measured positive downward from the top of the storage formation. Because the water and CO2 are under their own hydrostatic equilibrium, pw ðzÞ and pn ðzÞ are straight lines, the slopes of which are determined by the fluid densities qw and qn . The capillary pressure, defined as the pressure difference between CO2 and water, is zero at the bottom of the transition zone (z ¼ hc ), and below which only water is present. If the capillary-pressure curve employs the BC-type representation (as in Fig. 1(b)), the pressure fields of water and CO2 under gravity-capillary equilibrium are depicted in Fig. 9(b). The capillary pressure is equal to the capillary entry pressure (Pc;e ) at the bottom of the transition zone. Knowing Pc ðzÞ from Eq. (8) and the capillary-pressure–saturation relationship, the CO2 saturation distribution and column height under gravity-capillary equilibrium can be computed analytically, if the total amount of CO2 in the transition zone, V, is giR ven. The dimensionless variable V is defined as V ¼ H1 Sn ðzÞ dz, where H is the thickness of the aquifer. Analytical results and high-resolution numerical simulations were performed using different capillary entry-pressure representations. CO2 dissolution is not modeled here. The simulation domain is a 1D vertical column of 50 m (i.e., H ¼ 50 m). Each gridblock is 0.01 m in size, for which numerical diffusion effects are negligible. Two capillary-pressure curves, namely, the ‘Snt ¼ 0:005’ and the ‘BC-type’ cases in Fig. 5, are compared. The CO2 and water densities are qn ¼ 727 kg/m3 and qw ¼ 982 kg/m3, and the viscosities are ln ¼ 0:070 mPa s and lw ¼ 0:342 mPa s, respectively, all evaluated at the formation pressure and temperature described in Section 2. The permeability kz is 1013 m2 (100 md), the porosity is 0.15, and the relative-permeability curves are defined in Eq. (4). The top boundary is impermeable, and the initial in-situ pressure at the bottom boundary is held constant.

The total amount of CO2 in the transition zone, V, is fixed as 0.0128. Initially, the CO2 saturation in the top 0.8 m of the domain is Sn ¼ 1  Swi ¼ 0:8, below which Sn ¼ 0 (such that V ¼ 0:0128). This corresponds to the sharp-interface condition, which ignores the capillary-pressure effects. Since capillary pressure is modeled here, the CO2 column height is expected to expand from the initial condition. The major driving forces for fluid migration are gravity and capillary forces, which are the major forces during the postinjection period of CO2 storage. The simulation stops when gravity-capillary equilibrium is reached. The CO2 saturation distributions under gravity-capillary equilibrium for the two entry-pressure representations are shown in Fig. 10. The analytical and numerical solutions agree perfectly. While the CO2 column height in the ‘BC-type’ case is 2.96 m, the column height in the ‘Snt ¼ 0:005’ case is 10.91 m. We define ‘plume’ as the main body of the capillary transition zone, which is the shaded area in Fig. 10. The thickness of the plume is denoted by hplume . For the ‘BC-type’ case, the plume fills the entire transition zone, i.e., hc ¼ hplume . However, for the ‘Snt ¼ 0:005’ case, the transition zone is made of a plume and an additional CO2 column below the plume. This additional column is composed of very small amount of CO2, because Snt is only 0.005. For simplicity, hereafter we refer to this low-CO2-saturation column below the CO2 plume as a ‘saturation tail’. The length of the tail is denoted by htail . Hence, for the ‘Snt ¼ 0:005’ case, hc ¼ hplume þ htail . In this case, the saturation tail is much longer than the thickness of the plume. The saturation tail occurs only when the capillary-pressure curve employs a VG-type representation. When gravity-capillary equilibrium is reached in the tail:

h

i VG PVG c ðSw ¼ 1  Snt Þ  P c ðSw ¼ 1Þ ¼ ðqw  qn Þ g htail ;

the length of the tail can be evaluated as: VG htail ¼ ½ PVG c ðSw ¼ 1  Snt Þ  P c ðSw ¼ 1Þ  = ½ ðqw  qn Þ g ;

0

0 P

P

pw ( z )

pn ( z )

pw ( z )

pn ( z )

hc hc z

z

Pc ( S w = 1) = Pc ,e

Pc ( S w = 1) = 0

(a)

z (m)

0

0.1

Sn 0.2

Fig. 9. CO2 and water pressure distribution under gravity-capillary equilibrium: (a) VG-type representation; (b) BC-type representation. z is the depth measured positive downward from the top of the formation. hc is the column height. The shaded area indicates the capillary transition zone.

0.3

0.4

0

0

0

2

2

4 6

Pc ( S w = 1) = Pc , e Analytical

8

Simulation

10

(b)

ð9Þ

where the superscripts VG denotes VG-type representation. In this example, PVG c ðSw ¼ 1Þ ¼ 0. As long as the CO2 plume has enough CO2 to supply the saturation tail, the length of the tail is not influenced by the thickness of the plume. Because the tail is composed of very small amount of CO2, it can form even if the plume is thin. The tail is longer when the density contrast between the resident water and the supercritical CO2 is low, which is usually the case in deep storage formations. The influence of the tail becomes increasingly evident as the CO2 plume migrates under the caprock and becomes stretched and thinner over time. We reiterate that the saturation tail is the consequence of employing a VG-type capillary-pressure curve, and it is not an artifact of numerical diffusion – it is corroborated by analytical calculations and fine-grid simulations. Fig. 11 shows the growth of the CO2 column height over time. In the ‘BC-type’ case, the column height tends to stabilize quickly. In

12

z (m)

492

4 6

Sn 0.2

0.3

0.4

≈ 0.005

Pc ( S w = 1 − S nt ) Analytical

8 10 12

(a)

0.1

Simulation

Pc ( S w = 1) = 0

(b)

Fig. 10. CO2 saturation distribution within the transition zone under gravitycapillary equilibrium: (a) ‘BC-type’ case; (b) ‘Snt ¼ 0:005’ case. The shaded area is the CO2 plume. In (b), a ‘saturation tail’ is observed below the plume.

493

B. Li et al. / Advances in Water Resources 62 (2013) 488–498

12 10

hc (m)

8

Total hc ('entryhc (BC-type) plateau' case) hc of major CO2 body hplume (Snt = 0.005) ('Snt = 0.005' case) Total = 0.005' hc (Shc 0.005) nt =('Snt case)

6 4 2 0 0.01

0.1

1

10 100 Time (years)

1000

10000

Fig. 11. Development of the capillary transition zone for different entry-pressure representations (dissolution is not modeled). hc is the column height of the transition zone. hplume is the thickness of the plume. For the ‘BC-type’ case, hc ¼ hplume . For the ‘Snt ¼ 0:005’ case, the difference between hc and hplume is the length of the saturation tail.

the ‘Snt ¼ 0:005’ case, the thickness of the CO2 plume also stabilizes quickly, but the saturation tail keeps growing. The growth is fast at first, and then it slows down and stabilizes after about 1000 years, when gravity-capillary equilibrium is finally reached in the saturation tail. Although the column heights for the two cases are dramatically different, the mass of CO2 in the saturation tail is so small that it almost has no impact in simulations that do not account for dissolution, and it does not affect the extent and the migration speed of the CO2 plume. Fig. 7 illustrates this point.

4.1.3. Important factors The saturation tail is part of the capillary transition zone. The characteristic time scale for the development of the transition zone can be written as:

t ¼

kz ðDqgÞ2

7

Total hc ('entry-plateau' hc (BC-type) case)

6

hplume (Snt CO2 = 0.005) hc of major body ('Snt = 0.005' case)

3

5

hc (Shc 0.005) Total = 0.005' nt =('Snt case)

2.5

4

1 1 þ kw kn

   dPc ; dSn

ð10Þ

3.5

3 2

CBC-type ('entryplateau' case) CS('Snt = nt = 0.005 0.005' case)

2 1.5 1

1 0 0.01



/

where the superscript  denotes characteristic values, / is the porosity, kz is the permeability in the vertical direction, Dq is the density difference between water and CO2, and kw and kn are the mobilities of water and CO2, respectively. Detailed derivation of this expression can be found in [31]. As indicated by Eq. (10), the lower the density contrast between the resident water and the CO2, the longer it will take for the transition zone, and hence the saturation tail, to develop. On the other hand, the tail will grow longer eventually (see Eq. (9)). In addition, larger vertical permeability and larger CO2 mobility facilitate the development of the saturation tail. An important element in the CO2 relative-permeability curve is the

C ( × 10 –3)

hc (m)

4.1.2. CO2 dissolution and molecular diffusion However, for long-term CO2 sequestration simulations, where CO2 dissolution is modeled, the CO2 delivered to the long saturation tail quickly dissolves in the water. The tendency toward gravity-capillary equilibrium keeps delivering CO2 to sustain the tail, until the water surrounding the tail is fully saturated with the dissolved CO2. The longer the saturation tail tends to be, the larger area of water will it contact. Consequently, more CO2 will be delivered to the water and dissolve. As mentioned in Section 4.1.1, in deep storage formations, the tail can even be much longer than the thickness of the plume. Fig. 12 illustrates the evolution of the CO2 column subject to dissolution. The CO2 solubility (expressed as a mole fraction, xCO2 ) is 0.0254. Molecular diffusion is not modeled here. All other properties, model grid, and boundary conditions are unchanged. The initial condition is the same sharpinterface condition as before, and the water in the initial twophase region is fully saturated with dissolved CO2. Initially there is no dissolved CO2 in the single-phase region. The amount of dis-

solved CO2 is characterized by the dimensionless variable R C ¼ H1 xCO2 ðzÞ dz, where H ¼ 50 m. As indicated in Fig. 12(a), the growth of the CO2 column height for the ‘Snt ¼ 0:005’ case is not as quick as that in Fig. 11. This is because the saturation tail dissolves as it grows. However, CO2 continues to be transferred from the CO2 plume to the underlying water through the tail, and the dissolution causes the plume to shrink. Eventually, much more CO2 is dissolved in the ‘Snt ¼ 0:005’ case compared with the ‘BCtype’ case, as observed in Fig. 12(b). Over long time periods, molecular diffusion is an important mass-transfer process in CO2 solubility trapping, which delivers CO2 from the two-phase region to the underlying water. If CO2 diffusion is strong, the impact of the saturation tail may be obscured to some extent. The results of simulations modeling both dissolution and diffusion are shown in Fig. 13. A typical CO2 diffusivity in water in a porous medium [24], 109 m2/s, is used in the simulations. Compared with the ‘BC-type’ case, the presence of the saturation tail for the ‘Snt ¼ 0:005’ and ‘Snt ¼ 0:01’ cases is still prominent, and it causes the CO2 plume to shrink more quickly (see Fig. 13(a)) due to faster CO2 dissolution (see Fig. 13(b)). The larger the Snt value, the faster is the dissolution. Therefore, while molecular diffusion tends to dampen the importance of the saturation tail, which leads to accelerated dissolution, the phenomenon persists nevertheless. So far, all of the simulations covered in this section have been in 1D, where convective mixing of CO2-rich and CO2–free water – a multi-dimensional process – is suppressed. Discussion of convective mixing is covered in Section 4.2.

0.1

1

10 100 Time (years)

(a)

1000

10000

0.5 0.01

0.1

1

10 100 Time (years)

1000

10000

(b)

Fig. 12. Development of the capillary transition zone for different entry-pressure representations (dissolution is modeled): (a) CO2 column height (hc ); (b) amount of dissolved CO2 (C).

494

B. Li et al. / Advances in Water Resources 62 (2013) 488–498

Total hc ('BC-type' case) hc (BC-type)

hc of major ('Snt = hplume (SCO2 0.01) nt = body 0.01' case) Total hcnt('Snt = 0.01' case) = 0.01) hc (S

4 3 2

4 3.5 C ( × 10 3 )

5

hc (m)

4.5

hc of major ('Snt = hplume (SCO2 0.005) nt = body 0.005' case) Total hcnt('Snt = 0.005' case) hc (S = 0.005)

3 2.5 2

CBC-type ('BC-type' case)

1.5

CS('Snt = 0.005' case) nt = 0.005

1

1

CS('Snt = 0.01' case) nt = 0.01

0.5

0

0 0

200

400 600 Time (years)

800

1000

(a)

0

200

400 600 Time (years)

800

1000

(b)

Fig. 13. Development of the capillary transition zone for different entry-pressure representations (dissolution and diffusion are modeled): (a) CO2 column height (hc ); (b) amount of dissolved CO2 (C). Complete dissolution of the CO2 occurs at 1000 years for the ‘Snt ¼ 0:005’ case, and 500 years for the ‘Snt ¼ 0:01’ case.

1

0.8

H2O CO2

0.6 0.4

Snc ≈ 0

0.2

Relative permeability

Relative permeability

1

0.8

H2O

0.6 0.4 0.2

CO2

Snc ≈ 0

0

0 0

0.2 0.4 0.6 0.8 Sw

(a)

1

0

0.2 0.4 0.6 0.8 Sw

1

(b)

Fig. 14. CO2–water drainage relative-permeability curves measured in the laboratory: (a) Bennion and Bachu [32]; (b) Perrin and Benson [33].

critical CO2 saturation (Snc , see Eq. (4)). It is defined as the minimum CO2 saturation that allows the CO2 to flow in the porous medium for a drainage process. The CO2 becomes mobile only if its saturation is above Snc . If Snc is larger than Snt (width of the entry-slope region, see Fig. 1(a)), the saturation tail will not form due to zero mobility of the CO2. Laboratory measurements of CO2–water drainage relative-permeability curves suggest that Snc can be close to zero. Fig. 14 shows two such examples [32,33]. It is difficult in practice to measure Snc accurately through a typical core-flood experiment. In addition, it remains unclear whether the Snc measured under a certain flow rate in the laboratory is representative under different flow rate conditions in the reservoir. More importantly, the CO2 relative permeability at Sw ¼ 1  Snt ¼ 0:995 is only Oð107 Þ (calculated from Eq. (4)) in the simulations described above. However, even with such very low CO2 mobility, the saturation tail is still able to grow and to influence the simulation results in the long term. Therefore, when the VG-type representation is used, the simulation result is also very sensitive to the value of Snc . 4.2. Density-driven convection Density-driven convection is an important process in CO2 solubility trapping. As the CO2 that accumulates below the caprock dissolves, the water becomes denser than the underlying CO2–free water. Convective mixing of the heavier CO2-rich and the lighter CO2–free water is expected to take place, and the mixing behavior can be influenced by the presence or absence of a saturation tail. The onset of convection and the critical wavelength have been studied by several authors using stability analysis [3,18,19], and high-resolution computations have been performed to model the

evolution of the fingers and to estimate the mass-transfer rates [18,20–22]. However, most of these and related papers ignored capillary pressure and assumed a sharp interface separating the two-phase (supercritical CO2 and CO2-rich water) and single-phase (CO2-rich and CO2–free water) regions, as sketched in Fig. 15(a). The sharp-interface condition assumes that the water in the twophase region is immobile. Convection is only allowed in the single-phase region, where diffusion of the CO2 into the water leads to an unstable density stratification. This treatment allows one to only analyze and simulate miscible convection in the single-phase region below the supercritical CO2 layer [3,18,20,21]. Considering a capillary transition zone in modeling the convective mixing differs from using a sharp-interface assumption. Fig. 15(b) sketches the mass-transfer processes when a BC-type capillary-pressure curve is used. In this model, water is mobile in the two-phase region, leading to convection in this region. Elenius et al. [23] studied the cross-flow between the two-phase and single-phase regions driven by convection, and its impact on the onset time of convection and the instability of the CO2–rich water layer. They concluded that the cross-flow reduces the onset time, increases the mass-transfer rate, and causes greater CO2 dissolution compared with the sharp-interface assumption. The domains of their study were limited in the single-phase region, and the cross-flow was modeled as a boundary condition. Capillary pressure was not directly modeled in their study. If a VG-type capillary-pressure curve is used, the tendency to form a saturation tail will transfer supercritical CO2 from the CO2 plume to the underlying water-saturated region, as illustrated in Fig. 15(c). The CO2 in the tail will quickly dissolve, and the tail will be resaturated by newly delivered CO2. Because convection mixes the CO2-rich and CO2–free water, the water surrounding the saturation tail tends to be undersaturated. As the saturation tail dissolves, it is continuously resupplied from the overlying plume, resulting in a greater cross-flow between the two-phase and single-phase regions, compared with the BC-type capillary-pressure curve. Thus, the dissolution rate of the CO2 plume is accelerated. Fine-grid, 2D, two-phase numerical simulations were performed to investigate the effect of the entry-pressure representation. The porous medium is 10 m  50 m in the x and z directions (i.e., L ¼ 10 m, and H ¼ 50 m), and each gridblock is 0.1 m  0.1 m. All the properties (including the CO2 solubility and diffusivity) and the top and bottom boundary conditions are the same as the simulations described in Section 4.1.2. The left and right boundaries are periodic, which amounts to ‘rolling up’ the domain about a vertical axis and connecting the two boundaries together. The initial saturation distribution is the same as in Section 4.1.2,

495

B. Li et al. / Advances in Water Resources 62 (2013) 488–498

Sn = 1 – Swi

Sn = 1 – Swi no gradient

Sn gradient

dissolved CO2

dissolved CO2

(a)

(b) diffusion

CO2 plume ‘Saturation tail’. The tail forms and dissolves, facilitating the convective mixing.

Sn = 1 – Swi

convection within single-phase region

Sn gradient

convection between two-phase and single-phase regions mass transfer of supercritical CO2 into the ‘saturation tail’

dissolved CO2

two-phase region, Sn= 1 –Swi two-phase region, Sn< 1 –Swi interface between two-phase and single-phase regions maximum depth of the ‘saturation tail’ interface between CO2-rich and CO2-free water

(c)

Fig. 15. Sketch of the density-driven convection process: (a) using the sharp-interface assumption, ignoring the capillary pressure; (b) using a BC-type capillary-pressure curve; (c) using a VG-type capillary-pressure curve.

RR 1 the dimensionless variable C ¼ LH xCO2 ðx; zÞ dxdz. The C’s for all the cases are shown in the figure. The VG-type model clearly yields more CO2 dissolution than the BC-type model. This is more evident for larger Snt values. Fig. 16 is produced using a first-order (in space and time discretizations), two-phase simulator, and it provides a rough qualitative comparison. While not as accurate as high-order simulations, this is adequate to demonstrate the importance of the entry-pressure representation on convective mixing. Performing high-order multiphase simulations is exceptionally challenging, and to our knowledge, all of the high-order simulations in the literature that model the convective mixing are for miscible singlephase flow (e.g., [18,22]).

i.e., the CO2 saturation in the top 0.8 m of the domain is Sn ¼ 1  Swi ¼ 0:8, below which Sn ¼ 0. Initially, the water in the two-phase region is fully saturated with dissolved CO2, and there is no dissolved CO2 in the single-phase region. Four cases employed with different capillary-pressure models are studied. The ‘BC-type’ case uses Eq. (5), and the remaining three VG-type cases use Eq. (6) with Snt values as 0.005, 0.01, and 0.02, respectively. All four capillary-pressure curves are identical except for their entry-pressure representations. The simulation results are compared in Fig. 16, which shows the development of convective fingers after 10 years from the initial condition. The amount of the dissolved CO2 is characterized by

0

BC-type C = 2.19×10-3

0

Snt = 0.005 C = 2.68×10-3

0

Snt = 0.01 C = 3.14×10-3

0

Snt = 0.02 C = 3.67×10-3

x CO 2 0.025 10

10

10

20

20

20

20

0.020

z (m)

10

0.015 0.010 0.005 50

50 0 2 4 6 8 10 x (m)

(a)

50 0 2 4 6 8 10 x (m)

(b)

0

50 0 2 4 6 8 10 x (m)

0 2 4 6 8 10 x (m)

(c)

(d)

Fig. 16. Impact of different entry-pressure representations on density-driven convection, simulated by a first-order two-phase simulator: (a) BC-type representation; (b)Snt ¼ 0:005; (c) Snt ¼ 0:01; (d) Snt ¼ 0:02. (b)–(d) are all VG-type representations. The total amounts of the dissolved CO2 (C) are also shown.

496

B. Li et al. / Advances in Water Resources 62 (2013) 488–498

and artificially constructed pore-networks, which can exclude surface-irregularity effects, have provided evidence that the entryslope region is a reflection of the pore-throat-size distribution [35,36]. However, it is unknown whether the importance of the entry-slope region varies with different rock types. In addition, whether the entry-slope region is still important as the sample size increases should be examined. Larger samples tend to exhibit more small-scale heterogeneities in, for example, porosity, permeability, and the capillary-pressure–saturation relationships [37–39]. Heterogeneity in the capillary-pressure–saturation relationships, also called capillary heterogeneity, creates local capillary barriers that may prevent the saturation tail and hence the accelerated dissolution from occurring. To study the effect of small-scale heterogeneities on CO2 solubility trapping, simulations were performed on a 2D domain, which is 0.2 m by 10 m in the x and z directions, and each gridblock is 0.01 m  0.01 m. Factors affecting the accuracy and efficiency of the simulations considering capillary heterogeneity were discussed in [40,41]. The fluid properties, boundary conditions, and the initial CO2 distribution are the same as those used in Section 4.2. The permeability field is randomly generated, following a uniform distribution between 6  1014 m2 (60 md) and 1013 m2 (100 md) (see Fig. 18(a)). The porosity is 0.15, and the relative-permeability curves are those in Eq. (4), and all of them are homogeneous throughout the domain. The capillary-pressure curves are heterogeneous and are assumed to follow the J-function scaling relation [42]:

5. Discussion In this section, we discuss the entry-slope region of the capillary-pressure curve measured in laboratory, and study the impact of small-scale heterogeneities. Finally, we elaborate some important issues about implementing the capillary-pressure models in a simulator, and offer our suggestions. 5.1. Laboratory measurements and small-scale heterogeneities The capillary-pressure curves measured in the laboratory using mercury porosimetry usually have an entry-slope region, as sketched in Fig. 17(a). There has been a wide discussion in the literature on the physical meaning and importance of the entry-slope region. Consensus, however, has not been reached. Many authors, e.g., Schowalter [8] and Jennings [9], suggested that the true entry pressure can be found on the extension of the plateau region, as illustrated by the dash line in Fig. 17(b). They stated that the entry-slope region is of little importance, and that it is mainly created by mercury conforming to the irregularities on the surface of the rock sample but without actually entering the pores. Katz and Thompson [10] argued that the capillary-pressure curve from Sw ¼ 1 to the major inflection point (see Fig. 17(a)) should be neglected. Instead of ignoring the entire entry-slope region, Nabawy et al. [11] used Pittman’s method [34] and located the entry pressure near the secondary inflection point, which often exists in laboratory-measured capillary-pressure curves (see Fig. 17(a)). Usually, the corrected curve still has an entry-slope region (see the solid line in Fig. 17(b)). The entry-slope region may be the result of a certain distribution of the pore-throat sizes. Numerical simulations on natural

Pc

rffiffiffiffi / Pc ¼ 1:633  10 J; k 7

Major inflection point

Extension of the plateau

102

101

101

10-1

10-1 0

0.5 Sw

1

0

0.5 Sw

(a)

1

ð12Þ

(b)

where Swi ¼ 0:2. For most part of the curve (Swi 6 Sw 6 1  Snt ), it is exactly equal to Eq. (11). At 1  Snt < Sw 6 1, however, it uses a steep slope to connect the plateau to the endpoint JðSw ¼ 1Þ ¼ 0:5. The Snt values for the four VG-type cases are 0.005, 0.01, 0.02, and

Fig. 17. Sketch of the capillary-pressure curve measured from mercury porosimetry: (a) the original curve; (b) the corrected curve using Pittman’s method [34]. The dash line is the extension of the plateau region.

k (×10-14 m2) 10

0

1 5

BC-type BC-type Snt == 0.005 0.005 Snt Snt == 0.01 0.01 Snt 0.02 Snt == 0.02 Snt Snt == 0.04 0.04 Snt

0.8 4

9

0.6 3

4 8 6 7

J

z (m)

2

0.4 2

0.3 1.5 0.2 1.0 0.1 0.5 0 0.95

1 0.2

8 0

0.1 0.2 x (m)

(a)

6

ð11Þ

8 0:5 > Sw Swi > if Swi 6 Sw 6 1  Snt < 1S wi  J ¼  0:5  > nt Swi w >  0:5 1S þ 0:5 if 1  Snt < Sw 6 1; : 1S 1Swi Snt

Entry-slope region

100

Secondary inflection point

J

100

 0:5 Sw  Swi ; 1  Swi

where J is a BC-type J-function. Swi is 0.2, and the unit of P c and k are bar and m2, respectively. If k is 1013 m2 (100 md), Eq. (11) is equivalent to Eq. (5). If a gridblock has a lower k than its neighbors, its capillary-pressure curve and entry pressure are greater. This serves as a local capillary barrier that prevents the nonwetting CO2 from entering the gridblock, unless the capillary pressure overcomes the entry pressure. Five cases with different entry-pressure representations of J are compared. The ‘BC-type’ case uses Eq. (11), where JðSw ¼ 1Þ ¼ 1. The remaining four cases use the VG-type J-function:

Pc

102

J¼

0.975

1

Sw

0 0.2

0.4

0.6

0.8

1

Sw

(b)

Fig. 18. Input for simulations accounting for capillary heterogeneity: (a) permeability field; (b) J-functions.

C ( × 10–3)

B. Li et al. / Advances in Water Resources 62 (2013) 488–498

16 14 12 10 8 6 4 2 0 BC-type Snt Snt == Snt Snt == 0.005 0.005 0.01 0.01

Snt Snt == 0.02 0.02

Snt Snt = 0.04 0.04

Fig. 19. Amount of the dissolved CO2 after 20 years, subject to capillary heterogeneity.

0.04, respectively. The J’s for all the five cases are plotted in Fig. 18(b). The amounts of the dissolved CO2 for all the cases after 20 years from the initial condition are compared in Fig. 19. The saturation distributions for the ‘BC-type’ and ‘Snt ¼ 0:04’ cases are shown in Fig. 20. While the accelerated dissolution is not significant for the ‘Snt ¼ 0:005’ case, it is still observed for all cases and can be quite important for larger Snt values. Therefore, the entry-pressure representation remains an important consideration in modeling CO2 solubility trapping even in the presence of small-scale capillary heterogeneity. Whether the entry-slope region of a VG-type capillary-pressure curve is real and important in natural reservoir rocks still requires more experimental observations. A new method for measuring capillary-pressure curves developed by Pini et al. [43] recently may be a viable approach. Compared to mercury porosimetry, which provides a single capillary-pressure curve for a small rock sample, this method can quantify the extent of capillary heterogeneity on much larger rock samples. Therefore, this method may improve our understanding of the capillary entry-pressure by taking small-scale heterogeneities into consideration. 5.2. Implementation issues in simulators Entry-pressure representation depends not only on the interpretation of laboratory measurements but how capillary entry pressure is treated numerically in the simulator. Typically, a simulator takes the input of a capillary-pressure curve in either functional form or tabular form. Regardless of how the capillarypressure curve is provided, the simulator can treat the entry pressure in several ways. This is particularly important for BC-type curves, whose entry pressures are nonzero. For example, a simulator may require Pc ¼ 0 in fully water-saturated regions, but allow the Pc to jump to the entry pressure at Sw ¼ 1  e, where e is a small threshold value. In other words, P c is allowed to be discontin-

497

uous. Some simulators however cannot handle discontinuous Pc , and BC-type curves are avoided altogether (e.g., TOUGH2 version 2.0 [44]). Other simulators use work-arounds to deal with this issue. For example, one approach is to use the ‘virtual P c ’ idea described by our previous work, where P c is assigned to the entry pressure even in the water-saturated regions [40,41]. The second way is to connect the plateau of the BC-type curve to Pc ðSw ¼ 1Þ ¼ 0 using a steep slope, which is equivalent to transforming the BC-type curve into VG-type. As we have shown, the amount of dissolution will be very different for these two types. Unfortunately, the treatment of the entry pressure is often not included in the user manual of a simulator. It is therefore important that the user consider the simulator’s treatment before running a CO2 sequestration simulation. 6. Concluding remarks The typical shape of a capillary-pressure curve is either convex (e.g., Brooks–Corey model) or S-shaped (e.g., van Genuchten model). Both models are widely used in simulating long-term CO2 sequestration. The two models represent the entry of the nonwetting phase into the porous medium differently. The Brooks– Corey-type (BC-type) model uses a plateau that ends with a nonzero capillary entry pressure, while the van-Genuchten-type (VG-type) model has an ‘entry-slope’ region, which connects the end-point (usually zero) to the plateau. In this paper, we simulate long-term CO2 sequestration using these two representations of the capillary entry pressure. Although the two entry-pressure representations produce similar results when CO2 dissolution is ignored, the VG-type representation significantly accelerates the rate of CO2 solubility trapping compared with the BC-type representation. For the representative examples we study, the acceleration can be so large as to have complete dissolution of the CO2 plume versus persistence of over 50% of the plume over a 5000-year period. We also observe that the simulation results are very sensitive to the steepness of the VG-type curve around the entry-slope region. The cause of such strong sensitivity to the entry-pressure representation is studied. Compared with the BC-type curve, the entryslope region of the VG-type curve yields a longer CO2 column height due to the interaction of gravity and capillary forces, leading to the development of a so-called ‘saturation tail’. The saturation tail grows longer in deep storage formations, where the density contrast between the resident water and the supercritical CO2 is low. The growth is facilitated if the vertical permeability is large and if the critical CO2 saturation (Snc ) value is small. The saturation tail consists of very small amount of CO2, which can dissolve quickly into the surrounding water. The tendency toward gravity-capillary equilibrium keeps transferring more CO2 from the CO2 plume to the underlying water through this saturation tail. Density-driven convection then carries the CO2–rich water away

Fig. 20. CO2 saturation distributions after 20 years, subject to capillary heterogeneity: (a) ‘BC-type’ case; (b) ‘Snt ¼ 0:04’ case.

498

B. Li et al. / Advances in Water Resources 62 (2013) 488–498

and brings CO2–free water close by, allowing more CO2 to dissolve. This vertical mass-transfer mechanism introduced by the VG-type representation leads to faster CO2 dissolution compared with the BC-type representation. We also show that this mechanism remains prominent even in the presence of small-scale capillary heterogeneity. Capillary-pressure curves measured in the laboratory tend to have an entry-slope region. Consensus has not been reached regarding the importance of such entry-slope region, neither is it clear whether the importance varies with rock types. Careful measurements, preferably on relatively large cores, are needed to study this issue. Note that the intention of this work is not to argue whether such an entry-slope region is real in natural systems. The intention is, however, to point out that the VG-type capillary-pressure model, which is widely used in making predictions of long-term CO2 sequestration, will lead to higher levels of CO2 dissolution compared with the BC-type model. This has not been reported in the literature. We should also note that the entry-pressure representation depends not only on the capillary-pressure model, but also on the simulator design. Details about implementing capillary-pressure models in simulators can be found in Section 5.2. We recommend that the simulator user be aware of how the capillary entry pressure is handled in the simulator, in order to understand how this affects the rate of dissolution. Acknowledgments The authors thank Global Climate and Energy Project (GCEP) for supporting this research, and acknowledge Dr. Gary Li, Dr. Yaqing Fan, and Ruslan Iskhakov for helpful discussions. References [1] Benson SM, Cook P, Anderson J, Bachu S, Nimir HB, Basu B. Underground geological storage. In: IPCC Special Report on Carbon Dioxide Capture and Storage, Chapter 5. Intergovernmental Panel on Climate Change. Cambridge, UK: Cambridge University Press; 2005. [2] Lindeberg E, Wessel-Berg D. Vertical convection in an aquifer column under a gas cap of CO2. Energy Convers Manage 1997;38:S229–34. http://dx.doi.org/ 10.1016/S0196-8904(96)00274-9. [3] Ennis-King JP, Paterson L. Role of convective mixing in the long-term storage of carbon dioxide in deep saline formations. SPE J 2005;10(3):349–56. http:// dx.doi.org/10.2118/84344-PA. [4] Kopp A, Class H, Helmig R. Investigations on CO2 storage capacity in saline aquifers. Part 1: Dimensional analysis of flow processes and reservoir characteristics. Int J Greenhouse Gas Control 2009;3(3):263–76. http:// dx.doi.org/10.1016/j.ijggc.2008.10.002. [5] Bear J. Dynamics of fluids in porous media. Courier Dover Publications; 1972. [6] van Genuchten MT. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci Soc Am J 1980;44(5):892–8. [7] Brooks RH, Corey AT. Hydraulic properties of porous media. Fort Collins, CO, USA: Colorado State University; 1964. [8] Schowalter TT. Mechanics of secondary hydrocarbon migration and entrapment. AAPG Bull 1979;63(5):723–60. [9] Jennings JB. Capillary pressure techniques: application to exploration and development geology. AAPG Bull 1987;71(10):1196–209. [10] Katz AJ, Thompson AH. Quantitative prediction of permeability in porous rock. Phys Rev B 1986;34(11):8179–81. [11] Nabawy BS, Geraud Y, Rochette P, Bur N. Pore-throat characterization in highly porous and permeable sandstones. AAPG Bull 2009;93(6):719–39. http:// dx.doi.org/10.1306/03160908131. [12] Ennis-King J, Paterson L. Engineering aspects of geological sequestration of carbon dioxide. Paper SPE 77809 presented at the Asia Pacific oil and gas conference and exhibition, Melbourne, Australia; 8–10 October 2002. http:// dx.doi.org/10.2118/77809-MS. [13] Doughty C, Pruess K. Modeling supercritical carbon dioxide injection in heterogeneous porous media. Vadose Zone J 2004;3(3):837–47. [14] Pruess K, Nordbotten J. Numerical simulation studies of the long-term evolution of a CO2 plume in a saline aquifer with a sloping caprock. Transp Porous Med 2011;90(1):135–51. http://dx.doi.org/10.1007/s11242-011-9729-6. [15] Class H, Ebigbo A, Helmig R, Dahle HK, Nordbotten JM, Celia MA, et al. A benchmark study on problems related to CO2 storage in geologic formations. Comput Geosci 2009;13(4):409–34. http://dx.doi.org/10.1007/s10596-0099146-x.

[16] Kumar A, Ozah R, Noh M, Pope GA, Bryant S, Sepehrnoori K, et al. Reservoir simulation of CO2 storage in deep saline aquifers. SPE J 2005;10(3):336–48. http://dx.doi.org/10.2118/89343-PA. [17] Pruess K, Oldenburg C, Moridis G. TOUGH2 User’s Guide, Version 2.0. Technical report LBNL-43134, Lawrence Berkeley National Laboratory, Berkeley, California; 1999. [18] Riaz A, Hesse M, Tchelepi HA, Orr FM. Onset of convection in a gravitationally unstable diffusive boundary layer in porous media. J Fluid Mech 2006;548(1):87–111. http://dx.doi.org/10.1017/S0022112005007494. [19] Mahidjiba A, Robillard L, Vasseur P. Onset of penetrative convection of cold water in a porous layer under mixed boundary conditions. Int J Heat Mass Transfer 2006;49:2820–8. http://dx.doi.org/10.1016/j.ijheatmasstransfer. 2006.02.019. [20] Pruess K, Zhang K. Numerical modeling studies of the dissolution–diffusion– convection process during CO2 storage in saline aquifers. Technical report LBNL-1243E, Lawrence Berkeley National Laboratory, California; 2008. [21] Pau GSH, Bell JB, Pruess K, Almgren AS, Lijewski MJ, Zhang K. High-resolution simulation and characterization of density-driven flow in CO2 storage in saline aquifers. Adv Water Resour 2010;33(4):443–55. http://dx.doi.org/10.1016/ j.advwatres.2010.01.009. [22] Neufeld JA, Hesse MA, Riaz A, Hallworth MA, Tchelepi HA, Huppert HE. Convective dissolution of carbon dioxide in saline aquifers. Geophys Res Lett 2010;37(22):L22404. http://dx.doi.org/10.1029/2010GL044728. [23] Elenius MT, Nordbotten JM, Kalisch H. Effects of a capillary transition zone on the stability of a diffusive boundary layer. IMA J Appl Math 2012;77(6):771–87. http://dx.doi.org/10.1093/imamat/hxs054. [24] Dahle HK, Eigestad GT, Nordbotten JM, Pruess K. A model-oriented benchmark problem for CO2 storage. In: Workshop on modeling and risk of assessment of geological storage of CO2; 2009. [25] Cao H. Development of techniques for general purpose simulators. Ph.D. thesis, Stanford University, USA; 2002. [26] Jiang Y. Techniques for modeling complex reservoirs and advanced wells. Ph.D. thesis, Stanford University, USA; 2007. [27] Kumar A. A simulation study of carbon sequestration in deep saline aquifers. M.S. thesis, University of Texas at Austin, USA; 2004. [28] Aziz K, Settari A. Petroleum reservoir simulation. London, England: Applied Science Publishers; 1979. [29] Silin D, Patzek T, Benson SM. A model of buoyancy-driven two-phase countercurrent fluid flow. Transp Porous Med 2009;76(3):449–69. http:// dx.doi.org/10.1007/s11242-008-9257-1. [30] Fayers FJ, Sheldon JW. The effect of capillary pressure and gravity on twophase fluid flow in a porous medium. Trans AIME 1959;216:147–55. [31] Nordbotten JM, Dahle HK. Impact of the capillary fringe in vertically integrated models for CO2 storage. Water Resour Res 2011;47(2):W02537. http:// dx.doi.org/10.1029/2009WR008958. [32] Bennion DB, Bachu S. Dependence on temperature, pressure, and salinity of the IFT and relative permeability displacement characteristics of CO2 injected in deep saline aquifers. Paper SPE 102138 presented at the SPE annual technical conference and exhibition, San Antonio, TX, USA; 24–27 September 2006. http://dx.doi.org/10.2118/102138-MS. [33] Perrin JC, Benson SM. Relative permeability explorer, Benson Lab, Stanford University; 2012. Available from: . [34] Pittman ED. Relationship of porosity and permeability to various parameters derived from mercury injection-capillary pressure curves for sandstone. AAPG Bull 1992;76(2):191–8. [35] Silin D, Tomutsa L, Benson SM, Patzek TW. Microtomography and pore-scale modeling of two-phase fluid distribution. Transp Porous Med 2011;86(2):495–515. http://dx.doi.org/10.1007/s11242-010-9636-2. [36] Joekar-Niasar V, Hassanizadeh S, Leijnse A. Insights into the relationships among capillary pressure, saturation, interfacial area and relative permeability using pore-network modeling. Transp Porous Med 2008;74(2):201–19. http:// dx.doi.org/10.1007/s11242-007-9191-7. [37] Perrin JC, Benson SM. An experimental study on the influence of sub-core scale heterogeneities on CO2 distribution in reservoir rocks. Transp Porous Med 2010;82(1):93–109. http://dx.doi.org/10.1007/s11242-009-9426-x. [38] Shi JQ, Xue Z, Durucana S. Supercritical CO2 core flooding and imbibition in tako sandstone: influence of sub-core scale heterogeneity. Int J Greenhouse Gas Control 2011;5(1):75–87. http://dx.doi.org/10.1016/j.ijggc.2010.07.003. [39] Krevor SCM, Pini R, Li B, Benson SM. Capillary heterogeneity trapping of CO2 in a sandstone rock at reservoir conditions. Geophys Res Lett 2011;38(15):L15401. http://dx.doi.org/10.1029/2011GL048239. [40] Li B. Including fine-scale capillary heterogeneity in modeling the flow of CO2 and brine in reservoir cores. M.S. thesis, Stanford University, USA; 2011. [41] Li B, Benson SM, Tchelepi HA. Modeling fine-scale capillary heterogeneity in multiphase flow of CO2 and brine in sedimentary rocks. In: Proceedings of the XIX international conference on water resources University of Illinois at Urbana-Champaign, IL, USA; June 2012. p. 17–22. [42] Leverett MC. Capillary behavior in porous solids. Trans AIME 1941;142:152–69. [43] Pini R, Krevor SCM, Benson SM. Capillary pressure and heterogeneity for the CO2/water system in sandstone rocks at reservoir conditions. Adv Water Resour 2012;38:48–59. http://dx.doi.org/10.1016/j.advwatres.2011.12.007. [44] Personal communication with TOUGH2 developers; 2012.