University of Southern California Dissertations and Theses

The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media

(USC Thesis Other)

The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content UNIVERSITY OF SOUTHERN CALIFORNIA

The Role of Counter-Current Flow in the Modeling and
Simulation of Multi-Phase Flow in Porous Media

Mohammad Javaheri

PHD DISSERTATION

MORK FAMILY DEPARTMENT OF CHEMICAL ENGINEERING AND
MATERIALS SCIENCE
LOS ANGELES, CALIFORNIA
December 2013

© MOHAMMAD JAVAHERI 2013
ii

Abstract

Carbon Capture and Storage (CCS) is the capture of carbon dioxide (CO
2
) from large
source points, such as power plants, and storing it in geological formation. It is a potential
means of reducing CO
2
concentration in the atmosphere to mitigate global warming and
ocean acidification. Geological formations are considered to be the most promising
storage sites for carbon dioxide. Deep saline aquifers have the most storage capacity
compared to other geological sites.
After CO
2
is injected into an aquifer, a fraction of it, that can be large, is
immobilized in the form of residual phase. Residual entrapment is an important
component in successful storage of CO
2
in saline aquifers. Trapped CO
2
is eventually
dissolved into the brine. Estimation of the fraction of the injected CO
2
that is trapped
depends on CO
2
saturation profiles and saturation history. In the simulation of CO
2

injection into saline aquifers, conventional relative permeabilities functions are
commonly used and it is assumed that viscous coupling is negligible. We show that the
vertical migration of CO
2
in saline aquifers is often dominated by counter-current flow.
Experimental and simulation studies of two- and three-phase flow in porous media show
that the counter-current relative permeability is less than co-current relative permeability,
that are commonly used (current industry standard) in the simulation of multiphase flow.
This study focuses on including a velocity- -dependent relative permeability model
in the simulation of multiphase flow in porous media to account for the flow dynamics
(co- vs. counter-current flow) on relative permeability and fluid mobility. We use both
co-current and counter-current relative permeability in the simulation of CO
2
injection
iii

from a well into an aquifer and show that the plume saturation profile is influenced by the
dynamic relative permeability changes (transitions between co-current and counter-
current flow). The fraction of the injected CO
2
that is trapped and the time-scale of
vertical migration of CO
2
both increase when the transitions from co- to counter-current
relative permeability are included. We conclude that using one set of relative
permeabilities is not sufficient in the simulation of CO
2
injection into saline aquifers. We
also demonstrate that a velocity-dependent mobility does not introduce numerical
instability and that the simulation time does not change significantly.
We have applied the velocity-dependent relative permeability in the simulation of
gas (CO
2
) injection into an oil reservoir in the context of enhanced oil recovery with three
mobile phases. Cumulative oil and gas production increases just slightly in the new
model, but water production decreases quite significantly. That is due to a better sweep
efficiency when counter-current relative permeability is included in the calculation of
phases’ mobilities.
Counter-current relative permeability can be measured directly using a controlled
experimental design to eliminate the effect of boundaries on saturation distribution. It can
also be estimated by pore-scale simulations, e.g. Lattice Boltzman simulations. Both
methods have been used in the literature and both show a reduction in the relative
permeability in counter-current settings. We have tested an indirect procedure to
conclude a reduction of relative permeability in counter-current flows.
The numerical simulations are done using IMPES and fully-implicit methods. No
numerical instabilities were observed in any of the two methods by including counter-
current relative permeability either in a discrete scheme or by using a velocity-dependent
iv

scheme. We conclude that the new model can predict CO
2
saturation profile more
accurately than the conventional standard model.

v

Dedicated to My Parents,
Ahmad and Maryam

vi

Acknowledgements
I would like to thank my advisor, Kristian Jessen, whose ideas contributed to make my
Ph.D. productive. I am very thankful to him for the time we spent to discuss various
aspects of my research and for his passion and enthusiasm throughout the course of this
dissertation.
I also thank the committee members of my dissertation, Dr. Iraj Ershaghi and Dr.
Julian Domaradzki for their valuable comments. I especially thank Dr. Ershaghi as the
head of the petroleum engineering program in our department who always supported me
during my doctoral studies at the University of Southern California.
The years that I spent at the University of Southern California could not be
enjoyable without accompany of good friends. I would like to thank all of them,
especially my friends at our department that we spent good times together. The friendship
that we started as graduate students will definitely continue in future. I would also like to
thank the support of the administrative assistants in our department who were always
very helpful.
I thank the Graduate School at the University of Southern California for granting
the Provost Fellowship in the first two years of my Ph.D. studies. I also thank the Global
Climate and Energy Project (GCEP) at Stanford University for the financial support of
my project. I acknowledge SUPRI-B research group at Stanford University for providing
the General Purpose Research Simulator (GPRS) code.
Lastly, I would like to thank my family, especially my parents, for their constant
love and encouragements. They have always supported me in any way they could, even if
they are on the other side of the globe. Thank you.
vii

Table of Contents
Abstract ............................................................................................................................... ii
Acknowledgements ............................................................................................................ vi
List of Tables .......................................................................................................................x
List of Figures .................................................................................................................... xi
CHAPTER ONE: INTRODUCTION ..................................................................................1
1.1 Background ................................................................................................................1
1.2 Motivation ..................................................................................................................4
1.3 Research Objectives ...................................................................................................4
1.4 Organization of the Manuscript .................................................................................6
CHAPTER TWO: LITERATURE REVIEW ......................................................................7
2.1 Background ................................................................................................................7
2.2 Viscous Coupling in Multiphase Flow in Porous Media ...........................................7
2.3 Counter-Current Relative Permeability .....................................................................8
2.4 Relative Permeability Anisotropy ............................................................................12
2.5 Three-phase Relative Permeability Models .............................................................13
CHAPTER THREE: NON-WETTING PHASE PLUME MIGRATION IN
VERTICAL TWO-PHASE COUNTER-CURRENT FLOW ..................................14
3.1 Introduction ..............................................................................................................14
3.2 Experimental Approach ...........................................................................................18
3.2.1 Apparatus and characterization .......................................................................18
3.2.2 Segregation Experiments .................................................................................23
3.3 Simulation of Segregation Experiments ..................................................................26
3.4 Results: Prediction and Optimization ......................................................................28
3.4.1 Simulation with co-current model input ..........................................................28
3.4.2 Adjustment of model parameters .....................................................................31
3.5 Discussion and Conclusion ......................................................................................38
CHAPTER FOUR: DYNAMIC MOBILITY IN THE BRINE/CO
2
SYSTEM (IMPES
FORMULATION) ....................................................................................................46
4.1 Introduction ..............................................................................................................46
4.2 Modeling and Simulation.........................................................................................47
4.2.1 Hysteresis in Relative Permeability .................................................................48
4.2.2 Hysteresis in Capillary Pressure ......................................................................51
4.2.3 Flow Transitions ..............................................................................................52
4.3 Results ......................................................................................................................52
4.3.1 Homogeneous Aquifer .....................................................................................53
4.3.2 Heterogeneous Aquifer (Short Correlation) ....................................................60
4.3.3 Heterogeneous Aquifer (Medium Correlation) ...............................................62
4.3.4 Heterogeneous Aquifer (Long Correlation) ....................................................65
4.4 Discussion of the Results .........................................................................................68
4.5 Concluding Remarks ................................................................................................72
viii

CHAPTER FIVE: RESIDUAL ENTRAPMENT DURING SIMULTANEOUS
INJECTION OF CO
2
AND BRINE .........................................................................73
5.1 Introduction ..............................................................................................................73
5.1.1 Simultaneous injection of CO
2
and brine ........................................................74
5.2 Statement of the Problem .........................................................................................74
5.3 Modeling Approach .................................................................................................76
5.3.1 Simulation model .............................................................................................76
5.3.2 Scaling groups .................................................................................................79
5.3.3 Capillary pressure and relative permeability ...................................................81
5.4 Impact of Dimensionless Numbers on Residual Entrapment of CO
2
......................82
5.4.1 CO
2
saturation and residual entrapment (co-current model) ...........................82
5.4.2 Gravity number and residual entrapment of CO
2
............................................92
5.4.3 Segregation distance of brine and CO
2
............................................................94
5.4.4 Rate of brine injection .....................................................................................94
5.5 Co-Current vs. Mixed Model ...................................................................................95
5.5.1 Segregation distance of the fluids ....................................................................98
5.5.2 Heterogeneous permeability fields ..................................................................98
5.6 Concluding Remarks ..............................................................................................101
CHAPTER SIX: TRANSITIONS BETWEEN CO-CURRENT AND COUNTER-
CURRENT MOBILITY IN TWO-PHASE FLOWS .............................................102
6.1 Introduction ............................................................................................................102
6.2 A New Mobility Model ..........................................................................................103
6.3 Calculation Examples ............................................................................................106
6.4 Calculation Results ................................................................................................110
6.4.1 Example 1 ......................................................................................................110
6.4.2 Example 2 ......................................................................................................116
6.4.3 Effect of Vertical Permeability ......................................................................121
6.4.4 Interpolation Scheme .....................................................................................122
6.5 Concluding Remarks ..............................................................................................124
CHAPTER SEVEN: TRANSITIONS BETWEEN CO-CURRENT AND COUNTER-
CURRENT MOBILITY IN THREE-PHASE FLOWS ..........................................126
7.1 Introduction ............................................................................................................126
7.2 Three-Phase Relative Permeability Calculation ....................................................126
7.2.1 Stone I model .................................................................................................127
7.2.1.1 Aziz and Settari’s modification of the Stone I model ..........................128
7.2.1.2 Residual oil saturation (S
om
) ................................................................128
7.2.2 Counter-current Relative Permeability in Three-Phase Flows ......................129
7.2.2.1 Velocity-dependent relative permeability in three-phase flows ..........129
7.3 Reservoir Model and Fluid Properties ...................................................................132
7.3.1 Fluid Properties .............................................................................................132
7.4 Calculation Example 1: Simultaneous Injection of Water and Gas (Two Wells) .134
7.5 Example Calculation 2: Simultaneous Injection of Water and Gas (Five-Spot
Pattern) .................................................................................................................141
ix

7.6 Concluding Remarks ..............................................................................................146
CHAPTER EIGHT: RECOMMENDATIONS FOR FUTURE STUDIES .....................147
REFERENCES ................................................................................................................150

x

List of Tables
Table 3-1. Comparison of observed and adjusted input parameters ................................. 37
Table 5-1. Properties of the four aquifer settings (100,000 ppm) ..................................... 78
Table 5-2. Gravity and capillary numbers for the four aquifer settings ............................ 93
Table 6-1. Specifications for the 2D aquifer model ........................................................ 108
Table 6-2. CO
2
properties ............................................................................................... 108
Table 6-3. Brine properties ............................................................................................. 109
Table 6-4. Specifications for the 3D aquifer model ........................................................ 109
Table 6-5. Numerical comparison of the three mobility models for Example 1 ............ 116
Table 6-6. Numerical comparison of the three mobility models for Example 2 ............ 121
Table 7-1. Properties of the gas phase ............................................................................ 133
Table 7-2. Properties of the oil phase ............................................................................. 133
Table 7-3. Properties of the water phase ......................................................................... 134
Table 7-4. Comparison of numerical performance in the co-current and mixed models 140
Table 7-5. Numerical comparison between co-current and mixed models .................... 145

xi

List of Figures
Figure 1.1. Deep underground storage of CO
2
(IPCC, 2005) ............................................. 3
Figure 3.1. Dimensions and design of the Borosilicate glass column for resistivity
measurement ............................................................................................................. 19
Figure 3.2. Steady-state relative permeability functions for drainage and imbibition
processes ................................................................................................................... 22
Figure 3.3. Normalized capillary pressure data extracted from Dawe et al. (1992) ......... 23
Figure 3.4. Schematic diagram of experimental setup for segregation experiment .......... 24
Figure 3.5. Saturation profiles of iC
8
in the column at various times; left: Experiment
A and right experiment B .......................................................................................... 26
Figure 3.6. Saturation of iC
8
in the column observed in the experiments and predicted
based on co-current model parameters (upper panel: Experiment A, lower panel:
Experiment B) ........................................................................................................... 30
Figure 3.7. Comparison of experiment and simulation for experiment A (upper panel)
and experiment B (lower panel) when endpoint non-wetting phase relative
permeability is kept constant; Simulation input was adjusted based on
experiment A only ..................................................................................................... 35
Figure 3.8. Comparison of observations and simulation for experiment A (upper
panel) and experiment B (lower panel) when endpoint non-wetting phase relative
permeability is kept constant; Simulation input was adjusted based on both
experiments ............................................................................................................... 36
Figure 3.9. Measured and adjusted relative permeability for wetting (brine) and non-
wetting (iC
8
) phases .................................................................................................. 38
Figure 3.10. The sensitivity of iC
8
migration time (to the top of the column) to
relative permeability multiplier ................................................................................. 39
Figure 3.11. Viscous coupling during two-phase iC8-brine flow in a cylindrical pore .... 43
Figure 3.12. Viscous coupling during two-phase CO
2
-brine flow in a cylindrical pore .. 43
Figure 4.1. Aquifer dimensions and the location of injector ............................................ 48
Figure 4.2. Logarithm of four permeability fields with different correlations in
heterogeneity a) homogeneous field, b) short cor
r
elation length, c) medium
correlation length, and d) long correlation length ..................................................... 50
xii

Figure 4.3. Co-current and counter-current relative permeabilities of brine and CO
2
..... 51
Figure 4.4. Residual entrapment of gas as a function of time in a homogeneous
aquifer ....................................................................................................................... 55
Figure 4.5. CO
2
saturation after injection for the co-current model (left) and the
mixed model (right) .................................................................................................. 58
Figure 4.6. CO
2
saturation at different times (left) and corresponding regions where
counter-current flow occurs (right) for a heterogeneous aquifer with short
correlation length and k
v
/k
h
=0.3; Regions with counter-current flow are shown in
black. Injection stops at 2400 days ........................................................................... 59
Figure 4.7. Residual entrapment of gas as a function of time in a heterogeneous
aquifer (short correlation length) .............................................................................. 60
Figure 4.8. a) CO
2
saturation from co-current model, b) Trapped saturation from co-
current model, c) CO
2
saturation from mixed model, d) Trapped saturation from
mixed model, all after 15000 days for a heterogeneous aquifer with medium
correlation length and permeability ratio of 0.5 ........................................................ 61
Figure 4.9. CO
2
saturation at different times (left) and corresponding regions where
counter-current flow occurs (right) for a heterogeneous aquifer with a medium
correlation length and a permeability ratio of 0.3; Counter-current flow occurs in
regions with black color on the right panels. Injection stops at 2400 days............... 63
Figure 4.10. Residual entrapment of gas as a function of time in a heterogeneous
aquifer (medium correlation length) ......................................................................... 65
Figure 4.11. Saturation distributions (left column), regions with counter-current flow
in the vertical direction (middle column) and regions with counter-current flow
in the horizontal direction (right column) for a heterogeneous aquifer with long
correlation length and k
v
/k
h
=0.3 ............................................................................... 67
Figure 4.12. Residual entrapment of gas as a function of time in a heterogeneous
aquifer (long correlation length) ............................................................................... 68
Figure 4.13. Relative change in the calculated residual entrapment of CO
2
between
the co-current flow model and mixed-flow model ................................................... 69
Figure 5.1. Schematic of simulation model ...................................................................... 77
Figure 5.2. Heterogeneous permeability field with a medium correlation length ............ 79
Figure 5.3. a) Saturation profiles and b) trapped saturation profiles of CO
2
for the four
aquifer settings in the homogeneous field with a horizontal permeability of 500
mD and k
v
/k
h
=0.01 .................................................................................................... 85
xiii

Figure 5.4. Residual entrapment of CO
2
vs. time for the a) homogeneous aquifer, b)
heterogeneous aquifer with k
v
/k
h
=0.01 .................................................................... 86
Figure 5.5. a) Saturation profiles and b) trapped saturation profiles of CO
2
for the four
aquifer settings in the homogeneous field with a horizontal permeability of 500
mD and k
v
/k
h
=0.1 ...................................................................................................... 88
Figure 5.6. a) Saturation profiles and b) trapped saturation profiles of CO
2
for the four
aquifer settings in the heterogeneous field with a horizontal permeability of 500
mD and k
v
/k
h
=0.1 ...................................................................................................... 89
Figure 5.7. Residual entrapment of CO
2
vs. time for the a) homogeneous aquifer, b)
heterogeneous aquifer with k
v
/k
h
=0.1 ...................................................................... 90
Figure 5.8. a) Saturation profiles and b) trapped saturation profiles of CO
2
for the four
aquifer settings in the homogeneous field with a horizontal permeability of 500
mD and k
v
/k
h
=1 ......................................................................................................... 91
Figure 5.9. Residual entrapment of CO
2
vs. time for the homogeneous aquifer with
k
v
/k
h
=1 ...................................................................................................................... 92
Figure 5.10. Gravity number and residual entrapment ..................................................... 93
Figure 5.11. Segregation distance vs. permeability ratio .................................................. 94
Figure 5.12. Trapped fraction of CO
2
in the co-current and mixed flow models with
k
v
/k
h
=0.1 ................................................................................................................... 96
Figure 5.13. a) Saturation profiles and b) trapped saturation profiles of CO
2
in the co-
current and mixed modes at the end of injection and end of simulation for the
homogeneous deep aquifer with k
v
/k
h
=0.1 .............................................................. 97
Figure 5.14. Trapped fraction of CO
2
in the co-current and mixed flow models with
k
v
/k
h
=0.01 ................................................................................................................. 98
Figure 5.15. a) Saturation profiles and b) trapped saturation profiles of CO
2
in the co-
current and mixed modes at the end of injection and end of simulation for the
heterogeneous deep aquifer with k
v
/k
h
=0.1 ............................................................. 99
Figure 5.16. Residual entrapment of CO
2
in the homogeneous and heterogeneous
deep aquifer (k
v
/k
h
=0.1) with the co-current and mixed flow models .................... 100
Figure 6.1. Co-current and counter-current relative permeability curves ....................... 104
Figure 6.2. Capillary pressure function used in the simulation studies of this chapter .. 105
xiv

Figure 6.3. CO
2
saturation distribution as predicted by the co-current (left column),
linear interpolation mixed (middle column) and discrete mixed model (right
column) at 2500 days (top row), 4000 days (middle row), and 10000 days
(bottom row) ........................................................................................................... 111
Figure 6.4. The difference in CO
2
saturation between the co-current and the mixed
models (left: linear interpolation; right: discrete) at 2500 days (top row), 4000
days (middle row), and 10000 days (bottom row) .................................................. 112
Figure 6.5. Trapped CO
2
saturation profiles in the three models at 2500 days (top
row), 4000 days (middle row), and 10000 days (bottom row) ............................... 113
Figure 6.6. Fraction of the injected CO
2
that is trapped as a function of time ............... 114
Figure 6.7. Relative permeability variations in a cell (#546) located 75m from the top
of the domain and half way along the length of the aquifer as observed in the
three mobility models ............................................................................................. 115
Figure 6.8. CO
2
saturation profile in the top layer in the co-current (left column),
linear interpolated mixed (middle column) and discrete mixed (right column)
models at 2500 days (top row), 4000 days (middle row), and 10000 days (bottom
row) ......................................................................................................................... 117
Figure 6.9. CO
2
saturation profile in the three models in the 8th layer (from the top) in
the three models at 2500 days (top row), 4000 days (middle row), and 10000
days (bottom row) ................................................................................................... 117
Figure 6.10. Trapped CO
2
saturation as predicted by the three mobility models; left
column: co-current; middle column: mixed (linear interpolation); right column:
mixed (discrete), at 3000 days (1
st
row), 4000 days (2
nd
row), 6000 days (3
rd

row) and 12000 days (4
th
row) ................................................................................ 118
Figure 6.11. Fraction of the injected CO
2
that is trapped as a function of time ............. 119
Figure 6.12. Relative permeability variations in a cell (#16273) located just above the
injection zone as observed in the three mobility models ........................................ 120
Figure 6.13. Relative difference in the fraction of CO
2
that is trapped between the co-
current and mixed models after 20000 days ........................................................... 122
Figure 6.14. CO
2
saturation profile in the top layer in the co-current (left column),
linear interpolated mixed (middle column) and quadratic mixed (right column)
models at 2500 days (top row), 4000 days (middle row), and 10000 days (bottom
row) ......................................................................................................................... 123
Figure 6.15. Fraction of the CO
2
trapped in the co-current and mixed models .............. 124
xv

Figure 7.1. Relative permeability of the oil-water system .............................................. 131
Figure 7.2. Relative permeability of the gas-oil system ................................................. 131
Figure 7.3. Methane saturation profiles as predicted by co-current model (left), mixed
model (middle) and their difference (right) in the 3
rd
layer from the top ............... 135
Figure 7.4. Oil saturation profiles as predicted by co-current model (left), mixed
model (middle) and their difference (right) in the 3
rd
layer from the top ............... 136
Figure 7.5. Water saturation profiles as predicted by co-current model (left), mixed
model (middle) and their difference (right) in the 3
rd
layer from the top ............... 137
Figure 7.6. Water saturation profiles as predicted by co-current model (left), mixed
model (middle) and their difference (right) in the 5
th
layer from the top ............... 138
Figure 7.7. Comparison of injection and production rates between co-current and
mixed models in layer 2 from top ........................................................................... 139
Figure 7.8. Bottomhole pressure in one of the injection wells (injecting water) ............ 140
Figure 7.9. Gas saturation profiles as predicted by co-current model (left), mixed
model (middle) and their difference (right) in the 1
st
layer from the top ................ 142
Figure 7.10. Oil saturation profiles as predicted by co-current model (left), mixed
model (middle) and their difference (right) in the 1
st
layer from the top ................ 143
Figure 7.11. Water saturation profiles as predicted by co-current model (left), mixed
model (middle) and their difference (right) in the 8
th
layer from the top ............... 144
Figure 7.12. Comparison of injection and production rates between co-current and
mixed models .......................................................................................................... 145

1

Chapter One: Introduction

1.1 Background
Carbon dioxide (CO
2
) is one of the main components of greenhouse gases (GHG) which absorb
the heat radiation from the surface of the earth that eventually causes the earth’s temperature to
rise. This increase in the temperature is known as the greenhouse effect. The atmospheric
concentration of CO
2
rose from 280 ppm in 1800 to 370 ppm in 2000 mainly due to the
consumption of fossil fuels (Fankhauser et al., 2001). The CO
2
concentration may continue to
increase to between 500 and 1000 ppm by the year 2100 (Houghton et al., 2001). In the last two
decades, there has been worldwide concern regarding global warming and its effect on the
planet. Some may dispute the validity of these findings, but among the scientific community
there is a consensus about this issue (Oreskes, 2004) and that the increase in temperature in the
last 50 years is attributed to the increased concentration of greenhouse gases (Inter-governmental
Panel on Climate Change, 2001).
There are some ways to reduce emission of CO
2
into the atmosphere: (1) improved/alternate
energy uses, (2) CO
2
capture and utilization, and (3) CO
2
capture and storage (CCS) or long-term
disposal (Bachu, 1994). Carbon capture and storage is an approach to mitigate global warming
by capturing CO
2
from large point sources and storing it instead of releasing it into the
atmosphere.
Reduction of CO
2
concentration in the atmosphere requires either abandoning fossil fuel
energy or storing most of the CO
2
currently released to the atmosphere. Fossil fuels will remain
the mainstay of energy production in the 21
st
century (U.S. Department of Energy (DOE), 2004)
until renewable energy resources are developed and applied commercially. However, reducing
2

the stream of CO
2
entering the atmosphere is not a simple task. Fossil fuels, which provide about
85 percent of the world’s energy (the rest is made up of renewable, nuclear, and hydro-electricity
energy), are made of hydrocarbons, and burning them releases large quantities of CO
2
(Herzog
and Golomb, 2004). Even as renewable energy sources emerge, fossil-fuel burning will remain
substantial for years. The fossil fuel in greatest supply, coal, emits the most CO
2
per unit of
energy produced. A great challenge for the engineers of the 21
st
century will be to develop
systems for capturing the CO
2
produced by burning fossil fuels and sequestering it safely away
from the atmosphere (National Academy of Engineering (NAE), 2008). CO
2
can also be
captured from power plants and be utilized in a variety of industries. However, this is considered
to only delay CO
2
emissions to the atmosphere and is not considered a solution for the long-term
fate of CO
2
. Furthermore, the energy currently needed to capture and transport CO
2
ultimately
leads to the production of more CO
2
(Bachu et al., 1994). Therefore, long-term CO
2
disposal is
the only remaining option for reducing CO
2
concentration in the atmosphere.
Sequestration is the capture and safe storage of CO
2
that could potentially be emitted to the
atmosphere (Tang, 2006). Sequestration refers to any storage plan that can keep CO
2
out of the
atmosphere (Park, 2005). In general, CO
2
storage sites can be divided into two categories:
geological sites and marine sites. Oceanic sequestration needs more research in areas such as
developing technology, reducing the cost and investigating unknown physical and chemical
processes. Geological sequestration has been identified as one of the most promising and viable
sequestration technologies. Geological sequestration of CO
2
is the capture of CO
2
from major
sources, transporting it usually by pipeline, and injecting it into underground formations such as
oil and gas reservoirs, saline aquifers, and unmineable coal seams for a geologically significant
period of time (Bachu, 2002; U.S. DOE, 2004). Figure 1.1 shows different kinds of underground
3

sites for a geological sequestration of CO
2
. Any viable system for sequestering CO
2
must meet
the following requirements: (1) effective and cost competitive, (2) stable for long-term storage,
and (3) environmentally friendly (Park, 2005; U.S. DOE, 2004).

Figure 1.1. Deep underground storage of CO
2
(IPCC, 2005)

CO
2
sequestration in deep geological formations has been suggested as a way of reducing
GHG emissions. The capacity of oil reservoirs for CO
2
disposal is limited (Bachu, 1994; Seo,
2004). Injection of CO
2
into gas fields can be accomplished once giant gas fields are fully
depleted in the next decades. Unlike coal bed methane reserves and oil reservoirs, sequestration
of CO
2
in deep saline aquifers does not produce value-adding by-products, but it has other
advantages. They are generally unused and are available in many parts of the world (Seo, 2004).
4

It has been estimated that deep saline formations in the United States could potentially store up
to 500 gigatonnes (Gt) of CO
2
(U.S. DOE, 2004). It is estimated that deep saline formations are
likely to have a capacity of at least 1000 Gt worldwide (IPCC, 2005), while some studies show a
capacity of an order of magnitude greater (Herzog and Golomb, 2004; Bachu, 2002). Currently,
30 Gt per year of CO
2
is emitted due to human activities (Marland, et al., 2006). Most existing
large CO
2
point sources are also within easy access of a saline formation. Therefore,
sequestration in saline formations is compatible with a strategy of transforming large portions of
the existing U.S. energy and industrial assets to near-zero carbon emissions via low-cost carbon
sequestration retrofits (U.S. DOE, 2004). However, it is important to investigate the behaviour of
CO
2
injected into aquifers for effective and safe use of storage.

1.2 Motivation
The use of reservoir simulation in multiphase flow in porous media is very mature. Existing
reservoir simulation softwares are used for the displacement of brine by CO
2
in saline aquifers.
Currently, the effect of counter-current flow on the mobility of phases is neglected in simulation
packages. However, as our study shows, it can have a great impact on the vertical migration of
CO
2
in saline aquifers. We try to include a general mobility function in the simulation of
multiphase flow in porous media which accounts for neglecting coupling terms in the
generalized Darcy equation.

1.3 Research Objectives
The objective of this study is to investigate the effect of a dynamic relative permeability
representation in the simulation of two- and three-phase flows in porous media. In conventional
5

simulation of multiphase flow in porous media, it is assumed that the driving force (pressure
gradient and gravity) on one phase does not affect the flow of another phase. Experimental
(Lelievre, 1966; Bentsen and Manai, 1993; Dullien and Dong, 1996) and theoretical studies
(Whitaker, 1986; Kalaydjian, 1987) have shown that the viscous coupling occurs in porous
media and that sometimes this effect cannot be neglected (Dullien and Dong, 1996).
It is assumed that viscous coupling is negligible in the simulation of CO
2
injection into saline
aquifers. We show that the vertical migration of CO
2
in saline aquifers, due to buoyancy, is
dominated by counter-current flow. Experimental (Lelievre, 1966; Bentsen and Manai, 1993)
and Simulation (Bourbiaux and Kalaydjian, 1990; Li et al., 2005) studies show that the counter-
current relative permeability is less than co-current relative permeability. The difference between
the two is believed to be the result of viscous coupling between the phases (Bourbiaux and
Kalaydjian, 1990). The degree of the coupling defines the differences between the two modes of
flow (co-current vs. Counter-current). In this study we investigate the impact of integrating this
effect in the simulation of CO
2
injection into saline aquifers, and particularly on residual
entrapment of CO
2
. We start by studying two-phase flow problems and extend the scope of the
study to three phase flow. In the two-phase flow modeling we use the implicit pressure and
explicit saturation (IMPES) and fully-implicit formulations while in the three-phase flow
modeling fully-implicit method is used.
The final chapter of the manuscript tries to answer whether capillary pressure plays a
significant role in the simulation of CO
2
injection into saline aquifers or not.

6

1.4 Organization of the Manuscript
The manuscript is organized as follows. We start (Chapter 2) with a literature review of counter-
current relative permeability. Chapter three presents the results of a gravity segregation
experiment in a synthetic porous material. The results show that counter-current relative
permeabilities are lower than co-current values. Chapter four discusses the significance of
including the dynamic relative permeability (mobility) in the simulation of two-phase flow of
CO
2
and brine. Chapter five describes the impact of considering counter-current relative
permeability in the simultaneous injection of CO
2
and brine (as a methodology to increase
storage capacity of saline aquifers) in four aquifer settings reported by Kopp et al. (2009).
Chapter six includes velocity-dependent relative permeability in the two-phase flows. Chapter
seven demonstrates the potential impact of velocity-dependent relative permeability in three-
phase flows. The effect of capillary pressure and its representation is presented in chapter eight
while chapter nine provides a summary of the dissertation and recommendations for future
research are discussed.

7

Chapter Two: Literature Review

2.1 Background
Multiphase flow in porous media is extended from the single phase flow with the addition of
relative permeability to Darcy’s equation. Relative permeability is a function of phase
saturations, saturation history and flow velocity (non-Darcy flow). It is assumed that the flow of
one phase does not affect the flow of another phase (there is no viscous coupling between the
phases). The validity of this assumption is questioned by some experimental and numerical
studies.
In this chapter we start by reviewing studies on viscous coupling and counter-current
relative permeability in porous media. Then studies related to the anisotropy of relative
permeability (mostly in bedding and layered systems) are presented.

2.2 Viscous Coupling in Multiphase Flow in Porous Media
Muskat (1946) proposed the effective permeability to be used in the Darcy’s equation for
multiphase flow in porous media. Later experimental (Odeh, 1959; Lelievre, 1966; Lefebvre du
Prey, 1978; Bourbiaux and Kalaydjian, 1988; Maini et al., 1989) and theoretical studies (de la
Cruz and Spanos, 1983, 1988; Whitaker, 1986; Kalaydjian, 1987) showed that the use of Muskat
relative permeability may not be accurate.
Some researchers (de la Cruz and Spanos, 1983, 1988; Whitaker, 1986; Kalaydjian, 1987)
have used the volume averaging of the Navier-Stokes equation to arrive at the two-phase flow
equation of immiscible fluids in porous media:
8

(2.1)

(2.2)
In equations (2.1) and (2.2), k
11
, k
12
, k
21
and k
22
are the four independent generalized permeability
coefficients. The coefficients k
12
and k
21
are related to the viscous drag between the two phases.
In fact, there are two types of drag; one due to the flow of the fluid over the solid particles and
the other one is the momentum transfer between the two fluid phases. The four permeability
coefficients are related to the traditional effective permeability of the two phases, but there is a
coupling between the two equations. Unlike the standard relative permeability measurements,
two experiments are needed to determine the four generalized permeability coefficients. Viscous
coupling is supposed to be more significant in three-phase flow (Dehghanpour et al., 2011), but
most of the research in this area is limited to the two-phase displacements.

2.3 Counter-Current Relative Permeability
The first experiments to compare co-current and counter-current relative permeabilities based on
Muskat's equation were presented by Lelievre (1966). He measured co-current and counter-
current relative permeabilities for different viscosity ratios and found that the relative
permeability in counter-current flow settings is always less than the relative permeability in co-
current flow. The observed difference between co- and counter-current flows was larger for the
wetting phase and in some cases the wetting phase relative permeability was reduced by as much
as 50%.
9

Bourbiaux and Kalaydjian (1990) studied the dynamics of spontaneous imbibition in co-
current and counter-current experiments as related to oil recovery from a porous media. They
observed that the oil recovery rate was lower in counter-current flow settings than in co-current
flow settings. They subsequently simulated their experiments and found a significant difference
between experimental and simulated recoveries for counter-current flow settings. By decreasing
the oil relative permeabilities by 60% from the co-current values, or the water relative
permeabilities by 45% from their co-current values, or by decreasing both oil and water relative
permeabilities by 30% from their co-current values, they managed to reproduce the observed oil
recovery more accurately. They attributed this reduction of relative permeabilities in counter-
current flow to viscous coupling that acts differently in co-current and counter-current flow
settings.
Eastwood and Spanos (1991) studied two-phase flow in porous media and found that the
relative permeability in counter-current flow is less than the relative permeability in co-current
flow. Based on the two-phase flow model of de la Cruze and Spanos (1983), they derived that
the Muskat relative permeabilities for oil in counter-current flow settings are less than those for
co-current flow settings.
Bentsen and Manai (1993) measured relative permeabilities in co-current and counter-
current steady state experiments to find the four independent generalized permeability
coefficients needed for modeling two-phase flow in porous media. They used oil and water in
their experiments with water as the wetting phase. They used a microwave approach to measure
saturations in-situ. The results of their experiments show that counter-current relative
permeabilities are always less than co-current relative permeabilities for both phases and for all
saturations studied in the experiments. The counter-current relative permeability for water
10

(wetting phase) was always at least 25% less than the co-current relative permeability. Relative
permeabilities of oil (non-wetting phase) were always at least 20% less in counter-current flow
compared to co-current flow. From the experiments, they concluded that the generalized
coefficients lie midway between co-current and counter-current relative permeability curves.
Bentsen (1998) studied the effect of viscous coupling on the effective mobility of the
wetting and non-wetting phases in two-phase flow based on experimental observations and
theoretical analysis. He measured co-current and counter-current relative permeabilities of both
phases to validate his work. In his work, the counter-current relative permeabilities of both
phases were found to be approximately 70% of those for co-current flow. This conclusion is
consistent with the results of Bourbiaux and Kalaydjian (1990). Based on the experiments, he
also concluded that the cross mobilities are not equal.
Al-Wadahi et al. (2000) studied counter-current three-phase flow in a porous media driven
by gravity and capillary forces. They conducted three-phase flow experiments and used a
commercial simulator to interpret their observations. Multiphase flow was achieved using water
(at irreducible saturation), benzyl alcohol, and decane. Benzyl alcohol and decane were the
mobile phases in the counter-current experiment. They applied artificial neural networks to
match the experimental observations and to derive a set of relative permeability and capillary
pressure functions. Based on their results, the relative permeabilities of both decane and benzyl
alcohol in co-current setting were decreased significantly to arrive at a good agreement between
numerical calculations and experimental observations for the counter-current experiment. The
required reductions in the three-phase relative permeabilities were more significant than those
reported for two-phase experiments: Benzyl alcohol relative permeability was reduced by a
11

factor of 2 while decane relative permeability was reduced by a factor of 5 compared to co-
current flow settings.
Langaas and Papatzacos (2001) investigated variations in co-current and counter-current
steady state relative permeability as a function of saturation, wettability, driving force and
viscosity ratio based on displacement calculations in a uniform pore structure. Their results
demonstrate that the relative permeabilities in counter-current settings are always less than those
in co-current flow settings. They attributed this difference to viscous coupling, and to different
levels of capillary forces.
Li et al. (2005) used the same procedure of Al-Wadahi et al. (2000) and verified the
observations. They included hysteresis in the capillary pressure and relative permeabilities in
their numerical calculations, and concluded that counter-current relative permeabilities are less
than co-current relative permeabilities.
Shad and Gates (2008) studied segregated multiphase flow in a single fracture for co-
current and counter-current flow regimes. They use the Stoke’s equation (capillary is ignored)
and conclude that the relative permeability is not only a function of phase saturations but also
flow pattern and the relative permeability curves used for flow through fractures in co-current
regime cannot directly be applied to counter-current regime.
The experimental and numerical studies summarized above all suggest that relative
permeabilities can be significantly different between co-current and counter-current flow
settings. This general observation can be attributed to the fact that in co-current flows the
pressure gradients of the flowing phases are aligned and the flow of one phase assists the flow of
the other phase. In counter-current flow, on the other hand, pressure gradients are in opposite
directions and the flow of one phase may interfere with the flow of the other phase. In other
12

words, the viscous coupling between the two (or more) flowing phases may have a negative
impact on the overall mobility in counter-current flows and a neutral to positive impact on co-
current flow settings.

2.4 Relative Permeability Anisotropy
Regardless of viscous coupling, relative permeability can be anisotropic even in a homogeneous
and isotropic porous medium (Prats and Lake, 2008). In the industry, the same relative
permeability is used in horizontal and vertical displacements, except for gravity drainage. Corey
and Rathjens (1956) investigated the effect of stratification on relative permeability for flow
parallel and perpendicular to bedding. A course and fine grain stratified core samples were
combined and oil and gas were allowed to flow parallel and perpendicular to the core interface
and the resulting relative permeability were measured. It was found out that the gas relative
permeability varied dramatically between the two settings and was much lower in series
composites. Honarpour et al. (1994) reports the dependency of relative permeability on small-
scale rock laminations while Iverson et al. (1996) reported differences in relative permeability
between flow parallel and perpendicular to bedding planes in Tansleep sandstones and the
consequent effects on oil recovery. Crotti et al. (1998) studied vertical and horizontal floods on
core sand samples and showed that the measured relative permeability in horizontal and vertical
floods were different. Directional relative permeability is included in Schlumberger’s Eclipse
black-oil simulator (Eclipse Manual, 2010).

13

2.5 Three-phase Relative Permeability Models
Many models have been introduced in the literature for the calculation of three-phase relative
permeabilities (Corey, 1956; Stone, 1970; Stone 1973; Dietrich and Bonder, 1976; Aziz and
Settari, 1979; Fayers and Matthews, 1984; Aleman and Slattery, 1988; Baker, 1988; Delshad and
Pope, 1989; Kokal and Maini, 1989; Balbinski et al., 1997; Moulu et al., 1999; Blunt, 1999).
These models assume that the wetting and non-wetting relative permeabilites depend only on
their own saturations, k
rw
=f(S
w
) and k
rg
=f(S
g
), in a water-wet system (Larsen and Skauge, 1995,
Baker, 1995). However, the oil relative permeability depends on the saturation of two phases,
and the difference between these models is in the prediction the oil relative permeability.

14

Chapter Three: Non-Wetting Phase Plume Migration in Vertical Two-Phase Counter-
Current Flow
1

3.1 Introduction
In this chapter an experimental approach to study two-phase counter-current flow settings in
porous media is discussed. The modeling results are used to infer the necessity to incorporate
counter-current relative permeability in the simulation of multi-phase flows where gravity plays
a significant role. The experiments were performed at small scale. However, the results are used
for large-scale problems in the context of CO
2
sequestration.
The migration of CO
2
after injection into an aquifer is governed by multiphase flow
phenomena in porous media. Reservoir parameters (such as pressure, temperature and absolute
permeability) are not the only factors that affect the paths and time scale for CO
2
migration.
Saturation functions, such as relative permeability and capillarity, are also important in dictating
feasible injection rates and subsequent plume migration (Juanes et al., 2006; Burton et al., 2009;
Saadatpoor et al., 2010).
Relative permeability and the parameters that affect it have been researched extensively,
especially in the petroleum industry (Honarpour et al., 1986). In the context of CO
2
/brine
systems, Flett et al. (2004) and Juanes et al. (2006) demonstrated that hysteresis in relevant
saturation functions can have a significant impact on the long-term immobilization of CO
2
.

Bennion & Bachu (2008, 2010) summarized the characteristics of co-current drainage and

1
The results of this chapter is published in the International Journal of Greenhouse Gas Control:
Javaheri, M., Nattwongasem, D., Jessen, K., 2013. Relative permeability and non-wetting phase plume migration in
vertical counter-current flow settings, International Journal of Greenhouse Gas Control, 12, pp. 168-180,
doi:10.1016/j.ijggc.2012.10.006.
15

imbibition relative permeabilities of CO
2
/brine and H
2
S/brine systems for various rock types.
Krevor et al. (2012) measured the relative permeability of four rock samples and confirmed that
the residual/trapped CO
2
saturation of all four samples depend on the maximum CO
2
saturation
before imbibition.
Numerical simulators that are currently used in the modeling of CO
2
injection into saline
aquifers assume, implicitly, that saturation functions (e.g. relative permeability) do not depend
on flow direction. However, the migration of the injected CO
2
in an aquifer occurs in both co-
and counter-current flow settings. The lateral flow of CO
2
and brine is in the same direction
away from the injection zone (co-current flow), while the upward flow of CO
2
due to buoyancy
and the corresponding downward flow of brine are in opposite directions (counter-current flow).
The vertical counter-current flow is especially dominant when the viscous force is smaller than
gravity (e.g. after injection). The difference between co- and counter-current flow and its impact
on CO
2
migration and entrapment was addressed through numerical calculations by Javaheri and
Jessen (2011). They applied the dependency of flow directions on relative permeability and
demonstrated that a reduction in relative permeability functions due to transition from co-current
to counter-current flow (in the vertical direction) increases the residual entrapment of CO
2
and
retards its upward migration.
The assumption that the relative permeability is the same in co- and counter-current flows
is not supported by other experimental and theoretical studies that have demonstrated a reduction
in the relative permeability functions of counter-current flows relative to the co-current values
(Lelievre, 1966; Bentsen and Manai, 1993). In the previous chapter a summary of counter-
current relative permeability was presented that shows almost all direct and indirect
16

measurements show a reduction of relative permeability in counter-current flows. It was shown
that the general form of Darcy’s equation can be written as (Li et al., 2005b):

(3.1)
where K is the absolute permeability, k r,ij are the generalized relative permeabilities, is the
viscosity, p is the pressure, ρ is the density and subscripts w and n refer to wetting and non-
wetting phases.
Although the effects of viscous coupling in two-phase flow are well documented in the
litterature, commercial and research simulators do commonly not allow for the use of a relative
permeability tensor as dictated by Equation (3.1). Furthermore, standardized experimental
protocols for measurement of the generalized relative permeabilities have to our best knowledge
not been established. Accordingly, most flow simulations in subsurfance settings are based on
implementations of Muskat’s extension of Darcy’s equation that do not destinguish between co-
current and counter-current flow.
To validate the applicability and accuracy of co-current relative permeability in settings
where fluid displacements occur in a counter-current setting, as is the case for CO
2
injection into
saline aquifers, we have designed experiments that attempt to minimize the effects of boundary
conditions: An initial fluid distribution is established and the wetting-phase saturation is
measured indirectly throughout the segregation process based on the observed changes in
resistivity along the porous medium.
The objective of this experiment is to test the accuracy of saturation functions that are
measured in co-current flow settings, in predicting the segregation of two fluids in counter-
current buoyancy-driven settings. Four-electrode resistivity measurements were used to estimate
17

the wetting-phase saturation along the porous column as a function of time. To represent the
dynamics of CO
2
migration in saline aquifers the experiments were performed at ambient
condition using analog fluids, brine/iso-octane (iC
8
), in a glass-bead pack. Analog fluids have
also been used in the study of mixing and mass transfer in the context of miscible displacements
(Shojaei et al., 2012) iC
8
has a density that is similar to that of supercritical CO
2
, but its viscosity
is higher and iC
8
represents as such not the perfect analog. However, the higher viscosity of iC
8

eliminates a potentially unstable non-wetting phase front that would otherwise complicate and
potentially prevent the interpretation of relative permeability effects from the experimental
observations. Prior to the segregation experiments, the packed column was characterized in terms
of porosity, permeability, co-current relative permeability and static capillary pressure.
We use numerical calculations to interpret our experimental observations and adjust the
input parameters to match the saturations observed in the experiments. The results and
observations of this study are then placed in the context of large-scale CO
2
migration in saline
aquifers and we deduce that it is necessary to incorporate the effect of counter-current flow of
CO
2
and brine for accurate simulation of plume distribution and migration time scales.
In the following sections, we present a detailed description of the experimental setup
including characterization of fluid and porous medium properties with emphasis on model input
parameters for numerical calculations. The experimental observations and a comparison with
numerical calculations are presented subsequently, followed by a discussion of the key findings
from our analysis of segregation experiments and modeling efforts.

18

3.2 Experimental Approach
3.2.1 Apparatus and characterization
A glass column packed with glass beads served as porous medium in this work with emphasis on
water-wet systems. The porous medium consists of a borosilicate glass column with adjustable
plungers at both ends. The glass column has an inner diameter of 50 mm and a height of 500
mm. Stainless steel frits (type 316) with an average pore size of 10 m were attached to the end
of each plunger. The frits were used as current electrodes and to create a uniform flow
distribution of injected fluids. Holes were drilled along the plungers to house electrical wires that
connected the metal frits. This allowed us to connect a power source without disturbing the flow.
The glass column was further modified by drilling holes along the length to house the potential
electrodes. Three holes were drilled on the perimeter of the column, at selected locations, with a
spacing of 120 degrees. A total of 36 holes were drilled along the glass column to establish a
total of 11 equally spaced sections. Stainless steel wires (type 316) were used as electrode
material to be consistent with the frits. At each level along the column, the three wires were
connected together and used as the potential electrodes. Figure 3.1 shows the design of the glass
column used in our segregation experiments including the modifications needed for our four-
electrode resistivity measurements.
To represent and study the migration of a supercritical CO
2
plume in a saline aquifer, we
performed two segregation experiments at low-pressure using the synthetic porous material and
analog fluids. We used BT13 (170 Mesh) glass beads with an average particle diameter of 88 μm
to represent the porous medium, while the immiscible two-phase brine/iC
8
fluid system was used
to represent brine/supercritical CO
2
at reservoir conditions. The density of iC
8
is similar to
supercritical CO
2
at high pressure, but its viscosity is higher. The higher value of the non-wetting
19

phase viscosity provides for a more favorable mobility ratio. This, in turn, inhibits the onset of
viscous instabilities and facilitates the interpretation of the experimental observations.

Figure 3.1. Dimensions and design of the Borosilicate glass column for resistivity measurement

Brine was prepared from deionized water and NaCl at a concentration of 20,000 ppm
resulting in a density of 1013 kg/m
3
at 21C. At similar conditions, the non-wetting phase (iC
8
)
has a density of 692 kg/m
3
that is comparable to high-pressure CO
2
(CO
2
density is about 700
kg/m
3
in the pressure range of 10 MPa (at 35 C) and 20 MPa (at 65 C). To visualize the
propagation of the non-wetting phase in the column we used an oil soluble dye (Sudan Red 7B).
The viscosity of brine and iC
8
are 1.0 and 0.48 mPa.s, respectively, as reported by Cinar et al.
(2006). The viscosity of iC
8
is higher than that of supercritical CO
2
(e.g. approximately 0.06
mPa.s at 10 MPa and 35 C) and one may argue that iC
8
is not a perfect analog. However, the
more favorable viscosity ratio of approximately 2 for brine/iC
8
relative to approximately 16 for
20

brine/CO
2
promotes a stable non-wetting phase front and suppresses viscous instabilities in the
packed column. This in turn enables us to study relative permeability effects in counter-current
flow settings that would otherwise be unfeasible. We address the higher mobility ratio of
brine/CO
2
systems in the context of our findings in the discussion section.
Pendant drop measurements were performed to determine the interfacial tension of the
fluid system and a value of 47.3 mN/m was observed from repeated measurements. This IFT
value corresponds to the IFT of a 200,000 ppm brine/CO
2
system at 75 C and 7.5MPa according
to Bachu and Bennion, (2009).
We used a dry-packing approach where glass beads were loaded into the glass column with
one plunger in place while resting on a shaker to ensure an effective packing of the column. The
top plunger was then inserted to provide a good support for the beads. The column was
subsequently evacuated from the top and brine was loaded from the bottom to fully saturate the
packed column under continuous evacuation. The porosity of the packed column was measured
to 38.6% from gravimetric considerations while the permeability was measured to 4.8 Darcy
from the stabilized pressure drop at a fixed injection rate. At fully saturated conditions, the
porosity at each section of the column was calculated from resistivity measurements and found to
be in good agreement with the overall porosity.
Dynamic measurement of in-situ brine saturations can be done at laboratory scale by using
four-electrode resistivity measurements (Rust, 1952; Lewis et al., 1988; Knackstedt et al., 2007).
In this approach, changes in the electrical potential across a given section of the porous media
are recorded. From Archie’s equation (Archie, 1942) or the resistivity index (RI), changes in
saturation can then be computed as a function of time. Nakatsuka et al. (2010) used a four-
electrode resistivity approach to measure the saturation of CO
2
in consolidated cores. Their work
21

focused primarily on the presence of clay and how the clay content influences resistivity
measurements. They found that Archie’s equation is sufficient for estimating CO
2
saturation in
clean and homogeneous sands whereas the resistivity index (RI) provides for a more accurate
estimation of CO
2
saturation in heterogeneous porous materials.
Co-current bounding relative permeability curves for primary drainage and imbibition
processes were measured by conventional steady-state flow experiments (Honarpour et al.,
1986). Figure 3.2 reports the co-current relative permeability measurements and a Corey-type
representation for the drainage process. As observed from Figure 3.2, hysteresis in the relative
permeability is significant and was hence included in subsequent numerical calculations. The
relative permeability was measured in a co-current steady-state setting and provides for initial
input to our numerical calculations and interpretation of the experimental observations. From our
steady-state relative permeability measurements we observe an irreducible wetting phase
saturation of 0.2 and a maximum trapped non-wetting phase saturation of 0.16. The endpoint
relative permeability of iC
8
is 0.15. This may seem low, but laboratory measurement of relative
permeability of CO
2
/brine systems on reservoir cores show an inverse proportionality between
the endpoint relative permeability of CO
2
and the absolute permeability of the core (Bennion and
Bachu, 2005, 2010). In another study, Bachu and Bennion (2008) reported that the maximum
relative permeability of CO
2
is about one-fifth to one-seventh of permeability to brine (at 100%
saturation) for sandstone, carbonate and shale formations from central Alberta in western
Canada. Furthermore, the low endpoint relative permeability of the non-wetting phase, observed
in our experiments, is consistent with our visual observations that flow instabilities did not occur
during the relative permeability measurements.
22

Capillary pressure functions were adapted from Dawe et al. (1992), who measured
capillary pressure for drainage and imbibition processes in glass-bead packs. As their
measurements were performed for different glass bead sizes and fluid systems, their capillary
pressure observations were initially scaled using Leverett J-scaling (Leverett, 1941). The
primary drainage and imbibition capillary pressures corresponding to the glass bead used in our
segregation experiments were then constructed from the J-functions extracted from Dawe et al.
(1992). van Genuchten (1980) capillary pressure functions were used to smooth the capillary
pressure data for both drainage and imbibition processes. To match the irreducible wetting phase
saturation and maximum trapped non-wetting saturation obtained from our steady-state relative
permeability experiments, we re-normalized the phase saturations. In addition, the drainage
capillary pressure was measured using the porous-plate technique (Figure 3.3). The bounding J-
functions used in our subsequent numerical calculations are compared with our measurements in
Figure 3.3.

Figure 3.2. Steady-state relative permeability functions for drainage and imbibition processes
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
S
w
k
r

k
r,w
, Drainage
k
r,iC
8
, Drainage
k
r,w
, Imbibition
k
r,iC
8
, Imbibition
k
r,w
, Corey
Drainage
k
r,iC
8
, Corey
Drainage
23

Figure 3.3. Normalized capillary pressure data extracted from Dawe et al. (1992)

3.2.2 Segregation Experiments
To capture the dynamic changes in the electric potential across the individual sections of the
column during our segregation experiments, two data acquisition modules from National
Instruments (NI-DAQ) were used to provide and acquire signals. A 0.5 V AC signal with
frequency of 1 kHz was generated from the analogue channel of one NI-DAQ unit. A shunt
resistor (2 kΩ) was connected in series to the top of the packed column to prevent current
overflow (and accelerated corrosion). The wires from the packed column were connected to the
analogue input channels of the second NI-DAQ module. Both NI-DAQ modules were connected
to a PC to control the units and to collect data. The potential difference across each section of the
packed column and across the shunt resistor was continuously recorded using LabView software.
Figure 3.4 shows a schematic diagram of the dynamic segregation experiment including
connections to the source and acquisition systems. Once the column was fully saturated with
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
S
w
Leverett J-function

Drainage, Dawe et al. (1992)
Imbibition, Dawe et al. (1992)
Drainage, van Genuchten fit
Imbibition, van Genuchten fit
Drainage, Measured
24

brine, the resistivity of each section of the packed column was determined to validate the average
porosity measurement as noted above.

Figure 3.4. Schematic diagram of experimental setup for segregation experiment

A Teledyne 260D syringe pump was used to inject iC
8
at the top of the column at a rate
below the critical velocity as calculated from the petrophysical properties discussed earlier
(Lake, 1989). The dyed iC
8
was then injected to create an initial non-wetting saturation
distribution (plume) as uniform as possible. After the desired amount of iC
8
was injected, the
packed column was inverted and the segregation experiment initiated (time zero). The potential
differences at each section of the packed column were continuously recorded allowing for
subsequent calculation of the wetting phase saturation as discussed below.
The electric current through the column at any given time was calculated from the potential
drop across the shunt resistor. From the potential differences across the column sections, the
resistances of the sections of the column were then calculated. The brine saturation in each
25

section of the packed column was then calculated based on the resistivity index (RI) with a
saturation exponent, n, using:

(3.2)
where S
w,j
is the brine saturation in section j of the column.

is the resistivity for section j
with 100% brine saturation and

is the resistivity of section of the packed column with two
phases present. The saturation exponent (n) was determined by minimizing the volume error for
the injected non-wetting phase using the objective function ( ):

(3.3)

(3.4)
where V
j
is the pore volume of section j of the column as calculated from the porosity (ϕ), r is the
inner diameter of the column, and l
j
is the length of section j. The minimization of f was
performed subject to the constraint:

(3.5)
This approach allows us to determine the brine saturations at each section of the packed column
as a function of time.
To repeat the segregation experiment, two pore-volumes of isopropanol were injected to
remove iC
8
and brine from the column. Deionized water was subsequently used to flush the
column and remove the isopropanol and finally, brine was loaded from the bottom to displace the
deionized water.
Two gravity segregation experiments were performed with the setup. In the first
experiment (experiment A) 90 cm
3
of iC
8
(~ 0.23 PV) was injected into the column while in the
second experiment (experiment B) 120 cm
3
of iC
8
(~ 0.31 PV) was used. Figure 3.5 reports the
26

saturation distribution in the packed column as a function of time for the two segregation
experiments. In experiment A, we observe that the non-wetting phase does not reach the top of
the column after 24 hours. This indicates that the fluids in the column reach a gravity-capillary
equilibrium before the plume reaches the top. To ensure that the plume reaches the top of the
column, experiment B was performed with a larger initial volume (plume) of the non-wetting
phase. From Figure 3.5 (right) we see that the plume in experiment B reaches the top of the
column after approximately 24 hours and continues to migrate upwards at later times until the
experiment was stopped after 72 hours.

Figure 3.5. Saturation profiles of iC
8
in the column at various times; left: Experiment A and
right experiment B

3.3 Simulation of Segregation Experiments
To simulate and interpret the experimental observations, we assume 1D fluid flow and the
physical properties of brine and iC
8
(viscosity and density) are assumed to be constants, as
determined in the experiments. An IMPES formulation was used to solve the incompressible

27

two-phase flow problem. In the numerical calculations, we use a spatially refined model to
reduce numerical diffusion and upscale (average) the calculation results to compare with the iC
8

saturations observed from the experiments. We note that a grid refinement of the model by a
factor of five is sufficient to converge the numerical calculations. The permeability is assumed to
be constant at 4.8 Darcy in all cells, while the porosity in each cell of the refined model is set
equal to the porosity of relevant section in the experimental setup. The initial distribution of the
non-wetting phase saturation in the refined calculation was constructed from a piece-wise linear
interpretation of the initial non-wetting phase saturation of the experiment, preserving the overall
pore volume of the non-wetting phase.
Corey-type relative permeability functions were used in the simulation with:

(3.6)

(3.7)
where

and

are the wetting and non-wetting relative permeability, respectively, and
is the scaled saturation defined as

. The exponents, and , and the
endpoint non-wetting phase relative permeability,

, were measured directly via steady-state
experiments. Killough’s model (Killough, 1976) was used to represent hysteresis in the relative
permeability of both phases. Residual entrapment of the non-wetting phase was calculated from
Land’s model (Land, 1968) with Land’s coefficient, C, obtained from the maximum residual
non-wetting phase saturation (

) as observed in the steady-state
relative permeability measurements.
The bounding drainage and imbibition capillary pressure functions were obtained from the
J-functions reported in Figure 3.3. The scanning curves of the capillary pressure were obtained
28

from Killough’s model (Killough, 1976) by appropriate interpolation between the drainage and
imbibition curves.
In order to include the effect of viscous drag in simulation studies with both co-current and
counter-current flows, one can use the generalized form of Darcy’s equations (Equation (3.1)
with four permeability coefficients. Another approach is to use the conventional Darcy equations
(ignoring off-diagonal terms in Equation (3.1), but to make a distinction between co-current and
counter-current flows in the calculation of the relative permeability functions (Lelievre, 1966;
Bentsen and Manai, 1993). We have chosen the second approach, since the transport equations
remain the same and introduce a relative permeability multiplier to differentiate between the
flow directions.

3.4 Results: Prediction and Optimization
We start by presenting a comparison between experimental observations and numerical
calculations based on co-current input parameters and proceed with the adjustment of the input
parameters subsequently.

3.4.1 Simulation with co-current model input
In the first step of our modeling efforts, we perform numerical calculations in a purely predictive
mode. That is, numerical calculations were performed based on measured input parameters
including bounding relative permeability measurements and capillary pressure data. We note that
the relative permeability measurements were performed in co-current flow settings. The
numerical calculations presented in this section are hence based on the co-current measurements
and assume that there is no difference between co-current and counter-current saturation
29

functions. The expected change in the wetting phase saturations were then calculated and
compared with our experimental observations. Figure 3.6 compares the saturation profiles of the
non-wetting phase as observed in the experiments with the calculated saturation profiles based
on co-current model parameters. It is clear, from Figure 3.6, that the non-wetting phase plume is
predicted to move significantly faster by the numerical calculations than what is observed in the
experiments. In the numerical calculation, the non-wetting phase is predicted to reach the top of
the column after only five hours for experiment A and after three hours in experiment B, while
this is not observed in the experiment before one day (24 hours for experiment B). Therefore, we
conclude that the numerical calculations based on co-current saturation functions fail to replicate
the observed data. This is particularly evident for experiment B, where the arrival time of the
non-wetting phase plume at the top of the column is approximately an order of magnitude in
error.
30

Figure 3.6. Saturation of iC
8
in the column observed in the experiments and predicted based on
co-current model parameters (upper panel: Experiment A, lower panel: Experiment B)

0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
Height (cm)
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw

0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
Height (cm)
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw

Experiment
Simulation
t=0 min t=10 min t=60 min t=300 min t=1440 min
t=0 min t=60 min t=300 min t=1440 min t=4320 min
31

3.4.2 Adjustment of model parameters
To obtain an improved agreement between our experimental observations and numerical
calculations, select model input parameters were adjusted by optimization. We attempt to match
our experimental observations by adjusting three parameters related to the relative permeability
of wetting and non-wetting phases. These include the saturation exponent of brine and iC
8
, and
the endpoint relative permeability of iC
8
. The selected parameters were adjusted by minimizing
the error (f) between the calculated and observed saturations:

(3.8)
where S sim is the upscaled non-wetting phase saturation obtained from the refined numerical
calculation, S exp is the corresponding saturation observed in the experiment, and is the
weight factor. Four discrete times (other than initial time) were included in the optimization of
each experiment (marked in Figure 3.6 - Figure 3.8) to honor the time-dependent saturation
profiles, with the weight factor successively increased (from one to four) to favor later times.
The end time for each experiment (one day for experiment A and three days for experiment B) is
represented by t end. We have included in the denominator, because the non-wetting phase
saturation is zero in some upper segments of the column at early times of the experiment. also
serves as a weight factor in the optimization to avoid overrepresentation of low saturation values.
We used a value of 0.2 for in the optimization. Additional discussion of the optimization
approach is provided below.
Phase relative permeabilities were measured in a co-current setting prior to the segregation
experiments. However, the segregation experiments with brine and iC
8
occur in a counter-current
flow setting. We assume here that the model parameters other than the relative permeability (e.g.
32

capillary pressure) do not depend on the flow direction and focus on modifying the relative
permeabilities to match the observations of the counter-current flow experiments.
As discussed previously, several studies (Lelievre, 1966; Bourbiaux and Kalaydjian, 1990;
Bentsen and Manai, 1993;) have demonstrated a reduction in the relative permeability functions
for counter-current flow relative to co-current values, if the conventional Darcy formulation is
used in the transport equations. In the modification of the model input parameters, two
approaches were used to mimic counter-current relative permeability:
1) Fix the endpoint values of the relative permeability for both phases and adjust the saturation
exponents and (see Equations (3.6) and (3.7) – two adjustable parameters.
2) Adjust the endpoint relative permeability for the non-wetting phase,

, as well as the
saturation exponent for both phases ( and – three adjustable parameters.
The second approach is based on the experimental measurement of counter-current relative
permeability in the brine/oil system (Manai, 1991; Lelievre, 1966) where a distinct reduction in
the endpoint relative permeability of the non-wetting phases was observed. We note that it is not
documented whether this reduction is related to the boundary conditions of the experiments or
other effects, e.g. a change in fluid distribution in the counter-current displacement compared to
that of the co-current displacement experiment.
We start by following the first approach (P1) and keep the endpoint relative permeability
constant while the difference between calculations and observations were minimized using
Equation (3.8) for both experiments with two adjustable parameters ( , ). It is noticed that the
exponents of relative permeability of both phases are increased significantly in the optimization.
Figure 3.7 compares the experimental and calculated non-wetting phase saturations using this
approach and shows a considerable improvement in agreement between calculations and
33

experiments. The upper and lower panels of Figure 3.7 show the results for experiment A and
experiment B, respectively. The saturations are shown at five different times (initial time and
four later times used in the optimization) for each experiment. The experimental observations at
these times were used in Equation (3.8).
In the second approach (P2), the endpoint relative permeability of the non-wetting phase
was included in the parameter estimation. In this approach, three parameters were used to
improve the agreement between calculations and experimental observations ( , , and k re).
Figure 3.8 compares the experimental saturation profiles for experiment A and B with the
calculated saturations after parameter adjustment (P2). We observe that the use of an additional
adjustable parameter does not improve, with any significance, the agreement between
calculations and experiments relative to the previous approaches.
From the first approach (P1) we found that the exponents of the relative permeability
increase if the endpoint relative permeability of the non-wetting phase is kept constant and that
the increase is more significant for the non-wetting phase. In the second approach (P2), the
endpoint is decreased and the exponents of the relative permeability increase relative to the co-
current inputs. The adjusted parameters from the second approach are consistent with previous
studies where the endpoint relative permeability of the non-wetting phase is reduced during
counter-current flow. If the non-wetting phase endpoint relative permeability is kept constant, the
saturation exponents increase further in the parameter estimation. In the second approach, the
reduction of the relative permeability of the non-wetting phase during counter-current flow may
be a consequence of factors beyond viscous drag: A new configuration of the two phases in the
pore space can e.g. result in additional drag. The main contribution to the reduction in relative
34

permeability cannot be isolated from the presented experiments and additional work is warranted
on this topic.
35

Figure 3.7. Comparison of experiment and simulation for experiment A (upper panel) and
experiment B (lower panel) when endpoint non-wetting phase relative permeability is kept
constant; Simulation input was adjusted based on experiment A only

0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
Height (cm)
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw

0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
Height (cm)
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw

Experiment
Simulation
t=0 min t=10 min t=60 min t=300 min t=1440 min
t=0 min t=60 min t=300 min t=1440 min t=4320 min
36

Figure 3.8. Comparison of observations and simulation for experiment A (upper panel) and
experiment B (lower panel) when endpoint non-wetting phase relative permeability is kept
constant; Simulation input was adjusted based on both experiments

0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
Height (cm)
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw

0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
Height (cm)
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw
0 0.5
0
5
10
15
20
25
30
35
40
45
50
S
nw

Experiment
Simulation
t=0 min t=10 min t=60 min t=300 min t=1440 min
t=0 min t=60 min t=300 min t=1440 min t=4320 min
37

Table 3-1 summarizes the adjustment of the parameters and shows that the sum of relative
errors (f) is reduced significantly relative to the initial input by any of the two approaches used in
matching the experimental observations.

Table 3-1. Comparison of observed and adjusted input parameters
Experiments\Parameters

f
Co-current 3.8 1.7 0.15 25.4
Adjustment (P1) 4.3 2.7 0.15 11.5
Adjustment (P2) 4.4 2.4 0.11 10.9

Finally, we compare the measured and adjusted relative permeabilities (from P1 and P2) in
Figure 3.9. The difference between measured and adjusted relative permeabilities clearly
demonstrates that the relative permeability measured in a co-current flow regime cannot be used
to predict the dynamics in a counter-current flow setting accurately.
38

Figure 3.9. Measured and adjusted relative permeability for wetting (brine) and non-wetting
(iC
8
) phases

3.5 Discussion and Conclusion
Direct measurement of counter-current relative permeability is not trivial. In this study, we used
an alternate approach to show that relative permeabilities that are obtained in a co-current flow
setting cannot be applied directly in the prediction of dynamics in counter-current flow settings.
A comparison of counter-current measurements and calculations based on co-current
observations (relative permeability) show that the migration time of a non-wetting phase plume
may be misrepresented by approximately an order of magnitude. In the context of CO
2
storage in
aquifers, a mismatch between calculations and actual migration time scales of this order is
clearly not acceptable. Figure 3.10 demonstrates the sensitivity of the time it takes the non-
wetting phase to reach the top of the column to the relative permeability of the non-wetting
phase. In this series of numerical calculations, the non-wetting phase relative permeability was
39

gradually reduced from its co-current value by a multiplier. The relative permeability multiplier
does not change the curvature of the relative permeability, and accordingly, the coefficient at
which the error in migration time is minimum, differs from the adjusted values reported in
previous section. The error bars shown in Figure 3.10 include contributions from temporal errors
introduced by the sampling frequency as well as saturation errors introduced by the four-
electrode measurements.

Figure 3.10. The sensitivity of iC
8
migration time (to the top of the column) to relative
permeability multiplier

While our parameter estimation does not result in a unique set of model parameters (P1
and P2 both substantially match the experimental observations), the approach clearly confirms
that there is a considerable difference between the input parameters based on co-current
observations and the adjusted parameters that provide for an improved agreement between
calculations and experimental observations of the segregation process. Furthermore, the
40

agreement between calculations and experiments does not appear to depend strongly on the
endpoint relative permeability of the non-wetting phase (which is reduced in P2). This is, in part,
due to the range of non-wetting phase saturations that is observed during the experiments: The
brine saturation never approaches the irreducible saturation.
The reduction in the relative permeability that is required to match the experimental
observations can be explained by the interfacial effects (viscous and capillary) at the brine/iC
8

interfaces. This viscous effect (coupling) changes the velocity profile during counter-current
flow (as compared to co-current flow), and therefore, the mobility of both phases is reduced
relative to co-current flow setting, where the relative permeabilities were measured. In previous
experimental studies of counter-current relative permeability (Lelievre, 1966; Manai, 1991), the
endpoint relative permeability of the non-wetting phase was invariably reduced. This observation
may be a direct result of the boundary conditions that can affect the pressure gradients and fluid
distributions in the experiment. It may also be due to new interface configurations in the pore
spaces that can increase the contact area and related viscous force. A consistent explanation of
this phenomenon has, to the best of our knowledge, not yet been presented.
Although the viscosity of the non-wetting phase in our experimental efforts is higher than
the viscosity of supercritical CO
2
, it can be shown that adequate representation of changes in
flow direction will be important for brine/CO
2
systems as well. From Navier-Stokes equations
for a single cylindrical pore with radius R, one can express the radial variation of the velocity in
the z direction as:

(3.9)

41

where v represents the fluid velocity in the bulk flow direction (z), µ is the fluid viscosity, r is
the radial direction and p is the fluid pressure. Equation (3.9) is commonly used in deriving
Hagen-Poiseuille equation and assumes steady state flow, that the radial component of the
velocity is zero and that the flow is axisymmetric. In two-phase flow settings, assuming that
brine is the strongly wetting phase (located in the interval

, where

is the radius of
fluid interfaces), it can be shown that the flow rates of non-wetting (Q
nw
) and wetting (Q
w
)
phases are given by (Bacri et al., 1990; Dullien and Dong, 1996):

, (3.10)
where the coefficients

are given by:

(3.11)

(3.12)

(3.13)
Similarly to Equation (3.1) that arises from volume averaging of porous media, Equation
(3.10) shows the dependency of the flow rate of each phase on the pressure gradients of both
phases. We note that the variable a in Equations (3.11) - (3.13) represents the square root of the
non-wetting phase saturation. In addition, we observe that only one of the coefficients (c
11
)
depends on both phase viscosities. Equations (3.10) - (3.13) allow us to investigate the impact of
the viscosity ratio on the coupling between the flowing phases. In other words, are the off-
diagonal (coupling) terms in Equation (3.10 significant in settings where the viscosity ratio is
high as is the case for the CO
2
/brine system? If that is the case for a capillary tube, we would
also expect to see a significant difference between co-current and counter-current relative
42

permeabilities in realistic porous materials in the context of CO
2
/brine systems. Figure 3.11
reports the variation of the coefficients (c
ij
), normalized with the maximum value of c
11

(corresponding to single phase flow of the non-wetting phase) as a function of the non-wetting
phase saturation. In Figure 3.11, the viscosity ratio is equal to that of the analog system used in
our experiments (M ~ 2). From Figure 3.11, we observe that the off-diagonal terms dominate
over c
11
for saturations of the non-wetting phase less than 0.2. This indicates that the flow of the
wetting phase facilitates (drives) the flow of the non-wetting phase at low saturations via the
viscous coupling. Similarly, we observe that the off-diagonal terms dominate over c
22
at non-
wetting phase saturations greater than 0.4. This indicates that the viscous coupling facilitates the
flow of the wetting phase at wetting phase saturations below 0.6. For this viscosity ratio, it is
hence clear that a change in the direction of the pressure gradient in one phase would result in
counter-current flow rates that differ significantly from the co-current ones for a fixed magnitude
of the pressure gradients.
43

Figure 3.11. Viscous coupling during two-phase iC8-brine flow in a cylindrical pore

Figure 3.12. Viscous coupling during two-phase CO
2
-brine flow in a cylindrical pore
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Non-wetting phase saturation
Normalized c
ij

c
11
c
22
c
12
= c
21
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.05
0.1
0.15
0.2
0.25
Non-wetting phase saturation
Normalized c
ij

c
11
c
22
c
12
= c
21
44

Figure 3.12 reports the variation of the coefficients (c
ij
) as a function of the non-wetting
phase saturation for a case where the viscosity ratio (M = 10) is representative of a CO
2
/brine
system at 20MPa, 65 C and a salinity of 70,000 ppm (µ
brine
~0.5 and µ
CO2
~0.05). We observe
that the diagonal coefficient, for the wetting phase (c
22
), as well as the off-diagonal coefficient
(c
12
) is reduced when compared to the analog system, while the non-wetting phase coefficient
(c
11
) is increased. However, the observations from the analog system (M~2) holds true: The
viscous coupling promoted (and dominates) the flow of the non-wetting phase for saturations
less than 0.1 and the coupling promotes the flow of the wetting phase for saturations less than
0.6. In the range of non-wetting phase saturations where the off-diagonal term is smaller than the
diagonal terms (0.1-0.4), the off-diagonal term still accounts for ~50% of the values of the
diagonal terms. Accordingly, the viscous coupling will have effect in the context of fractional
flow over the entire range of saturations even when the mobility ratio is high. Therefore, we
would also expect significant differences between co-current and counter-current flow in settings
that are relevant to CO
2
/brine systems.
During CO
2
injection into saline aquifers, the vertical migration of CO
2
occurs in a
counter-current flow setting where brine is displaced downwards and in the horizontal direction
(Javaheri and Jessen, 2011). The associated reduction of the CO
2
relative permeability during its
upward migration is not the same in all cells during numerical calculations and depends on the
phase velocities. Therefore, in order to represent the effect of counter-current flow in numerical
calculations, the coefficients of the transport equations (Darcy’s law) should depend on
velocities (and their relative directions), which have to be determined iteratively. A different
approach for incorporation the reduction in relative permeability during counter-current flow was
proposed by Javaheri and Jessen (2011), and includes the use of a separate set of relative
45

permeability for each phase in counter-current flows (or the use of relative permeability
multiplier). The relationship between relative permeability and phase velocities or the magnitude
of the relative permeability multiplier is not universal and has to be determined for each setting.
We have previously tested this idea via numerical calculations for CO
2
migration in saline
aquifers at different reservoir conditions and in nearly all cases we observe that counter-current
flow is the dominant mode of flow in the upward migration of CO
2,
particularly after injection
stops. The associated reduction of CO
2
relative permeability, over the co-current values,
increases the time-scale of its migration to the cap-rock. The increase in migration time can be
quite significant and accordingly, the effect must be included in estimates of storage capacity of
saline aquifers.
The reduction of relative permeability of both phases, observed in this work, is similar to
previous studies (Bourbiaux and Kalaydjian, 1990) and shows that the flow direction has a
significant impact on fluid mobility. In the context of CO
2
injection into saline aquifers, where
flow settings are different in vertical and horizontal directions (Javaheri and Jessen, 2011), the
application of a more sophisticated representation of relative permeability is warranted.
Therefore, we find it is necessary to use a displacement-dependent representation of relative
permeability in the simulation of CO
2
injection into saline aquifers to ensure that the migration
of CO
2
is represented with sufficient accuracy.

46

Chapter Four: Dynamic Mobility in the Brine/CO
2
System (IMPES Formulation)
1

4.1 Introduction
In the previous chapter we presented a study of relative permeability measurements and
modeling during counter-current flow in porous a media. In the study, relative permeability in
counter-current flow was observed to be less than relative permeability in co-current flow. In this
chapter we use the observation in the context of two-phase flow of CO
2
and brine and show that
an integration of the reduction in relative permeability as flow transitions between co-current and
counter-current settings can impact the predicted dynamics significantly.
We use a synthetic case study to show the impact of using dynamic relative permeability
(mobility) on residual trapping of CO
2
. In conventional simulation of multiphase flow in porous
media, the relative permeability of each phase depends only on the saturation (and saturation
history) of the existing phases. Here we consider two-phase flow of CO
2
and brine (salt
precipitation is not considered). This study is focused on a single injection well, but can be
applied in field-scale simulation of CO
2
injection into saline aquifers.
First, the modeling and simulation approach is presented. Then the results are reported and
the difference between including the dynamic mobility and ignoring it is highlighted. The
chapter is then completed by a discussion of the results and a chapter summary.

1
The results of this chapter is published in the International Journal of Greenhouse Gas Control:
Javaheri, M., Jessen, K., 2011. Integration of counter-current relative permeability in the simulation of CO
2
injection
into saline aquifers, International Journal of Greenhouse Gas Control, 5 (5), pp. 1272-1283.
doi:10.1016/j.ijggc.2011.05.015.

47

4.2 Modeling and Simulation
We use a two-dimensional aquifer model to study the impact of transitions between co- and
counter-current flow settings. The dimensions of the synthetic model are 2000 m by 300 m in the
horizontal and vertical directions, respectively. A volume of CO
2
accounting for 5% of the pore
volume is injected at a constant rate at the bottom of the aquifer over a 2400 day period. We
subsequently study the migration pattern for up to 15000 days. The top and bottom boundaries
are closed and the side boundaries are open. We use production wells to mimic the open
boundaries of a section of a larger aquifer. To do this, each cell at the vertical boundaries has a
completion that is operated at a constant bottom-hole pressure equal to the initial hydrostatic
pressure in the formation. The pressure at the top boundary is initially set at 20 MPa and the
temperature is set at 330 K. The properties of brine are assumed to be constant with values of
1030 kg/m
3
and 0.5 cp for density and viscosity, respectively. The density of CO
2
is calculated
from the Peng-Robinson equation of state, and the CO
2
viscosity is calculated from the
correlation of Fenghour et al. (1998).
The domain and its petro-physical properties are summarized in Figure 4.1. The domain
is represented by 41*41 cells (selected after an initial grid-refinement study). We study four
different permeability fields, with an average permeability of 125 mD in the horizontal direction.
The four permeability maps include one homogeneous field and three spatially heterogeneous
maps with different correlation lengths. The four settings used in this study are shown in Figure
4.2. Four different vertical-to-horizontal permeability ratios are considered for each aquifer
setting. We use an IMPES finite-volume approach to simulate the migration dynamics.
We include hysteresis in relative permeability of both phases as well as in capillary pressure. The
scanning curves of the relative permeabilities are assumed to be reversible between drainage and
48

imbibition processes. Killough’s hysteresis model (Killough, 1976) is used for both relative
permeability and capillary pressure as discussed below.

Figure 4.1. Aquifer dimensions and the location of injector

4.2.1 Hysteresis in Relative Permeability
In this work we use Land’s model (Land, 1968) with a trapping coefficient of 2 to account for
residual trapping. Hysteresis in relative permeability is based on the Killough’s model via sets of
bounding relative permeability curves; one set for co-current flow and another set for counter-
current flow. The endpoint relative permeability of CO
2
is reduced from 0.6 in co-current flows
to 0.4 in counter-current flows, a reduction of one-third that is consistent with experimental
observations (Bourbiaux and Kalaydjian, 1990; Bentsen and Manai, 1993; Manai, 1991). The
value of residual gas saturation is assumed to be the same for both co-current and counter-current
relative permeabilities during imbibition processes: There are no experimental observations
available in the literature to support whether the residual saturation depends on the mode of flow
300 m
2000 m
Injection
P
top
=20 MPa
T
top
=330 K
φ=0.25
k
h,avg
=125 md
k
v
/k
h
=0.1, 0.3, 0.5, 1
49

or not, and we assume, for simplicity, that it is independent of flow mode. The scanning curves
of the co-current and counter-current relative permeabilities of CO
2
start from the same
saturation on the drainage curves with different values of relative permeabilities (depending on
whether the flow is co-current or counter-current), and end at the same point at the residual
saturation of CO
2
corresponding to zero relative permeability. For the wetting phase, the
endpoints of the relative permeabilities of the two flow modes are the same (relative permeability
of the wetting phase is one when the saturation of the wetting phase is one, and zero at the
irreducible wetting phase saturation). We adjust the exponent of the Corey-type relative
permeability function from 2 for co-current flow to 3 for counter-current flow. The bounding
relative permeability curves that are used throughout our numerical calculations are shown in
Figure 4.3.
50

Figure 4.2. Logarithm of four permeability fields with different correlations in heterogeneity a)
homogeneous field, b) short cor
r
elation length, c) medium correlation length, and d) long
correlation length

The flowing fraction of gas is calculated from:

(4.1)
where S
f
is the flowing gas saturation, S
gc
is the critical gas saturation (assumed to be zero in this
work), S
g
is the average cell saturation, C is the Land trapping coefficient, and S
t,max
is the
maximum residual gas saturation in the block and is obtained from Land’s equation (Land,
1968).

a b
c d

5 10 15 20 25 30 35 40
5
10
15
20
25
30
35
40
1.8
2
2.2
2.4
2.6
2.8
3
51

Figure 4.3. Co-current and counter-current relative permeabilities of brine and CO
2

4.2.2 Hysteresis in Capillary Pressure
Killough’s model (1976) is used to include hysteresis in capillary pressure in our displacement
calculations. The drainage capillary pressure is obtained from:

(4.2)
where J is the Leverett J-function and S
g
is the CO
2
saturation. The imbibition capillary pressure
curve is obtained by scaling Equation (4.2) with respect to the residual saturation of CO
2
, to
ensure that the capillary pressure and relative permeabilities are consistent (identical residual
CO
2
saturation). Scanning curves for capillary pressure are obtained by interpolating between
drainage and imbibition curves using:

(4.3)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
S
w
k
r

I
Imbibition
Drainage
krw-co
krg-co
krw-counter
krg-counter
52

for transitions from drainage to imbibition, and:

(4.4)
for transitions from imbibition to drainage. When a flow reversal occurs, the relationship:

(4.5)
is used to find the unknown saturation parameters.

4.2.3 Flow Transitions
We track the individual phase velocities throughout the numerical calculations to identify
transitions between co- and counter-current flows. This is the procedure that is normally used to
ensure the selection of appropriate upwind directions for fluid mobilities in the pressure
equation. When a transition in flow direction is observed, we have implemented an option to
transition between the input co-current and counter-current relative permeability functions. In
our current implementation, we allow for an abrupt jump from e.g. the co-current relative
permeability curve to the counter-current relative permeability curve at a given saturation when a
change in flow direction occurs. The impact of our implementation on the numerical stability and
overall performance of the numerical calculations are discussed in more detail in section 4.4.

4.3 Results
We study the CO
2
migration dynamics in four different permeability fields derived from Ide et
al. (2007). In each case, three flow models are considered:
1. A co-current flow model using a single set of relative permeability functions.
2. A pure counter-current flow model using a single set of relative permeability functions.
53

3. A mixed-flow model allowing for transitions between co- and counter-current flows.
In the co-current flow model, the relative permeability for co-current flow is used in both
horizontal and vertical directions (via a single set of bounding relative permeability functions).
This model assumes that transitions from co-current to counter-current flow and vice versa do
not have any impact on the phase mobilities. In the counter-current flow model, the migration of
CO
2
in the vertical direction is forced to trace the counter-current relative permeability curves in
both directions. The impact of using co-current or counter-current relative permeabilities in the
horizontal direction has been studied as discussed later. Finally, in the mixed-flow model both
co-current and counter-current flows are considered: That is, co-current relative permeability is
used wherever the flow is co-current, and counter-current relative permeability is used wherever
the flow is counter-current. In the following sections we report and compare our findings for the
individual flow models and aquifer settings.

4.3.1 Homogeneous Aquifer
We start by considering a homogeneous aquifer with a constant permeability of 125 mD in the
horizontal direction. Four different vertical to horizontal permeability ratios of 0.1, 0.3, 0.5 and 1
are investigated. Figure 4.4 reports the fraction of the injected gas that is trapped by imbibition
processes as a function of time for the four different vertical to horizontal permeability ratios and
the three flow models. The main component of the residual trapping starts after the injection has
ended and imbibition processes initiate in the regions where brine was displaced by CO
2
. In
some cases, we observe that the fraction of the gas that is trapped increases and then decreases
during the injection period. This is caused by our calculation of the fractional trapping based on
the total gas that has been injected at a given time. Imbibition processes during the injection
54

period are observed to occur at some distance from the injection well where pressure gradients
are sufficiently low that gravity, capillary and viscous forces start to balance. As injection
continues, the ratio of the trapped gas to the total injected gas may decrease if the region under
imbibition remains unchanged.
For this aquifer setting, we observe that the integration of the counter-current relative
permeability model in our displacement calculations results in an increase in the residual
trapping. This is particularly evident as the vertical permeability increases. In contrast, for a
vertical to horizontal permeability ratio of 0.1, we observe only a modest difference between
calculations that utilize the co-current flow model and the mixed-flow model in terms of the
fraction of the gas that is ultimately trapped. However, when the permeability ratio (vertical to
horizontal) increases, the difference in the predicted trapping between the two flow models also
increases. We observe also, as expected, that as the vertical permeability increases, the ultimate
amount of gas that is trapped decreases. This is caused by a shift in the balance between
capillary, viscous and gravity forces towards a more gravity dominated flow setting. In terms of
the predicted residual entrapment, we observe that the use of a co-current flow model instead of
a mixed-flow model may underestimate the residual trapping by more than 20% (relative).

55

Figure 4.4. Residual entrapment of gas as a function of time in a homogeneous aquifer

After the injection ends, brine starts to imbibe into the plume. Therefore, a sharp increase
in the trapped gas saturation is observed at this time. The rate at which the gas is trapped
depends, in part, on the vertical permeability as discussed below. As the vertical permeability
increases from 12.5 mD to 125 mD, the resistance to imbibition in the vertical direction is
reduced (due to an increase in the vertical gas mobility) and the initial rate of residual trapping
increases. Higher vertical permeability also eases the vertical migration of the injected CO
2
and
the initial plume after injection is not as wide as when the vertical permeability is low. In
addition, more free gas accumulates at the top of the aquifer for higher vertical permeabilities,
and the final amount of trapped gas decreases.
The trapped gas saturation approaches a plateau when brine has invaded the entire initial
CO
2
plume. Any gas that has reached the top boundary at this point can be displaced further in
0 3000 6000 9000 12000 15000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=0.1

0 3000 6000 9000 12000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=0.3

0 3000 6000 9000 12000 15000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=0.5

0 3000 6000 9000 12000 15000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=1

co-current
counter-current
mixed
co-current
counter-current
mixed
co-current
counter-current
mixed
co-current
counter-current
mixed
56

the horizontal direction by capillary forces that gradually allow for brine to re-invade the top
portion of the aquifer. The reinvasion of brine into the upper part of the aquifer and the
associated horizontal displacement of CO
2
is a much slower process than the invasion of the
initial plume. Therefore, in settings with a relatively high vertical permeability, an abrupt rise is
observed in the fractional trapping of the injected gas, followed by a gradual increase until the
final plateau is reached. This general trend is evident from Figure 4.4 where the initial rate of
entrapment increases as the permeability ratio increases from 0.1 to 1.
By comparing the four panels in Figure 4.4, we observe that the initial rate of residual
entrapment after injection is predicted to be higher for the co-current model than for the counter-
current and mixed models. Furthermore, the difference in the initial rate of entrapment is
observed to depend on the vertical permeability: The lower the vertical permeability, the larger
the difference.
To explain the difference in the initial rate of entrapment, we report the fluid distributions
in terms of gas saturations after injection as predicted by the co-current model and the mixed
model in Figure 4.5. The rate of entrapment is controlled by two factors: a) the magnitude and
distribution of the gas saturation in the plume after injection and b) the resistance to imbibition of
brine into the plume. Based on conventional initial-residual considerations, a higher initial
saturation of CO
2
will result in a larger fraction of residual entrapment during an imbibition
process. The rate of imbibition into the plume is set by gravity and capillary forces as well as the
vertical permeability as discussed by Nattwongasem and Jessen (2009): The characteristic time
for vertical plume migration is proportional to the inverse of the vertical permeability for gravity
dominated flows.
57

For k
v
/k
h
= 0.1, we observe that the initial plume after injection is very similar for the two
flow models in terms of extent and saturation (Fig. 5). Accordingly, we expect the ultimate
residual entrapment to be similar for the two flow models (initial-residual consideration).
However, as the vertical fluid mobility is higher for the co-current model than for the mixed (and
counter-current) flow model, the initial rate of entrapment will be higher for the co-current
model.
As the vertical permeability is increased to k
v
/k
h
= 0.3, we observe differences in the initial
plumes as predicted by the co-current and mixed-flow models (Fig. 5). For the co-current model,
the plume is now shifted upwards relative to the plume of the mixed model with lower CO
2

saturations in the lower portion of the plume. Although the resistance to flow in the vertical
direction is still lower in the co-current model, less CO
2
will now be trapped per volume of
imbibition. This results in an initial rate of entrapment that is more similar to the mixed and
counter-current models.
A similar argument is valid for k
v
/k
h
of 0.5 and 1: The rate of imbibition is higher for the
co-current model (volume per time). However, the entrapment per volume is less than for the
mixed model (initial-residual argument) and the effective initial rate of entrapment as predicted
by co-current and mixed models becomes similar.

58

Figure 4.5. CO
2
saturation after injection for the co-current model (left) and the mixed model
(right)

To further investigate the differences in the residual entrapment of CO
2
observed from
Figure 4.4, we turn to the maps of the saturation distribution. Figure 4.6 reports the saturations
(left) and regions of counter-current flow (right) as predicted by the mixed-flow model at various
times for a permeability ratio of 0.3. After 200 days of injection, most of the region invaded by
CO
2
(except at the very tip of the plume and in a small region near the injector) is experiencing
counter-current flow. Regions with counter-current flow are shown in dark color. A similar
pattern is observed after 500 and 1000 days of injection. During the early time, only a very small
region in the bottom of the plume experiences imbibition. After 2400 days, when the injection is
stopped, the majority of the CO
2
plume is in a counter-current flow mode. Therefore, during the

5 10 15 20 25 30 35 40
10
20
30
40
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
k
v
/k
h
=0.1
k
v
/k
h
=0.3
k
v
/k
h
=0.5
k
v
/k
h
=1
59

injection period, the plume experiences primarily counter-current flow in its upward motion.
After the end of injection, brine starts to imbibe into the plume from the sides and bottom, and
this circulation initiates co-current flow at the bottom of the plume towards its center and
upwards. The same trend in development of counter-current flow during the plume migration is
observed for the heterogeneous aquifer settings that are discussed in the following sections.

Figure 4.6. CO
2
saturation at different times (left) and corresponding regions where counter-
current flow occurs (right) for a heterogeneous aquifer with short correlation length and
k
v
/k
h
=0.3; Regions with counter-current flow are shown in black. Injection stops at 2400 days
200 days
500 days
1000 days
2400 days
3000 days

5 10 15 20 25 30 35 40
10
20
30
40
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
60

4.3.2 Heterogeneous Aquifer (Short Correlation)
Next we consider a heterogeneous aquifer with a relatively short-ranged correlation in the
permeability distribution. A comparison of the residual gas vs. time as predicted by the three
flow models is shown in Figure 4.7. Similar to the homogeneous aquifer, the fraction of the
residual gas increases when counter-current relative permeability is included in the simulation of
the gas migration. The relative difference between the co-current and the mixed model is again
more significant as the vertical permeability increases. Residual trapping in the heterogeneous
aquifers is predicted to be more significant than for the homogeneous aquifer settings.

Figure 4.7. Residual entrapment of gas as a function of time in a heterogeneous aquifer (short
correlation length)

0 3000 6000 9000 12000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=0.1

0 3000 6000 9000 12000 15000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=0.3

0 3000 6000 9000 12000 15000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=0.5

0 3000 6000 9000 12000 15000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=1

co-current
counter-current
mixed
co-current
counter-current
mixed
co-current
counter-current
mixed
co-current
counter-current
mixed
61

To better illustrate this difference, we report the saturation of CO
2
as well as the residual
saturation of CO
2
after 15000 days in Figure 4.8 for a permeability ratio of 0.5. We observe that
after 15000 days, the injected CO2 is either trapped as a residual phase or mobile as a CO2-rich
phase at the top of the aquifer. Figure 4.8a shows the CO2 saturation after 15000 days as
predicted from the co-current model, and Figure 4.8c shows the CO2 saturation at the same time
as predicted by the mixed model. We observe that the plume is predicted to be larger by the
mixed model and that there is less free gas relative to what is predicted by the co-current model.
Figure 4.8b and Figure 4.8d show the corresponding trapped gas saturations as predicted from
the co-current and mixed models, respectively.

Figure 4.8. a) CO
2
saturation from co-current model, b) Trapped saturation from co-current
model, c) CO
2
saturation from mixed model, d) Trapped saturation from mixed model, all after
15000 days for a heterogeneous aquifer with medium correlation length and permeability ratio of
0.5

0
0.1
0.2
0.3
0.4
0.5
0.6
0
0.05
0.1
0.15
0.2
0.25
0.3

0
0.1
0.2
0.3
0.4
0.5
0.6
0
0.05
0.1
0.15
0.2
0.25
0.3
a
c
b
d
62

During injection, the plume is, again, almost entirely in a counter-current flow mode, and
the fluid mobility in the vertical direction is decreased. This causes more spreading of the plume
in the horizontal direction that, in turn, results in additional residual trapping. In addition, the
maximum gas saturation is higher in the mixed model due to the reduction of the vertical gas
mobility and, therefore, based on the Land’s model, the trapped gas increases in the mixed model
relative to the co-current flow model.
Residual trapping is more dominant near the injection well and below the top of the
aquifer, since the maximum saturation of CO
2
at those locations is higher than in the rest of the
plume. At the top of the aquifer, most of the CO
2
is still mobile and the level of residual trapping
is low. For this aquifer setting, the residual trapping is predicted to be higher in the mixed flow
model relative to the co-current model and the difference is observed to depend on the vertical to
horizontal permeability ratios. For permeability ratios of 0.5 and 1, the residual gas saturation is
increased by approximately 0.1 (16-18% relative change) when the mixed model is used in the
calculations over the results of the co-current model.

4.3.3 Heterogeneous Aquifer (Medium Correlation)
In a third example, we study a heterogeneous aquifer with a medium correlation in the
permeability distribution. The predicted displacement dynamics for this heterogeneous aquifer
setting are similar to the previous two cases. The development of counter-current flow is shown
in Figure 4.9. Due to the larger correlation in heterogeneity, the plume is less symmetrical as
compared to previous cases. However, during the injection period the plume is still controlled by
counter-current flow in its vertical migration, except for a small region near the injection well.
63

Figure 4.9. CO
2
saturation at different times (left) and corresponding regions where counter-
current flow occurs (right) for a heterogeneous aquifer with a medium correlation length and a
permeability ratio of 0.3; Counter-current flow occurs in regions with black color on the right
panels. Injection stops at 2400 days

200 days
500 days
1000 days
2400 days
3000 days

5 10 15 20 25 30 35 40
10
20
30
40
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
64

The velocity profiles reveal that after the injection is stopped, brine invades the center of
the plume mostly from the upper part of the plume. Therefore, a co-current flow (with counter-
current flow along the sides that allows more imbibition of brine into the plume) starts from the
bottom of the plume and invades more of the plume with time. A comparison between the
predicted entrapments of CO
2
for the three flow models and four permeability ratios are
presented in Figure 4.10. As the permeability ratio increases, the final entrapment of gas
decreases, the difference in the trapped gas between the co-current and mixed mode increases, a
sharp rise proportional to the permeability ratio is observed in the slope of the trapped fraction of
gas after injection, and a long tail in the development of the trapped gas for high vertical
permeabilities is noticed. In this example, the difference in the ultimate gas entrapment between
the co-current model and the mixed model is greater than 0.1 (~20% relative change) for
permeability ratios of 0.5 and 1.
65

Figure 4.10. Residual entrapment of gas as a function of time in a heterogeneous aquifer
(medium correlation length)

4.3.4 Heterogeneous Aquifer (Long Correlation)
In our last calculation example, a heterogeneous aquifer with a longer range of correlation in the
permeability distribution is investigated. For this case, we also focus on the magnitude of
counter-current flow in the horizontal direction. Saturation profiles as well as regions of the
plume where counter-current flow occurs in vertical and horizontal directions are shown in
Figure 4.11. Counter-current flow is dominantly observed in the vertical direction during the
injection period. This observation points to a potential problem when using a strictly counter-
current relative permeability model. It is observed that almost the entire plume experiences co-
current flow in the horizontal direction during the injection period, and that counter-current flow
in the horizontal direction only arise in a small region after the injection is stopped. Therefore,
0 3000 6000 9000 12000 15000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=0.1

0 3000 6000 9000 12000 15000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=0.3

0 3000 6000 9000 12000 15000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=0.5

0 3000 6000 9000 12000 15000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=1

co-current
counter-current
mixed
co-current
counter-current
mixed
co-current
counter-current
mixed
co-current
counter-current
mixed
66

using a single set of counter-current relative permeability functions reduces the mobility of the
fluids in the horizontal direction although the dominant flow mode is co-current. This may not be
critical for the vertical migration of the plume, but it is important in its lateral movement after
reaching the top of the aquifer. Therefore, it is recommended to use both sets of relative
permeability functions (mixed model). A comparison of the fraction of the injected CO
2
that is
trapped as a function of time is reported in Figure 4.12 for the four permeability ratios. In this
example, the difference between the fractions of the gas that is ultimately trapped is
underestimated more than 0.1 (17-21% relative change) for permeability ratios of 0.5 and 1 if the
impact of counter-current flow is not included in the displacement calculations.

67

Figure 4.11. Saturation distributions (left column), regions with counter-current flow in the
vertical direction (middle column) and regions with counter-current flow in the horizontal
direction (right column) for a heterogeneous aquifer with long correlation length and k
v
/k
h
=0.3

200 days
500 days
1000 days
2400 days
3000 days

5 10 15 20 25 30 35 40
5
10
15
20
25
30
35
40
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
68

Figure 4.12. Residual entrapment of gas as a function of time in a heterogeneous aquifer (long
correlation length)

4.4 Discussion of the Results
In the previous sections, we have demonstrated the importance of considering counter-current
relative permeability in the simulation of CO
2
injection into saline aquifers.
During the injection period, the leading edge of the plume is in a co-current flow mode,
while the rest of the plume experiences counter-current flow (except near the injection well).
When the plume reaches the top of the aquifer, nearly the entire region occupied by the plume is
under counter-current flow conditions. After injection stops, a co-current flow region starts to
develop from the bottom of the plume towards the top of the aquifer. This transition is caused by
0 3000 6000 9000 12000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=0.1

0 3000 6000 9000 12000 15000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=0.3

0 3000 6000 9000 12000 15000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=0.5

0 3000 6000 9000 12000 15000 15000
0
0.2
0.4
0.6
0.8
1
time (days)
fraction of gas trapped
k
v
/k
h
=1

co-current
counter-current
mixed
co-current
counter-current
mixed
co-current
counter-current
mixed
co-current
counter-current
mixed
69

a circular movement of brine that is displaced from the top of the aquifer by the less dense CO
2
-
rich phase.
Figure 4.13 summarizes our example calculations in terms of the relative difference in
fractional gas entrapment as predicted by the co-current flow model and the proposed mixed-
flow model. The relative differences in the range of 15-20% clearly demonstrate the need for
including the impact of counter-current flow in the simulation of CO
2
injection into saline
aquifers.

The current literature on counter-current relative permeability is largely limited to two-
phase oil/water systems. To our best knowledge, no such experimental observations are available
for the system of brine and CO
2
. Hence, there is a need for new experimental evidence with
emphasis on this important application.

Figure 4.13. Relative change in the calculated residual entrapment of CO
2
between the co-
current flow model and mixed-flow model
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
5
10
15
20
25
vertical to horizontal permeability
relative increase in trapped gas (%)

permeability field 1
permeability field 2
permeability field 3
permeability field 4
70

Our calculated results that are based on tracking the velocities of the flowing phases show
that the CO
2
plume is mainly in a counter-current flow mode during its vertical migration, and
therefore, it appears that the use of counter-current relative permeability in the vertical direction
is appropriate. In the horizontal direction, co-current flow is the predominant flow mode based
on the boundary conditions used.
From the results presented in this study, the use of the counter-current relative permeability
model may appear to be a good choice for representing the fluid mobilities during plume
migration. However, this observation cannot be generalized to settings where a part of the plume
is in a counter-current flow mode (rising plume) and another part is in a co-current flow mode
(CO
2
migrating horizontally below the cap-rock). The use of a counter-current relative
permeability model can in such settings lead to erroneous results, as the lateral movement of CO
2

at the top boundary will be underestimated. Accordingly, we recommend the use of a mixed-flow
model.
We have investigated other aquifer settings (boundary conditions) e.g. one side open and
one side closed and arrive at substantially identical results. It should be noted that if the
boundary conditions are changed, counter-current flow in horizontal direction can be more
significant and may not be ignored. This case is, however, readily handled by using a mixed-flow
model.
In addition, a separate set of numerical calculations with constant pressure boundaries
(constant pressure along the vertical column) was performed. Again, similar results were
obtained for the homogeneous and heterogeneous domains used in this study. This suggests that
the boundary conditions, especially when the boundaries are located far from the injection
71

location, do not have a significant impact on the dynamics of the plume migration and the rate of
residual entrapment.
We have neglected dissolution of CO
2
into brine in this study but do not expect that this
will change the observations presented in this paper significantly. Solubility of CO
2
in brine and
the associated (modest) increase in the brine density may play a more significant role at later
times when the plume spreads along the top of an aquifer and gravity currents may start to form.
An important issue is related to the transitions between the two sets of relative
permeability functions when the mode of flow is changed. We have used a simple rule for
transitioning between the curves that is based on the velocities of the two phases; whenever the
mode of flow changes, a discontinuous change in relative permeability is introduced. The
transitioning between relative permeability curves should not (and did not in our calculation
examples) introduce any numerical instability as similar discontinuities in the coefficients of the
pressure equation caused by phase transitions is routine in e.g. compositional simulation of gas
injection processes. In addition, a change in flow direction can only occur when the phase
potential passes through a value of zero.
The introduction of counter-current relative permeabilities (and transitions) in our flow
calculations had only a minor effect on the pressure solver in terms of the number of required
iterations.
It may be argued that the transition between the two sets of relative permeability functions
should follow a smooth path (similar to scanning curves in e.g. capillary pressure). However, as
noted above, any change in the flow direction of a given phase is caused by a variation in the
driving force that must pass through zero. Therefore, any discontinuous change in the value of
the phase relative permeability will not be problematic due to the vanishing driving force.
72

In our IMPES implementation, the flow directions/modes at each time step are recognized
from the previous time step. In adaptive-implicit or fully-implicit implementations, additional
caution may be warranted to ensure adequate stability/performance. We note, however, that this
study is a first attempt to incorporate the use of counter-current relative permeability and to
demonstrate the importance of identifying flow modes in the simulation of CO
2
injection into
saline aquifers. Additional experimental and theoretical studies are needed in this area to further
prove the validity of the concepts and its implementation in large-scale simulation.

4.5 Concluding Remarks
Based on the examples and analysis presented in this chapter, we arrive at the following
conclusions:
1. Counter-current flow and associated relative permeability is found to play an important role
in determining the migration dynamics of CO
2
that is injected into a saline aquifer, especially
when the average vertical permeability of the aquifer is high.
2. Counter-current flow can be a dominant mode of flow during the injection period and
influences both the plume evolution and the rate and amount of subsequent residual
entrapment.
3. Numerical calculations including transitions between co-current and counter-current relative
permeability functions are found to predict a higher residual saturation of CO
2
relative to
calculations with a single set of co-current relative permeability functions.
4. Simulation of CO
2
injection into saline aquifers based on a single set of relative permeability
functions may accordingly not be sufficiently accurate for predicting plume evolution and
long term fate of the injected CO
2
.
73

Chapter Five: Residual Entrapment during Simultaneous Injection of CO
2
and Brine

5.1 Introduction
CO
2
is commonly injected as a supercritical phase into saline aquifers in order to reduce the
injected volume. If injected in the bottom part of an aquifer (e.g. Kumar et al., 2005) the CO
2

will subsequently migrates towards the top due to buoyancy. The injected CO
2
can be trapped in
four different ways: residual, dissolution, hydrodynamic and mineral trapping. Mineral trapping
occurs in the very long term (hundreds to thousands of years) and can be neglected in the short to
medium time scales. Hydrodynamic trapping is different from other trapping mechanisms in that
the CO
2
is mobile and is trapped only because of a sealing barrier that hinders further upward
migration of CO
2
. This mobile CO
2
-rich phase can potentially leak to other formations wherever
the integrity of the cap-rock is weak or if pathways exist (e.g. abandoned wells) for CO
2
to
escape from the structural trap. Therefore, it is important to apply injection strategies to
minimize the fraction of the mobile CO
2
and to increasing the effective capacity of other
trapping mechanisms.
The reminder of this chapter is organized as follows: First we state the problem statement.
Then present our modeling approach including a discussion of relevant scaling analysis. We
proceed with a presentation and discussion of the results of our simulation study, and conclude
the chapter with a summary of our key findings.

74

5.1.1 Simultaneous injection of CO
2
and brine
Different strategies have been proposed to accelerate either dissolution or residual trapping of
CO
2
in aquifers. One such strategy includes co-injection of brine to accelerate CO
2
entrapment as
a residual phase (Hassanzadeh et al., 2009, Leonenko et al., 2006, Leonenko and Keith, 2008,
Anchliya and Ehlig-Economides, 2009, Akinnikawe et al., 2010). In these schemes, CO
2
is
injected in the bottom part of the aquifer, while brine is injected at or near the same location of
the injection well in the upper part of the aquifer. Brine can be injected in the same well as CO
2

or by using a separate well. In this setting, brine production is required from a well (or multiple
wells) at some distance to the injection well. The purpose of brine injection is to force the
injected CO
2
to travel further in the horizontal direction and to reduce the gravity override. This,
in turn, increases the volume of the aquifer that is contacted by the injected CO
2
. Brine injection
will enhance both dissolution and residual trapping of CO
2
, and can be continued after the
injection of CO
2
has ended.

5.2 Statement of the Problem
Kopp et al. (2009) have estimated the storage capacity of saline aquifers based on analysis of a
reservoir parameter database including a large number of oil reservoirs. By dimensional analysis
(via gravity and capillary numbers) they demonstrate that it is possible to classify aquifer settings
by the relevant driving forces, assuming that statistical geologic parameters of saline aquifers do
not differ significantly from oil and gas reservoirs. In addition, they demonstrate that the phase
relative permeabilities have a great impact on the storage capacity of a given aquifer.
In this study, we ignore the effect of dissolution trapping and focus entirely on residual
entrapment. Dissolution of CO
2
into the brine starts during the injection phase and may continue
75

into the very long term via gravity currents (Ennis-King and Paterson, 2003). Any engineering of
sequestration processes that aim at enhancing the residual trapping of CO
2
will also enhance
dissolution, as the magnitude of both depends on the sweep efficiency of the displacement
process (Shamshiri and Jafarpour, 2010). At time scales of the order of the injection period,
residual entrapment is more significant than dissolution (IPCC, 2005, Gasda et al., 2010).
However, dissolution trapping can be a significant mechanism over longer periods of time.
Dissolution trapping is significant especially when gravity currents are initiated in an aquifer
(Ennis-King and Paterson, 2003), and convective mixing increases the rate of mass transfer of
CO
2
into the brine. Such currents are, however, commonly initiated a long time after the
injection has ended (Hassanzadeh et al., 2006, Javaheri et al., 2009). At the relatively short time-
scales investigated in this work, trapping by dissolution may be neglected, to provide
conservative estimates of the trapped CO
2
. Residual CO
2
(and hydrodynamically-trapped CO
2
)
will eventually dissolve into brine.
In this work we aim at studying a wide range of relevant aquifer settings and utilize the
categories proposed by Kopp et al. (2009). Saline aquifers are categorized into five categories
(cold, warm, shallow, deep, and median) based on a statistical analysis of a huge number of
reservoirs (more than 1200 reservoirs). The properties of the brine and CO
2
(especially CO
2

properties) are different in each category due to differences in representative temperature and
pressure (Correlations of density and viscosity of CO
2
used for the four aquifer settings are listed
in Appendix 2). Therefore, the relative magnitudes of the driving forces that are acting on CO
2
in
the aquifer (gravity, viscous and capillary) are also different for each aquifer setting. The relative
magnitude of the driving forces on CO
2
during its migration from the injection well towards the
top of an aquifer plays a central role in the fraction of the CO
2
that is trapped. We compare the
76

significance of residual trapping for four aquifer settings (shallow aquifer excluded) and discuss
the observations in terms of on the gravity number (relative importance of gravity to viscous
forces) and the capillary number (capillary to viscous forces). The significance of relative
permeability on residual trapping is also addressed. It is addressed how CO
2
saturation and
residual entrapment of CO
2
can vary when the mixed flow model instead of the co-current model
is used.

5.3 Modeling Approach
In this section we present our modeling and simulation approach in terms of aquifer dimensions,
boundary conditions and the relevant saturation functions and fluid properties. To facilitate our
analysis, appropriate scaling groups are presented and discussed in the context of residual
entrapment.

5.3.1 Simulation model
The use of 3D models to simulate and analyze CO
2
injection into saline aquifers is accurate when
grid sizes are not large. As the grid sizes increase, accuracy of the results can be reduced
significantly. At the same time, if the number of the cells in a 3D model is increased sufficiently
to capture small-scale processes of importance to aquifer sequestration (Gasda et al., 2010),
simulation times may become prohibitively long. In this study we use a 2D domain with a
relatively high resolution. The results obtained from this 2D study may not be directly applicable
to 3D settings, but the analysis presented in this study can readily be extended to 3D settings.
We consider a 2D vertical cross section of an aquifer with a height of 200 m and a length
of 4000 m (see Figure 5.1). The lower and upper boundaries of the model are closed to flow.
77

CO
2
is injected from the lower-half of the left boundary and brine is injected from its upper-half.
The right boundary of the model is an open boundary that produces brine at a constant rate to
maintain the pressure stabilized. After an initial investigation of the sensitivity to grid
refinement, a computational grid of 100×50 cells was selected for all subsequent numerical
calculations. The initial pressure and temperature at the top boundary of each aquifer and the
average fluid properties are summarized in Table 5-1. Viscosity and density of the fluids listed in
Table 1 are calculated at the mean temperature and pressure of each aquifer setting.

Figure 5.1. Schematic of simulation model

200 m
4000 m
CO
2
Injection
φ=0.25
k
h,avg
=500 mD
Brine
Injection
Brine
Production
78

Table 5-1. Properties of the four aquifer settings (100,000 ppm)
Type
Depth
(m)
Pressure
(kPa)
Temperature
(°C)
ρ
b

(unsaturated)
(kg/m
3
)
µ
CO2

(cp)
ρ
CO2

(kg/m
3
)
IFT
(mN/m)
Cold 1524 15510 37.43 1069.36 0.0725 808.38 28.43
Deep 3495 35310 115.0 1030.11 0.0556 669.04 28.59
Median 1524 15470 55.13 1060.22 0.0548 683.40 28.65
Warm 1524 14950 104.49 1029.28 0.0285 339.55 32.85

We consider two permeability fields in this work: A homogeneous and a heterogeneous
field. The homogeneous setting has a horizontal permeability of 500 mD, while the vertical
permeability is varied in order to investigate a range of gravity numbers. Since the aquifer length
is four kilometres, horizontal permeability is selected somewhat high to facilitate the co-injection
of brine and CO
2
, and reflects that co-injection may not be suitable for low permeable aquifer
settings. The heterogeneous setting has an average permeability of 500 mD in the horizontal
direction. The heterogeneous field has a medium correlation length and is shown in Figure 5.2.
CO
2
is injected for 8 years and at the end of the injection period, ~ 15% of the pore volume
is occupied by CO
2
. Brine injection and production is initially considered during the CO
2

injection period. If brine injection and production is continued beyond the CO
2
injection phase,
CO
2
will eventually arrive at the brine production well. Simulations are continued until no
further changes in saturations occur.

79

Figure 5.2. Heterogeneous permeability field with a medium correlation length

Fluid properties (especially CO
2
properties) are greatly influenced by the differences in
pressure and temperature in the four aquifer settings considered in this work.

5.3.2 Scaling groups
In this work, we address the role of residual entrapment of CO
2
in the context of dimensionless
groups introduced by Zhou et al. (1994). Starting from continuity equation and the definition of
capillary pressure, Zhou et al. (1994) derived dimensionless groups that describe the two-
dimensional displacement of incompressible two-phase flow in porous media. The dimensionless
groups considering gravity, capillary and viscous forces are defined to determine the dominant
regime of flow in the relevant two-phase displacements:

(5.1)
log(k
x
)

10 20 30 40 50 60 70 80
10
20
30
40
50
2 2.2 2.4 2.6 2.8 3 3.2
80

(5.2)
where N
gv
and N
cv

are gravity and capillary numbers, respectively, representing the ratio of
gravity and capillary forces in the transverse direction to the viscous force in the horizontal
direction, u
t
is the total horizontal velocity, k
av
is the average horizontal permeability, and H and
L are characteristic lengths in the transverse and horizontal directions, respectively. p
c
*

represents the characteristic capillary pressure of the porous media and is defined by:

(5.3)
where S
wc
and S
gr
are irreducible water and residual gas saturations, respectively. The
interpretation of the dimensionless for the two-phase brine/CO
2
system and the definitions differ
slightly from the work of Zhou et al. (1994). We represent the capillary pressure as:

(5.4)
where J is the dimensionless capillary pressure (J-function), and R
v
represents a dimensionless
vertical permeability, defined as vertical permeability divided by average vertical permeability.
To further parameterize the relationship between the flow in horizontal and vertical directions we
use the shape factor defined as:

(5.5)
k
av
and k
ah
are average permeabilities in the vertical and horizontal directions, respectively.
The dimensional groups provide a criterion to determine the dominant driving force for
fluid movement and we demonstrate that CO
2
entrapment can be related to these dimensionless
groups. Over the range of aquifer pressures and temperature in a given aquifer, except for
regions near the wellbore where pressures are high, the density, and hence, compressibility of
81

CO
2
does not vary significantly, and the assumption of incompressible flow (used to arrive at the
dimensionless groups) seems to be reasonable.

5.3.3 Capillary pressure and relative permeability
We included hysteresis in both capillary pressure and relative permeability throughout our
numerical calculations and we note that it has been shown previously that hysteresis in capillary
pressure and relative permeability tends to reduce the dissolution of CO
2
into brine (Altundas et
al., 2010).
Hysteresis in capillary pressure is included via Killough’s model (Killough, 1976) where
scanning curves for capillary pressure are constructed by interpolation between bounding
drainage and imbibition capillary pressure curves. The residual saturation of each scanning curve
is calculated by Land’s formula (Land, 1968):

(5.6)
where S
gr
is the residual gas saturation and S
gi
is the initial gas saturation in the imbibitions
curve. C is the Land trapping coefficient and is a constant for all aquifer settings. The maximum
residual gas saturation is set to 0.3 (corresponding to C=2). As CO
2
saturation reduces from S
gi

to S
gr
along a scanning curve, the trapped CO
2
saturation increases from zero to S
gr
.
The drainage capillary pressure (in kPa) curve is calculated from:

(5.7)
with a maximum value of 36.8 kPa at the connate water saturation. Accordingly, the capillary
entry pressure is ~ 1 psi. It is assumed that the capillary pressure is the same for all aquifer
settings. This assumption may not be entirely realistic, and usually as the depth of the aquifer
82

increases, capillary increases because the size of pore throats (permeability) decreases. However,
because of the lack of data and also to reduce the number of independent parameters, we assume
that all aquifers have the same maximum capillary pressure.
Corey-type relative permeabilities are used in this study. Corey exponent for both phases is
set to 3. We use a connate water saturation of 0.25, and critical gas saturation of 0. It is assumed
that the scanning curves in relative permeability are reversible. We use the same relative
permeability functions for all aquifers settings. Again, as the depth of an aquifer setting
increases, the residual saturation of both phases and also the relative permeabilities may change
due to the changes in pore structure. This may even be true within a given aquifer. The relative
permeability functions are known to have a significant impact on the final residual entrapment
and as relative permeability decreases, depending on the vertical to horizontal permeability, the
fraction of residual CO
2
that is trapped increases (Javaheri and Jessen, 2011). However, in this
study we investigate the effect of the relative significance of driving forces on residual trapping
for different aquifers and assume that the bounding relative permeability curves are similar in all
aquifer settings.

5.4 Impact of Dimensionless Numbers on Residual Entrapment of CO
2

In this section, a reservoir simulator (CMG, 2010) is used for simulation. In the next sections an
in-house simulator is used.

5.4.1 CO
2
saturation and residual entrapment (co-current model)
For each aquifer setting, we performed four simulations: two for a homogeneous permeability
field, and two for the heterogeneous field. For each permeability field, three vertical
83

permeabilities are considered: 1, 0.1 and 0.01. According to the properties of CO
2
and brine the
mobility ratio is much higher than unity. Therefore, according to the definitions of the dominant
flow regime by Zhou et al. (1994), when the permeability ratio is 0.01, the system is viscous-
dominated, for the permeability ratio of 0.1 there is a balance between gravity and viscous
forces, and for the permeability ratio of 1, gravity is the dominant force. For all aquifer settings,
capillary force does not play a central role, because the capillary number is very small. Capillary
number for these settings is an order of magnitude less than gravity number.
The density of CO
2
in the warm aquifer setting is much lower than the density of CO
2
in
the other aquifer settings, and therefore, the density difference between brine and CO
2
is higher
for this setting. The viscosity of brine is also the lower for the warm aquifer. Therefore, the
gravity number for a warm aquifer (Equation (5.1 is larger compared to other aquifer settings for
the same total velocity. Gravity numbers are smaller for cold aquifers. We will demonstrate that
there is a significant difference between the four aquifer settings in terms of storage capacity for
fixed geological properties.
Figure 5.3 shows the distribution of CO
2
saturation (a) and trapped CO
2
saturation (b) for
the four aquifer settings in a homogeneous permeability field with k
v
/k
h
=0.01. All simulations
are continued for 200 years, because when the vertical permeability is low, it takes a long time
for the fluid redistribution to end. The fluid distribution is reported at three different times; after
8 years when injection and production stops, after 30 years when CO
2
is still in transition
between free and residual phases, and finally after 200 years. In all of the four aquifers settings,
CO
2
does not reach to the top boundary after 8 years. In this case, viscous forces are strong
enough to retard vertical migration of CO
2
to a later time (viscous-dominated flow). Saturation
profiles and trapped saturations of the four aquifers are fairly similar after 8 years and 200 years,
84

but differ at 30 years. While most of the CO
2
is trapped in the warm aquifer after 8 years, it is
mostly entirely in a mobile phase in the cold aquifer.
Figure 5.4a compares the amount of the injected CO
2
that is trapped as a function of time.
It is interesting to notice that ~ 10% of the injected CO
2
is trapped after the injection period in all
aquifer settings. The ultimate fraction of CO
2
that is trapped is similar for the four aquifer
settings while the dynamics of arriving at this level of entrapment differ. When the length of the
aquifer is used as a characteristic length in the horizontal direction in the definition of the gravity
number, average gravity numbers between 0.03 (for the cold aquifer) and 0.17 (for the warm
aquifer) are obtained, which confirms that viscous force dominate in all aquifer settings during
injection.
It can be shown from Zhou et al. (1994) that the flow potential gradient in the vertical
direction is a function of gravity and capillary numbers. Since the capillary number is much
higher than gravity number for these settings, it is the gravity number that controls the migration
of CO
2
in the vertical direction during/after injection. This is why we observe a significant
difference between the cold and warm aquifers in the transition from low to high levels of
residual entrapment. Figure 5.4b reports the fraction of CO
2
that is trapped as a function of time
for the heterogeneous aquifer. We observe that the trend is very similar to the homogeneous
aquifer setting.

85

Figure 5.3. a) Saturation profiles and b) trapped saturation profiles of CO
2
for the four aquifer
settings in the homogeneous field with a horizontal permeability of 500 mD and k
v
/k
h
=0.01
Cold Deep Median Warm
b) Trapped saturation of CO
2
a) Saturation of CO
2
8 years
30 years
200 years
8 years
30 years
200 years

0 0.2 0.4 0.6 0.8
0 0.1 0.2 0.3

0 0.2 0.4 0.6
86

Figure 5.4. Residual entrapment of CO
2
vs. time for the a) homogeneous aquifer, b)
heterogeneous aquifer with k
v
/k
h
=0.01

Next, we compare the saturation and trapped gas saturation in the four aquifer settings (for
both homogeneous and heterogeneous permeability fields) with k
v
/k
h
=0.1. These are reported in
Figure 5.5 for the homogeneous field and Figure 5.6 for the heterogeneous field. Unlike the
previous example calculations, there is a substantial difference in saturation profiles between the
four aquifer settings. After injection stops, the injected CO
2
has not reached to the top boundary
of the cold aquifer, while it has reached to the top boundary and has travelled the length of the
warm aquifer. The downward movement of CO
2
at the end of the warm aquifer is due to conning
0 20 40 60 80 100 120 140 160 180 200
0
0.2
0.4
0.6
0.8
fraction of gas trapped
0 20 40 60 80 100 120 140 160 180 200
0
0.2
0.4
0.6
0.8
time (years)
fraction of gas trapped

Cold
Deep
Median
Warm
a
b
87

near the brine production well. After 200 years, most of the CO
2
in the cold aquifer is trapped as
residual gas and only a small fraction of CO
2
is located at the top boundary of the aquifer. In
contrast, most of the CO
2
in the warm aquifer is accumulated as a free phase at the top of the
aquifer. A very small fraction of CO
2
is trapped in the warm aquifer after 200 years, and even,
some of the trapped CO
2
is due to conning when CO
2
is displaced downwards causing some
additional residual CO
2
. The saturation profiles of the cold aquifer in Figure 5.3 and Figure 5.5
are very similar. The average gravity number is about 0.3 (for the aquifer with k
v
/k
h
=0.3), which
shows that viscous forces are dominant in both cases. The median and deep aquifers have similar
features in terms of the distribution of saturation and trapped saturation. In addition, the
saturation profiles vary little after 30 years.
Figure 5.8 compares the saturation and trapped gas saturation in the four aquifer settings
for a homogeneous permeability field with k
v
/k
h
=1. In this gravity-dominated flow, the
segregation distance between CO
2
and brine is very small compared to the other two cases.
Saturation profiles of the four aquifers do not show significant differences, since all of them are
under gravity flow. Figure 5.9 shows the fraction of CO
2
that is trapped. The residual entrapment
of CO
2
does not vary significantly between the four aquifer settings and is a small fraction of the
injected CO
2
, compared with the other two permeability ratios. Except for the cold aquifer, the
segregation distance between the injected fluids is very close to the injectors.

88

Figure 5.5. a) Saturation profiles and b) trapped saturation profiles of CO
2
for the four aquifer
settings in the homogeneous field with a horizontal permeability of 500 mD and k
v
/k
h
=0.1

Cold Deep Median Warm
a) Saturation of CO
2
30 years
200 years
8 years
30 years
200 years

0 0.2 0.4 0.6 0.8
0 0.1 0.2 0.3
b) Trapped saturation of CO
2
8 years

0 0.2 0.4 0.6
89

Figure 5.6. a) Saturation profiles and b) trapped saturation profiles of CO
2
for the four aquifer
settings in the heterogeneous field with a horizontal permeability of 500 mD and k
v
/k
h
=0.1

Cold Deep Median Warm
30 years
200 years
8 years
30 years
200 years
8 years

0 0.2 0.4 0.6 0.8
0 0.1 0.2 0.3
a) Saturation of CO
2
b) Trapped saturation of CO
2

0 0.2 0.4 0.6
90

Figure 5.7. Residual entrapment of CO
2
vs. time for the a) homogeneous aquifer, b)
heterogeneous aquifer with k
v
/k
h
=0.1

0 20 40 60 80 100 120 140 160 180 200
0
0.2
0.4
0.6
0.8
fraction of gas trapped

0 20 40 60 80 100 120 140 160 180 200
0
0.2
0.4
0.6
0.8
time (years)
fraction of gas trapped
Cold
Deep
Median
Warm
a
b
91

Figure 5.8. a) Saturation profiles and b) trapped saturation profiles of CO
2
for the four aquifer
settings in the homogeneous field with a horizontal permeability of 500 mD and k
v
/k
h
=1

30 years
200 years
8 years
30 years
200 years
8 years
Cold Deep Median Warm
a) Saturation of CO
2

0 0.2 0.4 0.6 0.8
0 0.1 0.2 0.3
b) Trapped saturation of CO
2

0 0.2 0.4 0.6
92

Figure 5.9. Residual entrapment of CO
2
vs. time for the homogeneous aquifer with k
v
/k
h
=1

5.4.2 Gravity number and residual entrapment of CO
2

It was shown in the previous section that the vertical permeability plays an important role
in the migration of CO
2
in the aquifer. Based on Equations (5.1 and (5.2, gravity and capillary
numbers of the four aquifers are listed in Table 5-2. Fraction of the trapped gas is plotted as a
function of the gravity number in Figure 5.10. For the viscous-dominated flow characterized by
(Zhou et al., 1994):

(5.8)
and knowing that M >> 1 and N
cv
<< N
gv
in these cases, it is seen that the flow regime is
controlled by the gravity number. For N
gv
<< 1, the fraction of the trapped gas is independent of
the aquifer properties. In the gravity-dominated flow (N
gv
>> N
cv
and

), characterized
in these problems by N
gv
>> 1, fraction of the trapped gas is also independent of the aquifer
properties and is very low. It is only in the transition zone that the trapped gas depends on the
0 20 40 60 80 100 120 140 160 180 200
0
0.2
0.4
0.6
0.8
time (years)
fraction of gas trapped

Cold
Deep
Median
Warm
93

relative magnitudes of the gravity and viscous forces. In the transition zone, unlike the gravity-
and viscous-dominated regimes, the fraction of the trapped gas is significantly different between
the four aquifers. Capillary number is an order of magnitude less than gravity number and does
not play a significant role in this study (for 3D simulations it can be more important due to the
reduced velocity away from the injection well).

Table 5-2. Gravity and capillary numbers for the four aquifer settings
Type
Depth
(m)
Pressure
(kPa)
Temperature
(°C)
ρ
b

(unsaturated)
(kg/m
3
)
µ
CO2

(cp)
ρ
CO2

(kg/m
3
)
IFT
(mN/m)
N
gv
×
k
h
/k
v

N
cv
×
k
h
/k
v

Cold 1524 15510 37.43 1069.36 0.0725 808.38 28.43 2.74 0.59
Deep 3495 35310 115.0 1030.11 0.0556 669.04 28.59 9.72 1.51
Median 1524 15470 55.13 1060.22 0.0548 683.40 28.65 5.25 0.78
Warm 1524 14950 104.49 1029.28 0.0285 339.55 32.85 17.20 1.40

Figure 5.10. Gravity number and residual entrapment
10
-2
10
-1
10
0
10
1
10
2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
gravity number (N
gv
)
fraction of gas trapped
94

5.4.3 Segregation distance of brine and CO
2

The segregated distance from the wellbore depends strongly on the vertical permeability. It is
inversely proportional to the vertical permeability (Rossen, 2006), which is shown in Figure
5.11. The segregation distance in the warm aquifer is about an order of magnitude less than the
cold aquifer in the range of permeabilities studied.

Figure 5.11. Segregation distance vs. permeability ratio

5.4.4 Rate of brine injection
The impact of brine circulation rates on residual entrapment of CO
2
was considered in this study.
Three brine injection rates are considered: One case injects CO
2
only (zero brine rate) to observe
how effective brine can increase the sweep efficiency of the displacement. In the other two
examples, brine is injected at a volumetric rate equal to the CO
2
injection rate and at a twice the
rate of CO
2
. It was shown that brine injection has a significant impact on amounts and timing of
10
-2
10
-1
10
0
10
1
10
2
10
3
10
4
10
5
permeability ratio (k
v
/k
h
)
segregation distance (m)

Cold
Deep
Median
Warm
95

residual entrapment. This is consistent with another study (Hassanzadeh et al., 2009) that shows
brine injection increases the immobilization of CO
2
considerably. Brine circulation increases the
sweep efficiency of the displacement and results in less CO
2
accumulation at the top of the
aquifer.

5.5 Co-Current vs. Mixed Model
Transitions between co-current and counter-current flow and the impact on relative permeability
during migration of CO
2
in saline aquifers were addressed by Javaheri and Jessen (2011). In this
section we briefly demonstrate that a transition between flow modes and associated changes in
relative permeability increases the fraction of CO
2
that is trapped. In this example we utilized an
in-house simulator. The dimensions of the model are the same as the previous section. Two
relative permeability models are included in these displacement calculations: a co-current
relative permeability model that use a single set of relative permeability curves for both CO
2
and
brine, and a mixed model with two sets of relative permeability functions for each phase to
represent transitions between co-current and counter-current flow (relative permeability curves
are the same as Figure 3.3). Relative permeability in counter-current flow is decreased by 30%
from the corresponding co-current values. In the mixed model, the relative permeability is
decreased to the counter-current relative permeability wherever counter-current flow is observed.
We inject ~15% pore volume of CO
2
over 8 years into an aquifer with the same parameters as
shown in Figure 5.1.
The properties of each aquifer (initial pressure and temperature) are set and CO
2
density
and viscosity are correlated with pressure (from previous time step) and temperature which is
assumed to be fixed (see Appendix 2).
96

In the previous section we saw that a permeability ratio of 0.1 is in the transition from a
viscous regime to gravity regime with increased vertical permeability. In this section we use the
mixed model for this permeability ratio and show how different CO
2
saturation and trapped
saturations can be from the co-current mode. Figure 5.12 shows the comparison of the co-current
flow model and mixed flow model for the four aquifer settings. The difference between the two
is the least in the warm aquifer (gravity dominates the flow) and the most in the deep aquifer (in
the transition regime with N
gv
=0.97). In the other two aquifers viscous force is more than gravity
and there is a modest difference between the two regimes.

Figure 5.12. Trapped fraction of CO
2
in the co-current and mixed flow models with k
v
/k
h
=0.1
CO
2
saturation and trapped saturation profiles of the co-current and mixed models for the deep
aquifer with k
v
/k
h
=0.1 are shown in Figure 5.13. The profiles are shown at two times; end of
injection period (8 years) and end of simulation (60 years). The increases fraction of the trapped
CO
2
in the mixed model is related to a better sweep efficiency in this model.
0 10 20 30 40 50 60
0
0.2
0.4
0.6
time (years)
fraction of gas trapped
Cold

co-current
mixed
0 10 20 30 40 50 60
0
0.2
0.4
0.6
time (years)
fraction of gas trapped
Deep

co-current
mixed
0 10 20 30 40 50 60
0
0.2
0.4
0.6
time (years)
fraction of gas trapped
Median

co-current
mixed
0 10 20 30 40 50 60
0
0.2
0.4
0.6
time (years)
fraction of gas trapped
Warm

co-current
mixed
97

Figure 5.13. a) Saturation profiles and b) trapped saturation profiles of CO
2
in the co-current and
mixed modes at the end of injection and end of simulation for the homogeneous deep aquifer
with k
v
/k
h
=0.1
In the viscous-dominated regimes (k
v
/k
h
=0.01), the final gas entrapment is the same, but the
transient entrapment is different (Figure 5.14).
a) Co-current (upper panel) and mixed (lower panel) flow models (8 years)
b) Co-current (upper panel) and mixed (lower panel) flow models (60 years)
Saturation Trapped Saturation

0 0.2 0.4 0.6 0.8

0 0.05 0.1 0.15 0.2 0.25 0.3

0 0.2 0.4 0.6 0.8

0 0.1 0.2 0.3
98

Figure 5.14. Trapped fraction of CO
2
in the co-current and mixed flow models with k
v
/k
h
=0.01

5.5.1 Segregation distance of the fluids
The segregation distance of the CO
2
and brine in the mixed mode increases relative to the co-
current mode. In the transition regime, as Figure 5.13 shows, the increase in the segregation
distance is the most (60% increase). When one of the forces is very dominant, the difference in
the segregation between the two flow models is not significant.

5.5.2 Heterogeneous permeability fields
We investigated the profiles of CO
2
saturation for the deep aquifer (k
v
/k
h
=0.1) in a
heterogeneous permeability field with a large correlation length for both co-current and mixed
flow models. Saturation profiles of CO
2
are plotted in Figure 5.15 at two different times; end of
injection (8 years) and end of simulation (60 years).
0 10 20 30 40 50 60
0
0.2
0.4
0.6
0.8
1
time (years)
fraction of gas trapped
Cold

co-current
mixed
0 10 20 30 40 50 60
0
0.2
0.4
0.6
0.8
1
time (years)
fraction of gas trapped
Deep

co-current
mixed
0 10 20 30 40 50 60
0
0.2
0.4
0.6
0.8
1
time (years)
fraction of gas trapped
Median

co-current
mixed
0 10 20 30 40 50 60
0
0.2
0.4
0.6
0.8
1
time (years)
fraction of gas trapped
Warm

co-current
mixed
99

Figure 5.15. a) Saturation profiles and b) trapped saturation profiles of CO
2
in the co-current and
mixed modes at the end of injection and end of simulation for the heterogeneous deep aquifer
with k
v
/k
h
=0.1
a) Co-current (upper panel) and mixed (lower panel) flow models (8 years)
b) Co-current (upper panel) and mixed (lower panel) flow models (60 years)
Saturation Trapped Saturation

0 0.2 0.4 0.6 0.8

0 0.1 0.2 0.3

0 0.2 0.4 0.6 0.8

0 0.1 0.2 0.3
100

In the mixed model (Figure 5.15a lower left panel) there is a region (in the lower part of
the field near the injection source) beneath a low-permeable layer that has a high CO
2
saturation.
This feature is not seen in the co-current model. These zones can occur in settings with high
degrees of heterogeneity, but does not happen in a homogeneous setting.
Figure 5.16 shows the comparison of the residual entrapment of CO
2
in the homogeneous
and heterogeneous fields for the co-current and mixed modes. Layers of low permeability in the
heterogeneous field contribute to a higher sweep efficiency, and hence, more entrapment of CO
2
.

Figure 5.16. Residual entrapment of CO
2
in the homogeneous and heterogeneous deep aquifer
(k
v
/k
h
=0.1) with the co-current and mixed flow models

0 10 20 30 40 50 60
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
time (years)
fraction of gas trapped

co-current, homogeneous
co-current, heterogeneous
mixed, homogeneous
mixed, heterogeneous
101

5.6 Concluding Remarks
Based on the results and analysis presented in the previous sections, we arrive at the following
conclusions:
1. Brine injection increases the residual trapping of CO
2
significantly in all aquifer settings.
Brine injection retards the vertical migration of CO
2
by decreasing the gravity number,
making the gravity force less significant relative to the viscous force.
2. There is a significant difference between the four investigated aquifer settings in terms of
the significance of residual entrapment and the storage capacity. Cold aquifers provide
for the largest storage capacity, whereas warm aquifers have the least storage capacity.
The storage capacity of median and deep aquifers is between the warm and cold aquifers
and they behave similarly. Even for cases where the final fraction of the injected CO
2
that
is trapped is the same in all aquifers, the transient behavior is significantly different.
3. Residual entrapment of CO
2
depends on the gravity number and decreases when gravity
number increases.
4. Inclusion of a counter-current relative permeability model in the simulation of CO
2

injection in saline aquifers increases the estimates of the trapped saturation. When brine
is co-injected with CO
2
, counter-current flow is still a dominant flow mode in the vertical
direction while flow in the horizontal direction is predominantly co-current. The mixed
model is mostly important in the transition regime (between gravity and viscous forces).

102

Chapter Six: Transitions between Co-Current and Counter-Current Mobility in Two-
Phase Flows

6.1 Introduction
In the previous chapters, we addressed the difference between co-current and counter-current
relative permeabilities and their impact on simulation results of the CO
2
injection into saline
aquifers. It was shown that the vertical migration of CO
2
in the aquifers is dominantly in the
counter-current mode. Therefore, the challenge is how to define relative permeability at any
saturation and flow direction given the input co-current and counter-current relative
permeabilities.
In this chapter we introduce an approach to transition between co-current and counter-
current flows in the numerical simulation of CO
2
injection into saline aquifers. The difference
between co- and counter-current flow and its impact on CO
2
migration and entrapment was
addressed in Chapter four. In Chapter two we applied the dependency of flow directions on
relative permeability and demonstrated that a reduction in relative permeability functions due to
transitions from co-current to counter-current flow (in the vertical direction) increases the
residual entrapment of CO
2
and retards its upward migration. An IMPES scheme was used in all
previous numerical calculations. In this chapter, we use the fully-implicit black oil formulation
implemented in the GPRS code (made available from Stanford University). We focus on the
transition of flow from co-current to counter-current (and vice versa) and the subsequent changes
in fluid mobilities. Numerical stability and differences between the new model and the co-current
model are demonstrated based on two calculation examples.

103

6.2 A New Mobility Model
In this section we propose a method to calculate the phase mobilities in each cell by considering
the flow direction. In order to consider flow direction in fully implicit numerical calculations,
two methodologies can be used. In the first approach Equation (3.1) can be used in the model
formulation and a full tensorial relative permeability representation must be employed. In this
approach the dependency of the off-diagonal terms with repsect to saturation must be
determined/predicted. Appropriate changes in the Jacobian matrix must accordingly be
introduced. While the adjustment of the the Jacobian may be tedious, the measurement of the
tensorial relative permeability will represent a significantly larger effort that may not be feasible.
Another methodology is to use the same transport equations (Darcy’s equations), and track
co-current and counter-current flows. This way, existing simulators can easily be adjusted to
allow for a representation of transitions between co-current and counter-current flow: As the
saturation functions are different in co-current and counter-current flows, an additional set of
bounding saturation functions representing counter-current flow is needed. In this work, we
focus entirely on the impact of flow direction on relative permeability and defer its impact in the
context of capillary pressure for later studies. Accordingly, relative permeability functions for
pure counter-current flow are needed as input to flow simulations. Figure 6.1 shows a
representative set relative permeability functions for co- and counter-current flows that mimics
the experimental observations of Lelievre (1966) and Bentsen and Manai (1993). One then has to
define the transition between co-current and counter-current flow as observed by tracing the
phase velocities: In each time step, the velocity of each phase at the center of each cell is
calculated from the interface velocities. In the case of co-current flow, co-current relative
permeability for each phase is used. If counter-current flow is observed in a cell, counter-current
104

relative permeability is used for each phase in that cell. The transition from co-current to
counter-current can be implemented in two ways:
1. Introduction of a discrete (abrupt) transition in relative permeability from co-current to
counter-current at a fixed saturation.
2. Introduction of an interpolation scheme that provides for a continuous transition between
co-current and counter-current values at a fixed saturation based on the magnitude of the
phase velocities.
In this work, we have tested both transition approaches and compare their behavior with the
traditional co-current model below.

Figure 6.1. Co-current and counter-current relative permeability curves

The capillary pressure curve that is used in the simulation studies of this chapter is plotted in
Figure 6.2. The hysteresis in capillary pressure is ignored in this chapter.
105

Figure 6.2. Capillary pressure function used in the simulation studies of this chapter
If we assume a continum model for the fluids in porous media, the velocity and shear stress
should be equal in the interface of the two fluids. These constratints determine the velocity
profile and the associated relative permeability values at the pore scale. Therefore, an
interpolation scheme between co-current and counter-current relative permeability should
consider both velocity and viscosity of the two fluids. In brine/CO
2
systems, the viscosity of both
brine and CO
2
do not change significantly over the pressure and temperature ranges of the
formation and we assume that the interpolation scheme does not depend on viscosity.
Accordingly, the relative permeability is interpolated between co-current and cunter-current
relative permeabilities based on the velocity of both phases reconstructed at the cell centers. If
co-current flow is observed in a given cell, the input co-current relative permeability functions
(including hysteresis) are used. However, when a transition from co-current to counter-current
flow is observed in a given cell, we employ the following interpolation scheme to calculate the
phase relative permeabilities:
106

(6.1)

(6.2)
where F is kept equal to unity for co-current flow and is calculated from the following equation
for counter-current flow:

(6.3)
Here, v
min
and v
max
are the minimum and maximum values of the vertical component of the
velocity of the two phases in a cell, respectively. The coefficient n is set to 1 in this section
corresponding to a linear interpolation between the co-current and counter-current relative
permeability. The interpolation parameter (F) is accordingly always in the interval of 0 to 1.
When the phase relative permeability is interpolated between the co-current and counter-current
values, Equation (6.3), the elements of the Jacobian must be adjusted accordingly. A simpler
approach that is also tested in this work is to use the discrete transition where F is set to 1 for
cells with co-current flow and set to 0 for cells with counter-current flow.

6.3 Calculation Examples
In this section we present two example calculations where CO
2
is injected into a saline aquifer.
The example calculations are performed with three different mobility models:
a) co-current relative permeabilities are used throughout the calculations
b) interpolation of relative permeability between co-current and counter-current values
(mixed/interpolation)
c) transition (discrete) of relative permeability between co-current and counter-current
values (mixed/discrete)
107

The mixed models (b+c) use both co-current and counter-current relative permeability curves via
the interpolation scheme given in Equation (6.3) or via a discrete and abrupt transition as
discussed above. The interpolation in the mixed model can be either linear or quadratic.
In the first calculation example an aquifer is represented by a 2D homogeneous
computational domain with specifications listed in Table 6-1. Brine and CO
2
properties (listed in
Table 6-2 and Table 6-3, respectively) are represented via black-oil tables and correspond to a
formation temperature of 60 °C and a salinity of 50000 ppm. Figure 6.1 shows the input relative
permeability functions for CO
2
and brine. We use Killough’s model to represent hysteresis in the
non-wetting phase relative permeability and Land’s model (Land, 1968) is used to calculate the
trapped saturation of CO
2
with a Land’s trapping coefficient of 2.5. CO
2
is injected from the
lower part of the left edge of the aquifer (lower 200 ft) and we mimic an open boundary on the
right edge of the domain by placing production wells operated at a constant bottomhole pressure
equal to the hydrostatic pressure of the formation. The pressure at the top of the aquifer is
initially set to 3000 psi (~20.7 kPa). CO
2
injected for a period of 2500 days after which, we
monitor the migration and immobilization of the injected CO
2
until 10000 days.
In the second calculation example, CO
2
is injected into a 3D heterogeneous aquifer
model. The average horizontal permeability is 90 mD and the average vertical permeability is 9
mD. Additional aquifer properties are listed in Table 6-4. We use the same fluid properties and
saturation functions as in the first calculation example. To mimic open boundaries, production
wells were placed along the four edges of the aquifer and operated at a constant bottomhole
pressure equal to the hydrostatic pressure of the aquifer. The model was initialized with a
pressure of 3000 psi at the top of the aquifer and CO
2
is injected for a period of 2500 days in the
108

lower portion (bottom three layers) of the reservoir at the center of the domain. The subsequent
migration and entrapment of CO
2
was then simulated till 15000 days.

Table 6-1. Specifications for the 2D aquifer model
L
x
(ft) 5000
L
z
(ft) 500
N
x
100
N
z
40
k
h
(mD) 50
k
v
(mD) 20
Porosity 0.25
Temperature (°C) 60
CO
2
injection rate 0.015 Mt/yr/100 ft

Table 6-2. CO
2
properties
P (psi) B
g
(Rb/STB)  (cp)
308.00 0.052674 0.016701
616.01 0.024032 0.016846
924.01 0.014297 0.017059
1232.01 0.009234 0.017362
1540.01 0.006017 0.017784
1848.02 0.004036 0.018372
2156.02 0.003214 0.019196
2464.02 0.002868 0.02038
2772.02 0.002675 0.022136
3080.03 0.002548 0.024823
3388.03 0.002455 0.028851
3696.03 0.002383 0.034048
4004.03 0.002324 0.039317
4312.04 0.002275 0.043918
4620.04 0.002233 0.047813
4928.04 0.002196 0.051159
5236.05 0.002164 0.054099
5544.05 0.002135 0.056732
109

Table 6-3. Brine properties
P (psi) B
o
(Rb/STB)  (cp)
308.00 1.017198 0.513568
616.01 1.016314 0.514148
924.01 1.01544 0.514728
1232.01 1.014574 0.515308
1540.01 1.013717 0.515888
1848.02 1.012869 0.516468
2156.02 1.01203 0.517048
2464.02 1.0112 0.517628
2772.02 1.010379 0.518208
3080.03 1.009566 0.518788
3388.03 1.008763 0.519368
3696.03 1.007968 0.519948
4004.03 1.007183 0.520529
4312.04 1.006406 0.521109
4620.04 1.005638 0.521689
4928.04 1.004879 0.522269
5236.05 1.004129 0.522849
5544.05 1.003388 0.523429

Table 6-4. Specifications for the 3D aquifer model
L
x
(ft) 5000
L
y
(ft) 5000
L
z
(ft) 500
N
x
50
N
y
50
N
z
40
k
h,avg
(mD) 90
k
v,avg
(mD) 9
Porosity 0.25
Temperature (°C) 60
CO
2
injection rate 0.26 Mt/yr
110

6.4 Calculation Results
6.4.1 Example 1
First we consider the 2D calculation example. Figure 6.3 reports the CO
2
saturation distribution
as predicted by the three mobility models at three different times. After 2500 days, when
injection stops, most of the CO
2
has migrated to the top of the aquifer in the co-current model,
while a significant fraction of the CO
2
is still migrating upwards in the mixed models. In order to
highlight the differences in CO
2
saturation distribution as predicted by the three models more
clearly, we report the differences in CO
2
saturation distribution between the co-current model
and the two mixed models in Figure 6.4. From Figure 6.3 and Figure 6.4, we observe that the
injected CO
2
is predicted to travel further into the formation prior to full phase segregation by the
mixed models relative to the co-current model. In particular, the mixed/discrete model predicts
migration dynamics that differs significantly from the co-current model: In some portions of the
aquifer the saturation difference between the models are as high as 0.25. Furthermore, the lateral
spreading of CO
2
below the caprock is more pronounced in the co-current model compared to
the mixed models: The injected CO
2
is predicted to travel approximately twice the distance under
the caprock as compared to the mixed/discrete model after 2500 days. This causes more CO
2
to
be trapped in the mixed models during the post-injection imbibition processes as shown in Figure
6.5. Figure 6.5 also demonstrates that the co-current and mixed models predict entrapment of
CO
2
in very different locations of the aquifer.
111

Figure 6.3. CO
2
saturation distribution as predicted by the co-current (left column), linear
interpolation mixed (middle column) and discrete mixed model (right column) at 2500 days (top
row), 4000 days (middle row), and 10000 days (bottom row)

0 0.4 0.78
112

Figure 6.4. The difference in CO
2
saturation between the co-current and the mixed models (left:
linear interpolation; right: discrete) at 2500 days (top row), 4000 days (middle row), and 10000
days (bottom row)

-0.25 0 0.25
S
co-current
– S
mixed (interpolation)
S
co-current
– S
mixed (discrete)
113

Figure 6.5. Trapped CO
2
saturation profiles in the three models at 2500 days (top row), 4000
days (middle row), and 10000 days (bottom row)

Since imbibition processes are predicted to initiate faster in the co-current model (faster
segregation), the rate of entrapment of CO
2
is initially faster in this model. However, the final
entrapment of CO
2
in the mixed models are higher due to longer travel distance before full phase
segregation as observed in Figure 6.3 (less free gas under the cap rock). Figure 6.6 reports the
fraction of the CO
2
that is trapped versus time. The final entrapment of CO
2
is quite different
between the three models: CO
2
entrapment predicted by the mixed/discrete and the
mixed/interpolated models are approximately 50% and 25% higher than what is predicted by the
co-current model. We observe also that the saturation in the plume as predicted by the
mixed/interpolated model falls in between the co-current model and the mixed/discrete model.

0 0.15 0.3
114

Figure 6.6. Fraction of the injected CO
2
that is trapped as a function of time

Additional information can be extracted by tracking the relative permeability of selected
cells during the example calculations. Figure 6.7 reports the observed relationship between the
CO
2
relative permeability and CO
2
saturation for the three mobility models in a specific cell
located 75m from the top of the domain and approximately half way along length of the domain.
For the co-current mobility model, the drainage co-current relative permeability is initially traced
followed by an imbibition scanning curve that is obtained from the maximum CO
2
saturation and
Land’s trapping coefficient. For the mixed/discrete model we observe an early departure from
the co-current drainage curve to the counter-current drainage curve. As the CO
2
saturation
increases, we see a single transition back to the co-current drainage curve before the counter-
current drainage curve is again traced until the maximum saturation is reached and the cell traces
115

an imbibition scanning curve until complete entrapment has occurred. CO
2
relative permeability
observed from the mixed/interpolated model is, in general, smooth and falls between the co-
current and mixed/discrete models with the exception of late in the imbibition process where the
mixed/interpolated relative permeability coincides with the co-current imbibition scanning curve.

Figure 6.7. Relative permeability variations in a cell (#546) located 75m from the top of the
domain and half way along the length of the aquifer as observed in the three mobility models

Table 6-5 summarizes the computational aspects of the first example in terms of CPU time,
average time step size and the total number of Newton iteration. Despite the additional
nonlinearities that are introduced by the mixed mobility models, we observe only marginal
differences between the traditional co-current representation and the proposed mixed models.
The marginal changes that are observed in the total number of Newton iterations suggests that
the proposed models for phase mobility including transitions between co-current and counter-
116

current flow, do not introduce any numerical difficulties relative to the standard treatment of
relative permeability in fully implicit simulations.

Table 6-5. Numerical comparison of the three mobility models for Example 1

co-current
mixed
(interpolation)
mixed
(discrete)
Simulation time
(sec)
869 794 823
Number of time
steps
112 117 110
Average time step
size (days)
89.3 85.5 90.9
Newton iterations 1007 1011 988

6.4.2 Example 2
Next, we consider displacement calculations in 3D aquifer settings based on the three mobility
models. Figure 6.8 compares the CO
2
saturation in the top layer as predicted by the three
mobility models at different times. From Figure 6.8, we observe that the footprint of the CO
2
plume at the top of the formation differs significantly between the three models. The largest
difference in the footprint is observed between the co-current model and the mixed/discrete
model. This is due to the more aggressive reduction in the vertical mobility of the CO
2

introduced by the mixed/discrete model as flow transitions between co-current and counter-
current settings. The smaller footprint at the top of the aquifers (reduced vertical mobility) is
balanced by a larger lateral extend (horizontal mobility) of the CO
2
plume as shown in Figure 6.9
that reports the plume evolution in the 8
th
layer from the top of the aquifer.
117

Figure 6.8. CO
2
saturation profile in the top layer in the co-current (left column), linear
interpolated mixed (middle column) and discrete mixed (right column) models at 2500 days (top
row), 4000 days (middle row), and 10000 days (bottom row)

Figure 6.9. CO
2
saturation profile in the three models in the 8th layer (from the top) in the three
models at 2500 days (top row), 4000 days (middle row), and 10000 days (bottom row)

0 0.2 0.4 0.6

0 0.2 0.4 0.6
118

Figure 6.10. Trapped CO
2
saturation as predicted by the three mobility models; left column: co-
current; middle column: mixed (linear interpolation); right column: mixed (discrete), at 3000
days (1
st
row), 4000 days (2
nd
row), 6000 days (3
rd
row) and 12000 days (4
th
row)

The difference in the vertical mobility of CO
2
will impact the rate of entrapment as a
reduced vertical mobility will obstruct the imbibition processes at the bottom of the plume after
injection has ended. This is clearly demonstrated in Figure 6.10 that shows the distribution of
entrapped CO
2
as a function of time for the three mobility models. However, the ultimate amount
of CO
2
that is entrapped is expected to be larger for the mixed models where less CO
2
is allowed
to reach the top of the aquifer at early times. Figure 6.11 reports the amount of CO
2
that is
trapped as a function of time as predicted by the three models. Similar to the 2D model, we
observe significant differences between the behaviors of the three models: The co-current and
mixed/discrete models have large differences in the dynamics of CO
2
entrapment, while the

0 0.1 0.2 0.3

0 0.1 0.2 0.3
119

mixed/interpolated model falls in between. In contrast to the 2D model, the main difference in
entrapment dynamics for this example is the rate: The difference in time to reach 50%
entrapment, as predicted from the co-current model and the mixed/discrete model, is
approximately 2000 days. However, after 15000 days we observe only moderate differences in
the entrapped amount as predicted by all three mobility models.

Figure 6.11. Fraction of the injected CO
2
that is trapped as a function of time

Figure 6.12 demonstrates in more detail the difference in the migration dynamics between
the three mobility models in terms of the CO
2
relative permeability versus CO
2
saturation for a
given cell located above the injection zone. During the injection period the cell is primarily in
co-current flow settings except for a few transitions. After the injection is stopped (at the
maximum CO
2
saturation for the given cell) the flow transitions to counter-current settings until
~3000 days where it transitions back to co-current settings. The mixed/interpolated model shows
120

a more smooth transition pattern as compared to the mixed/discrete model. Since Land’s
coefficient is assumed to be the same in the three mobility models and the maximum CO
2

saturation predicted to be similar in the three models, the residual CO
2
saturation will also be
similar in all models for the given cell.

Figure 6.12. Relative permeability variations in a cell (#16273) located just above the injection
zone as observed in the three mobility models

Table 6-6 reports the CPU time, average time step size and the total number of Newton
iteration in the three models of the second example. The number of the time steps is the same in
the three examples. The total simulation time and Newton iteration are slightly higher in the
mixed models compared to the co-current model. However, they do not introduce any difficulties
in the numerical simulation.

121

Table 6-6. Numerical comparison of the three mobility models for Example 2

co-
current
mixed
(interpolation)
mixed
(discrete)
Simulation
time (sec)
666
9
6896 6902
Number of
time steps
166 166 166
Average time
step size (days)
90.4 90.4 90.4
Newton
iterations
156
6
1586 1583

6.4.3 Effect of Vertical Permeability
In the previous section, the results of the two examples were presented. In both examples,
the vertical permeability was 0.1 of the horizontal permeability. In this section we investigate the
changes of the vertical permeability relative to horizontal permeability, and its effect on the
migration of CO
2
in the aquifer in the counter-current model.
We revised Example 2 with various k
v
/k
h
ratio. The relative difference of the trapped CO
2

between the co-current and counter-current models (linear and squared interpolations) are plotted
in Figure 6.13. It is clear that higher vertical permeability increases the difference in the trapped
CO
2
between the co-current and mixed models.
122

Figure 6.13. Relative difference in the fraction of CO
2
that is trapped between the co-current and
mixed models after 20000 days

6.4.4 Interpolation Scheme
In the previous section we used a linear interpolation between the co-current and counter-
current relative permeability curves in the counter-current flow. In this section, we use another
interpolation formula (quadratic) and examine how the results will change with a modification in
the interpolation parameter F. The power n in Equation (6.3) is set to 0.5 for the squared
interpolation in this section. k
v
/k
h
=1 for the examples in this section.
Figure 6.14 compares CO
2
saturation in the top layer of the aquifer in the co-current and
mixed models. A decrease in the power n shifts the interpolation parameter towards the counter-
current curve and CO
2
saturation profile becomes closer to the discrete mixed model. This is
clearer in Figure 6.15 where fraction of the trapped CO
2
that is injected is plotted with time for
123

the co-current and mixed models. Comparing Figure 6.15 and Figure 6.11, it is concluded that an
increase in the vertical permeability increases the contrast between the co-current and mixed
models.

Figure 6.14. CO
2
saturation profile in the top layer in the co-current (left column), linear
interpolated mixed (middle column) and quadratic mixed (right column) models at 2500 days
(top row), 4000 days (middle row), and 10000 days (bottom row)

0 0.35 0.7
124

Figure 6.15. Fraction of the CO
2
trapped in the co-current and mixed models

6.5 Concluding Remarks
Throughout the calculation examples, we have used two sets of bounding relative
permeability curves: One set for co-current flow and one set for counter-current flow. While the
reduction in the bounding counter-current relative permeability relative to the co-current values
is consistent with experimental observations, it is not clear that the maximum residual CO
2

saturation will be identical for co-current and counter-current imbibition processes. Accordingly,
additional experimental evidence is needed to further understand the residual entrapment in
counter-current flow settings.

An additional challenge related to the modeling of transitions between co-current and
counter-current flow arises from the need for counter-current relative permeability data. There is
currently no standard approach for performing such measurements. However, the advent of
125

digital rocks and displacement calculations on rock images may serve as a substitute for
additional experimental efforts.
In the previous sections, we have presented, implemented and tested a new mobility model
that differentiates between co-current and counter-current flow. The new mobility model was
tested in a fully implicit simulator in the context of CO
2
injection into saline aquifers. 2D and 3D
displacement calculations clearly demonstrate that reductions in relative permeability due to
transitions between co-current and counter-current flow settings can impact rates and amounts of
CO
2
entrapment significantly. Based on our examples and analysis, we arrive at the following
conclusions:
1. A simple approach for representing the impact of transitions between co-current and counter-
current flow on fluid mobilities has been introduced.
2. The mixed mobility models allow us to represent reductions in relative permeability due to
flow transitions without the need for a full tensor representation of relative permeability
during multiphase flow.
3. The proposed methodology was tested in a fully implicit simulator without interfering with
the stability and general performance of the numerical calculations.
4. The representation of flow transitions was demonstrated to have a significant impact on both
amounts of residual entrapment and time scales for entrapment of CO
2
during and after
injection into an aquifer.

126

Chapter Seven: Transitions between Co-Current and Counter-Current Mobility in Three-
Phase Flows

7.1 Introduction
In the previous chapter, the transition between co- and counter-current relative permeability was
studied for two-phase flow settings. We applied a velocity-dependent transition coefficient to
calculate relative permeability of each phase by tracing relative permeability between the input
co- and counter-current relative permeabilities in the simulation of CO
2
migration in saline
aquifers.
In this chapter we apply the modification of phase mobilities to three-phase flow settings.
We present the simulation results of injecting gas (CO
2
) and water into an oil reservoir, where
both oil and water are mobile, and compare the results of oil, water and gas production rates with
the standard simulations where the effect of counter-current flow is not considered. The chapter
organization is as follows: First the modifications of the relative permeability calculations in the
simulation code are described. Then the example model is described. Finally the model is used
for simulation of enhanced oil recovery processes (simultaneous injection of water and gas and
water-alternating-gas processes) and a comparison of simulation results with and without the
modification for transitions between co- and counter-current flows is presented.

7.2 Three-Phase Relative Permeability Calculation
In this work, the modified Stone I model by Aziz and Settari (1979) is used in the calculation of
phase relative permeabilities.

127

Two sets of data are needed to calculate the relative permeabilities in a three-phase flow:
oil-water relative permeability data and gas-oil relative permeability data. The water relative
permeability is assumed to be a function of water only. The gas relative permeability is assumed
to be a function of gas only, although this assumption is challenged in some studies (Moulu et
al., 1998). However, the oil relative permeability depends on two saturations.

7.2.1 Stone I model
Stone (1970) introduced his first three-phase relative permeability model based on fluid flow
through channels. The oil relative permeability is calculated by:

(7.1)
with:

(7.2)

(7.3)
and:

(7.4)

(7.5)

(7.6)
In the above definitions,

is the relative permeability to oil taken from oil-water data,

is
the relative permeability to oil taken from gas-oil data,

is the connate water saturation, and

is the three-phase residual oil saturation. S
om
is different from the residual oil saturation (in
two-phase flow), because it can be further reduced by gas injection.
128

7.2.1.1 Aziz and Settari’s modification of the Stone I model
The Stone I model should reduce to the two-phase relative permeability data in the extreme cases
of

and

. However, this is true only when:

(7.7)
The model that is used in this work is the modified model by Aziz and Settari (1979):

(7.8)
where:

(7.9)

(7.10)
and

.
In the above definitions

, and the limit of (7.7) is removed.

7.2.1.2 Residual oil saturation (S
om
)
Residual oil saturation is needed for the calculation of three-phase relative permeability in the
Stone I model. The parameter S
om
in the Stone I model is a free parameter to adjust the residual
oil measurements. Fayers and Matthews (1984) proposed:

(7.11)
with:

(7.12)
where S
orw
is the residual oil saturation in oil-water data and S
org
is the residual oil saturation in
gas-oil data.
129

7.2.2 Counter-current Relative Permeability in Three-Phase Flows
The calculation of three-phase relative permeability is considered to be an interpolation of the
two two-phase relative permeabilities. Therefore, the velocity-dependent relative permeabilities
that were used for the calculation of phase mobilities in two-phase flows are extended for three-
phase flows as discussed in the following subsection.

7.2.2.1 Velocity-dependent relative permeability in three-phase flows
In the calculation examples presented in this chapter, it is assumed that water is the wetting
phase. Therefore, the relative permeability of water in counter-current flows is calculated by
considering the oil and gas phases as the non-wetting phase relative to water. The water relative
permeability is obtained using the input co-current and counter-current relative permeabilities of
water and oil. Since water is in contact with oil only, its relative permeability is related to the
average velocities of oil and water:

(7.13)
where

is an interpolating parameter that is calculated by:

(7.14)
where

and

. If a cell is in a co-current mode in a
specific direction (e.g. x, y and z), co-current relative permeability is used for that cell in that
specific direction. Equations (7.13) and (7.14) are used only for cells with counter-current flows.
Since directional relative permeabilities are used in the calculation of phase mobilities, it is
possible that the input co-current relative permeability is used in one direction (e.g. x) and that an
interpolation of the co-current and counter-current inputs is used in another direction (e.g. z).
130

The gas phase is the non-wetting phase and water + oil combined is assumed to be the
wetting phase. In the case of counter-current flow between gas and oil, the gas relative
permeability is accordingly interpolated between the input co- and counter-current input relative
permeabilities for gas:

(7.15)
where

is the interpolating parameter that is calculated by:

(7.16)
where

and

.
The relative permeability of oil is then calculated by:

(7.17)
It is noticed that the oil is in contact with both water and gas and, therefore, the interpolation for
this phase is done using the interpolating parameter of both other phases.
The input relative permeabilities used in subsequent calculation examples are shown in
Figure 7.1 and Figure 7.2. The derivatives of the relative permeabilities are calculated as:

(7.18)
because oil and gas are considered to be the non-wetting phase relative to water, and

(7.19)

and

are calculated numerically.

131

Figure 7.1. Relative permeability of the oil-water system

Figure 7.2. Relative permeability of the gas-oil system
132

7.3 Reservoir Model and Fluid Properties
The reservoir model is a 3D model with a homogeneous permeability field. In the first example,
the dimensions of the reservoir are 4000*4000*240 ft represented by 600 cells (10*10*6). The
horizontal permeability is 200 mD and the vertical permeability is 50 mD. Gas (CO
2
) is injected
from an injection well in one corner of the reservoir (in the bottom layer) with a rate of 70000
MSCF/D. Water is injected from the same corner on top of the injection well (top two layers). A
single production well is placed in the opposite corner of the reservoir (completed in all layers) at
a pressure below the hydrostatic pressure: the initial reservoir pressure is 4200 psi and the
production well is at operated at a constant BHP of 3500 psi. In the second example, the
reservoir has the dimensions of 2400*2400*240 ft with a homogeneous but non-isotropic
permeability (k
x
=k
y
=100 mD, k
z
=20 mD). There is a production well in the middle of the
reservoir that operates at a constant bottom-hole pressure of 3200 psi and it is completed
throughout the whole vertical direction. Four injections, located on the corners of the reservoir,
are injecting gas from the bottom two layers at a constant rate of 10000 MSCF/day each and
water from the top two layers at a constant rate of 2000 bbl/day each.

7.3.1 Fluid Properties
The water, oil and gas properties are listed in Table 7-1,

Table 7-2 and Table 7-3. In the base case, the initial water saturation is 0.16 (irreducible
saturation) and the initial oil saturation is 0.84.

133

Table 7-1. Properties of the gas phase
P B
O
VISC RGO
psi RB/STB Cp SCF/STB
1232.011 0.009234 0.017362 0
1540.013 0.006017 0.017784 0
1848.016 0.004036 0.018372 0
2156.019 0.003214 0.019196 0
2464.021 0.002868 0.02038 0
2772.024 0.002675 0.022136 0
3080.027 0.002548 0.024823 0
3388.029 0.002455 0.028851 0
3696.032 0.002383 0.034048 0
4004.035 0.002324 0.039317 0
4312.037 0.002275 0.043918 0
4620.04 0.002233 0.047813 0

Table 7-2. Properties of the oil phase
P B
O
VISC RGO
psi RB/STB cp SCF/STB
14.7 1.062 1.04 10
264.7 1.15 0.975 90.5
514.7 1.207 0.91 180
1014.7 1.295 0.83 370
2014.7 1.435 0.695 630
2514.7 1.5 0.641 770
3014.7 1.565 0.594 930
4014.7 1.695 0.51 1270
5014.7 1.671 0.549 1270

134

Table 7-3. Properties of the water phase
P B
O
VISC RGO
psi RB/STB cp SCF/STB
14.7 1.041 0.31 0
264.7 1.0403 0.31 0
514.7 1.0395 0.31 0
1014.7 1.038 0.31 0
2014.7 1.035 0.31 0
2514.7 1.0335 0.31 0
3014.7 1.032 0.31 0
4014.7 1.029 0.31 0
9014.7 1.013 0.31 0

7.4 Calculation Example 1: Simultaneous Injection of Water and Gas (Two Wells)
In this section we present and discuss the results of including the effect of counter-current flow
in three-phase flow simulations. The vertical permeability that is used for the base simulation
case is 0.25 of the horizontal permeability (k
x
=200mD). Methane is injected from one corner (in
the bottom layer) of the reservoir while water is injected into the reservoir from the top two
layers. The water injection rate is constant at 8000 bbl/d for 4000 days, and is then reduced to
2000 bbl/d for the next 1000 days (5000 days from initial time) and is then stopped. The methane
injection rate is initially set to zero until 4000 days. Then 5000 MSCF/D is injected until 5000
days and then increased to 15000 MSCF/D till the end of the simulation (10000 days). This
example includes a production well in the opposite corner of the reservoir that is operated at a
constant bottom hole pressure of 3500 psi (till 4000 days), followed by 3300 psi (till 5000 days).
The BHP is then further reduced to 3000 psi until the end of the simulation. The initial reservoir
pressure is 4200 psi, and the initial oil saturation in the reservoir is 0.84. Simulations were
performed with co-current and mixed models and the results are presented below. The total
135

injected water is 0.2 of the reservoir pore volume and the total injected gas is 0.37 of the
reservoir pore volume.
The gas (methane) saturation profiles are shown in Figure 7.3. The difference between the
two models is on the edge of the gas plume indicating that a better areal sweep is predicted by
the mixed model.

Figure 7.3. Methane saturation profiles as predicted by co-current model (left), mixed model
(middle) and their difference (right) in the 3
rd
layer from the top

0 0.35 0.7

0 0.35 0.7

0 0.05
1000 days
6000 days
10000 days
136

Figure 7.4 shows the oil saturation profiles. At 1000 days only water is injected. The
counter-current flow between oil and water reduces the mobility of both phases. At 6000 days
there is counter-current flow between oil and gas and it continues until the end of simulation.

Figure 7.4. Oil saturation profiles as predicted by co-current model (left), mixed model (middle)
and their difference (right) in the 3
rd
layer from the top

Water saturation profiles are shown in Figure 7.5 and Figure 7.6 and we note that the effect
of counter-current flow is seen primarily during the water injection period.
1000 days
6000 days
10000 days

0 0.35 0.7

0 0.35 0.7

-0.05 0 0.05
137

Figure 7.5. Water saturation profiles as predicted by co-current model (left), mixed model
(middle) and their difference (right) in the 3
rd
layer from the top

0 0.35 0.7

0 0.35 0.7

0 0.05
1000 days
6000 days
10000 days
138

Figure 7.6. Water saturation profiles as predicted by co-current model (left), mixed model
(middle) and their difference (right) in the 5
th
layer from the top

Figure 7.7 compares the injection and production rates as predicted by the co-current and
mixed models. The methane injection is exactly the same between the two models. The gas and
oil production rates are predicted to be very similar by the co-current and mixed models.
However, the water production rate is predicted to be higher in the co-current model. So, why is
the oil production rate predicted to be the same from the two models whereas the water
production rate is predicted to be different? A closer look at the production wells reveals that
completions in the lower half of the reservoir (in particular the completion at the bottom of the

0 0.35 0.7

0 0.35 0.7

0 0.05
1000 days
6000 days
10000 days
139

production well) are responsible for the water production. This can be verified from Figure 7.6
where regions of higher water saturation are observed in the bottom of the reservoir. The main
difference in water production between the co-current and mixed models is the production from
the most bottom production wells: The completions in the bottom layer of the reservoir account
for more than 90% of the water production rate. This is the reason that in the lower right graph of
Figure 7.7 there is a significant difference in water production and yet the oil and gas production
rates are similar.

Figure 7.7. Comparison of injection and production rates between co-current and mixed models
in layer 2 from top

Figure 7.8 compares the pressure in the water injection well. Since the vertical mobility is
reduced in some cells in the mixed model (located below this well), the injection pressure is
higher in this model. The maximum difference in pressure is about 26 psi at 1200 days.
0 2000 4000 6000 8000 10000
0
2
4
6
8
10
x 10
7
time (years)
cummulative injected
methane (MSCF)

co-current
mixed
0 2000 4000 6000 8000 10000
0
5
10
15
x 10
7
time (years)
cummulative gas
production (MSCF)
0 2000 4000 6000 8000 10000
0
1
2
3
4
x 10
7
time (years)
cummulative oil
production (STB)
0 2000 4000 6000 8000 10000
0
1
2
3
x 10
6
time (years)
cummulative water
production (STB)
140

Figure 7.8. Bottomhole pressure in one of the injection wells (injecting water)

Table 7-4 compares the numerical performance of the simulations with co-current and mixed
models, and we observe that the simulation time of the mixed model increases significantly (by a
factor of 7) compared to the co-current model.

Table 7-4. Comparison of numerical performance in the co-current and mixed models

co-current mixed
Simulation time (sec) 108 201
Number of time steps 217 214
Average time step size (days) 46.1 46.7
Newton iterations 528 1333

141

7.5 Example Calculation 2: Simultaneous Injection of Water and Gas (Five-Spot Pattern)
In this section the effect of counter-current flow on the simultaneous injection of water and gas is
studied. We observe that counter-current flow occurs mostly between the oil and water but it can
also occur between the injected gas and the oil in places where the water is near its residual
saturation.
Figure 7.9 shows the gas saturation profiles obtained from the co-current and mixed
models. Since gas is injected from the bottom of the reservoir and water is initially at its
irreducible saturation, counter-current flow occurs between gas and oil. The oil saturation
profiles exhibit more pronounced differences at the edges of the water front where counter-
current flow with water makes a difference in the segregation distance between the two models
(Figure 7.10). Water saturation profiles for layer 8 are shown in Figure 7.11.
142

Figure 7.9. Gas saturation profiles as predicted by co-current model (left), mixed model (middle)
and their difference (right) in the 1
st
layer from the top

0 0.35 0.7

0 0.35 0.7

-0.05 0 0.05
1000 days
4000 days
7300 days
143

Figure 7.10. Oil saturation profiles as predicted by co-current model (left), mixed model
(middle) and their difference (right) in the 1
st
layer from the top

Figure 7.12 presents a comparison of the production rate of gas, oil and water as predicted
by the co-current and the mixed models. Figure 7.12 confirms that the gas injection rate is the
same for the two models. Gas and oil productions are predicted to be similar from the two
models. However, the water production is predicted to be higher in the co-current model (~20%).
Based on Figure 7.9-Figure 7.11 we observe that there is little difference between the two
models except for the water production.

0 0.35 0.7

0 0.35 0.7

-0.05 0 0.05
1000 days
4000 days
7300 days
144

Figure 7.11. Water saturation profiles as predicted by co-current model (left), mixed model
(middle) and their difference (right) in the 8
th
layer from the top

The simulations were done on the General Purpose Research Simulator (GPRS) from
Stanford University. An implicit formulation was used in this example. No instability in the
numerical calculations was observed. However, the simulation time of the mixed model was
found to be higher. The average time step size in the co-current model is larger than in the mixed
model. Table 7-5 summarizes the performance of the simulator for the two simulations.
Although the number of time steps is higher in the co-current model, it is faster compared to the

0 0.35 0.7

0 0.35 0.7

0 0.05
1000 days
4000 days
7300 days
145

mixed model because the number of Newton iterations used in the co-current model is less than
for the mixed model.

Figure 7.12. Comparison of injection and production rates between co-current and mixed models

Table 7-5. Numerical comparison between co-current and mixed models

co-current mixed
Simulation time (sec) 3306 4797
Number of time steps 118 125
Average time step size (days) 61.9 58.4
Newton iterations 1678 1817

0 2000 4000 6000 8000
0
2
4
6
x 10
7
time (years)
cummulative injected
CO
2
(MSCF)

co-current
mixed
0 2000 4000 6000 8000
0
1
2
3
4
x 10
8
time (years)
cummulative gas
production (MSCF)
0 2000 4000 6000 8000
0
0.5
1
1.5
2
2.5
x 10
7
time (years)
cummulative oil
production (STB)
0 2000 4000 6000 8000
0
2
4
6
8
x 10
7
time (years)
cummulative water
production (STB)
146

7.6 Concluding Remarks
In this chapter, our representation of transitions between co- and counter-current flows was
extended to three-phase displacement problems. Velocity-dependent relative permeabilities were
used for both gas/oil and oil/water phases. We tested the new model in the simulation of
simultaneous gas and water injection into an oil reservoir, and from the presented results, we
arrive at the following observations/conclusions:
1. We observe that counter-current flow between gas and oil is modest and that it may be
reasonable to neglect the flow transitions for this pair. However, counter-current flow
between oil and water is more significant, primarily, in the vertical direction. Regions of
counter-current flow (between oil and water) are not static but vary with time.
2. Numerical simulation of three-phase flow using the proposed velocity-dependent relative
permeability model does not result in numerical instabilities (in the presented examples).
However, the simulation time increases in the mixed model. This result is in contrast to the
two-phase flow results where no difference between co-current and mixed model was
observed in terms of numerical performance.
3. Including counter-current flow in the numerical simulation of three-phase flows does not
interfere with the production of oil and gas (at least for the examples presented). The main
difference between the two models is observed for the water production rate which is lower
in the mixed model. The pressure of the injection wells is slightly higher in the mixed model
in the examples presented.

147

Chapter Eight: Recommendations for Future Studies

In this work we have emphasized the importance of including counter-current flow in the
numerical simulation of multi-phase flow in porous media. We have shown that there is a
difference between co-current and counter-current relative permeability and the counter-current
relative permeability should be included in simulation studies. We input counter-current relative
permeability as an additional input (in addition to co-current relative permeability) and use an
interpolation scheme (or discrete) to transition between the two sets. We have applied this
analysis in the numerical simulation of CO
2
injection into saline aquifers and also in the context
of enhanced oil recovery (in water-alternating-gas (WAG) and also simultaneous injection of
water and gas processes).
Although we have applied the counter-current flow in several processes, however, there are
some questioned that needs further studies:
1- In our numerical simulations we have assumed that the counter-current relative permeability
is given. However, it is necessary to obtain the counter-current relative permeability either
with experimental studies or pore-scale simulations. Therefore, pore-scale modeling of
counter-current flow is necessary. Pore-scale modeling (e.g. Lattice Boltzman simulation) is
needed to quantify the effect of counter-current flow on fluid mobility. In very small scales,
molecular dynamics can be used. If any changes of relative permeability in the pore-scale
simulation are observed then it can be argued how to upscale and include these changes in
larger scales.
2- In order to integrate counter-current flow in fluid flow simulations, one can use the
generalized Darcy equation. The generalized Darcy equation has non-zero off-diagonal terms
148

in the matrix that contains fluid mobilities. These terms can be obtained using pore-scale
simulations. Another method that was used in this study is to use the conventional Darcy
equations, but account for off-diagonal terms by inputting another set of relative permeability
(for counter-current flows). It would be interesting to use both methods and compare them in
terms of simulation results (pressures and saturations) and also numerical performance.
3- If two sets of relative permeabilities are used in the simulation of multi-phase flows, the
interpolation scheme between the input co-current and counter-current relative permeabilities
needs further studies. Two methods have been used in our studies: discrete method and
velocity-dependent method. The velocity-dependent method is smoother. However, further
studies needs to determine which parameters should be considered in the interpolation.
4- In three-phase flows the problem becomes more complicated. We considered a three-phase
system as a combination of two two-phase systems and used the same procedure that was
used in the two-phase flows. Given the fact that there is no general three-phase flow model
for relative permeability, accounting for counter-current flow in a three-phase system
becomes even more complicated. Any modification of phases’ mobilities in three-phase
flows impact the optimization of enhanced oil recovery mechanisms where three mobile
phased are present.
5- The nonlinear solver did not show any problem in its performance in the two-phase flows for
various examples that were used. However, some convergence problems arose in some three-
phase flow examples. A mathematical analysis to guarantee the convergence of the solver,
especially in cases where a jump in relative permeability occurs (e.g. in the discrete mixed
model), is very helpful. In two-phase flows we have tried various examples with varying
149

degrees of gravity, viscous and capillary forces. Stability was not an issue in any of them.
However, a mathematical framework in this area is needed.
6- Counter-current flow was observed in numerical simulation of CO
2
sequestration in saline
aquifers in various scales. With an increase or decrease of the grid numbers by an order of
magnitude in the examples that were used in this study, the pattern of counter-current flow
remains the same in the vertical migration of CO
2
. However, if a velocity-dependent relative
permeability is used in the interpolation scheme to calculate fluid mobilities, the degree of
decrease in relative permeability is different (even if all the cells in a region experience
counter-current flow). If the flow properties are needed to be upscaled in a region where both
co-current and counter-current flows are present, or even in a pure counter-current flow with
varying degree of counter-current effect, the upscaling scheme becomes more complicated.
7- The focus of this study has been on the effect of counter-current flow on fluid mobility.
However, other parameters, such as capillary pressure, can be affected by counter-current
flow. Since capillary pressure and relative permeability are related to each other, a change in
relative permeability can affect capillary pressure. If it is believed that fluid distribution in
counter-current flow is different from co-current flow, then it is not surprising that capillary
pressure can be changed. The change of capillary pressure is important in CO
2
migration in
saline aquifers, especially as the degree of heterogeneity increases in an aquifer. Also, it is
possible that in counter-current flow the residual saturation of CO
2
changes.
8- Our examples were tested mostly in black-oil formulations. It would be essential to test it in
compositional simulation studies in enhanced oil recovery methods.

150

References
Akinnikawe, O., Chaudhary, A., Vasquez, O., 2010. Increasing CO
2
-storage efficiency through a
CO
2
brine-displacement approach. In SPE International Conference on CO
2
Capture,
Storage, and Utilization. New Orleans, Louisiana.
Aleman, M.A., Slattery, J.C., 1988. Estimation of three-phase relative permeabilities, 3(2), 111-
131.
Al-Wadahi, M., Grader, A.S., Ertekin, T., 2000. An investigation of three-phase counter-current
flow using X-ray computerized tomography and neuro-simulation modeling. In Proceedings
of SPE Annual Technical Conference and Exhibition. Dallas, Texas.
Altundas, Y.B., Ramakrishnan, T.S., Chugunov, N., de Loubens, R., 2010. Retardation of CO
2

due to capillary pressure hysteresis : a new CO
2
trapping mechanism. In SPE International
Conference on CO
2
Capture, Storage, and Utilization. New Orleans, Louisiana.
Anchliya, A. & Ehlig-Economides, C., 2009. Aquifer management to accelerate CO
2
dissolution
and trapping. SPE 126688, In Proceedings of SPE International Conference on CO
2

Capture, Storage, and Utilization. San Diego, California.
Aziz, K., Settari, A., 1979. Petroleum Reservoir Simulation, Applied Science Publishers,
London.
Bachu, S, 2002. Sequestration of CO
2
in geological media in response to climate change: road
map for site selection using the transform of the geological space into the CO
2
phase space.
Energy Conversion and Management, 43(1), 87-102.
Bachu, Stephan, Gunter, W.D., Perkins, E.H., 1994. Aquifer disposal of CO
2
: hydrodynamic and
mineral trapping. Energy Conversion and Management, 35(4), 269-279.
Baker, L.E., 1988. Three-phase relative permeability correlations, SPE 17369, In Proceedings of
SPE Enhanced Oil Recovery Symposium, Tulsa, Oklahoma.
151

Baker, L.E., 1995, Three-phase relative permeability of water-wet, intermeadiate-wet, and oil-
wet sandstone, New Developments in Improved Oil Recovery, Geological Society Special
Publication, 84, 51-61.
Balbinski, E.F., Fishlock, T.P., Goodyear, S.G., Jones, P.I.R., 1999. Key characteristics of three-
phase oil relative permeability formulations for improved oil recovery predictions,
Petroleum Geoscience,
Bentsen, R., 1998. Effect of momentum transfer between fluid phases on effective mobility.
Journal of Petroleum Science and Engineering, 21(1-2), 27-42.
Bentsen, R.G., Manai, A.A., 1993. On the use of conventional cocurrent and countercurrent
effective permeabilities to estimate the four generalized permeability coefficients which
arise in coupled, two-phase flow. Transport in Porous Media, 11(3), 243-262.
Bourbiaux, B.J., Kalaydjian, F.J., 1990. Experimental study of cocurrent and countercurrent
flows in natural porous media. SPE Reservoir Engineering, 5(3), 361-368.
Blunt, M., 2000. An empirical model for three-phase relative permeability, In Proceedings of
SPE Annual Technical Conference and Exhibition. Houston, Texas.
CMG, 2010. Computer Modeling Group Ltd.
Corey, A.T., Rathjens, C.H., 1956. Effect of stratification on relative permeability. Journal of
Petroleum Technology, 8(12), 69-71.
Corey, A.T., Rathjens, C.H., Henderson, J.H., Wyllie, M.R.J., 1956. Three-phase relative
permeability. Journal of Petroleum Technology, 8(11), 63-65.
Crotti, M.A., Inlab, S.A., Rosbaco, J.A., 1998. Relative permeability curves: the influence of
flow direction and heterogeneities - Dependence of end point saturations on displacement
mechanisms. SPE 39657, In SPE/DOE Improved Oil Recovery Symposium. Tulsa,
Oklahoma.
152

de la Cruz, V., Spanos, T.J.T., 1983. Mobilization of oil ganglia. AIChE Journal, 29(5), 854-858.
Dehghanpour, H., Aminzadeh, B., Mirzaei, M., DiCarlo, D.A., 2011. Flow coupling during
three-phase gravity drainage. Physical Review E, 83(6), 1-4.
Delshad, M., Pope, G.A., 1989. Comparison of the three-phase oil relative permeability models.
Transport in Porous Media, 4(1), 59-83.
Dengen, Z., Fayers, F., Orr, F., 1997. Scaling of multiphase flow in simple heterogeneous porous
media. SPE Reservoir Engineering, 12(3), 173–178.
Dietrich, J.K., Bondor, P.L., 1976. Three-phase oil relative permeability models, SPE 6044, In
Proceedings of SPE Annual Fall Technical Conference and Exhibition, New Orleans,
Louisiana.
Dullien, F.A.L., Dong, M., 1996. Experimental determination of the flow transport coefficients
in the coupled equations of two-phase flow in porous media. Transport in Porous Media,
25(1), 97-120.
Eastwood, J., Spanos, T.J.T., 1991. Steady-state countercurrent flow in one dimension.
Transport in Porous Media, 6(2), 173-182.
Eclipse Manual, 2010. Schlumberger.
Ennis-King, J., Paterson, L., 2005. Role of convective mixing in the long-term storage of carbon
dioxide in deep saline formations. SPE Journal, 10(3), 5-8.
Fankhauser, S. et al., 2001. Vulnerability to climate change and reasons for concern: A synthesis.
In Climate Change 2001: Impacts, Adaption, and Vulnerability. 941-967.
Fayers, F.J., Matthews, J.D., 1984. Evaluation of normalized Stone’s methods for estimationg
three-phase relative permeabilities. SPE Journal, 224-232.
153

Fenghour, A., Wakeham, W.A., Vesovic, V., 1998. The viscosity of carbon dioxide. Journal of
Physical and Chemical Reference Data, 27(1), 31–44.
Gasda, S.E., Nordbotten, J.M., Celia, M.A, 2011. The impact of local-scale processes on large-
scale CO
2
migration and immobilization. Energy Procedia, 4, 3896-3903.
Hassanzadeh, H., Pooladi-Darvish, M., Keith, D.W., 2006. Stability of a fluid in a horizontal
saturated porous layer: effect of non-linear concentration profile, initial, and boundary
conditions. Transport in Porous Media, 65(2), 193-211.
Hassanzadeh, H., Pooladi-Darvish, M., Keith, D.W., 2009. Accelerating CO
2
dissolution in
saline aquifers for geological storage - mechanistic and sensitivity studies. Energy & Fuels,
23(6), 3328-3336.
Herzog, H., Golomb, D., 2004. Carbon capture and storage from fossil fuel use. Encyclopedia of
Energy, 1(6562), 277-287.
Honarpour, M.M., Cullick, A.S., Saad, N., 1994. Influence of small-scale rock laminations on
core plug oil/water relative permeability and capillary pressure. SPE 27968, In University of
Tulsa Centennial Petroleum Engineering Symposium. Tulsa, Oklahoma.
Houghton, J.T. et al., 2001. Climate change 2001: The scientific basis - contribution of Working
Group I to the third assessment report of the Intergovernmental Panel on Climate Change.
IPCC, 2005. Carbon dioxide capture and storage,
Ide, T.S., Jessen, K., Orr, F.M., 2007. Storage of CO
2
in saline aquifers: Effects of gravity,
viscous, and capillary forces on amount and timing of trapping. International Journal of
Greenhouse Gas Control, 1(4), 481-491.
Iverson, W.P., Dunn, T.L., Ajdari, I., 1996. Relative permeability anisotropy measurements in
Tensleep sandstones. SPE 35435, In SPE/DOE Improved Oil Recovery Symposium. Tulsa,
Oklahoma.
154

Javaheri, M., Jessen, K., 2011. Integration of counter-current relative permeability in the
simulation of CO
2
injection into saline aquifers. International Journal of Greenhouse Gas
Control, 5(5), 1272-1283.
Javaheri, M., Abedi, J., Hassanzadeh, H., 2009. Linear stability analysis of double-diffusive
convection in porous media, with application to geological storage of CO
2
. Transport in
Porous Media, 84(2), 441-456.
Kalaydjian, F., 1987. A macroscopic description of multiphase flow in porous media involving
spacetime evolution of fluid/fluid interface. Transport in Porous Media, 2(6), 537-552.
Kalaydjian, F., 1990. Origin and quantification of coupling between relative permeabilities for
two-phase flows in porous media. Transport in Porous Media, 5(3), 215-229.
Killough, J.E., 1976. Reservoir simulation with history-dependent saturation functions. SPE
Journal, 16(1), 37-48.
Kokal, S.L., Maini, B.B., 1990. An improved model for estimating three-phase oil-water-gas
relative permeabilities from two-phase oil-water and oil-gas data. Journal of Canadian
Petroleum Technology, 29(2), 105-114.
Kopp, A., Class, H., Helmig, R., 2009a. Investigations on CO
2
storage capacity in saline
aquifers. International Journal of Greenhouse Gas Control, 3(3), 263-276.
Kopp, A., Class, H., Helmig, R., 2009b. Investigations on CO
2
storage capacity in saline
aquifers—Part 2: Estimation of storage capacity coefficients. International Journal of
Greenhouse Gas Control, 3(3), 277-287.
Kumar, A., Noh, M., Pope, G.A., Sepehrnoori, K., Bryant, S., Lake, L.W., 2004. Reservoir
simulation of CO
2
storage in deep saline aquifers. SPE 89343, In SPE/DOE Improved Oil
Recovery Symposium. Tulsa, Oklahoma.
155

Land, C., 1968. Calculation of imbibition relative permeability for two- and three-phase flow
from rock properties. SPE Journal, 8(2), 149-156.
Langaas, 2001. Numerical investigations of the steady state relative permeability of a simplified
porous medium. Transport in Porous Media, 45, 241-266.
Larsen, J.A., Skauge, A., 1995. Comparing hysteresis models for relative permeability in WAG
studies, Society of Core Analysts Conference Paper 9506, San Francisco, California.
Lefebvre, E.P., 1978. Gravity and capillarity effects on imbibition in porous media. SPE Journal,
18(3), 195–206.
Lelievre, R., 1966. Étude d’écoulements diphasiques permanents à contre-courants en milieu
poreux: comparaison avec les écoulements de même sens. University of Toulouse, France.
Leonenko, Y., Keith, D.W, 2008. Reservoir engineering to accelerate the dissolution of CO
2

stored in aquifers. Environmental science & technology, 42(8), 2742-7.
Leonenko, Y., Pooladi-Darvish, M., Keith, D.W., Hassanzadeh, H., 2006. Accelerating the
dissolution of CO
2
in aquifers. In GHGT-8, 8th International Conference on Greenhouse
Gas Control Technologies. Trondheim, Norway.
Li, G., Karpyn, Z.T., Halleck, P.M., Grader, A.S., 2005. Numerical simulation of a CT-scanned
counter-current flow experiment. Transport in Porous Media, 60(2), 225-240.
Maini, B., Coskuner, G., Jha, K., 1990. A comparison of steady-state and unsteady-state relative
permeabilities of viscous oil and water in Ottawa sand. Journal of Canadian Petroleum
Technology, 29(2), 72–77.
Manai, A., 1991. The measurement of cocurrent and countercurrent relative permeabilities and
their use to estimate generalized relative permeabilities. University of Alberta.
156

Marland, G., Boden, T.A., Andres, R.J., 2006. Global, regional, and national fossil-fuel CO
2

emissions, Technical Report, Oak Ridge National Laboratory, doi
10.3334/CDIAC/00001_V2010.
Moulu, J.C., Vizika, O., Egermann, P., Kalaydjian, F., 1999. A new three-phase relative
permeability model for various wettability conditions, In Proceedings of SPE Annual
Technical Conference and Exhibition, Houston, Texas.
Moulu, J.C., Vizika, O., Kalaydjian, F., Duquerroix, J.P., 1998. A new model to calculate three-
phase relative permeabilities: Application and validation for a sandstone, Revue de l’institut
francais du petrole, 53(4), July-August issue.
Muskat, M., 1937. The Flow of Homogeneous Fluids Through Porous Media, McGraw-Hill
Book Company.
National Academy of Engineering, 2008. Available at: http://www.engineeringchallenges.org/.
National Research Council, 2001. Climate change science: An analysis of some key questions,
Washington, D.C.: National Academy Press.
Nattwongasem, D., Jessen, K., 2009. Residual trapping of CO
2
in aquifers during the counter-
current Flow. SPE 125029, In Proceedings of SPE Annual Technical Conference and
Exhibition. New Orleans, Louisiana.
Odeh, A., 1959. Effect of viscosity ratio on relative permeability. Trans. AIME, 216, 346–352.
Oreskes, N., 2004. Beyond the ivory tower. The scientific consensus on climate change. Science,
306(5702), 1686.
Park, A.H.A., 2005. Carbon dioxide sequestration: chemical and physical activation of aqueous
carbonation of Mg-bearing minerals and pH swing process. The Ohio State University.
Prats, M., Lake, L.W., 2008. The anisotropy of relative permeability. JPT Technical Notes,
60(3), 99.
157

Rossen, W., van Duijn, C.J., Nguyen, Q.P., Vikingstad,, A.K., 2006. Injection strategies to
overcome gravity segregation in simultaneous gas and liquid injection into homogeneous
reservoirs. SPE 99794, In Proceedings of SPE/DOE Symposium on Improved Oil Recovery.
Tulsa, Oklahoma.
Shad, S., Gates, I.D., 2010. Multiphase flow in fractures: co-current and counter-current flow in
a fracture. Journal of Canadian Petroleum Technology, 49(2), 48-55.
Shamshiri, H., Jafarpour, B., 2010. Optimization of geologic CO
2
storage in heterogeneous
aquifers through improved sweep efficiency. SPE 139643, In SPE International Conference
on CO
2
Capture, Storage, and Utilization. New Orleans, Louisiana.
Shojaei, H., Rastegar, R., Jessen, K., 2012. Mixing and mass transfer in multicontact miscible
displacements. Transport in Porous Media, 94(3), 837-857.
Stone, H.L., 1970. Probability model for estimating three-phase relative permeability. Journal of
Petroleum Technology, 22(2), 214-218.
Stone, H.L., 1973. Estimation of three-phase relative permeability and residual oil data. Journal
of Canadian Petroleum Technology, 12(4), 53-61.
Tang, X., 2006. Numerical modeling of deformations caused by carbon dioxide sequestration in
coal seams. West Virginia University.
US DOE, 2004. Carbon sequestration R&D. Available at:
http://www.fe.doe.gov/programs/sequestration/index.html.
Whitaker, S., 1986. Flow in porous media II: The governing equations for immiscible, two-phase
flow. Transport in Porous Media, 1(2), 105-125.

Abstract (if available)

Abstract Carbon Capture and Storage (CCS) is the capture of carbon dioxide (CO₂) from large source points, such as power plants, and storing it in geological formation. It is a potential means of reducing CO₂ concentration in the atmosphere to mitigate global warming and ocean acidification. Geological formations are considered to be the most promising storage sites for carbon dioxide. Deep saline aquifers have the most storage capacity compared to other geological sites. ❧ After CO₂ is injected into an aquifer, a fraction of it, that can be large, is immobilized in the form of residual phase. Residual entrapment is an important component in successful storage of CO₂ in saline aquifers. Trapped CO₂ is eventually dissolved into the brine. Estimation of the fraction of the injected CO₂ that is trapped depends on CO₂ saturation profiles and saturation history. In the simulation of CO₂ injection into saline aquifers, conventional relative permeabilities functions are commonly used and it is assumed that viscous coupling is negligible. We show that the vertical migration of CO₂ in saline aquifers is often dominated by counter-current flow. Experimental and simulation studies of two- and three-phase flow in porous media show that the counter-current relative permeability is less than co-current relative permeability, that are commonly used (current industry standard) in the simulation of multiphase flow. ❧ This study focuses on including a velocity- -dependent relative permeability model in the simulation of multiphase flow in porous media to account for the flow dynamics (co- vs. counter-current flow) on relative permeability and fluid mobility. We use both co-current and counter-current relative permeability in the simulation of CO₂ injection from a well into an aquifer and show that the plume saturation profile is influenced by the dynamic relative permeability changes (transitions between co-current and counter-current flow). The fraction of the injected CO₂ that is trapped and the time-scale of vertical migration of CO₂ both increase when the transitions from co- to counter-current relative permeability are included. We conclude that using one set of relative permeabilities is not sufficient in the simulation of CO₂ injection into saline aquifers. We also demonstrate that a velocity-dependent mobility does not introduce numerical instability and that the simulation time does not change significantly. ❧ We have applied the velocity-dependent relative permeability in the simulation of gas (CO₂) injection into an oil reservoir in the context of enhanced oil recovery with three mobile phases. Cumulative oil and gas production increases just slightly in the new model, but water production decreases quite significantly. That is due to a better sweep efficiency when counter-current relative permeability is included in the calculation of phases' mobilities. ❧ Counter-current relative permeability can be measured directly using a controlled experimental design to eliminate the effect of boundaries on saturation distribution. It can also be estimated by pore-scale simulations, e.g. Lattice Boltzman simulations. Both methods have been used in the literature and both show a reduction in the relative permeability in counter-current settings. We have tested an indirect procedure to conclude a reduction of relative permeability in counter-current flows. ❧ The numerical simulations are done using IMPES and fully-implicit methods. No numerical instabilities were observed in any of the two methods by including counter-current relative permeability either in a discrete scheme or by using a velocity-dependent scheme. We conclude that the new model can predict CO₂ saturation profile more accurately than the conventional standard model.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Modeling and simulation of complex recovery processes

$Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells$

PDF

Coarse-scale simulation of heterogeneous reservoirs and multi-fractured horizontal wells

PDF

Asphaltene deposition during co₂ injection in enhanced oil recovery applications

PDF

Mass transfer during enhanced hydrocarbon recovery by gas injection processes

PDF

A study of dispersive mixing and flow based lumping/delumping in gas injection processes

PDF

Pore-scale characterization and fluid-flow simulation: application of the Lattice Boltzmann method to conventional and unconventional rocks

PDF

Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection

PDF

The study of CO₂ mass transfer in brine and in brine-saturated Mt. Simon sandstone and the CO₂/brine induced evolution of its transport and mechanical properties

PDF

Efficient simulation of flow and transport in complex images of porous materials and media using curvelet transformation

$Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model$

PDF

Optimization of CO2 storage efficiency under geomechanical risks using coupled flow-geomechanics-fracturing model

PDF

Chemical and mechanical deformation of porous media and materials during adsorption and fluid flow

$Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing$

PDF

Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing

PDF

A study of diffusive mass transfer in tight dual-porosity systems (unconventional)

PDF

Effective flow and transport properties of deforming porous media and materials: theoretical modeling and comparison with experimental data

PDF

Integration of multi-physics data into dynamic reservoir model with ensemble Kalman filter

PDF

A coupled multiphase flow-geomechanics framework for modeling material, kinematic and contact nonlinearities in faulted reservoirs

PDF

Reactivation of multiple faults in oilfields with injection and production

PDF

Modeling and simulation of multicomponent mass transfer in tight dual-porosity systems (unconventional)

PDF

Optimization of coupled CO₂ sequestration and enhanced oil recovery

PDF

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

Asset Metadata

Creator Javaheri, Mohammad (author)

Core Title The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Publication Date 12/06/2013

Defense Date 05/31/2013

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag CO₂ sequestration,counter-current flow,multi-phase flow,OAI-PMH Harvest,porous media,relative permeability,reservoir simulation

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Jessen, Kristian (committee chair), Domaradzki, Julian A. (committee member), Ershaghi, Iraj (committee member)

Creator Email javaheri62@gmail.com,mohammaj@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-356870

Unique identifier UC11295233

Identifier etd-JavaheriMo-2211.pdf (filename),usctheses-c3-356870 (legacy record id)

Legacy Identifier etd-JavaheriMo-2211.pdf

Dmrecord 356870

Document Type Dissertation

Format application/pdf (imt)

Rights Javaheri, Mohammad

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA