University of Southern California Dissertations and Theses

Modeling and simulation of complex recovery processes

(USC Thesis Other)

Modeling and simulation of complex recovery processes

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

MODELING AND SIMULATION OF COMPLEX RECOVERY PROCESSES

by
Marjan Sherafati

A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirement for the Degree
DOCTOR OF PHILOSOPHY
(PETROLEUM ENGINEERING)

August 2018

Copyright 2018 Marjan Sherafati

ii

Dedication

To my parents, for their unconditional love and support

iii

Acknowledgement
I would like to express my sincere appreciation to anyone who has encouraged and supported me
throughout my doctoral studies.
In particular I would like to thank Professor Kristian Jessen, my PhD advisor who possesses every
quality one could wish to find in an advisor. Besides his deep knowledge in various aspects of
reservoir engineering, his curiosity and passion to solve challenging problems have constantly
inspired me during the course of my PhD at the University of Southern California. He gives his
students a lot of freedom in their research, which in turn nurtures self-esteem and creativity in
them. He also has great values and ethics. There are no words to describe my gratitude to him for
being incredibly understanding when I was going through a difficult time.
I have been privileged to work with two other world-renowned academics, Professor Iraj Ershaghi
and Professor Julian Domaradzki, as my dissertation committee members. Their valuable
comments and recommendations have significantly improved the quality of my doctoral
dissertation.
I would like to extend my gratitude to the Graduate School at USC for supporting my doctoral
studies through the Chevron’s PhD Fellowship. I also owe my gratitude to the Petroleum
Engineering program and the Graduate Student Government at USC for providing me with travel
grants to participate in professional conferences and exhibitions, and to MHA Petroleum
Consultants, Microseismic Inc., and California Resources Corporation for the internship
opportunities.
iv

I am extremely grateful to my friends and peers in Dr. Jessen’s research group. They helped me
greatly with their guidance and their support. My special thanks go to Dr. Hasan Shojaei, Dr.
Mohammad Javaheri, Dr. Shahram Farhadi and Dr. Mohammad Evazi for their valuable insights
and kind assistance in research, in helping me finding internships and employment, and in all else.
I am thankful to my amazing friends with whom I have shared wonderful memories. My special
thanks go to Mahshad Samnejad, Roya Ermagan, Shahrzad Hadian, Sona Sehat and Cyrus
Ashayeri. I am grateful to all of them, and to many more friends not mentioned here, for their
constant support and encouragement.
My deepest gratitude goes to my family; my lovely parents, Hamid and Mohtaram, who have loved
me unconditionally and guided me through life, and my brilliant siblings, Mandana, Arash and
Ardalan, who set a great example for me of success, kindness and strength. If it wasn’t for their
love, encouragement and support, I would not have been able to achieve the highest academic
degree.
I would like to especially thank my brother, Dr. Ardalan Sherafati, who has always been my idol
and my inspiration throughout life and my sister-in-law, Dr. Elham Foroozesh. Living halfway
across the globe from Iran, they are my only family in the US. I cannot thank them enough for
their kindness, their guidance, and for making me feel like I have a second home here in the US.
v

Abstract
Compositional reservoir simulation is a powerful tool commonly used to estimate the potential
hydrocarbon recovery from Enhanced Oil Recovery methods. Successful design and
implementation of EOR methods are highly dependent on the accuracy by which simulation tools
can represent the physics that govern the displacement behavior in a reservoir. In this work, we
investigate, and aim to improve the accuracy of some physical models that are frequently used to
describe and interpret multiphase flows via compositional reservoir simulation.
In the simulation of gas injection EOR methods, it is common to use relative permeability values
obtained from co-current experiments and it is assumed that the viscous coupling between flowing
phases is negligible. In this work we show that Water Alternating Gas (WAG) and Simultaneous
Water and Gas (SWAG) injection processes are dominated by counter-current flow in the vertical
direction. Previous experimental and simulation studies on counter-current flow have shown that
counter-current relative permeability is always less than the co-current relative permeability. Most
of the previous studies on the subject of counter-current flow have focused on two-phase flow. In
this work we focus on counter-current flow in a three-phase flow system such as WAG and SWAG
injection processes.
We use both co-current and counter-current relative permeability functions in the simulation of
these gas injection EOR processes (co-current flow in the horizontal direction and counter-current
flow in the vertical direction) and show that the saturation profiles, oil recoveries and gas oil ratios
are significantly influenced by the transition of flow from co-current to counter-current mode.
Accounting for counter-current relative permeability results in higher oil recovery and lower GOR
in all cases. The extent to which the conventional flow model and the flow model that accounts for
vi

counter-current relative permeability differ, relies on the properties of the reservoir (i.e. vertical to
horizontal permeability ratio) and the history of production and flooding of the reservoir (e.g.
initial fluid saturation in the reservoir at the beginning of the WAG/SWAG flood).
Next, we focus on the phase behavior of fluids in unconventional reservoirs. Despite the large
potential of unconventional resources, many unknowns still exist regarding the physics controlling
the extraction processes in these settings. These include accurate representation of phase
equilibrium in tight formations and effective implementation of relevant models in simulation
tools. When a fluid is confined in pore spaces of nanometer size, significant interfacial curvatures
may occur that can result in large capillary pressures between the liquid and vapor phases. The
pressure difference between the two phases will likely affect the vapor-liquid equilibrium state.
Previous efforts have shown that this effect is negligible for conventional reservoirs (with pores in
the micron range) and current commercial reservoir simulators commonly ignore the effect of
capillary pressure in the VLE calculations. However, experimental and modeling efforts have
shown that ignoring capillary pressure in the VLE calculations will not be a valid approximation
for unconventional (tight) reservoirs. We analyze the numerical aspects of including capillarity
phenomena in VLE calculations, in an effort, to arrive at robust and efficient algorithms for
stability analysis that can be used in compositional modeling/simulation of unconventional
reservoirs.
Next, we aim our focus on molecular diffusion and mass transfer during gas injection in fractured
reservoirs. Molecular diffusion can play a significant role in oil recovery during gas injection in
fractured reservoirs, especially when matrix permeability is low and fracture intensity is high.
Diffusion of gas components from a fracture into the matrix extracts oil components from matrix
vii

and delays, to some extent, the gas breakthrough. This in turn increases both sweep and
displacement efficiencies.
In current simulation models, molecular diffusion is commonly modeled using a classical Fick’s
law approach with constant diffusion coefficients. In the classical Fick’s law approach, the
dragging effects (off-diagonal diffusion coefficients) are neglected. In addition, the gas-oil
diffusion at the fracture-matrix interface is normally modeled by assuming an average composition
at the interface which does not have a sound physical basis.
We present a dual-porosity model in which the generalized Fick’s law is used to describe molecular
diffusion to account for dragging effects. Gas-oil diffusion across the fracture-matrix interface is
modeled based on film theory in which the gas in fracture and oil in the matrix are assumed to be
at equilibrium at the interface. A novel and consistent matrix-fracture (M-F) transfer function is
introduced for gas-oil diffusion based on film theory. Diffusion coefficients are calculated using
the Maxwell-Stefan approach and are pressure, temperature and composition dependent.
Finally, we will discuss a common problem in coarse scale modeling of molecular diffusion in
naturally fractured reservoirs. In commercial reservoir simulation tools, transfer rates between the
flowing and the stagnant domain in a dual-porosity model are modeled using the Warren and Roots
model, in which no attempt is made to solve diffusion equations within each block and quasi-
steady state transfer is assumed between the two domains. Previous researchers have shown that
such a simplification does not correctly represent the mass transfer between the domains in a
dynamic displacement.
A more accurate approach to modeling transfer rates in such reservoirs would be to further
discretize the matrix blocks but such an approach lacks computational efficiency for the purposes
viii

of large scale reservoir simulation, because orders of magnitude more computational cells are
required for this approach.
In the final chapter, we present a semi-analytical approach to model multicomponent molecular
diffusion within each matrix block in a coarse-scale simulation model and develop equations for
time-dependent transfer functions between the flowing and stagnant domains. The time-
dependency of the transfer functions are characterized based on the composition of domains.
Generalized Fick’s law is used to describe the diffusive fluxes to account for dragging effects
between components. Analytical and semi-analytical solutions to the multicomponent mass
transfer equations are obtained using the linearized theory of Toor and eigenvalue decomposition
of the diffusion coefficient matrix.
To demonstrate the accuracy of the proposed semi-analytical approach, various examples are
presented. In these examples the results of the suggested approach are compared to the Warren-
root model, analytical solution and high-resolution fine-scale models for 1-D linear and spherical
geometries for a variety of fluid compositions. We observe that our approach (coarse-scale model
using time-dependent transfer functions) accurately represents the analytical solution, with far
lower CPU time requirement compared to high-resolution fine-scale models. Furthermore, we
investigate experimental measurements of mass transfer rates in binary mixtures and further
validate the precision of the proposed model by comparison of the results obtained from our model
with the experimental data.
The novel and consistent matrix-fracture transfer function for coarse-scale models, based on semi-
analytical solutions to multicomponent mass transfer, overcome the inaccuracies of conventional
treatment of diffusion in dual porosity models, while eliminating the need for further discretization
ix

of the matrix blocks. This allows for accurate, yet efficient simulation of mass transfer in fractured
reservoirs.

x

Table of Contents
Dedication ....................................................................................................................................... ii
Acknowledgement ......................................................................................................................... iii
Abstract ........................................................................................................................................... v
List of Figures .............................................................................................................................. xiii
List of Tables ............................................................................................................................... xvi
Chapter 1 Introduction .................................................................................................................... 1
1.1. Background ...................................................................................................................... 1
1.2. Motivation ........................................................................................................................ 3
1.3. Organization of the Manuscript ........................................................................................ 5
Chapter 2 Literature Review ........................................................................................................... 7
2.1. Mathematical Modeling ................................................................................................... 7
2.2. Viscous Coupling in Multiphase Flow in Porous Media ................................................. 8
2.3. Counter-Current Relative Permeability ............................................................................ 9
2.4. Phase Behavior in Unconventional Reservoirs .............................................................. 12
2.5. Mass Transfer in Dual Porosity Systems ....................................................................... 15
Chapter 3 Simultaneous Water and Gas Injection Processes ........................................................ 18
3.1. Introduction .................................................................................................................... 18
3.2. Methodology .................................................................................................................. 20
3.2.1. Two-phase Relative Permeability ........................................................................... 20
3.2.2. Counter-Current Flow Calculations ........................................................................ 21
3.2.3. Three-Phase Relative Permeability ......................................................................... 21
3.2.4. Hysteresis ................................................................................................................ 23
3.2.5. IFT Scaling.............................................................................................................. 25
3.3. 2-D Model Setup ............................................................................................................ 26
3.4. Results of the 2-D Model ............................................................................................... 28
3.4.1. Example 1: Post Waterflood EOR .......................................................................... 29
3.4.2. Example 2: Reservoir Initially at Connate Water Saturation .................................. 32
3.5. 3-D Reservoir Model ...................................................................................................... 34
3.6. Results of the 3-D Model ............................................................................................... 35
xi

3.6.1. Example 3: Post Waterflood EOR .......................................................................... 35
3.6.2. Example 4: Reservoir Initially at Connate Water Saturation .................................. 37
3.7. A Different Three Phase Relative Permeability Model ................................................. 41
3.7.1. Example 5: Post Waterflood EOR, Interpolation ................................................... 42
3.7.2. Example 6: Reservoir at Connate Water Saturation, Interpolation ......................... 44
3.8. Discussion ...................................................................................................................... 48
3.9. Conclusions .................................................................................................................... 50
Chapter 4 Water Alternating Gas Injection Processes .................................................................. 52
4.1. Introduction .................................................................................................................... 52
4.2. Methodology .................................................................................................................. 54
4.3. Model Setup ................................................................................................................... 54
4.4. Results ............................................................................................................................ 56
4.4.1. Effect of Slug Size .................................................................................................. 60
4.4.2. Effect of Initial Water Saturation ............................................................................ 61
4.4.3. Impact of Counter-Current Relative Permeability Curves ...................................... 66
4.4.4. Use of Static Anisotropic Relative Permeability to Represent Counter-Current
Flow Effects ........................................................................................................................... 67
4.4.5. Impact of Capillary Pressure ................................................................................... 69
4.5. Summary and Discussion ............................................................................................... 71
4.6. Conclusions .................................................................................................................... 75
Chapter 5 Stability Analysis for Multicomponent Mixtures including Capillary Pressure .......... 76
5.1. Introduction .................................................................................................................... 76
5.2. Methodology .................................................................................................................. 78
5.2.1. Formulation of the stability analysis problem ........................................................ 80
5.2.2. Solution Methods .................................................................................................... 82
5.3. Results ............................................................................................................................ 86
5.3.1. Fluid Description .................................................................................................... 86
5.3.2. The Phase Envelope ................................................................................................ 87
5.3.3. Stability Testing ...................................................................................................... 88
5.3.4. Convergence Behavior ............................................................................................ 91
5.3.5. Effect of Pore Radius .............................................................................................. 95
5.3.6. Effect of Composition ............................................................................................. 98
5.4. Summary ...................................................................................................................... 101
xii

Chapter 6 Diffusion and Matrix-Fracture Interactions During Gas Injection in Fractured
Reservoirs ................................................................................................................................... 102
6.1. Introduction .................................................................................................................. 102
6.2. Mathematical Model .................................................................................................... 104
6.3. Molecular Diffusion ..................................................................................................... 107
6.3.1. Gas-Oil Diffusion .................................................................................................. 107
6.3.2. Gas-Gas and Oil-Oil Diffusion ............................................................................. 111
6.4. Transfer Functions........................................................................................................ 112
6.5. Calculation Examples ................................................................................................... 114
6.5.1. Example 1 ............................................................................................................. 117
6.5.2. Example 2 ............................................................................................................. 120
6.6. Discussion .................................................................................................................... 123
6.7. A Different Approach to Find Interface Compositions ................................................ 126
6.8. Conclusions .................................................................................................................. 128
Chapter 7 Coarse Scale Modeling of Multicomponent Diffusive Mass Transfer in Dual-porosity
Models......................................................................................................................................... 129
7.1. Introduction .................................................................................................................. 129
7.2. Methodology ................................................................................................................ 132
7.2.1. Binary Fluid System ............................................................................................. 133
7.2.2. General Formulation for Multicomponent Fluid System...................................... 136
7.2.3. Extension to Other Geometries ............................................................................. 139
7.3. Results and Discussion ................................................................................................. 139
7.3.1. Example 1 – Binary Fluid System in Spherical Geometry ................................... 140
7.3.2. Example 2 – Binary Fluid System in Linear Geometry ........................................ 142
7.3.3. Example 3 – Experimental Data from a CO2 Diffusion Experiments .................. 144
7.3.4. Example 4 – Multicomponent Fluid System in Spherical Geometry .................. 146
7.3.5. Implementation in Reservoir Simulation .............................................................. 149
7.4. Conclusions .................................................................................................................. 150
Chapter 8 Summary and Future Research Directions ................................................................. 151
8.1. Summary ...................................................................................................................... 151
8.2. Future Research ............................................................................................................ 154
Nomenclature .............................................................................................................................. 156
References ................................................................................................................................... 159
xiii

List of Figures
Figure 1.1 Classification of oil recovery methods .......................................................................... 2
Figure 3.1 Two-phase relative permeability curves for co-current and counter-current flow ...... 21
Figure 3.2 Relative permeability and IFT scaling ........................................................................ 26
Figure 3.3 2-D reservoir model for SWAG injections .................................................................. 28
Figure 3.4 Gas saturation distribution after 600 days of injection for different vertical
permeabilities in post water-flood tests for 2-D reservoir model, SWAG injection process ....... 30
Figure 3.5 Oil recovery vs. time for different vertical permeabilities in post water-flood tests for
2-D reservoir model, SWAG injection processes ......................................................................... 30
Figure 3.6 Gas saturation distribution after 600 days of injection for different vertical
permeabilities in reservoir with connate water saturation for 2-D reservoir model, SWAG
injection process ............................................................................................................................ 33
Figure 3.7 Oil recovery vs. time for different vertical permeabilities in res. with connate water
saturation for 2-D reservoir model, SWAG injection processes ................................................... 33
Figure 3.8 3-D reservoir model (left) and log of permeability distribution (right) for SWAG
injections ....................................................................................................................................... 34
Figure 3.9 Gas, oil and water saturation distribution after 800 days (1.5 PVI), post waterflood, 3-
D reservoir model, kv = 0.1 kh, SWAG injection process. .......................................................... 35
Figure 3.10 Gas, oil and water saturation distribution after 800 days (1.5 PVI), post waterflood,
3-D reservoir model, kv = kh, SWAG injection process. ............................................................. 37
Figure 3.11 Gas, oil and water saturation distribution after 800 days (1.5 PVI), connate water
saturation, 3-D reservoir model, kv = 0.1 kh, SWAG injection process. ..................................... 38
Figure 3.12 Gas, oil and water saturation distribution after 800 days (1.5 PVI), connate water
saturation, 3-D reservoir model, kv = kh, SWAG injection process. ........................................... 39
Figure 3.13 Recovery vs. PVI for all cases, 3-D reservoir model, SWAG injection process ...... 40
Figure 3.14 GOR (Sm3/Sm3) vs. PVI for all cases, 3-D reservoir model, SWAG injection
process........................................................................................................................................... 41
Figure 3.15 Gas, oil and water saturation distribution after 800 days (1.5 PVI), post water flood,
3-D reservoir model, kv = 0.1 kh, interpolation model for rel perm, SWAG injection process ... 43
Figure 3.16 Gas, oil and water saturation distribution after 800 days (1.5 PVI), post water flood,
3-D reservoir model, kv = kh, interpolation model for rel perm, SWAG injection process ......... 44
Figure 3.17 Gas, oil and water saturation distribution after 800 days, connate water saturation, 3-
D reservoir model, kv = 0.1 kh, interpolation model for rel perm, SWAG injection process ...... 45
Figure 3.18 Gas, oil and water saturation distribution after 800 days (1.5 PVI), connate water
saturation, 3-D reservoir model, kv = kh, interpolation model for rel perm, SWAG injection
process........................................................................................................................................... 46
Figure 3.19 Recovery vs. PVI for all cases, 3-D reservoir model, interpolation model for rel
perm, SWAG injection process ..................................................................................................... 47
Figure 3.20 GOR vs. PVI for all cases, 3-D reservoir model, interpolation model for rel perm,
SWAG injection process ............................................................................................................... 48
xiv

Figure 4.1 3-D reservoir model (left) and log of permeability distribution (right), WAG
injections ....................................................................................................................................... 55
Figure 4.2 Gas, oil and water saturation distribution for 0.1 PV slug size and Swi=0.7 after 800
days (1.5 PVI) of injection, Top layer .......................................................................................... 58
Figure 4.3 Gas, oil and water saturation distribution for 0.1 PV slug size and Swi=0.7 after 800
days (1.5 PVI) of injection, Middle layer ..................................................................................... 59
Figure 4.4 Gas, oil and water saturation distribution for 0.1 PV slug size and Swi=0.7 after 800
days (1.5 PVI) of injection, Bottom layer ..................................................................................... 59
Figure 4.5 Gas, oil and water saturation distribution for 0.1 PV slug size and Swi=0.2 after 800
days of injection, Top layer .......................................................................................................... 62
Figure 4.6 Gas, oil and water saturation distribution for 0.1 PV slug size and Swi=0.2 after 800
days of injection, Middle layer ..................................................................................................... 63
Figure 4.7 Gas oil and water saturation distribution for 0.1 PV slug size and Swi=0.2 after 800
days of injection, Bottom layer ..................................................................................................... 63
Figure 4.8 Recovery vs. PVI for all calculation examples ........................................................... 64
Figure 4.9 GOR vs. PVI for all cases, 3-D reservoir model, WAG injection process.................. 65
Figure 4.10 Results for two different counter-current relative permeability models, left: recovery
vs. PVI, right: GOR vs. PVI ......................................................................................................... 67
Figure 4.11 Comparison of anisotropic and dynamic relative permeability modeling ................. 68
Figure 4.12 J-function vs. saturation for primary drainage .......................................................... 70
Figure 4.13 Effect of capillary pressure on WAG simulation ...................................................... 70
Figure 5.1 Phase envelope of the 15-component oil with and without capillary pressure (rc =
10nm). Triangle = critical point. ................................................................................................... 88
Figure 5.2 Reduced tangent plane distance vs. pressure at constant temperatures of T=383.15K
(top), T=700K (middle) and T=780K (bottom). ........................................................................... 90
Figure 5.3 Convergence behavior of direct substitution, accelerated substitution and
minimization with lagged pw updates (Minimization 1) and minimization with the full Jacobian
(Minimization 2) at a) T=383.15 K, p = 83 atm and b) T=700 K, p = 119 atm for rc=10 nm ..... 92
Figure 5.4 Convergence of the capillary pressure for direct substitution, accelerated substitution,
minimization with lagged pc updates (Minimization 1) and minimization with the full Jacobian
(Minimization 2) at T=383.15K, p=83atm and T=700K, p=119atm for rc=10 n ......................... 94
Figure 5.5 Detection of phase boundaries by stability analysis at different capillary radii
(rc≤10nm). Triangle = critical point. ............................................................................................ 96
Figure 5.6 Detection of phase boundaries via stability analysis at different capillary radii
(rc>10nm). Triangle = critical point. ............................................................................................ 97
Figure 5.7 Calculated saturation pressure vs capillary radius at T=383.15 K .............................. 98
Figure 5.8 Detection of phase boundaries via stability analysis for a condensate system (rc=10
nm). ............................................................................................................................................... 99
Figure 5.9 Reduced tangent plane distance vs. pressure at constant temperature for T=250K
(top), T=300K (middle) and T=365K (bottom) for condensate fluid ......................................... 100
Figure 6.1 Schematic of two neighboring gridblocks containing oil and gas phases. ................ 107
Figure 6.2 Composition profiles near the interface during inter-phase mass transfer, after Taylor
and Krishna (1993). .................................................................................................................... 108
xv

Figure 6.3 Oil recovery (left) and producing GOR (right) versus time for example 1 when
diffusion is not considered (red), when diffusion is modeled using the conventional approach
(black) and when diffusion is modeled using our approach (blue). ............................................ 118
Figure 6.4 Oil recovery (left) and producing GOR (right) versus time for example 1 when
diffusion is not considered (solid red), when our approach is used (solid blue) and when our
approach is used and the off-diagonal diffusion coefficients are neglected (dashed blue). ....... 119
Figure 6.5 The CO2 composition (mole %) in the stagnant domain after 10 years from
calculations with diffusion based on the proposed approach (left) and absolute difference (mole
%) between the proposed approach and the conventional approach (right). .............................. 120
Figure 6.6 Oil recovery (left) and producing GOR (right) versus time for example 2 when
diffusion is not considered (red), when diffusion is modeled using the conventional approach
(black) and when diffusion is modeled using the proposed approach (blue). ............................ 121
Figure 6.7 Oil recovery (left) and producing GOR (right) versus time for example 2 when
diffusion is not considered (solid red), when our approach is used (solid blue) and when our
approach is used but the off-diagonal diffusion coefficients are neglected (dashed blue). ........ 122
Figure 6.8 Saturation variation in a matrix block (dotted line) and its corresponding fracture
(solid line) near the injection well for example 1 as calculated from the proposed approach (top)
and the conventional approach (bottom). .................................................................................... 123
Figure 6.9 Composition paths (mass fraction) of a matrix block near the injection well for
example 1 during the first 10 years as calculated using the conventional approach (solid line) and
the proposed approach (dashed line). .......................................................................................... 125
Figure 6.10 Oil recovery (left) and producing GOR (right) versus time for example 1 when
diffusion is not considered (solid red), when our approach presented in sec. 6.3.1 is used (solid
blue) and when our approach presented in sec. 6.7 is used (dashed blue) .................................. 127
Figure 7.1 Average concentration vs time (left) and correction factor vs time (right) for binary
fluid system in spherical geometry ............................................................................................. 142
Figure 7.2 Average concentration vs time (left) and correction factor vs time (right) for binary
fluid system in spherical geometry ............................................................................................. 143
Figure 7.3 Schematic of the assumed 1-D diffusion process ...................................................... 145
Figure 7.4 Amount of CO2 in vapor phase vs time ..................................................................... 146
Figure 7.5 Average composition of CO2 (top), nC4 (middle) and nC12 (bottom) in matrix block
vs time ......................................................................................................................................... 148

xvi

List of Tables
Table 3.1: Fluid description and compositions of oil and injection gas for SWAG injections .... 27
Table 5.1 Fluid properties and compositions ................................................................................ 86
Table 5.2 Non-zero binary interaction coefficients ...................................................................... 87
Table 6.1 Parameters for Corey-type relative permeabilities ..................................................... 115
Table 6.2 Fluid properties ........................................................................................................... 115
Table 7.1 Fluid description ......................................................................................................... 147

1

Chapter 1 Introduction
1.1. Background
Over the last few decades, the demand for oil and gas has been increasing while the discovery of
new oil and gas resources is declining. As a result, improved oil recovery (IOR) techniques, and
production from unconventional reservoirs play an important role in meeting the energy needs, in
the years to come.
The IOR techniques include all effort to produce oil after the primary recovery of oil, e.g.
intelligent well monitoring and control, secondary recovery methods and enhanced oil recovery
(EOR) processes. Secondary recovery methods include injection of water in an associated aquifer,
or injection gas in the gas cap to maintain the reservoir pressure. EOR processes include gas
injection (hydrocarbon, N2 and CO2), thermal methods (steam flooding and in-situ combustion)
and chemical flooding (polymer, alkaline and surfactant). Fig. 1.1 depicts the classification of oil
production techniques.
Gas injection and more particularly, CO2 injection is one of the most appealing EOR techniques
and has gained tremendous attention over the past decades. Availability of CO2 and easy
achievement of miscible conditions with the reservoir fluid are the main reasons of the popularity
of CO2 injection processes in the U.S and worldwide.
The main problem with gas injection processes is the low sweep efficiency experienced during
injection due to the unfavorable mobility ratio. Water is usually injected either at the same time as
gas or in alternating intervals in order to improve the sweep efficiency of the gas injection process.
2

Water Alternating Gas (WAG) injection combines the improved displacement efficiency of the
gas flooding with an improved macroscopic sweep by the injection of water (Christensen et al.,
2001): WAG injection targets oil that was not contacted (unswept zones) by previous oil recovery
methods.

Figure 1.1 Classification of oil recovery methods, adapted from the Oil and Gas Journal, Apr. 23, 1990
WAG floods are classified into miscible and immiscible floods. Miscible WAG floods are more
common since better recoveries are observed/reported in floods that approach miscible or near
miscible conditions between the injected gas and the reservoir oil. Usually, an initial slug of gas is
injected which mixes with the oil to some extent, creating near-miscible conditions. This slug is
then followed by a number of slugs of water and gas. If the initial slug of gas is particularly large
3

and the following slugs are small in size, the process is called hybrid WAG. Another tested process
consists of injecting gas and water simultaneously, which is called Simultaneous Water and Gas
Injection (SWAG) (Christensen et al. 2001).
Tight oil and shale reservoirs contribute significantly to oil and gas production in the US.
Currently, unconventional resources account for over 4 MMBOE per day of the US fossil fuel
production. Despite the large potential of unconventional resources, many unknowns still exist
regarding the physics controlling the extraction processes in these settings. These include accurate
representation of phase equilibrium in tight formations and effective implementation of relevant
models in simulation tools.
1.2. Motivation
The first step in any EOR project is the estimation of potential incremental recovery from the
project. In order to obtain an estimate of the hydrocarbon recovery, compositional reservoir
simulation is commonly used. With the help of compositional simulation, one can perform an
economic evaluation of the project. Moreover, compositional simulators are used to design EOR
projects before implementation in the field. Successful design of EOR projects is, accordingly,
dependent on the accuracy by which the available simulation tools can represent the physics that
govern the displacement behavior in a reservoir.
Over the past decades, a tremendous research effort has been dedicated to better understand the
physics of WAG injection processes and to better represent it in simulators (Larsen and Skauge
1999, Li et al. 2003, Yan et al. 2004, Bhambri and Mohanty 2005, Brown et al. 2013). However,
due to the high complexity of the nature of the multiphase, multicomponent displacement
processes experienced in a WAG process, the interaction of flow and phase behavior in porous
4

media and its proper modeling is not completely understood. This has motivated us to investigate
the accuracy of some physical models that are frequently used in compositional simulation of
WAG injection in an effort to enhance the simulation and design of these processes (Chapters 3
and 4).
When the pore sizes are in the submicron range (such as in the Bakken formation), the confined
fluid in the tight pores behave differently than unconfined bulk fluids in terms of phase equilibrium
and phase properties. Also, capillary pressure values in nano-size pores reaches high values that
ignoring its effect on vapor-liquid equilibrium calculations might result in inconsistent simulation
results. As a result, precise evaluation of phase behavior and the effects of capillary pressure on
vapor-liquid equilibrium in unconventional reservoirs is expected to result in more accurate
simulation and development of such resources. This has motivated us to develop algorithms that
can efficiently and robustly predict the behavior of reservoir fluids in confined spaces. Such
algorithms can be easily implemented in reservoir simulators and improve the accuracy of
production and GOR predictions (Chapter 5).
In current simulation models, molecular diffusion which is an important recovery mechanism
during gas injection in fractured reservoirs, is represented by simplified models. A classical Fick’s
law approach (which neglects interactions among components) with constant diffusion coefficients
is commonly used for molecular diffusion. In addition, the gas-oil diffusion at the fracture-matrix
interface is normally modeled by assuming an average composition at the interface which does not
have a sound physical basis. This has motivated us to develop and implement a dual-porosity
model that uses sophisticated physical models for molecular diffusion and matrix-fracture
interactions during gas injection in fractured reservoirs (Chapter 6).

5

1.3. Organization of the Manuscript
The remainder of this manuscript is organized as follows. In the next chapter, a brief literature
review of compositional simulation is presented.
Chapter 3 discusses the significance of including the dynamic relative permeability (mobility) in
the simulation of simultaneous water and gas injection processes. A number of examples in both
2-D and 3-D settings are presented to evaluate the effect of key reservoir parameters on prediction
results.
Chapter 4 describes the impact of considering relative permeability in water alternating gas
injection processes. An anisotropic 3-D reservoir model is used and several cases are presented to
demonstrate the impact of initial water saturation, slug size, different interpretations of counter-
current flow, and capillary pressure.
In Chapter 5, we start from fundamental thermodynamic definitions (e.g. Helmholtz energy) and
develop formulae for two-phase equilibrium in the presence of high capillary pressure values. We
then present new algorithms for stability analysis in unconventional reservoirs and demonstrate
the accuracy of the proposed algorithms by comparing the obtained results with separate phase
boundary calculations.
In Chapter 6, we present a dual-porosity model for gas injection in fractured reservoirs in which
the generalized Fick’s law is used for molecular diffusion; and gas-oil diffusion at fracture-matrix
interface is accounted for using the film theory. Using field-scale numerical examples, our
approach is compared with the conventional approach in which the molecular diffusion is modeled
using the classical Fick’s law with constant diffusion coefficients and the gas-oil diffusion is
modeled by assuming average composition at the interface.
6

Chapter 7 explores a common problem in using dual-porosity models for modeling of molecular
diffusion in naturally fractured reservoirs. The mass transfer between matrix and fracture in a dual-
porosity model is usually modeled using Warren and Roots formulation, which assumes that a
pseudo-steady state flow is developed. Starting from an exact analytical solution to the molecular
diffusion, we present an approximation of the analytical solution that does not rely on the Warren
and Root assumption of pseudo-steady state flow. A time-dependent transfer function is presented
that accounts for both early- and late-time behavior of matrix-fracture interactions. A correction
factor to the Warren and Root model is proposed, and calculations for binary fluid component
systems and multicomponent systems are presented. Numerical and experimental examples are
presented to establish the accuracy of the suggested correction factor.
7

Chapter 2 Literature Review
2.1. Mathematical Modeling
The description of flow and transport in a porous medium is typically stated by writing the mass
conservation equation (also referred to as continuity equation). The multicomponent multiphase
continuity equation can be written as
𝜕 𝜕𝑡
𝜙 ∑𝑥 𝑖𝑗
𝜌 𝑗 𝑆 𝑗 𝑛𝑝
𝑗 =1
+𝛻 .∑𝑥 𝑖𝑗
𝜌 𝑗 𝑣⃗
𝑗 𝑛𝑝
𝑗 =1
=0, 𝑖 =1,…,𝑛𝑐 , (2.1)
where xij is the mole fraction of component i in phase j, 𝜙 is the porosity, Sj is the saturation of
phase j, ρj is the molar density of phase j, 𝑣⃗
𝑗 is the velocity of phase j, t represents time, np represents
the number of phases, and nc is the number of components.
It should be noted that the effects of dispersion, chemical reaction, adsorption and temperature
variation have been neglected in derivation of Eq. (2.1) (Lake, 1989; Orr, 2007). The extended
Darcy equation is commonly used to calculate phase velocities in the above equation (replaces the
momentum balance):
𝑣⃗
𝑗 =−
𝑘 𝑘 𝑟𝑗
𝜇 𝑗 (𝛻 𝑃 𝑗 +𝜌 𝑚𝑗
𝑔⃗ℎ) , (2.2)
where 𝑘 𝑟𝑗
is the relative permeability of phase j, 𝜇 𝑗 is the viscosity of phase j, 𝜌 𝑚𝑗
is the mass
density of phase j, 𝑃 𝑗 is the pressure of phase j, 𝜌 𝑚𝑗
𝑔 ℎ is the flow potential due to gravity and k
represents the absolute permeability. The relationship between the phase pressures is described by
capillary pressure.
8

𝑃 𝑗 −𝑃 𝑘 =𝑃 𝑐𝑘𝑗 , 𝑗 =1,…,𝑛𝑝 , 𝑘 =1,…,𝑛𝑝 , 𝑘 ≠𝑗 (2.3)
The phases are assumed to be in thermodynamic equilibrium which can be represented by equality
of chemical potentials
𝜇 𝑖𝑗
=𝜇 𝑖𝑘
, 𝑗 =1,…,𝑛𝑝 , 𝑘 =1,…,𝑛𝑝 , 𝑘 ≠𝑗 (2.4)
It can be shown that the equality of chemical potentials leads to the equality of fugacities, with the
fugacities calculated using an equation of state (EOS) for reservoir fluids (Danesh, 1998; Orr,
2007). The phase saturations and mole fractions in each phase must sum to one
∑𝑆 𝑗 𝑛𝑝
𝑗 =1
=1 (2.5)
∑𝑥 𝑖𝑗
𝑛𝑐
𝑖 =1
=1 , 𝑗 =1,…,𝑛𝑝 (2.6)
The above equations provide us with enough information to obtain the solution to a flow/transport
problem in a porous medium (Orr, 2007).

2.2. Viscous Coupling in Multiphase Flow in Porous Media
As mentioned in the previous section, commercial simulation techniques use the extended version
of Darcy’s equation, proposed by Muskat (1937) to calculate phase velocities. Experimental
(Maini et al. 1990; Bourbiaux and Kalaydjian 1990; Lelievre 1966; Odeh 1959) and theoretical
studies (Kalaydjian 1987; Whitaker 1986; de la Cruz and Spanos 1983) have shown that the use
of Muskat’s relative permeability may not be accurate.
9

Kalaydjian (1987) and de la Cruz and Spanos (1983) use volume averaging of Navier-Stokes
equation to arrive at a multi-phase equation of fluid flow in porous media. For a three-phase
system, the volume averaging of Navier-Stokes equation can be presented as follows
𝑢 = −𝐾
[

𝑘 𝑟 (1,1)
𝜇 1
⁄
𝑘 𝑟 (1,2)
𝜇 2
⁄
𝑘 𝑟 (1,3)
𝜇 3
⁄
𝑘 𝑟 (2,1)
𝜇 1
⁄
𝑘 𝑟 (2,2)
𝜇 2
⁄
𝑘 𝑟 (2,3)
𝜇 3
⁄
𝑘 𝑟 (3,1)
𝜇 1
⁄
𝑘 𝑟 (3,2)
𝜇 2
⁄
𝑘 𝑟 (3,3)
𝜇 3
⁄
]

[

𝜕 𝑃 1
𝜕𝐿
𝜕 𝑃 2
𝜕𝐿
𝜕 𝑃 3
𝜕𝐿
]

. (2.7)
In this formulation, kr(i,i) is the relative permeability of each phase and kr(i,j) is the relative
permeability of each phase due to the flow of the other two phases. These coefficients represent
the viscous drag between the phases. In fact, there are two types of drag; a) one due to the flow of
fluid over solid particles/surfaces and b) the other is the momentum transfer between the fluid
phases.
In the conventional methods of modeling multi-phase flow, it is assumed that the phases act
independently and viscous coupling is neglected in many settings. During counter-current flow,
however, the viscous coupling is more significant and the non-zero off-diagonal terms in the
relative permeability tensor cannot be ignored.

2.3. Counter-Current Relative Permeability
Lelievre (1966) first conducted experiments to compare co-current and counter-current relative
permeabilites. He performed steady-state measurements of co- and counter-current relative
permeabilities for different viscosity ratios and found that in all cases, relative permeability in the
counter-current flow settings is less than the relative permeability in co-current settings. Much
10

larger differences were observed for the wetting phase than for the non-wetting phase and in some
cases, the reduction in the relative permeability of the wetting phase was as high as 50%.
Bourbiaux and Kalaydjian (1990) conducted experiments of spontaneous imbibition in co-current
and counter-current experiments and observed lower recoveries in the case of counter-current flow
settings. Simulation of their work later revealed that a minimum reduction of 30% of the relative
permeabilities of the co-current flow setting is needed to match the data from counter-current flow
experiments. This reduction was then attributed to the differences in the viscous coupling between
the two flow settings.
Eastwood and Spanos (1991) conducted two-phase flow experiments based on previous modeling
done by de la Cruz and Spanos (1983), and their experiments revealed that relative permeability
in counter-current flow is less than that of co-current flow. They concluded that relative
permeabilities for oil in counter-current flow is always less than those for co-current flow setting.
The work of Bensten and Manai (1993) was focused on the measurement of the independent
generalized permeability coefficients as they appear in equation (2.5). In order to do so, they
performed steady state co-current and counter-current relative permeability measurements. They
used a two-phase oil-water system in their experiments and used a microwave approach to measure
in-situ saturations. They observed that counter-current relative permeabilities are less than co-
current relative permeabilities for all experiments. Counter-current relative permeabilities were
lower by as much as 25% for water (wetting phase) and 20% for oil (non-wetting phase). They
arrived at the conclusion that generalized relative permeability coefficients lie somewhere between
co-current and counter-current relative permeabilities.
11

Al-Wadahi et al. (2000) studied counter-current flow in a three-phase system driven by gravity
and capillary forces. They used water (at irreducible saturation), benzyl alcohol and decane as their
three phases, but only two phases (benzyl alcohol and decane) were mobile in the counter-current
flow experiments. They used a commercial simulator to simulate their experiments and applied
artificial neural networks to derive a set of relative permeability and capillary pressure functions.
They observed that for the counter-current flow experiment, relative permeabilities for both of the
mobile phases were reduced significantly from their co-current flow values. The existence of a
third phase in the experiment caused the reductions in relative permeability values to be higher
than those reported for two-phase experiments. In their work, relative permeabilities were reduced
from their co-current values by a factor of 2 and 5, for benzyl alcohol and decane, respectively.
Langaas and Papatzacos (2001) studied the effect of saturation, wettability, driving force and
viscosity ratio on steady state co-current and counter-current relative permeability based on
displacement calculations in a uniform pore structure. Based on their results, counter-current
relative permeabilities are always lower than co-current relative permeabilities. They attributed
the difference to viscous coupling and different levels of capillary forces.
Li et al. (2005) continued the work of Al-Wadahi et al. (2000) and included hysteresis in relative
permeabilities and capillary pressure in their numerical simulations and arrived at the same
conclusions.
Javaheri and Jessen (2011) conducted numerical simulation of counter-current flow for CO2
sequestration in saline aquifers and concluded that counter-current flow is the dominant flow mode
during the injection of CO2 into the saline aquifer and hence that viscous coupling can play an
important role in determining the migration dynamics of CO2 in the subsurface.
12

Based on the above mentioned experimental and modeling studies, we observe that in the case of
counter-current flow setting, relative permeabilities are considerably less than in the case of co-
current flow. This can be attributed to the fact that pressure gradients in of co-current flow are in
the same direction for all of the flowing phases while the pressure gradients in counter-current
flow are in opposite directions. In other words, viscous coupling has a negative impact on overall
mobility in counter-current flow settings.
All the above-mentioned work is based on two-phase flow, even in the experiments of Al-Wadahi
et al. (2000) where only two of the three phases were mobile and the third phase (i.e. water) was
at an irreducible saturation and hence did not flow. To the best of our knowledge, no previous
studies exist that focus on the counter-current flow setting for three-phase flow. The next two
chapters of this manuscript will focus on the simulation of three-phase flow systems that exhibit
counter-current flow and the impact of accounting for the differences between co-current and
counter-current flow in compositional simulators that are used in the design of EOR projects.

2.4. Phase Behavior in Unconventional Reservoirs
Previous work has shown that fluid properties in a confined nanopore space differ from that of
bulk measurements. (Alharthy et al. 2013) In conventional methods of compositional simulation,
it is assumed that the thermodynamic properties of the fluid in the pore space are similar to that of
unconfined fluids. However, when the pore size is reduced to submicron-size, significant deviation
is observed from the behavior of unconfined fluids. Alteration of fluid properties and transport
processes in tight porous materials can be attributed, in part, to the limited number of molecules
13

present in these pores (Alharthy et al. 2013) rendering bulk (continuum) assumptions inaccurate.
Zarragoicoechea and Kuz (2003) showed that in order to properly account for the behavior of
confined fluids, the critical properties of the components should be altered as a function of the
molecular size to the pore size.
Singh et al. (2009) and Travalloni et al. (2010) have proposed models for the modifications in
critical fluid properties using trigonometric functions of the pore diameter. Singh et al. (2009),
proposed a model for evaluation of relative shift in critical pressure and temperature for methane,
n-butane, and n-octane in mica and graphite split-pore systems with a range of pore sizes from 0.8
nm to 5 nm.
Based on Singh’s (2009) work, Alharthy et al. (2013) proposed the below equations for calculation
of changes in critical pressure and temperature with respect to pore diameter. The cutoff pore
diameter in their work was set to 3 nm, meaning that for pore diameters of 3 nm or higher, critical
pressure and temperature of the confined fluid are assumed similar to those of an unconfined fluid.
∆𝑃 𝑐𝑟𝑖𝑡 = −0.4097𝑑 +1.2142 (2.8)
∆𝑇 𝑐𝑟𝑖𝑡 = 1.0983𝑒 (−0.929𝑑 )
(2.9)
In their work, Alharthy et al. (2013) use a classification similar to IUPAC pore size classification
to propose three thermodynamic phase behavior paths for fluids in oil shale samples.
1. Micropores, which have a mode of 2-3 nm pore sizes. This class account for about 3-5 %
of the pores in a sample and is dependent on mineralogy. Fluids in these pores will have a
confined pore phase behavior (as described by equations 2.8 and 2.9).
14

2. Mesopores with a mode of 25 nm account for 10-15% of the pores. Fluid behavior in this
class of pores display a phase behavior, which is the partially shifted phase behavior (e.q.
2.8 and 2.9 above).
3. Macropores have a mode of 11000 nm and represent the majority of the pore space. Fluid
in these pores is considered to be unconfined and has the properties of the bulk fluid.
Alharthy et al. (2013) implemented the above mentioned thermodynamic paths in a simulator and
concluded that including the effect of confined fluid phase behavior promotes more relevant flow
physics in unconventional shale reservoirs.
Many other researchers (e.g. Nojabaei et al., 2012, Firincioglu et al., 2012, Wang et al., 2013 and
Sandoval et al., 2016) have argued that the critical pressure and temperature remain unaffected
from capillarity. Around the critical point, the two phases become identical and the interfacial
tension (and hence the capillary pressure) approaches zero.
The pressure difference between the two phases, caused by the high interfacial curvatures in nano-
sized pores will likely affect the vapor-liquid equilibrium calculations. (Wang et al., 2013)
Researchers have shown that this effect is negligible for conventional reservoirs (Shapiro and
Stenby, 1997; Shapiro et al., 2000). All of the current commercial reservoir simulators ignore the
effect of capillary pressure on the vapor-liquid equilibrium. However, experimental and modeling
works (Du and Chu, 2012; Firincioglu et al., 2012; Wang et al., 2013; Nojabaei et al., 2012;
Nojabaei et al., 2014 and Parsa et al. 2014) have shown that ignoring capillary pressure in the
vapor-liquid equilibrium might not be a valid assumption for unconventional reservoirs.
15

Wang et al. (2013) have extensively studied the effect of high capillary pressure on two-phase split
calculations and have used the results of their work in a dual-porosity simulator. In their work,
they use Leverett-J function to represent the pore size distribution in the reservoir and calculate
capillary pressure values from J-function curves obtained from similar reservoir rocks.
Sandoval et al. (2016) proposed a new algorithm for phase envelope calculations in confined fluids
allowing for negative phase pressures.
In this section we presented several research papers on phase behavior of fluids in unconventional
reservoirs, but it should be noted that the list is not exhaustive. Despite a substantial amount of
research on the topic of phase behavior of confined fluids, the effect of capillarity in stability
analysis has not been extensively studied/documented to the best of our knowledge.

2.5. Mass Transfer in Dual Porosity Systems
In general, viscous, gravity, diffusion and capillary forces can contribute to mass transfer during
flow in porous media. Different investigators have studied mass transfer between high-
permeability fracture and low-permeability matrix in DP systems in the context of gas injection
processes.
da Silva and Belery (1989) showed by compositional simulations that Fick’s molecular diffusion
was very important in gas injection processes in fractured reservoirs because of the large contact
area for diffusion in these systems. They suggested that as the fracture spacing is reduced, the
molecular diffusion becomes more significant and may dominate over viscous forces. Hoteit and
Firoozabadi (2009) studied the role of diffusion during gas injection into fractured oil reservoirs
16

and demonstrated that diffusion can have a significant effect on oil recovery away from the
miscibility pressure.
In gas injection processes in fractured reservoirs low injection rates are desired to avoid early
breakthrough and hence improve the sweep efficiency. This in turn reduces the contributions from
viscous forces in such processes. Contributions from capillary forces will be marginal in miscible
and near-miscible gas injection processes. Findings by Burger and Mohanty (1997) confirm that
mass transfer from capillary forces is minimal in near-miscible displacements.
Li et al. (2000) performed experiments in which CO2 was injected at miscible conditions after
water flooding in artificially fractured cores containing dead oil. They found that gravity drainage
can significantly improve oil recovery after water flooding, and that gravity drainage declines as
initial water saturation increases and matrix permeability decreases. Asghari and Torabi (2008)
observed from their experiments that CO2 gravity drainage improved the oil recovery from a DP
system at miscible conditions. However, they could not match the experimental observations using
a commercial compositional simulator.
Darvish et al. (2006) performed a series of experiments to study the mass transfer between matrix
and fracture during miscible CO2 injection. They observed that diffusion played a more significant
role than gravity drainage in recovering oil from a long vertical tight matrix block when CO2 was
injected at a low rate. Using a commercial compositional simulator, they failed to accurately
reproduce the experimental observations mainly because the simulator did not take into account
the gas-oil diffusion between gas in the fracture and oil in the matrix.
To overcome this limitation of the commercial simulator, they placed a dummy two-phase layer
between the matrix (initially containing oil) and fracture (initially containing gas) in their
17

simulation model to initiate the diffusion between gas and oil. Then they manually changed the
composition of the dummy two-phase layer until a reasonable match between the simulation results
and experimental data was obtained.
However, one does not have to manually insert a dummy two-phase layer between gas and oil, and
adjust its composition in order to model the gas-oil diffusion. A more sophisticated model, known
as the film theory, has been used for gas-oil diffusion in the chemical engineering literature for
several years. According to the film theory, when gas and oil phases are placed in contact, a two-
phase interface forms between them in which the gas and oil are at equilibrium. Moreover, the
molar flux of each component across the interface must be continuous (Taylor and Krishna, 1993).
Moortgat and Firoozabadi (2013) used the film theory for gas-oil diffusion in their mixed finite
element, discontinuous Galerkin approach and obtained reasonable agreement between their
simulation results and experimental data from Darvish et al. (2006).
Vega et al. (2010) performed miscible CO2 injection experiments on a vertical low-permeability
siliceous shale core with an artificial fracture. They found that both diffusion and convection were
significant in their experiments, but failed to match the experimental observations using a
commercial compositional simulator.

18

Chapter 3 Simultaneous Water and Gas Injection Processes
1

3.1. Introduction
The Water Alternating Gas (WAG) injection process is a method of Enhanced Oil Recovery
originally proposed to improve the sweep of the injected gas, by using water to control the mobility
of the displacement and to stabilize the front (Caudle and Dyes, 1958). A similar effect can be
obtained by injection water and gas simultaneously in what is called a Simultaneous Water and
Gas Injection (SWAG) process. Because the microscopic displacement efficiency of oil by gas is
normally better than by water, the SWAG injection process combines an improved microscopic
displacement efficiency of the gas flood with an improved macroscopic sweep by water injection.
As a result, SWAG has been demonstrated to be an efficient EOR process in settings where gravity
segregation or viscous instabilities will result in a poor sweep for a regular gas flood. (Caudle and
Dyes, 1958)
SWAG injection has been tested in a few reservoirs in the United States and has shown promising
results (Christiansen et al., 1998). One of the first studies on modeling the SWAG process was
conducted by Stone (1982). Jenkins (1984) expanded the work of Stone (1982) and presented a
new model to estimate the gravity segregation for uniform and steady state injection of gas and
water. In his model, water and gas are uniformly injected along the thickness of the reservoir.
Rossen et al. (2010) further investigated the model presented by Jenkins and presented a new

1
This chapter was previously published as Sherafati M., Javaheri, M., & Jessen, K. (2014, April). The Role of
Counter-Current Flow in Simultaneous Water and Gas Injection Processes. In SPE Western North American and
Rocky Mountain Joint Meeting
19

approach for calculating the segregation length (distance to full segregation of injected fluids).
They also performed numerical simulations to support/validate their analytical work. They studied
different schemes for simultaneous injection and their effect on the segregation length and overall
displacement performance (recovery). They considered foam for additional mobility control of the
injected gas in their studies and concluded that injecting over the full height of the reservoir, with
injection of water on the top of gas provides for the optimum performance of a SWAG process.
SWAG processes result in complex saturation patterns because the saturation of the three phases
will follow a non-monotonic path as one trace the distance between injectors and producers. This
type of behavior calls for a rigorous description of the relative permeability for the three phases.
During the injection of water and gas, gas migrates to the top portion of the reservoir, while water
segregates to the bottom. The paths and rates of the migration are not only governed by reservoir
parameters such as pressure, temperature and permeability, but are also governed by saturation
functions such as relative permeability and capillary pressure, including hysteresis phenomena in
both.
The previous studies of counter-current flow settings, supporting that relative permeability in the
counter-current flow setting can be significantly different from co-current flow setting, indicate
that flow transitions during EOR operations such as SWAG (or WAG) can have a significant
impact on the displacement dynamics and hence on the optimal design of such EOR operations.
In this section, we investigate the effects of including a mobility reduction due to flow reversals in
displacement calculations, with emphasis on the impact on common EOR process design
parameters. The most important parameter to take into consideration is the oil recovery. Recovery
for SWAG processes is strongly dependent on the travel distance to full segregation of gas and
water. The further the segregation point is from the injection interval, the better the mobility
20

control of the gas in the reservoir, and the better the recovery. We perform simulations of SWAG
processes for a range of cases with and without mobility reductions due to flow reversal in the
vertical direction.
We start with a presentation of the relative permeability model used in the displacement
calculations followed by displacement calculations with a homogeneous 2-D simulation model.
We then compare the displacement dynamics of SWAG processes, with and without impact of
flow reversals, for a range of relevant conditions. We then extend the scope to a heterogeneous 3-
D model and report/compare the calculation results with and without the impact of flow reversals.
Finally, we test a different method for calculating three-phase relative permeabilities, and conclude
the chapter with an analysis and discussion of our results and observations.

3.2. Methodology
In this section, we present our approach to counter-current flow calculations, introduce the three-
phase relative-permeability model, and describe the methods used to represent hysteresis and IFT
scaling.
3.2.1. Two-phase Relative Permeability
Two-phase co-current and counter-current flow relative permeability curves were selected using
Corey type equations. In all cases, relative permeabilities are lower for counter-current flow than
for co-current flow settings. The endpoint values and Corey exponents used to represent counter-
current flow curves were based on previous experimental work (Bourbiaux and Kalaydjian 1990,
21

Al-Wadahi et al. 2000). Fig. 3.1 shows the two-phase relative permeability curves used in this
work.

Figure 3.1 Two-phase relative permeability curves for co-current and counter-current flow
3.2.2. Counter-Current Flow Calculations
In the process of the flow simulation, every grid block is checked for flow reversal. A detection of
flow reversals is done based on the value of the vertical velocity of the phases. We used a flow
direction flag for each phase which indicates if that phase is flowing upward or downward
(indicated by the vertical phase velocity). We use this flag to decide whether counter-current flow
is occurring for each two-phase set (oil-gas, water-oil). If flow reversals are detected, counter-
current flow relative permeabilities are used in the vertical direction. Note that the transition from
co-current flow relative permeability to counter-current relative permeability in our model is
discrete while an initial attempt to introduce continuous transitions was presented by Javaheri
(2014)
3.2.3. Three-Phase Relative Permeability
22

In order to calculate three-phase relative permeability, the modified Stone 1 model, proposed by
Aziz and Settari (1979) was used. This method interpolates between two-phase oil-gas and water-
oil relative permeability curves to generate three-phase relative permeabilities.
The Stone I model suggests that, in a water-wet medium, the relative permeabilities of gas and
water (non-wetting and wetting phases) are the same as their two-phase counterparts as observed
in gas-oil and oil-water displacements, respectively. Relative permeability of the intermediate
wetting phase, oil, depends on the phase configuration within the pore space, which in turn depends
on both water and gas saturations (Pejic and Maini, 2003). The modified Stone 1 model is defined
by the following set of equations:

𝑘 𝑟𝑜
=𝑘 𝑟𝑜𝑤𝑚𝑎𝑥 .𝑆 𝑜𝑛𝑜𝑟 .𝛽 𝑤 .𝛽 𝑔 (3.1)

𝛽 𝑤 =
𝑘 𝑟𝑜𝑤 𝑘 𝑟𝑜𝑤𝑚𝑎𝑥 (1−𝑆 𝑤𝑛𝑜𝑟 )
(3.2)

𝛽 𝑔 =
𝑘 𝑟𝑜𝑔 𝑘 𝑟𝑜𝑤𝑚𝑎𝑥 (1−𝑆 𝑔𝑛𝑜𝑟 )
(3.3)

𝑆 𝑜𝑛𝑜𝑟 =
𝑆 𝑜 −𝑆 𝑜𝑚
1−𝑆 𝑤𝑐
−𝑆 𝑜𝑚
(3.4)

𝑆 𝑤𝑛𝑜𝑟 =
𝑆 𝑤 −𝑆 𝑤𝑐
1−𝑆 𝑤𝑐
−𝑆 𝑜𝑚
(3.5)
23

𝑆 𝑔𝑛𝑜𝑟 =
𝑆 𝑔 1−𝑆 𝑤𝑐
−𝑆 𝑜𝑚
(3.6)

In the above equations, Som is the three-phase residual oil saturation, Swc is the connate water
saturation, krog and krow are the oil relative permeabilities in two-phase oil-gas and oil-water
displacement processes respectively, and krowmax is the maximum oil relative permeability in a two-
phase oil-water displacement process.
3.2.4. Hysteresis
In this work, we use the three-phase hysteresis model suggested by Blunt (2000). In this model,
assuming a water-wet porous medium, gas is trapped by both water and oil. Based on Land’s
method (Land, 1968), the amount of gas that can be trapped is related to the maximum saturation
of gas reached during the displacement process prior to the onset of imbibition:

𝑆 𝑔𝑟
=
𝑆 𝑔 𝑚𝑎𝑥
1+𝐶 𝑔 𝑆 𝑔 𝑚𝑎𝑥 , (3.7)

where Sgr is the amount of gas that can be trapped, Sg
max
is the historical maximum gas saturation,
and Cg is the Land trapping constant. Cg can be found from the residual gas saturation to water in
a two-phase gas-water displacement experiment. The flowing saturation of gas is then found from

24

𝑆 𝑔𝑓
=
1
2
[(𝑆 𝑔 −𝑆 𝑔𝑟
)+√(𝑆 𝑔 −𝑆 𝑔𝑟
)
2
+
4
𝐶 𝑔 (𝑆 𝑔 −𝑆 𝑔𝑟
)] . (3.8)

Hence, imbibition relative permeability of a given gas saturation is equal to the drainage relative
permeability of its corresponding flowing gas saturation.
Oil can be trapped by water in a water-wet system. To account for this, Blunt (2000) proposed to
consider the combined trapping of hydrocarbon (oil and gas) by water. The amount of oil that is
trapped is then given by the difference between the trapped hydrocarbon and the trapped gas
saturations. The residual hydrocarbon saturation is found from

𝑆 ℎ𝑟 =
𝑆 ℎ
𝑚𝑎𝑥
1+𝐶 ℎ
𝑆 ℎ
𝑚𝑎𝑥
, (3.9)

where Sh
max
is the maximum hydrocarbon saturation during the displacement, max(So+Sg), and Ch
is the smaller value of Cg and Co, i.e. min(Co, Cg). Co is found in a similar way to Cg, from two-
phase oil-water displacement. The flowing hydrocarbon saturation can then be obtained from

𝑆 ℎ𝑓 =
1
2
[(𝑆 ℎ
−𝑆 ℎ𝑟 )+√
(𝑆 ℎ
−𝑆 ℎ𝑟 )
2
+
4
𝐶 ℎ
(𝑆 ℎ
−𝑆 ℎ𝑟 )] (3.10)
The flowing oil saturation is finally calculated from Sof = max (Shf – Sgf, 0). The flowing oil and
gas saturations are then used in Eqs. 3.1 through 3.6 to calculate relative permeabilities in a three-
25

phase system.
3.2.5. IFT Scaling
In this work we consider injection of a gas composition that result in a near-miscible displacement
of the oil. As a result, we have to take into account the reduction in the interfacial tension between
oil and gas caused by the mixing of oil and gas during the displacement. A reduction in the
interfacial tension causes a reduction in the residual saturation of oil to gas flooding (Sorg), which
in turn causes the oil-gas relative permeability curves to change. Following Coats (1980), this
change in relative permeabilities can be accounted for by a scaling factor (F), which indicates how
close a mixture is to a critical point (Shojaei et al., 2012). The scaling factor (F) is calculated from

𝐹 = (
𝜎 𝜎 0
)
𝛽 , (3.11)

where σ0 is the reference interfacial tension, corresponding to the immiscible relative permeability
curves, and β is an adjustable scaling exponent, typically in the range of 0.1 to 0.25 (Coats 1980).
Here, we use a value of β equals to 0.15. Relative permeabilities of oil and gas are then calculated
from

𝑘 𝑟𝑔
=𝐹 ∗𝑘 𝑟𝑔
𝑛𝑜 −𝑖𝑓𝑡 +(1−𝐹 )∗𝑘 𝑟𝑔
𝑟𝑒𝑓 −𝑖𝑓𝑡 , (3.12)

𝑘 𝑟𝑜𝑔 =𝐹 ∗𝑘 𝑟𝑜𝑔 𝑛𝑜 −𝑖𝑓𝑡 +(1−𝐹 )∗𝑘 𝑟𝑜𝑔 𝑟𝑒𝑓 −𝑖𝑓 𝑡 , (3.13)

26

where krog
ref-ift
and krg
ref-ift
are two-phase oil and gas relative permeabilities calculated at a reference
IFT, respectively. krog
no-ift
and krg
no-ift
are oil and gas relative permeabilities when IFT tends to zero.
It is assumed that at the limit of zero interfacial tension, relative permeability curves will be linear
functions of saturation with unit slope and that there is no residual saturation of oil to gas flooding.
Fig. 3.2 shows oil-gas relative permeability curves when the model is at reference IFT and when
there is no IFT between the phases. Note that the gas relative permeability in the case of zero IFT
does not reach 1, because the gas-oil relative permeability curves are plotted at connate water
saturation.

Figure 3.2 Relative permeability and IFT scaling
3.3. 2-D Model Setup
In order to study the effect of counter-current flow on displacement performance, we start with a
two-dimensional homogeneous reservoir model that helps us identify the segregation distance and
related oil recovery. The 2-D model has dimensions of 120*3*20 m and is, for simulation purposes,
divided into 40*1*20 blocks, each of size 3*3*1m. The reservoir has a constant porosity of 25%
27

and horizontal permeability of 200 mD. The vertical permeability ranges from 50 mD to 200 mD
in different cases. Injectors are set on the left side of the model. Water is injected into the top 10
blocks at the rate of 4.5 Rm
3
/day, while gas is injected in the bottom 10 blocks with the rate of 4.5
Rm
3
/day. A production well is located on the right boundary of the reservoir and the producer is
operating at a constant pressure of 10 MPa (1450 psia). The initial pressure of the reservoir is 12
MPa (1740 psia).
The fluid description used to represent the oil/gas system is reported in Table 3.1 and we use the
Soave-Redlich-Kwong (Soave, 1972) EOS in our equilibrium calculations. Solubility of non-
hydrocarbon gases and hydrocarbons in the water phase is not considered in this work. At the
reservoir temperature of 336.15K, the injected gas is multicontact miscible with the oil at 14 MPa
and hence, the displacement examples will be near-miscible.
Fig. 3.3 represents the physical properties of the two-dimensional domain. We use an IMPEC
finite-difference compositional approach to simulate the flow dynamics.
Table 3.1: Fluid description and compositions of oil and injection gas for SWAG injections
Critical
Temperature(K)
Critical
Pressure(atm)
Component
Omega(ω)
Molecular
weight
Mole
fraction in
injection
gas
Mole
fraction in
reservoir
oil
CO
2
304.2 72.86 0.2236 44.01 0.8 0.0
Methane 190.6 45.39 0.0115 16.04 0.2 0.3
n-Butane 425.1 37.46 0.2002 58.12 0.0 0.4
n-Dodecane 658.0 17.96 0.5764 270.33 0.0 0.3

28

Figure 3.3 2-D reservoir model for SWAG injections
3.4. Results of the 2-D Model
In this section, we consider the simultaneous injection of water and gas and investigate the impact
of flow reversals on displacement performance by comparing two mobility models:
1. A pure co-current flow model using a single set of relative permeability curves (including their
hysteresis).
2. A mixed-flow model allowing for transition between co-current and counter-current relative
permeability.
In the co-current flow model, co-current relative permeability curves are used for both horizontal
and vertical directions. This model assumes that transitions between co-current and counter-
current flows do not affect phase mobilities. In the mixed flow model, both co-current and counter-
current relative permeability are considered. Co-current relative permeability curves are used
whenever flow is co-current and counter-current relative permeability curves are used whenever
flow is counter-current in the vertical direction. In the following section, we report and compare
our findings for two displacement calculation examples.

29

3.4.1. Example 1: Post Waterflood EOR
We start by considering a homogeneous reservoir after an initial water flood. We assume the
reservoir has been previously water-flooded and currently exists at residual oil saturation to water
(Sorw = 0.3). Gas injection is then performed to recover bypassed/trapped oil in the reservoir while
the simultaneous water injection aims to control the mobility of gas and to increase the sweep
efficiency of the EOR process. Three different vertical to horizontal permeability ratios were
considered: 0.25, 0.5, 1.0. Fig. 3.4 shows the gas saturation distribution for both co-current flow
and mixed flow models after 600 days of injection (1.5 pore volumes injected) while Fig. 3.5
reports the corresponding oil recovery.
30

Figure 3.4 Gas saturation distribution after 600 days of injection for different vertical permeabilities in post water-
flood tests for 2-D reservoir model, SWAG injection process

Figure 3.5 Oil recovery vs. time for different vertical permeabilities in post water-flood tests for 2-D reservoir
model, SWAG injection processes
31

For this setting, we observe that the integration of counter-current relative permeability in our
displacement calculations results in higher recoveries from the EOR process. This is particularly
evident as the vertical relative permeability increases. In contrast, for a vertical to horizontal
permeability ratio of 0.25, we observe only a modest difference between the co-current flow model
and mixed flow model in terms of recovery and saturation distribution. However, when the vertical
to horizontal permeability ratio increases, the difference between the results of the two flow models
also increases. It is also observed, as expected, that with the increase in vertical permeability, the
amount of trapped oil remaining in the reservoir increases: In reservoirs with high vertical
permeability, gravity segregation between gas and water occurs closer to the injection wells which
lead to a reduced sweep efficiency.
As seen in Fig. 3.4, full segregation does not occur over the length of the reservoir for a vertical to
horizontal permeability ratio of 0.25 (i.e. vertical permeability of 50 mD). As a result, we obtain a
better sweep in lower vertical to horizontal permeability ratios than for higher ratios.
Counter-current flow has significantly affected the segregation length in all cases. In the case of
100 mD vertical permeability, full segregation is opbserved at approximately 66 meters from the
injection well. Including counter-current permeability effects shifts this point to 90 meters from
the injection well. Also, in the case of 200 mD vertical permeability, the point of complete
segregation has shifted from 36 meters to 51 meters from the injectors. The shift of this point closer
to the production interval means that injected gas contacts a bigger portion of the reservoir before
fully segregating to the top of the reservoir, which leads to higher oil recovery. Incremental oil
recovery is reported as the incremental fraction of oil produced for each of the examples in Fig.
3.5. Including counter-current flow in the calculations is observed to increase the oil recovery and
that this increase is higher for higher vertical permeabilites. By comparing the curves, we can see
32

that most of the difference between recovery predicted by the co-current and mixed flow models
happens before breakthrough. After breakthrough, both curves follow a similar path until they
reach a plateau.
3.4.2. Example 2: Reservoir Initially at Connate Water Saturation
Next, we consider a case where simultaneous gas/water injection is implemented in a newly
discovered reservoir starting at the connate water saturation (Swc = 0.2). The predicted
displacement dynamics for this reservoir setting are similar to the previous case and a comparison
of gas saturation distributions is given in Fig. 3.6. Similar to the post water-flood case, the portion
of the reservoir that is not contacted by gas increases with increasing vertical permeability. When
the mixed flow model is used in the calculation, we observe a slight change in the saturation
distribution. In this example, the difference is not as evident as in the previous calculation example.
To better illustrate this difference, we report the oil recovery versus time for each flow model in
Fig. 3.7. We observe that the impact of including counter-current relative permeability is
significantly less than in the previous case. However, the impact of including counter-current
relative permeability is amplified as the vertical permeability increases.
33

Figure 3.6 Gas saturation distribution after 600 days of injection for different vertical permeabilities in reservoir
with connate water saturation for 2-D reservoir model, SWAG injection process

Figure 3.7 Oil recovery vs. time for different vertical permeabilities in res. with connate water saturation for 2-D
reservoir model, SWAG injection processes
34

3.5. 3-D Reservoir Model
Next, we use a three-dimensional reservoir model with a heterogeneous permeability distribution,
with a range of 7.26-9760 mD. The model has the dimensions of 150*150*20 m and is divided
into 50*50*10 grid blocks for simulation. The porosity of the model is constant at 30%. The
injection interval is on the left corner of the reservoir model and the production interval is located
at the right corner of the model, representing a quarter of a 5-spot injection pattern. Rossen et al.
(2010) report that the optimal scheme for WAG injection is to inject water at the top section of the
reservoir, and to inject gas at the bottom section. As a result, in our model we assume that, with
modern well technology, we have the ability to alternate between the perforations for gas and water
injection. The injector is constrained by a constant rate of 250 rm3/day, while the producer is kept
at constant pressure of 110 bar (1595 psia). The initial pressure at the top of the reservoir is 120
bar (1740 psia). Fig. 3.8 represents the physical properties of the three-dimensional domain and
the distribution of logarithm of horizontal absolute permeability in the model.

Figure 3.8 3-D reservoir model (left) and log of permeability distribution (right) for SWAG injections
35

3.6. Results of the 3-D Model
In this section we report the results of simulation for the 3-D model. As in the previous section,
two initial reservoir conditions are studied. In the first example, we study a reservoir that has
previously been water flooded and is now at residual oil saturation to water flooding (Sorw = 0.3).
The second example is a “fresh” reservoir at its connate water saturation (Swc = 0.2), with no
history of production operations. We also investigate 2 variations of vertical permeability. In the
first case, the vertical permeability is 10% of the horizontal permeability for each cell, and in the
second case, vertical and horizontal permeabilities are equal.
3.6.1. Example 3: Post Waterflood EOR
We start by considering a reservoir after a water-flood with Sorw = 0.3. Fig. 3.9 and Fig. 3.10 show
the saturation distribution for this case for kv = 0.1*kh and kv = kh, respectively.

Figure 3.9 Gas, oil and water saturation distribution after 800 days (1.5 PVI), post waterflood, 3-D reservoir model,
kv = 0.1 kh, SWAG injection process.
36

In these figures, the top row shows saturation distribution after 800 days of injection when we do
not include counter-current relative permeability (i.e. “conventional simulation”) and the bottom
row shows the saturation distribution for the simulations that account for counter-current relative
permeability in the vertical direction.
As can be seen, there is a noticeable difference between the saturation distributions resulting from
the two displacement calculations. We see that when accounting for counter-current relative
permeability, oil is drained from a larger fraction of the reservoir due to a reduced relative mobility
in the vertical direction for water and gas that force gas to contact a larger fraction of the reservoir
before full segregation occurs.
Another observation is that in the reservoir model with lower vertical permeability, segregation of
gas and water does not happen quickly. Because of the low vertical permeability, the effect of
gravity segregation has been reduced and oil is displaced from a larger portion of the reservoir. In
contrast, with a higher vertical permeability, we see that segregation happens close to the injection
37

point and oil is swept only from a small portion of the reservoir.

Figure 3.10 Gas, oil and water saturation distribution after 800 days (1.5 PVI), post waterflood, 3-D reservoir
model, kv = kh, SWAG injection process.
3.6.2. Example 4: Reservoir Initially at Connate Water Saturation
Next, we consider a reservoir at its connate water saturation with no previous history of floods.
Based on our relative permeability model, the connate water saturation for this system is 0.2. As a
result, the reservoir is initiated with 20% water and 80% oil in all cells. Figures 3.11 and 3.12 show
the saturation distribution after 800 days (1.5 PVI) of simultaneous injection of water and gas for
two different vertical permeability settings.
38

Figure 3.11 Gas, oil and water saturation distribution after 800 days (1.5 PVI), connate water saturation, 3-D
reservoir model, kv = 0.1 kh, SWAG injection process.
The results from these displacement calculations demonstrate, as before, that accounting for
counter-current flow in our relative mobility calculations affects the saturation distribution in the
system. By comparing these results with the 2-D calculation examples, we see that the impact of
flow reversal and associated reduction on fluid mobilities is less significant for a reservoir initially
at connate water saturation.
39

Figure 3.12 Gas, oil and water saturation distribution after 800 days (1.5 PVI), connate water saturation, 3-D
reservoir model, kv = kh, SWAG injection process.
Next, we report and compare the oil recovery and producing GOR as predicted by the above
displacement calculations. Fig. 3.13 shows these results for oil recovery vs. pore volumes injected
while Fig. 3.14 shows the producing GOR as a function of pore volumes injected.
The top row in both figures shows results for the post water flood example (Example 3) and the
bottom row shows results for the reservoir at connate water saturation (Example 4). The left
column shows results for kv = 0.1 kh while the right column displays results for kv = kh. In each
of the figures, we report a comparison between the recovery (and GOR) of the conventional model
(co-current flow) and the new model with counter-current relative permeability in the vertical
direction (mixed flow). We call this new model the mixed flow model, because we use co-current
relative permeability in the horizontal direction and transitions between co-current and counter-
current relative permeability in the vertical direction.
40

Figure 3.13 Recovery vs. PVI for all cases, 3-D reservoir model, SWAG injection process
Recovery factors are higher in all cases for the mixed flow model than for the co-current flow
model. This is contributed to the fact that a lower vertical mobility of phases during counter-current
flow, causes segregation of gas and water to occur closer to the production well and, hence, a larger
portion of the reservoir is swept before the injection fluids fully segregate.
Differences between the two flow models are more evident with higher absolute vertical
permeability. Also, the effect of accounting for counter-current relative permeability is larger for
post-water flood SWAG injection.

41

Figure 3.14 GOR (Sm3/Sm3) vs. PVI for all cases, 3-D reservoir model, SWAG injection process
Finally, in Fig. 3.14 we see that the producing gas oil ratio (Sm3/Sm3) is always lower for the
mixed-flow model than for the co-current flow model.
3.7. A Different Three Phase Relative Permeability Model
In order to gauge the impact of flow transitions and its coupling with the relative permeability
model, we apply another widely used three-phase relative permeability model, which is the
saturation weighted interpolation method. Baker (1988) used saturation-weighted interpolation
between the two-phase oil-water and gas-oil relative permeabilities for the oil phase to find the
42

three-phase oil relative permeability:
𝑘 𝑟𝑜
=
(𝑆 𝑤 −𝑆 𝑤𝑖
)𝑘 𝑟𝑜 (𝑤 )
+(𝑆 𝑔 −𝑆 𝑔𝑟
)𝑘 𝑟𝑜 (𝑔 )
(𝑆 𝑤 −𝑆 𝑤𝑖
)+(𝑆 𝑔 −𝑆 𝑔𝑟
)
. (3.14)
In the above equation, Sgr is the residual gas saturation in the two-phase oil/water experiment
(usually zero because oil-water experiments are performed in absence of gas). kro(w) is the oil
relative permeability in an oil-water two-phase experiment and kro(g) is the oil relative permeability
in a gas-oil two-phase experiment.
Baker (1988) compared three-phase relative permeability calculations using 12 different empirical
correlations and concluded that saturation-weighted interpolation provides the best agreement with
experimental data. Later, Delshad and Pope (1989) compared seven empirical three-phase relative
permeability models against three experiments. They also showed that saturation-weighted
interpolation was superior to Stone’s models (Blunt, 2000). Note that this treatment of the three-
phase relative permeability is the default calculation in Eclipse simulator.
The above formulation for 3-phase relative permeability calculation was implemented in our
IMPEC simulator. Except for oil 3-phase relative permeability calculations, the rest of the
formulations are the same as presented in the previous sections. We use the same 3-D reservoir
model to study this relative permeability model and we look at the two different initialization
examples explored above.
3.7.1. Example 5: Post Waterflood EOR, Interpolation
In the same manner as the previous sections, we start by considering a reservoir after a waterflood
with oil at its residual saturation to water. Figures 3.15 and 3.16 show the saturation distribution
after 800 days (1.5 PVI) of injection for two different vertical to horizontal absolute permeability
43

ratios of 0.1 and 1, respectively.

Figure 3.15 Gas, oil and water saturation distribution after 800 days (1.5 PVI), post water flood, 3-D reservoir
model, kv = 0.1 kh, interpolation model for rel perm, SWAG injection process

With the new model for calculating three-phase oil relative permeability, we observe the same
differences in the saturation distributions as we did using the modified Stone 1 model. In both
cases, segregation of gas and water happens closer to the production interval for the mixed flow
model, and with the increase in the vertical absolute permeability, segregation happens earlier (as
44

expected).

Figure 3.16 Gas, oil and water saturation distribution after 800 days (1.5 PVI), post water flood, 3-D reservoir
model, kv = kh, interpolation model for rel perm, SWAG injection process
3.7.2. Example 6: Reservoir at Connate Water Saturation, Interpolation
Next, we consider a case where the reservoir has not been flooded previously, and is at connate
water saturation (Swi = 0.2). Saturation distributions are reported in Figures 3.17 and 3.18. As
seen in the figures, a larger portion of the reservoir is contacted by the injected fluids when we
account for counter-current flow effects in the phase mobilities in the vertical direction. This
suggests that counter-current flow plays an important role in the dynamics of a system where we
45

are simultaneously injecting gas and water.

Figure 3.17 Gas, oil and water saturation distribution after 800 days, connate water saturation, 3-D reservoir model,
kv = 0.1 kh, interpolation model for rel perm, SWAG injection process
To better observe this effect, we report the recovery and gas oil ratio as a function of pore volumes
injected in Figures 3.19 and 3.20. The top row in both figures shows results for the post water
flood example (Example 5) and the bottom row shows results for the reservoir at connate water
saturation (Example 6). The left column shows results for kv = 0.1 kh and the right column shows
results for kv = kh. In each of the plots, we see a comparison between the recovery (and GOR) of
46

the conventional model (co-current flow) and the new model which represents counter-current
flow effects in the vertical direction (mixed flow).

Figure 3.18 Gas, oil and water saturation distribution after 800 days (1.5 PVI), connate water saturation, 3-D
reservoir model, kv = kh, interpolation model for rel perm, SWAG injection process
When counter-current relative permeability is accounted for in the vertical direction, phase
mobilities in the vertical direction are lower than they would be in the case of co-current flow.
This in turn translates into a delay of the gravity segregation between water and gas. As a result, a
larger percentage of the reservoir is swept by the injection fluids, and higher recovery factors are
observed.
In all of the figures, the dashed curve shows the results for the mixed flow model and the solid
curve shows the results for pure co-current flow model. The observed trends of the recovery factor
47

and the gas-oil ratio as a function of pore volumes injected is similar for both flow models. The
two different relative permeability models provide comparable results (Figs 3.19 and 3.20)

Figure 3.19 Recovery vs. PVI for all cases, 3-D reservoir model, interpolation model for rel perm, SWAG injection
process
Comparing the recovery results for the two different relative permeability models (modified Stone
1 and interpolation model), the two models reach at comparable results in all cases (See Figures
3.13 and 3.19). For the case of Kv = 0.1 Kh in the post water flood example, a large difference is
observed in terms of incremental oil recovery.
48

Figure 3.20 GOR vs. PVI for all cases, 3-D reservoir model, interpolation model for rel perm, SWAG injection
process
3.8. Discussion
In the previous sections, we have demonstrated the importance of considering counter-current
relative permeability in the simulation and design of simultaneous water and gas injection
processes. We observe, that the largest differences in saturation between the co-current flow model
and the mixed flow model is found at the gas/water fronts. We also studied the pressure distribution
in the system and observed that the largest differences in the pressure distribution in the reservoir
between the two flow models are observed near the producing well, but that the pressure difference
is relatively small in all cases (e.g. in the order of 0.001 Mpa).
The impact of including phase mobility reductions due to counter-current flow is more evident in
49

terms of oil recovery versus time. A relative difference of 50 percent (e.g. Fig. 3.5) clearly
demonstrates the need of understanding/including counter-current relative permeability models in
the simulation of simultaneous water and gas injection processes. The importance of counter-
current relative permeability is higher when lower oil saturations are present in the reservoir. This
is due to the fact that the density contrast between the injected gas and the reservoir fluid is higher
when oil is at residual saturation (i.e. the mobile fluids are mainly gas-water, with higher density
difference than gas-oil). Higher density contrast results in faster gravity segregation in the
reservoir. As a result, effects of counter-current relative permeability are more dominant in the
reservoir and affects the saturation distributions and recoveries more significantly.
The current literature on counter-current relative permeability is largely limited to two-phase
oil/water systems. To the best of our knowledge, no such study has been performed on three-phase
systems with all phases mobile. From the results presented in this study, the use of counter-current
relative permeability appears to be a good choice for representing fluid mobilities in the design of
post-waterflood SWAG processes.
We have also investigated other settings with different boundary conditions, e.g. injection of water
over the top 15 meters and gas over the bottom 5 meters of the injection interval, injection of water
over the top 5 meters and gas over the bottom 5 meters, with 10 meters of space between the
perforations, and injecting at a higher water rate. The results from the alternate boundary
conditions were similar to the cases presented here.
In the horizontal direction, flow of all phases is predominately co-current. As a result, there is no
need to include counter-current flow transitions in the horizontal direction. In the vertical direction,
however, counter-current flow is the dominant flow mode because of gravity segregation of gas
50

that is injected from the bottom, and water that is injected at the top. Furthermore, the introduction
of the counter-current relative permeability model did not significantly increase the simulation
time. Since we are using an IMPES simulator, flow transitions at each time step is identified from
the previous time step. In fully implicit simulation, additional studies will be required to ensure
the numerical stability and efficiency. This study is the first attempt to incorporate flow transitions
and to study its effect on the design of SWAG process. Further investigation is warranted to
validate of the proposed flow model implementation for three-phase flows, both from an
experimental point of view and in the context of fully implicit large-scale simulation.
3.9. Conclusions
Based on the examples and analysis presented in the previous sections, we arrive at the following
conclusions:
1. Counter-current flow and the associated reduction in phase relative permeabilities can play
an important role in the dynamics of simultaneous water and gas injection processes, especially
when the vertical permeability is relatively high.
2. Full segregation of the injection fluids (i.e. gas and water) occurs closer to the producing
well when counter-current relative permeability modelling is included in the displacement
calculations. This suggests that transitions between co-current and counter-current relative
permeability should be integrated in the design of SWAG processes to provide a more accurate
estimation of the segregation length
3. Simulation of simultaneous water and gas injection processes using a single set of
permeability curves may not provide sufficient physical detail to ensure accurate design of SWAG
51

processes.

52

Chapter 4 Water Alternating Gas Injection Processes
2

4.1. Introduction
In this chapter, we investigate the impact of mobility changes due to flow reversals from co-current
to counter-current flow on the displacement performance of Water Alternating Gas (WAG)
injection processes. In WAG processes, the injected gas will migrate towards the top of the
formation while the injected water will migrate towards the bottom of the formation. The
segregation of gas, oil and water phases will result in counter-current flow occurring in the vertical
direction in some portions of the reservoir during the displacement process.
Previous experimental and theoretical studies of counter-current flow have shown that the relative
mobility of each of the phases in a porous medium is considerably less when counter-current flow
prevails as compared to co-current flow settings. A reduction of the relative permeability in the
vertical direction results in a dynamic anisotropy in phase mobilities. This effect has, to the best
of our knowledge, not previously been considered in the modeling and simulation of WAG
processes.
A new flow model that accounts for flow reversals in the vertical direction has been implemented
and tested in a three-phase compositional reservoir simulator. In order to investigate the impact of
flow reversals, results from the new flow model are compared to cases where counter-current flow

2
This chapter was previously published in Transport in Porous Media journal: Sherafati M., Jessen K., Dynamic
Relative Permeability and Simulation of WAG Injection Processes. Transp Porous Med (2017) 117: 125
53

effects on the phase mobilities are ignored. A range of displacement settings, covering relevant
slug sizes, have been investigated to gauge the impact of mobility reductions due to flow reversals.
Significant differences, in terms of saturation distribution, producing GOR and oil recovery, are
observed between the conventional flow model (ignoring mobility reductions due to counter-
current flow) and the proposed new model that accounts for reductions in phase mobility during
counter-current flow. Accordingly, we recommend that an explicit representation of flow
transitions between co-current and counter-current flow (and the related impact on phase
mobilities) should be considered to ensure accurate and optimal design of WAG injection
processes.
Several authors have focused on the simulation of both miscible WAG and immiscible WAG
injection. Larsen and Skauge (1999) developed a methodology to simulate immiscible WAG
processes that accounts for three-phase flow effects and hysteresis, while Li et al. (2003) developed
a compositional model for WAG processes in viscous oil reservoirs. Bhambri and Mohanty (2005)
and Yan et al. (2004) developed streamline simulators to model WAG processes and take into
account a potential 4th phase (CO2-rich liquid phase). More recently, Al-Kobaisi et al. (2013)
studied CO2 trapping in WAG processes based on three-phase, three-dimensional displacement
calculations.
In this work, we investigate the effects of including a mobility reduction due to flow reversals in
displacement calculations, with emphasis on the impact on common EOR process design
parameters. Important design parameters to take into consideration include oil recovery, GOR, and
Net Present Value (NPV). In this chapter we investigate WAG processes for a range of cases with
54

and without mobility reductions due to flow reversals and investigate how the relevant design and
operational parameters differ when the impact of flow reversals is included.
We start with a discussion of relative permeability in counter-current flow settings followed by a
presentation of the relative permeability description (including hysteresis) used in the subsequent
displacement calculations. We then introduce the reservoir model and compare the displacement
dynamics of WAG processes, with and without considering the impact of flow reversals, for a
range of relevant conditions. We conclude the manuscript with an analysis and discussion of our
observations.
4.2. Methodology
Calculations for three-phase relative permeabilities are performed as described in the previous
chapter. The modified Stone 1 model is used for interpolating between two-phase relative
permeabilities to arrive at a three-phase relative permeability curve for the intermediate wetting
phase (assumed here to be oil). IFT scaling is performed using Coats’ model (Coats, 1980), and
hysteresis in the three-phase system is taken into account using Blunt’s model (Blunt 2000).
4.3. Model Setup
In this study, we use a three-dimensional reservoir model with a heterogeneous permeability
distribution. The model has the dimensions of 150*150*20 m and is divided into 50*50*10 grid
blocks for displacement calculations. The porosity of the model is constant at 30% resulting in a
model pore volume of 1.35·10
5
m
3
. An injection well is located on the left corner of the reservoir
model, while a production well is located on the right corner of the model, representing a quarter
of a 5-spot injection pattern. Rossen et al. (2010) demonstrated that the optimal scheme for WAG
55

injection is to inject water at the top portion of the reservoir, and to inject gas at the bottom section.
We adopt that injection strategy, assuming that we can alternate between the perforations for gas
and water injection. The injector is constrained at a constant rate of 250 rm
3
/day, while the
producer is operated at a constant bottom-hole pressure of 110 bar (1595 psia). The initial pressure
at the top of the reservoir is 120 bar (1740 psia). At the given injection rate, the injection of one
pore volume (1 PV) corresponds to ~ 540 days of injection.
Fig. 4.1 shows the physical properties of the three-dimensional domain and the distribution of the
logarithm of the absolute horizontal permeability in the model. The vertical permeability in the
reservoir is 0.1 times the horizontal permeability. We use our in-house IMPES finite-difference.
The 3-D reservoir model used is the same 3-D model used in chapter 3 for the simulation of SWAG
injection processes.

Figure 4.1 3-D reservoir model (left) and log of permeability distribution (right), WAG injections

56

Table 4.1 summarizes the compositional fluid description for the oil/gas system, while the binary
interaction coefficients are listed in Table 4.2.

Table 4.1 Fluid description and compositions of oil and injection gas, WAG injections
Critical
Temperature(K)
Critical
Pressure(atm)
Component
Omega(ω)
Molecular
weight
Mole fraction
in injection gas
Mole fraction
in reservoir oil
CO
2
304.2 72.86 0.2236 44.01 0.8 0.0
Methane 190.6 45.39 0.0115 16.04 0.2 0.3
n-Butane 425.1 37.46 0.2002 58.12 0.0 0.4
n-Dodecane 658.0 17.96 0.5764 270.33 0.0 0.3

The SRK equation of state (Soave, 1972) was used throughout our calculations. At the reservoir
temperature of 336.15K, the injected gas is multi-contact miscible with the oil at 140 bar and
hence, the following displacement examples are at near-miscible conditions. Apart from the four
components in the oil/gas system, water is also present in the reservoir. It is assumed in our
compositional simulator that there is no solubility between water and the oil/gas system.
Table 4.2 Binary Interaction Coefficients
CO
2
Methane n-Butane n-Dodecane
CO
2
0.000 0.105 0.120 0.115
Methane 0.105 0.000 0.000 0.000
n-Butane 0.120 0.000 0.000 0.000
n-Dodecane 0.115 0.000 0.000 0.000

4.4. Results
In order to investigate the effects of flow dependent phase mobility on the performance of WAG
injection schemes, we compare two different flow models. The first flow model uses co-current
relative permeability input only, and accounts for relative permeability hysteresis and IFT scaling.
This approach is the model commonly used in commercial simulation tools. The second flow
57

model allows for transitions between the co-current and counter-current relative permeability
curves. Since there is no evidence of counter-current flow occurring in the horizontal direction (in
these displacements), co-current relative-permeability curves are used exclusively in the horizontal
direction. In the vertical direction, the flow model allows for transitions between co-current and
counter-current relative permeability curves. We refer to this flow model as the “mixed flow
model”. The second model also includes hysteresis effects and IFT scaling as discussed previously.
WAG injection is often performed once a reservoir has been water flooded to target/recover the
residual and bypassed oil in the reservoir. Accordingly, we consider a reservoir that exists after an
initial water flood, with the remaining oil at its residual saturation to water (Sorw = 0.3). Initially,
we consider a slug size of 0.1 PV. The WAG ratio throughout our work is constant at unity and
the displacement calculations start with the injection of a gas slug, followed by subsequent slugs
of water and gas.
Saturation distributions for co-current and mixed flow models after 800 days of injection (1.5 PVI)
are reported in Figs. 4.2, 4.3 and 4.4. In these figures, the saturation distribution for gas, oil and
water are shown for the top, middle and bottom layers of the reservoir model. We observe that the
saturation distribution is predicted to be quite different between the two flow models. In particular,
we observe that the sweep of the middle layer (Fig. 4.3) is predicted to be higher by the mixed
flow model relative to the static flow model. This is a direct result of a reduction in fluid mobility
in the vertical direction due to counter-current flow settings: The travel distance to full segregation
of gas and water is predicted be longer by the mixed flow model.
To quantify the observed differences in saturation distributions, we calculate the root-mean-square
difference (RMSD) between the calculation results (at the final time) obtained by the co-current
flow model and the mixed flow model. For the first calculation example, the RMSD for gas, oil
58

and water saturations are 0.1160, 0.096 and 0.08, respectively. A comparison of producing GOR
and recovery factors is presented and discussed in more detail below.

Figure 4.2 Gas, oil and water saturation distribution for 0.1 PV slug size and Swi=0.7 after 800 days (1.5 PVI) of
injection, Top layer
59

Figure 4.3 Gas, oil and water saturation distribution for 0.1 PV slug size and Swi=0.7 after 800 days (1.5 PVI) of
injection, Middle layer

Figure 4.4 Gas, oil and water saturation distribution for 0.1 PV slug size and Swi=0.7 after 800 days (1.5 PVI) of
injection, Bottom layer
60

In the following subsections, we investigate the effects of slug size, initial water saturation, and
capillary pressure on the displacement dynamics and WAG performance. In each subsection that
follows, we analyze the effect of one parameter, while all the other input parameters for the
displacement calculations are kept constant as detailed above. Table 4.3 reports the summary of
the calculation examples that are considered in this study.
Table 4.3 Summary of tests for WAG injection process with different slug sizes
Example number Model initialization Injection time (days) Slug size (PV)
1
Post waterflood
(S
w,init
= 0.70)
800 0.10
2 800 0.25
3 800 0.50
4 Connate water saturation
(S
w,init
= 0.2)
800 0.10
5 800 0.25
6 800 0.50

4.4.1. Effect of Slug Size
The slug sizes of the injected gas used in field trials are mainly in the range of 0.1 to 3 PV
(Christensen et al., 2001). The optimal slug size for a field is governed by project economics, but
the most popular slug sizes used in WAG projects in the US is 0.1-0.5 PV (Christensen et al.,
2001). We considered 3 different slug sizes for WAG injection: 0.5, 0.25 and 0.1 PV. We present
the saturation distribution, oil recovery and gas-oil ratio (GOR) vs. pore volume injected for each
of these slug sizes separately.
We start by considering a slug size of 0.5 PV. With an injection rate of 250 rm
3
/day, it takes ~ 270
days to inject 0.5 pore volumes of fluid. As a result, over the time of the calculation example (800
days) we inject ~ 3 slugs. Significant differences are again observed in the distribution of fluid
saturations between the predictions of the two flow models with RMSD’s of 0.1675, 0.0946, and
61

0.1250 for gas, oil and water, respectively. We observe that all RMSD’s are larger for the 0.5 PV
slug example compared to the 0.1 PV slug example, and we attribute this increase in differences
between calculation results from the two flow models to the use of larger slugs that results in more
significant phase segregation and less trapping of gas.
Next, we consider a case where the slug size of the injection fluid is 0.25 of the pore volume of
the reservoir. It takes approximately 140 days to inject this slug size. The displacement dynamics
predicted for this WAG scenario are similar to the 0.5 PV slug example: The portion of the
reservoir that is contacted by the injected gas increases when counter-current relative permeability
is included in the displacement calculations. The RMSD’s for this example are 0.0461, 0.0812 and
0.0525 for gas, oil and water respectively. Here we observe that the RMSD’s for this slug size are
much lower than for the examples using 0.1 PV and 0.5 PV slug sizes. The calculated RMSD’s
are reported only at the final time of simulation, 800 days. At this instance, the reservoir is on a
gas injection cycle for the 0.1 PV and 0.5 PV slug sizes while water is injected for the 0.25 PV
slug example. During the water injection cycle, the density contrast between the injection fluid
and reservoir fluid is considerably less than during the gas injection cycle. As a result, we observe
less difference between the two flow models in terms of final saturations (lower RMSD’s for the
0.25PV slug example that is on water cycle at the final stage). Accordingly, we note that the RMSD
values would be more comparable if all three examples were going through a gas injection cycle
at the final stage.
4.4.2. Effect of Initial Water Saturation
The following examples consider the case where WAG injection is implemented in a newly
discovered reservoir starting at the connate water saturation (Swc = 0.2). The differences between
62

the predicted displacement dynamics, as obtained from the two flow models in this reservoir
setting, are similar to the previous examples and a comparison of gas, oil and water saturation
distributions is provided in Fig. 4.5, Fig. 4.6 and Fig. 4.7 for a 0.1 PV slug size WAG scheme.
When the counter-current relative permeability model is included in the displacement calculations,
we observe changes in the saturation distribution. The Root Mean Square of the differences
between the two flow models in terms of gas, oil, and water saturations are reported in Table 4.4
below.

Figure 4.5 Gas, oil and water saturation distribution for 0.1 PV slug size and Swi=0.2 after 800 days of injection,
Top layer
63

Figure 4.6 Gas, oil and water saturation distribution for 0.1 PV slug size and Swi=0.2 after 800 days of injection,
Middle layer

Figure 4.7 Gas oil and water saturation distribution for 0.1 PV slug size and Swi=0.2 after 800 days of injection,
Bottom layer
64

To better illustrate the difference between the two flow models, we report the oil recovery and
GOR versus pore volume injected (PVI) for each of the above example calculations in Fig. 4.8 and
Fig 4.9. The subplots in Fig. 4.8 report oil recovery vs. pore volume injected (PVI) for a different
slug size and initial water saturation.

Figure 4.8 Recovery vs. PVI for all calculation examples
The oil recovery is reported as the fraction of the oil in place for each of the examples. We observe
that accounting for counter-current relative permeability in the calculations results in significant
incremental oil recovery, as high as 10% (corresponding to a relative increase of 58%) for the 0.25
PV slug injection into a previously water-flooded reservoir. By comparing the recovery curves,
we observe that most of the differences between the results from the co-current and counter-current
flow models appear before breakthrough. For all displacements, the use of the mixed flow model
65

results in a later breakthrough of the injected phases. Once breakthrough has occurred, both models
predict a relatively similar rate of recovery.
Similar to Fig. 4.8, Fig. 4.9 reports the producing GOR from the calculation examples.

Figure 4.9 GOR vs. PVI for all cases, 3-D reservoir model, WAG injection process
From Fig. 4.9, we observe that the use of the counter-current relative permeability model has a
significant impact on the producing GOR. In all cases, the mixed flow model reduces the GOR vs.
time. Lower GORs can be explained by the differences in the gas saturation distribution between
the two flow models. The amount of gas that is injected in all models is the same. In the mixed
flow model, gas travels a longer distance and contacts a larger portion of the pore volume before
reaching the production interval. As a result, gas reaches the production well at a later time for the
mixed flow model in comparison to the co-current flow model. The differences between flow
66

models are more pronounced for the post-water-flood examples. This can be explained, in part, by
the density difference between the gas-water pair being larger than the density difference for the
gas-oil pair.
4.4.3. Impact of Counter-Current Relative Permeability Curves
As discussed above, co-current and counter-current relative permeability curves were generated
using Corey-type equations. The relative permeability curves used in our simulations were based
on previous experimental work on counter-current relative permeability (Bourbiaux and
Kalaydjian 1990, Al-Wadahi et al. 2000) that demonstrated that counter-current relative
permeabilities are lower than the co-current values. The differences between the two relative
permeability curves can be as high as 50% in terms of endpoint relative permeability.
In this section, we report calculation results using our base case model (co-current relative
permeability only) combined with a different set of relative permeability curves to represent
counter-current flow settings. A new set of counter-current relative permeability curves, that lies
between the co-current and counter-current relative permeability curves of the base model were
used. In this setting, the reduction in relative mobility (~ 25%) due to counter-current flow is less
than in the previous calculation examples (50% reduction). We refer to this model as Mixed Flow
B. Simulation of WAG injection was then repeated with an initial water saturation of 0.7 (post
water-flood initialization). The slug size was set to 0.1 PV and a total of 1.5 pore volumes were
injected (800 days of injection). The results, in terms of recovery and GOR, are reported in Figure
4.10.
67

Figure 4.10 Results for two different counter-current relative permeability models, left: recovery vs. PVI, right:
GOR vs. PVI
We observe that the calculated recovery (left panel) from mixed flow model B lies in between the
original mixed flow model (mixed flow A) and the co-current flow model. This is expected, since
the modified counter-current relative permeabilities are more similar to the co-current values. The
ultimate recoveries calculated from the co-current flow model, the mixed flow model A, and the
mixed flow model B, are 0.16, 0.26 and 0.20, respectively. While the impact of flow direction on
fluid mobility is reduced, an incremental recovery of 4% is till observed (corresponding to 25%).
Similar observations and trends are observed for the gas-oil ratio versus PVI in Fig. 4.10 (right
panel).
4.4.4. Use of Static Anisotropic Relative Permeability to Represent Counter-Current Flow
Effects
68

Commercial and research simulators do not commonly allow for transitions between co-current
and counter-current relative permeability values. One possible way to account for reductions in
relative permeability due to counter-current flow in commercial simulators is the use of anisotropic
relative permeability: co-current relative permeability curves can be applied in the horizontal
direction while counter-current curves can be applied in the vertical direction, regardless of the
flow mode in the vertical direction. Figure 4.11 compares this approach with the results obtained
from dynamic modeling of transitions between flow modes in the vertical direction.

Figure 4.11 Comparison of anisotropic and dynamic relative permeability modeling
In this figure, recovery and GOR versus PVI is reported for three different flow models: 1) only
co-current relative permeabilities are used (static flow model), 2) including transitions between
co-current and counter-current relative permeabilities in the vertical direction (mixed flow A) and
3) use of counter-current relative permeabilities in the vertical direction, with no dynamic
transition (mixed flow model C).
69

We observe that mixed flow model C predicts a very high recovery relative to the other two
models. Recovery factors for the static flow model, the mixed flow model A and the mixed flow
model C are 0.16, 0.26 and 0.47, respectively. In terms of gas-oil ratio, the static model predicts a
higher GOR than the dynamic flow model, and the differences between these two flow models
magnify at later times during the displacement process. The significant differences are primarily
due to the permanent reduction of fluid mobility in the vertical direction used in the static
representation of counter-current flow effects. A similar effect would be observed by reducing the
absolute permeability on the vertical direction combined with the use of co-current relative
permeability only. Accordingly, the use of such representation will not provide an accurate
approximation of the counter-current flow effects.
4.4.5. Impact of Capillary Pressure
A direct representation of capillary pressure was not included in the above calculation examples
assuming that viscous and gravity forces dominate the displacement dynamics. In order to evaluate
the effect of capillary pressure, we repeat our simulation for the co-current flow model using
Eclipse 300, and include capillarity via a primary drainage capillary pressure function. Hysteresis
in capillary pressure is ignored to represent an extreme case (pc for primary drainage is higher than
for secondary imbibition/drainage cycles). We apply Leverett J-scaling, which is based on a
reference capillary pressure curve:
𝐽 (𝑆 )=
𝑃 𝑐 (𝑆 )√𝑘 𝜑 ⁄
𝜎 cos𝜃 .
(4.1)
We use values for a Berea Sandstone, as extracted by Brooks and Corey (1964), in our calculation
example with Eclipse (see Fig. 4.12).
70

Figure 4.12 J-function vs. saturation for primary drainage
WAG displacement calculations were then performed twice: once including the J-function for the
calculation of capillary pressure and once without capillary pressure. Both simulations were
performed with the co-current relative permeability curves, including their hysteresis. Fig. 4.13
reports the results from the Eclipse model, with and without the inclusion of capillary pressure.

Figure 4.13 Effect of capillary pressure on WAG simulation
71

We observe, from Figure 4.13, that the effect of capillary pressure on the recovery factor is
minimal compared to the differences we observe when including the impact of flow reversals on
phase mobilities. The simulation without capillary pressure results in a slightly lower production.
By comparing Fig. 4.13 with Fig. 4.8 (recovery vs. PVI from our simulator), we observe that the
calculations performed with Eclipse 300 result in an earlier breakthrough. The differences between
the two simulators is attributed to the fact that Eclipse 300 uses Killough’s model (Killough, 1976)
for hysteresis in relative permeability, while Blunt’s model (Blunt, 2000) is used in our simulator.
Blunt’s model is not available in Eclipse 300 and hence the trapping of oil and gas will invariably
differ between the two approaches.
4.5. Summary and Discussion
In the previous sections, we have studied the impact of mobility reductions as a result of flow
transitions between co-current and counter-current flow. In particular, we have demonstrated the
potential impact of flow reversals on the WAG displacement dynamics.
𝑁 𝑔𝑣
=
∆ρ.g.L.𝑘 𝑎𝑣
𝐻 .𝑞 .𝜇 𝑜 . (4.2)
From Eq. 4.2, assuming an average value for the Darcy velocity throughout the reservoir, we can
calculate two transverse gravity numbers, one for the oil-water pair and one for gas-oil pair.
Combined with a calculation of the endpoint mobility ratio, M, we can identify the flow regime
for the displacement calculations.
The value for the dimensionless number Ngv.M/(1+M) is evaluated to 8.32 for the oil-water pair
and 59 for the gas-oil pair. According to Zhou et al. (1993), the oil-water pair is in the transition
72

zone while the gas-oil pair is in a gravity dominated regime. Since we are injecting both water and
gas, the actual gravity number for the displacement lies somewhere between the values for oil-
water and gas-oil pairs, which indicates a gravity dominated displacement throughout the
reservoir.
The results from our displacement calculations indicate that when a higher initial water saturation
exist in the reservoir, corresponding to a post water-flood setting, larger differences are observed
between the two flow models, in comparison to a reservoir that is initially at a connate water
saturation. This can, in part, be attributed to the density contrast between the injected fluid and the
reservoir fluid. When larger volumes of water exist in the reservoir, the average liquid phase has
a higher density than when the reservoir is saturated by oil. As a result, a previously water flooded
reservoir has a higher density contrast with the injection gas. This promotes gravity segregation in
the reservoir. Accordingly, the effect of accounting for mobility reductions during counter-current
flow is more evident for this setting and results in larger differences in terms of saturation
distribution, oil recovery and producing GOR.
We also observe that the differences become more evident for smaller slug sizes. The differences
in terms of calculated oil recovery for a 0.1 PV slug size are significantly larger than for 0.5 PV
slugs. We find a relative difference of 58% in oil recovery for 0.1 PV slug injection into a
previously water-flooded reservoir after 800 days of injection (from ~17% to ~27% absolute
recovery).
Selection of the optimal slug size for a WAG injection process depends on the economics of the
project. Usually, relatively small slug sizes are preferred over larger slug sizes, based on
73

observation from previous field trials (e.g. Song et al., 2014) where smaller slug sizes result in a
faster return of the investment (i.e. earlier incremental production of oil).
Accounting for mobility change due to viscous coupling (flow reversals) results, in some settings,
in a 58% relative increase in the calculated incremental oil recovery. Since the aim of EOR
simulation studies is to assess potential incremental recovery and to design the EOR process (e.g.
selection of slug sizes), the observed change in the calculated recovery between conventional flow
models and the proposed approach can have a significant impact on project economics and design:
The proposed flow model will provide additional input, in terms of expected range of recoveries,
and allow for more successful design of WAG EOR projects.
As observed from table 4.4 below, the RMSD’s for oil, gas and water saturations are larger for
larger slug sizes for both model initializations (post water-flood and connate water saturation).
This is not in conflict with the calculated recovery factors: The difference in terms of oil recovery
is larger for smaller slug sizes. The use of the counter-current flow model results in a different
arrangement of the phases for larger slugs but does not considerable affect the calculated recovery.
Table 4.4 RMSD of gas, oil and water saturations for different slug sizes and different initializations
Example number Slug size (Pore
Volume)
Model Initialization RMSD gas RMSD
oil
RMSD
water
1 0.1 Post Water-flood 0.1160 0.0906 0.0800
2 0.25 Post Water-flood 0.0461 0.0812 0.0525
3 0.5 Post Water-flood 0.1675 0.0946 0.1250
4 0.1 Connate Water
Saturation
0.1126 0.0784 0.0736
5 0.25 Connate Water
Saturation
0.0356 0.0668 0.0459
6 0.5 Connate Water
Saturation
0.1538 0.1138 0.1232

74

Given the significant differences observed for the fluid distribution in the reservoir (sweep), the
proposed flow model will also impact inverse modeling studies through a model that represents
the physics at play more accurately. Accordingly, additional research and investigation is needed
to evaluate the impact of using this flow model in history matching applications.
In the horizontal direction, the flow of all phases is predominately co-current. As a result, there is
no need to include counter-current flow transitions in the horizontal direction (although it can
readily be included following the proposed approach). In the vertical direction, however, counter-
current flow is the dominant flow mode because of gravity segregation when gas is injected at the
bottom of the formation, and water that is injected at the top. Furthermore, the introduction of the
counter-current relative permeability model was not observed to significantly increase the
simulation time. Since we are using an IMPES-type simulator, flow transitions at each time step
is identified from the previous time step. In fully implicit simulation, additional studies are
required to ensure the numerical stability and efficiency. This study is the first attempt to
incorporate flow transitions and to study the potential impact on the design of WAG injection
processes. Further investigation is warranted to validate of the proposed flow model and its
implementation for three-phase flows, both from an experimental point of view and in the context
of fully implicit large-scale simulation.
Accurate evaluation of counter-current relative permeability curves is essential in improving the
physical models used in WAG displacement calculations. Lab experiments on cores have been
used in the past (Lelievre 1966, Bourbiaux and Kalaydjian 1990, Eastwood and Spanos 1991,
Bensten and Manai 1993, Al-Wadahi et al., 2000) to measure counter-current relative permeability
data. However, it is important to note that the results of such core experiments can be affected by
the boundary conditions imposed by the setup. As a result, the counter-current relative
75

permeabilities obtained from laboratory tests may not be reliable. As an alternative, counter-
current relative permeabilities can be calculated via the use of pore-scale modeling (Mani and
Mouhanty 1998, Cruichshank et al., 2002) where the impact of the imposed boundary conditions
can be evaluated.
4.6. Conclusions
Based on the calculation examples and analysis presented in this chapter, we arrive at the following
set of conclusions:
1. A reduction in relative phase mobilities due to transitions between co-current and counter-
current flow in the vertical direction plays an important role in the dynamics of WAG
injection processes, especially when smaller slug sizes are used.
2. Using a single set of relative permeability curves for the simulation of WAG processes (i.e.
only co-current flow including hysteresis and IFT scaling) may not provide sufficient
physical detail for accurate design of such processes.
The absolute permeability values for the reservoir model, used in this study, are between 7 mD
and 9000 mD (highly permeable sandstone). Accordingly, capillary pressures will be moderate.
Inclusion of capillary pressure in WAG displacement calculations does not significantly affect the
results. Accordingly, the calculation examples presented in this chapter, which do not include
capillary pressure explicitly, are expected to provide a good representation of the physics at play
during WAG injection processes in relevant settings.
76

Chapter 5 Stability Analysis for Multicomponent Mixtures
including Capillary Pressure
3

5.1. Introduction
Recent advances in drilling/completion technology have enabled hydrocarbon extraction from
tight unconventional resources such as shale gas, liquid-rich shales and tight oil systems. These
resources currently contribute significantly to the total oil and gas production in the US.
Hydrocarbon in organic-rich shales is one of the most significant unconventional resources (Yan
et al., 2013). Despite the large potential and the advancement in production technology for these
resources, gaps still exist between the physical models used in currently available simulation tools
and the well-documented additional complexity of fluid flow and mass transfer in micro- and
meso-porous materials. A central challenge is related to the understanding and modeling of phase
equilibrium in confined spaces, e.g. in pores at the nanometer scale, and to develop appropriate
and accurate tools/algorithms to be used for estimation of reserves and for forecasting of
production (Sandoval et al., 2014). Conventional experimental PVT analysis is usually performed
based on the assumption that the porous medium does not influence the phase behavior (Du and
Chu, 2012). However, fluid properties and VLE behavior in confined spaces depart substantially
from the corresponding bulk measurements (PVT cell) where vapor-liquid interface curvature can
be neglected (Wang et al., 2013).

3
This chapter was previously published as “Sherafati, M., & Jessen, K. (2017). Stability analysis for
multicomponent mixtures including capillary pressure. Fluid Phase Equilibria, 433, 56-66."
77

When the fluid is confined in small pores, significant interfacial curvatures can arise that results
in large a capillary pressure between the equilibrium liquid and vapor phases: A pressure difference
that is known to impact the VLE. It has, however, been demonstrated previously that this effect is
negligible for conventional reservoirs with pore spaces in the micron range (Shapiro and Stenby,
1997, Shapiro et al., 2000).
Current commercial reservoir simulators commonly ignore the effect of capillary pressure in VLE
calculations. However, experimental and modeling work presented over the past few years (e.g.
Du and Chu, 2012, Wang et al., 2013, Firincioglu et al., 2012, Nojabaei et al., 2013, Nojabaei et
al., 2014, Wang et al., 2014) have shown that neglecting capillary pressure in VLE calculations is
not a good assumption for unconventional reservoirs. Furthermore, a standard and reliable
measurement of fluid properties and VLE for confined fluids in ultra-tight rocks is still not
available (Du and Chu, 2012).
Numerous researchers have worked on analyzing phase envelopes of confined fluids. Nojabaei et
al. (2013) coupled capillary pressure and phase equilibrium calculations for binary and
multicomponent mixtures in confined pores and integrated their model into a compositional
reservoir simulator. Based on example calculations, they conclude that better estimates of the
producing gas-oil ratio (GOR) are achieved when the capillary pressure is included in the
equilibrium calculations. Du and Chu (2012) proposed a thermodynamic model to calculate the
confined PVT properties for Bakken crude oil systems considering the effects of capillarity in the
porous medium. Wang et al. (2013) evaluated the effect of capillarity on VLE calculations using
a Leverett J-function approach to approximate a saturation dependent capillary pressure. Alharthy
et al. (2013) proposed a correlation to shift the phase envelope depending on pore sizes, dividing
78

the pore space into macro-, meso- and micro-porosity, while Sandoval et al. (2016) proposed a
new algorithm to calculate phase envelopes including the effect of capillary pressure.
Despite a substantial amount of research on the topic of phase behavior of confined fluids, the
effect of capillarity in stability analysis has not been extensively studied/documented to the best
of our knowledge. In this chapter, we will start from Gibbs energy and the requirement for phase
stability and introduce modifications needed to account for capillarity in stability analysis. We
propose four algorithms for stability testing: One based on direct/accelerated substitution and one
based on minimization. The convergence behavior of the proposed algorithms is illustrated for a
range of relevant fluid systems and proximity to critical points. The stability testing is
demonstrated to be consistent with the phase boundaries predicted by phase envelope calculations
including capillarity. Accordingly, the proposed algorithms provide for reliable tools to study the
effect of capillary pressure on phase stability over a wide range of compositions and capillary radii
as needed in compositional simulation of unconventional reservoirs/formations.
5.2. Methodology
In this section, we discuss the stability analysis of a confined fluid and the related impact of
capillary pressure. Starting from Helmholtz energy, one can derive the equilibrium condition, in
terms of the chemical potential, µ, for a system that includes the contribution from surface energy
(Firoozabadi, 1999)
𝜇 𝑖 𝑙 (𝑇 ,𝑝 𝑥 ,𝒙 )=𝜇 𝑖 𝑣 (𝑇 ,𝑝 𝑦 ,𝒚 ) , 𝑖 =1,..,𝑛 𝑐 ,
(5.1)
In Eq. 5.1, the phase pressures (p
l
and p
v
) are related through the Young-Laplace equation (written
for a spherical interface)
79

𝑝 𝑣 −𝑝 𝑙 =
2𝜎 𝑅 .
(5.2)
The interfacial tension, σ, is a given function of phase compositions (x and y), temperature, T,
phase pressures and the radius of curvature of the interface, R. In a capillary tube (pore) the
capillary pressure, Pc, is evaluated from
𝑃 𝑐 =𝑝 𝑣 −𝑝 𝑙 =
2𝜎 𝑐𝑜𝑠 𝜃 𝑟 𝑐 , (5.3)
where rc is pore radius and 𝜃 is the contact angle as measured through the denser phase. Although
porous rocks are made up by complex networks of non-cylindrical pores of variable radii, we use
here a single value (average) representation of the pore radius for simplicity. We also assume, by
default, that the rock is strongly liquid wet (𝜃 =0). We use the correlation of Macleod (1923) and
Sugden (1924) to calculate interfacial tension (IFT),
𝜎 =[∑𝜒 𝑖 (𝑥 𝑖 𝜌 𝐿 −
𝑛 𝑐 𝑖 =1
𝑦 𝑖 𝜌 𝑉 )]
𝐸 ,
(5.4)
where 𝜌 𝐿 and 𝜌 𝑉 are phase molar densities, 𝑥 𝑖 and 𝑦 𝑖 are phase compositions, 𝜒 𝑖 is the component
parachor and 𝐸 is an (adjustable) exponent. Researchers have suggested a variety of values for the
exponent E. Weignaug and Katz (1943) suggested that E=4, while Hough-Stegemeier (1961)

modified the exponent to E=3.66 for low IFT systems. In 1984, Lee and Chien (1984)

suggested
a value of E=3.91 using critical scaling theory of IFT to evaluate the dependency on temperature.
Danesh et al. (1991) presented a density-dependent model for exponent E. Schechter and Guo
(1998) recommended E= 3.88 based on experimental data for hydrocarbon systems. We use E =
80

4 (following Weignaug and Katz). From Eq. 5.4, one can calculate the interfacial tension, and
hence the capillary pressure, for a given set of phase compositions and densities of liquid and vapor
phases given a pore radius and a contact angle. This allows for a simple algorithm to update the
phase pressures by direct substitution as discussed in more detail later.
5.2.1. Formulation of the stability analysis problem
The tangent plane distance (TPD) of Gibbs has been used repeatedly in developing algorithms for
stability analysis of unconfined fluids (Michelsen, 1982). To apply this idea for confined fluids,
we make the assumption that the volume (or mass) contribution to changes in Gibbs free energy
is larger than the interface contribution. In terms of nucleation theory, this is equivalent to
assuming that the volume (radius) of any trial phase exceeds the radius of a critical nucleus. Based
on this assumption, we can write the TPD at isothermal conditions in the following form
𝑇𝑃𝐷 (𝒘 )= ∑𝑤 𝑖 (𝜇 𝑖 𝑤 (𝒘 ,𝑝 𝑤 )−𝜇 𝑖 𝑧 (𝒛 ,𝑝 𝑧 )) ,
𝑛 𝑐 𝑖 =1
(5.5)
where z is the feed composition existing at pressure p
z
and w is a trial phase composition existing
at p
w
. The assumption of neglecting the interfacial contribution to changes in Gibbs energy, as
outlined above, is consistent with the formulation of Sandoval et al. (2016) for tracing phase
boundaries. This is seen from Eq. 5.5 where the TPD will equal zero on the phase boundary where
the incipient phase (w) exists is in equilibrium with the bulk phase (z).
At this stage, we can follow the approach of Michelsen (1982), and introduce a modified TPD via
the variable transformation
81

𝑙𝑛𝑊 𝑖 =𝑙𝑛𝑤 𝑖 −𝑘 , (5.6)
where k is the reduced TPD at the stationary points of Eq. 5.5. The modified reduced TPD, tm, for
a system including capillary pressure can then be written in terms of fugacity coefficients
𝑡𝑚 (𝑾 )=1+∑𝑊 𝑖 (ln𝑊 𝑖 +ln𝜑̂
𝑖 𝑤 +ln𝑝 𝑤 −𝑑 𝑖 −1)
𝑁𝑐
𝑖 =1
, (5.7)
with
𝑑 𝑖 =ln𝑧 𝑖 +ln𝜑̂
𝑖 𝑧 +ln𝑝 𝑧 , 𝑖 =1,..,𝑁𝑐 (5.8)
and
𝑝 𝑤 −𝑝 𝑧 = ±𝑃 𝑐 . (5.9)
If the feed composition (z) is assumed to be a liquid and the trial phase (w) is vapor, the ± sign in
Eq. 5.9 is positive and if the feed composition is vapor and the trial phase is liquid, the sign
becomes negative (assuming that the liquid is the wetting phase). The focus of stability analysis is
to find the stationary points of Eq. 5.7. If the value of the tm is greater than zero for all stationary
point, the feed composition is stable, and if not, we have demonstrated that the feed is unstable
and will form at least two phases. In order to find the stationary points of Eq. 5.7, we need to solve
𝜕𝑡𝑚 𝜕 𝑊 𝑖 =0, 𝑖 =1,..,𝑛 𝑐 , (5.10)
corresponding to the following system of nonlinear equations
82

ln𝑊 𝑖 +ln𝜑̂
𝑖 𝑤 −𝑑 𝑖 +ln𝑝 𝑤 =0, 𝑖 =1,..,𝑛 𝑐 . (5.11)
5.2.2. Solution Methods
To find the stationary points of the modified tangent plane distance, Eq. 5.11 must to be solved for
the independent variables (W, p
w
). Due to its simplicity, direct substitution is an obvious choice
for solving the above equations where subsequent iterates of W are determined from
ln𝑊 𝑖 𝑡 +1
=𝑑 𝑖 −ln𝜑̂
𝑖 𝑤 ,𝑡 −ln𝑝 𝑤 ,𝑡 ,𝑖 =1,..,𝑛 𝑐 . (5.12)
The trial phase pressure (p
w
) is initialized as equal to the feed pressure (p
z
). In subsequent
iterations, the capillary pressure is calculated from Eq. 5.3 based on phase compositions and the
phase pressures. The trial phase pressure is then updated accordingly by direct substitution. It is
important to note that the capillary pressure is calculated at the end of the iteration, such that the
capillary pressure updates are lagging one iteration behind the direct substitution of the trial-phase
compositions.
In order to locate the stationary points of the modified tangent plane distance, it is necessary
to test multiple initial estimates. For two-phase VLE problems, we select two estimates: One
vapor-like initial estimate and one liquid-like. Wilson’s equation provides a reasonable
approximation for initial estimates
𝐾 𝑖 =
𝑃 𝑐 ,𝑖 𝑝 𝑧 exp[5.42(1 + ω
𝑖 )(1−
𝑇 𝑐 ,𝑖 𝑇 )] , (5.13)
where 𝑃 𝑐 ,𝑖 , 𝑇 𝑐 ,𝑖 and ω
𝑖 are the critical properties and acentric factor of the component i and 𝐾 𝑖 is
the equilibrium constant. The vapor-like and liquid-like initial estimates (W) are then calculated
from
83

𝑊 𝑖 =𝐾 𝑖 𝑧 𝑖 , (5.14)
and
𝑊 𝑖 =
𝑧 𝑖 𝐾 𝑖 ⁄ .
(5.15)
To study the convergence behavior, we define an error for the trial-phase composition as the norm
of the vector e
𝑒 𝑖 =ln𝑊 𝑖 +ln𝜑̂
𝑖 𝑤 −𝑑 𝑖 +ln𝑝 𝑤 , (5.16)
and consider the solution converged when the ||e|| is less than 10
−8
.
To enhance the convergence behavior of the direct substitution method, we utilize the acceleration
approach (Dominant Eigenvalue Method) of Orbach and Crowe (1971) with extrapolation
attempted after each five direct substitutions.
Direct substitution and accelerated direct substitution are rapidly convergent for most relevant
settings. However, near critical points the convergence rates can be very slow and a second order
method may be preferred. To formulate a second order convergent method, we use the variable
transformation suggested by Michelsen (1982)

and arrive at the following set of independent
variables
𝐯 =[2√𝑊 1
,2√𝑊 2
,..,2√𝑊 𝑛𝑐
]
𝑇 .
(5.17)
The gradients of the Eq. 5.7 with respect to vi is written as
84

𝑔 i
=
𝜕𝑡𝑚 𝜕 𝑣 𝑖 =√𝑊 𝑖 (ln𝑊 𝑖 +ln𝜑̂
𝑖 𝑤 −𝑑 𝑖 +ln𝑝 𝑤 ) , (5.18)
while the Hessian matrix is calculated as
𝐻 𝑖𝑗
=
𝜕 2
𝑡𝑚 𝜕𝑣
𝑖 𝜕𝑣
𝑗 =𝛿 𝑖𝑗
+√𝑊 𝑖 𝑊 𝑗 𝜕𝑙𝑛 𝜑̂
𝑖 𝑤 𝜕𝑊
𝑗 +
1
2
𝑔 𝑖 𝑣 𝑖 𝛿 𝑖𝑗
𝑖 ,𝑗 =1,..,𝑁 𝑐 . (5.19)
The stationary points of Eq. 5.7 can be found via the restricted step algorithm suggested by Hebden
(1973)
(𝑯 +𝛽 𝑰 )∆𝐯 +𝒈 =𝟎 , (5.20)
to handle cases where the Hessian matrix is indefinite (Michelsen, 1982). In Eq. 5.20, β is found
iteratively to arrive at ||∆𝒗 || equal to a specified maximum step length (see Hebden, 1973 for
additional discussion of step adjustments). In this formulation, the trial-phase pressure is not
included as an independent variable and must be updated separately via direct substitution (as
discussed above) or by a Newton iteration on the residual form of Eq. 5.9. Here, we use the Newton
update of the trial-phase pressure to promote convergence of the overall (sequential) scheme. We
note that the Newton update only requires the pressure derivative of the molar volume of the trial
phase (available from the calculation of fugacity coefficients).
Alternatively, the pressure of the trial phase can be treated as an independent variable by forming
a new variable vector
𝐬 =[2√𝑊 1
,2√𝑊 2
,..,2√𝑊 𝑛𝑐
,𝑝 𝑤 ]
𝑇 .
(5.21)
85

With this set of independent variables, we define a new set of equations to be solved as follows
𝐹 (𝒔 )={
𝑔 𝑖 𝑖 =1,..,𝑁 𝑐
𝑃 𝑐 ±(𝑝 𝑤 −𝑝 𝑧 ) 𝑖 =𝑁 𝑐 +1
(5.22)
The Jacobian of the above function can be defined as follows
𝐽 𝑖 ,𝑗 =
{

𝐻 𝑖𝑗
𝑖 =1,..,𝑁 𝑐 , 𝑗 =1,..,𝑁 𝑐

√𝑊 𝑖 (
𝜕𝑙𝑛 𝜑̂
𝑖 𝑤 ∂𝑝 𝑤 +
1
𝑝 𝑤 ) 𝑖 =1,..,𝑁 𝑐 , 𝑗 =𝑁 𝑐 +1

√𝑊 𝑗 𝜕 𝑃 𝑐 𝜕𝑊
𝑗 𝑖 =𝑁 𝑐 +1, 𝑗 =1,..,𝑁 𝑐

𝜕 𝑃 𝑐 ∂𝑝 𝑤 ±1 𝑖 =𝑁 𝑐 +1, 𝑗 =𝑁 𝑐 +1

(5.23)

Hebden’s restricted-step algorithm can again be used to solve the above set of equations
(𝑱 +𝛽 𝑰 )∆𝐬 +𝑭 =𝟎 . (5.24)
In the formulation of Eq. 5.22 and Eq. 5.23, both the tangent plane distance and the error in the
capillary pressure are minimized simultaneously. However, the Jacobian matrix (Eq. 5.23) is less
attractive to evaluate, due to the dependency of capillary pressure (via IFT) on phase compositions
and pressures. In addition, the factorization of the J is more CPU intensive than the factorization
of H (non-symmetric versus symmetric matrix). We discuss this further in the results section.
In the following section, we compare the performance of the 4 different algorithms discussed
above including: a) Direct substitution, b) Accelerated direct substitution, c) Sequential
86

minimization of TPD with separate Newton update of trial-phase pressure and d) Minimization
of the TPD with a full set of independent variables.

5.3. Results
First, we introduce the fluid description used in our example calculations and report the relevant
phase envelopes. We then report results from stability analysis calculations with emphasis on
convergence behavior near phase boundaries of oil and condensate systems. Finally, we consider
the effects of feed composition and pore radius on phase boundaries and stability analysis
calculations.
5.3.1. Fluid Description
In this study, we use a 15-component fluid description that is representative of a sample from a
realistic reservoir fluid. Table 5.1 lists the components, their critical properties, acentric factors,
molecular weights, volume shift parameters, parachors and composition of oil and condensate
systems. The last 5 components of the fluid composition are pseudo-components created during
characterization of the C7+ fraction of the reservoir fluid. The Soave-Redlich-Kwong (SRK)
equation of state was used in all calculations. Table 5.2 reports the non-zero binary interaction
coefficients for the components.
Table 5.1 Fluid properties and compositions
Comp
Tc
[K]
Pc
[atm]
Acentric
factor [-]
M
W

[g/mole]
Volume
shift [-]
Parachor
Oil
composition
Condensate
composition
N2 126.2 33.6 0.04 28.016 -0.1927 43.85 0.0018 0.0106
CO2 304.2 72.9 0.228 44.01 -0.0817 82.24 0.0082 0.0158
CH4 190.6 45.4 0.008 16.043 -0.1595 71.52 0.2292 0.7702
C2 305.4 48.2 0.098 30.069 -0.1134 113.71 0.0721 0.0939
87

C3 369.8 41.9 0.152 44.096 -0.0863 151.14 0.0737 0.0522
iC4 408.1 36.0 0.176 58.123 -0.0844 181.78 0.0158 0.0073
nC4 425.2 37.5 0.193 58.123 -0.0675 191.03 0.0523 0.0202
iC5 460.4 33.4 0.227 72.15 -0.0608 224.14 0.0225 0.0056
nC5 469.6 33.3 0.251 72.15 -0.039 231.73 0.036 0.0078
C6 507.4 29.3 0.296 86.177 -0.008 274.03 0.0484 0.0061
PS1 565.85 29.32 0.34612 120.312 -0.01648 331.01 0.196107 0.0099
PS2 683.44 19.84 0.57838 207.159 -0.04094 546.37 0.113893 0.0004
PS3 798.99 13.15 0.90175 354.274 -0.0616 884.73 0.066598 3.041E-6
PS4 899.68 9.57 1.19183 574.812 -0.05768 1391.97 0.041047 1.298E-8
PS5 1013.31 7.57 1.39383 1055.440 0.01401 2497.41 0.022355 1.987E-11

Table 5.2 Non-zero binary interaction coefficients
N2 CO2 CH4
N2
CO2 0.017
CH4 0.0311 0.12
C2 0.0515 0.12
C3 0.0852 0.12
iC4 0.1033 0.12
nC4 0.08 0.12
iC5 0.0922 0.12
nC5 0.1 0.12
C6 0.08 0.12
PS1 0.08 0.1 0.028349
PS2 0.08 0.1 0.044813
PS3 0.08 0.1 0.062256
PS4 0.08 0.1 0.077679
PS5 0.08 0.1 0.0952

5.3.2. The Phase Envelope
We first present the phase envelope for the 15-component fluid using the oil feed composition.
The phase envelope was initially constructed for two scenarios: a) without confinement and b)
with confinement corresponding to an average pore size of 10nm. The construction of phase
envelopes follows the work of Sandoval et al. (2016). Figure 5.1 shows the phase envelopes and
88

demonstrates that capillary pressure effects suppress the bubble point pressures for the entire
temperature range of the bubble point branch. Close to the critical point, the two curves converge,
since the two phases become identical (at the critical point) and the interfacial tension approaches
zero. In the dew point region, the saturation pressure is increased by capillary effects and the
cricondentherm is shifted to the right.

Figure 5.1 Phase envelope of the 15-component oil with and without capillary pressure (rc = 10nm). Triangle =
critical point.
5.3.3. Stability Testing
The phase envelope constructed in Figure 5.1 serves as the basis for testing the stability analysis
algorithms outlined above. The critical temperature and pressure of the fluid sample shown in
Figure 5.1 are 721.38 K and 113.16 atm, respectively. Stability analysis was performed along
constant temperature lines (383.15K, 700K and 780K) to identify the lowest pressure at which the
feed composition is stable at a given temperature (last positive/zero value of tm): This pressure
100 200 300 400 500 600 700 800 900
0
20
40
60
80
100
120
140
Temperature (K)
Pressure (atm)

normal envelope
P
Liquid
P
Gas
unstable feed
780 800 820 840
30
40
50
60

89

must be identical to the saturation pressure dictated by the phase envelope at that temperature.
Figure 5.1 reports the results of stability analysis for three different regions of the phase envelope.
At each temperature, the lowest pressure where the given feed is stable is marked with a black
square. First, stability analysis was performed on the bubble-point side of the PT phase envelope
(T=383.15K). The highest pressure corresponding to a negative tangent plane distance is in
agreement with the suppressed liquid phase saturation pressure at the given temperature. The
stability analysis indicates, as expected, an unstable vapor-like trial phase.
Next, we choose a point relatively close to the critical point (T=700K). In this region, the direct
substitution method takes considerably longer to converge compared to the minimization
approaches. The vapor-like trial phase converges to a negative value of TPD, while the liquid-like
trial phase arrives at the trivial solution. Finally, we perform stability analysis on the dew-point
branch of the envelope (T=780K), close to the cricondentherm. The highest pressure resulting in
a negative TPD is again in agreement with the relevant phase envelope. Here, the liquid-like trial
phase results in a negative TPD. Figure 5.2 reports the reduced TPD (tpd=TPD/RT) at the
stationary point as a function of pressure along constant temperatures of T=383.15 K, T=700 K
and T=780 K.
90

Figure 5.2 Reduced tangent plane distance vs. pressure at constant temperatures of T=383.15K (top), T=700K
(middle) and T=780K (bottom).
Let us first examine the tangent plane distance for T=383.15 K on the bubble point side, with liquid
as the reference phase and vapor as the trial phase. Starting at a pressure higher than the saturation
pressure, tpd is positive at the stationary point before reaching the saturation pressure (tpd = 0),
indicating a stable feed. As we cross the phase boundary, tpd at the stationary points becomes
80 81 82 83 84 85 86 87 88 89 90
-0.04
-0.02
0
0.02
0.04
T=383.15K
pressure(atm)
tpd
116 117 118 119 120 121 122 123 124 125 126
-4
-2
0
x 10
-4
T=700K
pressure(atm)
tpd
80 81 82 83 84 85 86 87 88 89 90
-6
-4
-2
0
2
x 10
-3
T=780K
pressure(atm)
tpd
91

negative, an indication of an unstable feed. By further reducing the pressure we find even lower
values of the tpd at the stationary points inside the two-phase region of the phase envelope.
Near the critical point (T=700K), and at pressure above the phase boundary, all solution
approaches arrive, as expected
23
, at the trivial solution (z = w). Right above the phase boundary,
small positive tpd values are observed, corresponding to a narrow “shadow region” (Rasmussen et
al., 2006) in the vicinity of the critical point. When the saturation pressure is reached, negative
values of tpd are observed by further reducing the pressure. On the dew point branch (T=780K),
the reference phase (z) and the trial phase (w) are assumed to be vapor and liquid, respectively. At
pressures well above the dew-point pressure, all algorithms arrive at the trivial solution (with
tpd=0) before reaching the so-called “shadow region” where positive values of tpd are found right
above the saturation pressure. Again, as the pressure is reduced further, tpd transitions to negative
values just below the saturation pressure.
The results presented above demonstrate that the proposed algorithms for stability testing
accurately identify the phase boundaries over the full range of relevant temperatures. For most
calculations (with the exception of the vicinity of the critical point) a modest number of iterations
is required to converge to the stationary points. A comparison of the convergence behavior of the
different algorithms is provided below.
5.3.4. Convergence Behavior
The calculation examples provided above illustrates the accuracy and robustness of the suggested
algorithms for stability analysis including effects of capillary pressure. All calculations to locate
the stationary points of Eq. 5.7 were performed four times: Once using direct substitution, once
using accelerated substitution (extrapolation based on the largest eigenvalue), once using the
92

minimization approach with separate Newton update of the trial-phase pressure (Eqs. 5.9, 5.17,
5.19) and once using the minimization approach including the trial-phase pressure as an
independent variable (Eqs. 5.21, 5.22, 5.23). All algorithms identify the same values of TPD at the
stationary points. The main difference between these approaches is the rate of convergence. The
left tab of Figure 5.3 provides a comparison of the convergence behavior of the four algorithms at
T= 383.15 K and 83 atm (just below the bubble point branch) and an average pore radius of 10
nm.

Figure 5.3 Convergence behavior of direct substitution, accelerated substitution and minimization with lagged pw
updates (Minimization 1) and minimization with the full Jacobian (Minimization 2) at a) T=383.15 K, p = 83 atm
and b) T=700 K, p = 119 atm for rc=10 nm
As can be seen, all algorithms perform well and converge after a moderate number of iterations.
The right tab of Figure 5.3 provides a comparison of the convergence rates at T=700 K (in the
vicinity of the critical point). Near the critical point, the minimization approaches require a
0 5 10 15
-25
-20
-15
-10
-5
0
5
iteration
log(||err||)
T=383.15 K
0 10 20 30
-30
-25
-20
-15
-10
-5
0
5
iteration
log(||err||)
T=700 K

Direct Substitution
Accelerated Direct Substitution
Minimization 1
Minimization 2
93

considerably lower number of iterations to converge than both the direct substitution approach and
the accelerated substitution approach. The accelerated substitution, as tested here with a single
eigenvalue, is not sufficient to promote the convergence of the direct substitution significantly
(since the eigenvalue approaches unity near the critical point). At these conditions, it is preferable
to use a minimization approach (also pointed out by Michelsen (1982)). In practical applications,
one may start with direct substitution and transition to a minimization approach if the convergence
is progressing slowly.
As discussed above, the update of the capillary pressure (and trial-phase pressure) is performed
separately from the update of trial-phase compositions in the first three algorithms. Near the critical
point, the algorithms proposed here will behave (become) identical to the original work of
Michelsen (1982)

and preserve the favorable convergence behavior of the minimization
approaches. When the capillary pressure (and trial-phase pressure) is updated outside the main
iterative scheme, one would expect the convergence of the trial-phase pressure (and capillary
pressure) to be slower than the convergence of the trial-phase compositions (particularly for the
minimization approach with lagging trial-phase pressure). To investigate the convergence rate of
the capillary pressure in more detail, we define the error in the capillary pressure (epc) relative to
the converged solution (Pc
*
) based on a tolerance of 1.1e-15.
𝑒 𝑝𝑐
=𝑎𝑏𝑠 (𝑃 𝑐 −𝑃 𝑐 ∗
) . (5.25)
Figure 5.4 depicts the convergence behavior as observed for direct substitution, accelerated
substitution and the two minimization approaches.
94

Calculations at the two temperatures are representative of the bubble point branch (T=383.15K,
p=83 atm) and the near-critical region (T=700K, p=119 atm) of the phase envelope, respectively.
A comparison of Figs. 5.3 and 5.4 demonstrates that the error in the capillary pressure is reduced
rapidly for the first three approaches as well: For both calculation examples (383.15K and 700K)
reported in Fig. 5.4 the capillary pressure is stabilized (small error compared to the converged
solution) at rates that are similar to the reduction of the errors in trial-phase compositions. As a
result, a separate update of the capillary pressure does not significantly affect the convergence of
the stability analysis calculations and the proposed algorithms can be implemented via modest
revision of existing codes for stability testing.

Figure 5.4 Convergence of the capillary pressure for direct substitution, accelerated substitution, minimization with
lagged pc updates (Minimization 1) and minimization with the full Jacobian (Minimization 2) at T=383.15K,
p=83atm and T=700K, p=119atm for rc=10 n
0 5 10 15
-30
-25
-20
-15
-10
-5
0
5
iteration
log(p
c
error)
T=383.15 K
0 10 20 30
-30
-25
-20
-15
-10
-5
0
5
iteration
log(p
c
error)
T=700 K

Direct Substitution
Accelerated Direct Substitution
Minimization 1
Minimization 2
95

A comparison of the convergence rates for the two minimization approaches in Figs. 5.3 and 5.4
demonstrates the expected quadratic convergence behavior for the algorithm using a full set of
independent variables, while a separate update (Newton) of the trial-phase pressure results in a
super linear convergence. Additional computational work (a factor of two) is required for
factorization of the matrix J given in Eq. 5.23 (non-symmetric matrix) relative to the factorization
of the symmetric matrix H given in Eq. 5.19. In addition, working with a full set of independent
variables requires the evaluation of compositional and pressure derivatives of the capillary
pressure in Eq. 5.23. Accordingly, the algorithm with separate update of the trial-phase pressure
(minimization 1) provides for the more efficient approach for stability testing when capillary
pressure is included.
5.3.5. Effect of Pore Radius
In the calculation examples presented above, the capillary pressure calculations, based on Eq. 5.3,
were performed for an average radius. In practice, however, a wide range of pore radii is observed
in tight formations. The porosity of the Middle Bakken region is approximately 6% and the
corresponding oil effective permeability is in the range of 30 to 100 nanodarcies. From the Kozeny-
Carman equation (Kozeny, 1927, Carman, 1997) Nojabaei et al. (2013) concluded that the
corresponding pore-throat radius should be in the range of 10-40 nm. In the calculation examples
provided above, we used a constant pore radius of 10 nm. Here, we examine the effect of pore
radius on the calculated phase boundaries to further test the proposed algorithms for stability
testing by repeating the example calculations for a range of pore radii of 3-50 nm. We present the
calculation results in 2 separate figures: one for pore radii smaller than 10 nm and one for pore
radii larger than 10 nm. Figure 5.5 reports the phase envelopes along with phase boundaries
detected by stability analysis at T=383.15 K for a capillary radius of 3, 5 and 10 nm.
96

Figure 5.5 Detection of phase boundaries by stability analysis at different capillary radii (rc≤10nm). Triangle =
critical point.
As the pore radius decreases, the depression of saturation pressures becomes more significant.
However, stability testing is invariably able to detect the phase boundary without any convergence
problems. Although we report results for calculation with pore radii down to 3 nm, continuum-
scale considerations may not be suitable at this scale and more detailed analysis using molecular
dynamics simulation may be more appropriate (Evans et al., 1986, Restagno et al., 2000,
Giovambattista et al., 2006, Jin and Firoozabadi, 2016). Next, we repeat the calculations at
T=383.15 K for capillary radii of 10, 20 and 50 nm as reported in Figure 5.6. As expected, when
the average pore radius is increased, the phase envelope including capillary pressure approaches
that of the bulk fluid. The differences between the two curves are already minimal for a capillary
radius of 50 nm.

100 200 300 400 500 600 700 800 900
0
20
40
60
80
100
120
140
Temperature (K)
Pressure (atm)

P
l
, r
c
=10nm
P
g
, r
c
=10nm
P
l
, r
c
=5nm
P
g
, r
c
=5nm
P
l
, r
c
=3nm
P
g
, r
c
=3nm
97

Figure 5.6 Detection of phase boundaries via stability analysis at different capillary radii (rc>10nm). Triangle =
critical point.
Figure 5.7 reports the calculated saturation pressures (detected via stability testing along a constant
temperature line in p-T space) for varying capillary radii at a constant temperature of T=383.15 K.
The saturation pressure of the bulk fluid at the given temperature is 99 atm. We observe that the
saturation pressure is reduced to less than half (46 atm) for a capillary radius of 3 nm. Changes in
the saturation pressure with respect to capillary radius is non-linear and increases exponentially
with decreasing capillary radius.
100 200 300 400 500 600 700 800 900
0
20
40
60
80
100
120
140
Temperature (K)
Pressure (atm)

P
l
, r
c
=50nm
P
g
, r
c
=50nm
P
l
, r
c
=20nm
P
g
, r
c
=20nm
P
l
, r
c
=10nm
P
g
, r
c
=10nm
360 380 400
80
85
90
95

98

Figure 5.7 Calculated saturation pressure vs capillary radius at T=383.15 K
5.3.6. Effect of Composition
Next, we consider a condensate reservoir with the fluid description and feed composition given in
Table 5.1, assuming an average pore radius of 10 nm. Figure 5.8 depicts the phase envelope and
the phase boundaries as identified by stability testing at 4 different temperatures. In comparison to
the heavier reservoir fluid considered above (oil), the changes in saturation pressures along the
bubble point branch are smaller for the condensate fluid description. This is mainly because the
interfacial tension is smaller due to smaller density difference between the two equilibrium phases.
As a result, smaller values of capillary pressure are observed that, in turn, result in smaller changes
in the saturation pressure relative to the bulk fluid. This is in agreement with previous results
reported for multicomponent systems based on the multicomponent Kelvin equation (Shapiro et
al., 2000). The critical point for this fluid composition is at T=247.63 K and P=111.94 atm. The 4
different temperatures along the phase envelope represent the bubble point branch, the near-critical
5 10 15 20 25 30 35 40 45 50
40
50
60
70
80
90
100
r
c
(nm)
P
sat
(atm)
99

region and the dew point branch including the cricondentherm. In all cases, the stability analysis
detects the phase boundaries accurately.

Figure 5.8 Detection of phase boundaries via stability analysis for a condensate system (rc=10 nm).
Figure 5.9 explores the reduced TPD at the stationary points as a function of pressure at constant
temperatures of T=250 K, T=300 K and T=365 K for the condensate system. Similar to Figure 5.2,
we observe, near the critical point, that the stability analysis finds the trivial solution until the
pressure is close, but still greater than the saturation pressure. Very close to the phase boundary
(0.005 atm above psat), a positive value of tpd (1.25e-8) is found. This is expected as the shadow
region is very narrow near the critical point. Finally, as we cross the phase boundary, we start
observing increasingly negative values of the tpd. For the temperatures on the dew point branch
(away from the critical point), we observe a change of sign in the reduced tangent plane distance
more clearly as we cross the phase boundary.
100 150 200 250 300 350 400
0
50
100
150
Temperature (K)
Pressure (atm)

normal envelope
critical point
P
Liquid
P
Gas
unstable feed
100

Figure 5.9 Reduced tangent plane distance vs. pressure at constant temperature for T=250K (top), T=300K (middle)
and T=365K (bottom) for condensate fluid
The proposed algorithms do not allow for negative values of vapor and liquid phase pressures. As
a result, special care must be taken in the lower portion of the dew point branch. Sandoval et al.
(2016) suggested an algorithm for phase envelope calculations allowing for negative phase
pressures, which can be incorporated in the algorithms for stability analysis proposed here.
110 111 112 113 114 115 116 117 118 119 120
-4
-2
0
x 10
-4
T=250K
pressure(atm)
tpd
145 146 147 148 149 150 151 152 153 154 155
-4
-2
0
2
4
x 10
-3
T=300K
pressure(atm)
tpd
75 76 77 78 79 80 81 82 83 84 85
-0.02
-0.01
0
0.01
0.02
T=365K
pressure(atm)
tpd
101

Negative phase pressures are, however, unlikely to be encountered during compositional
simulation at realistic reservoir conditions.
5.4. Summary
Accurate detection of phase stability is crucial in the modeling and simulation of production from
oil and gas reservoirs. Here we propose a set of robust, efficient and accurate algorithms for
stability testing, including capillary pressures, of relevance to compositional simulation of tight
oil/gas resources. The efficiency and accuracy of the proposed algorithms were demonstrated for
a range of relevant temperature and pressure conditions in oil and condensate systems. The
formulation and proposed solution strategies for stability testing provide a simple path to
incorporate effects of capillary pressure in existing codes for stability testing.
For large values of the capillary pressure, the saturation pressures deviate significantly from the
bulk fluid calculations for both oil and condensate systems. Along the bubble point branch of the
phase envelope, saturation pressures are suppressed while, along the dew point branch, the
calculated saturation pressures are increased compared to the bulk behavior (for both oil and
condensate systems). In the vicinity of the critical point, stability analysis calculations including
capillary pressure approach calculations for bulk-fluids, as the interfacial tension approaches zero.
The analysis presented here assumes a single average pore size to represent unconventional
resources. In practice, calculations will be more complicated as a wide range of pore sizes are
observed in reservoir rocks. However, a single average value of the pore radius will allow
engineers to gauge the potential impact of confinement on saturation pressures and related
production operations.
102

Chapter 6 Diffusion and Matrix-Fracture Interactions
During Gas Injection in Fractured Reservoirs
6.1. Introduction
Molecular diffusion, viscous displacement and gravity drainage contribute to recovery during gas
injection in naturally fractured reservoirs. The relative significance of these mechanisms depends
on several factors including matrix permeability, fracture intensity, fluids properties, injection rate
and reservoir pressure and temperature.
Viscous flow does not contribute directly to the recovery of oil because the injected gas channels
through high-permeability fractures that only represent a small fraction of the total pore volume.
Gravity drainage, which is driven by density difference between oil and gas, plays a significant
role when matrix permeability is high. Molecular diffusion may become the dominant recovery
mechanism in cases with low matrix permeability and high fracture intensity.
Contrary to conventional reservoirs where the impact of molecular diffusion is generally small,
molecular diffusion can play an important role in fractured reservoirs because of the large surface
area available for diffusion. Diffusion of gas components from a fracture into the matrix extracts
oil components from matrix and delays, to some extent, the gas breakthrough via high-permeability
fractures. As a result, both displacement and sweep efficiencies are increased. Different
investigators have demonstrated the role of molecular diffusion in fractured reservoirs (e.g. da
Silva and Belery, 1989; Hu and Whitson, 1991; Darvish et al., 2006; Hoteit and Firoozabadi, 2009;
Vega et al., 2010; Jamili et al., 2011).
103

Molecular diffusion is driven by gradients in concentrations (or more generally by gradients in the
chemical potential). Diffusion will oftentimes not play a significant role in regions of a porous
medium where fluids velocities are high (e.g. fractures) because the characteristic time for
diffusion is relatively large. However, in regions of the porous media where flow velocities are
low (e.g. matrix blocks) diffusion can dominate mass transfer over viscous flow (Perkins and
Johnston, 1963; Shojaei et al., 2012).
Single-phase multicomponent molecular diffusion can be modeled by different approaches
including the classical Fick’s law, the generalized Fick’s law and the Maxwell-Stefan (MS) model.
In the classical Fick’s law approach, the diffusive flux of each component is only a function of its
own concentration gradient. In other words, the interactions among different species (e.g. dragging
effects) are neglected. The generalized Fick’s law approach, in contrary, takes into account the
interactions between components; i.e. the diffusive flux of each component depends on the
concentration gradients of other components as well. In the MS model, diffusion is driven by a
gradient in the chemical potential and component interactions are accounted for using component
velocities (friction). It can be shown that under certain conditions the generalized Fick’s law and
MS model are equivalent (Taylor and Krishna, 1993).
In most reservoir simulation models (e.g. da Silva and Belery, 1989; Coats, 1989; Darvish et al.,
2006; Vega et al., 2010; Jamili et al., 2010), multicomponent molecular diffusion is modeled using
a classical Fick’s law approach with effective diffusion coefficients. This means the diffusive flux
of each component will always be in the opposite direction of its concentration gradient. However,
it has been shown that because of the dragging effects, diffusion may occur from a region of low
to high concentration, known as reverse diffusion (Duncan and Toor, 1962). Therefore,
104

considerable errors can be introduced in the simulation results if the dragging effects (off-diagonal
diffusion coefficients) are ignored.
The gas-oil diffusion at the fracture-matrix interface is usually modeled by assuming an average
composition at the interface (equipartition hypothesis). da Silva and Belery (1989) argue that for
practical purposes, this is a reasonable assumption for matrix blocks surrounded by interconnected
fractures. However, the equipartition hypothesis and associated calculation of gas-oil diffusion
coefficients do not have a rigorous physical basis.
In chemical engineering literature, the gas-oil diffusion has been modeled based on film theory in
which oil and gas are assumed to be in equilibrium at the interface, and component fluxes are
required to be continuous across the interface (e.g. Krishna and Standart, 1976; Taylor and
Krishna, 1993). This approach has also been used in petroleum engineering literature to model
gas-oil diffusion in laboratory scale experiments (e.g. Hu and Whitson, 1991; Irani et al., 2009;
Guo et al., 2009; Haugen and Firoozabadi, 2009; Hoteit, 2013; Moortgat and Firoozabadi, 2013).
In this chapter, we present a dual-porosity model for field-scale simulation of gas injection in
fractured reservoirs in which the generalized Fick’s law is used for diffusion; and gas-oil diffusion
at the fracture-matrix interface is modeled based on film theory. Diffusion coefficients are
calculated using the MS model and are pressure, temperature and composition dependent. Our
approach is compared with the conventional approach via a set of calculations examples. The
significance of dragging effects of M-F transfer function are also investigated.
6.2. Mathematical Model
We use a dual-porosity approach to model gas injection in fractured reservoirs. Dual-porosity
models assume that there are two communicating domains in a fractured reservoir: a flowing
105

domain (fracture) and a stagnant domain (matrix). Mass (or volume) balance equations are
formulated independently for these two domains. Transfer of fluids between these two overlapped
domains is accounted for using a source/sink term in the mass (or volume) balance equations
(Barenblatt et al., 1960; Warren and Root, 1963). We choose the dual-porosity formulation because
it provides a practical representation of fractured reservoirs with reasonable accuracy and
computational efficiency.
The mass conservation equation for a component i in the flowing domain (fracture) containing np
phases and nc components, where advection, diffusion and gravity are the main physical
mechanisms, can be written as (Jessen et al., 2008)
𝜙 𝜕 𝐶 𝑖 𝜕𝑡
+∇∙𝐹 𝑖 +∇∙𝐻 𝑖 −Γ
i
=𝑞 𝑖 , 𝑖 =1,…,𝑛 𝑐 , (6.1)
while the mass conservation equation for a component i in the stagnant domain (matrix) is written
as
𝜙 𝜕 𝐶 𝑖 𝜕𝑡
+Γ
i
=0, 𝑖 =1,…,𝑛 𝑐 , (6.2)
with the overall molar density Ci is given by
𝐶 𝑖 =∑ 𝑥 𝑖𝑗
𝜌 𝑗 𝑆 𝑗 𝑛 𝑝 𝑗 =1
, 𝑖 =1,…,𝑛 𝑐 . (6.3)
The overall molar advective flux Fi is given by
𝐹 𝑖 =∑ 𝑥 𝑖𝑗
𝜌 𝑗 𝑣 𝑗 𝑛 𝑝 𝑗 =1
, 𝑖 =1,…,𝑛 𝑐 . (6.4)
and the overall molar diffusive flux Hi is given by
106

𝐻 𝑖 =∑ (𝐻 𝑖 ,𝑗 +∑ 𝐻 𝑖 ,𝑗𝑘
𝑛 𝑝 𝑘 =𝑗 +1
)
𝑛 𝑝 𝑗 =1
, 𝑖 =1,…,𝑛 𝑐 . (6.5)
Hi,j denotes the molar diffusive flux of component i within phase j (gas-gas and oil-oil diffusion)
and Hi,jk represents the molar diffusive flux of component i at the interface between phases j and k
(gas-oil diffusion). More details on calculation of molar diffusive fluxes due to oil-oil, gas-gas and
oil-gas diffusion are given in the subsequent section.
In Eqs. (6.1) through (6.4), ϕ is the porosity; xij is the mole fraction of component i in phase j; ρj
and Sj represent the molar density and saturation of phase j, respectively; vj is the velocity vector
of phase j; qi is the source/sink term of component i; and Γi represents the transfer of component i
from matrix to fracture. The phase velocities in Eq. (6.1) are calculated using Darcy’s law, from
the pressure field which is obtained by solving the following volume-balance equation
(𝑉 𝑡 +∆𝑉 𝑡 )−𝑉 𝑐𝑒𝑙𝑙 =0,
1
𝑉 𝑡 ∆𝑉 𝑡 =𝑐 𝑡 ∆𝑝 +∆𝑡 ∑
𝜕 𝑉 𝑡 𝜕 𝑛 𝑖 𝑛 𝑐 𝑖 =1
(∇.𝐹 𝑖 −𝑞 𝑖 −Γ
i
) . (6.6)
where ct is the total fluid compressibility; p is the pressure; Vt is the total fluid volume in a given
control volume (cell); Vcell is the pore volume of the cell; and ni is the number of moles of
component i in the cell. Capillarity and gravity are not considered in this work where we focus on
diffusive mass transfer.
We use a finite volume IMPES (implicit pressure explicit saturation/composition) formulation
with a Cartesian grid in this chapter. The phase-equilibrium calculations are performed using the
Peng-Robinson (PR) equation of state (EOS) and the advective flux is calculated using the standard
single-point-upwind (SPU) scheme. The approach introduced by Michelsen (1982) for phase-split
and stability analysis calculations is utilized in the compositional simulator.
107

6.3. Molecular Diffusion
Molecular diffusion may occur within a single phase or between two different phases. For two
neighboring gridblocks (l and m) with interface area Alm that contain oil and gas phases (Fig. 6.1),
the cross-sectional areas that are available for gas-gas, oil-oil and gas-oil diffusion are
Ag=Alm×min(Sg,l,Sg,m), Ao=Alm×min(So,l,So,m) and Ac=Alm-Ag-Ao respectively.

Figure 6.1 Schematic of two neighboring gridblocks containing oil and gas phases.
In the following subsections we discuss how to calculate molecular diffusion within a phase or
between two different phases based on generalized Fick’s law and film theory.
6.3.1. Gas-Oil Diffusion
First, we consider two different phases which are placed in contact as shown in Fig. 6.2. According
to the film theory, thermodynamic equilibrium will prevail at the interface between the two phases.
In addition, the interface is assumed to offer no resistance to mass transfer. In other words, there
will be no accumulation at the interface, provided there are no interface chemical reactions
(continuity of molar fluxes). In Fig. 6.2, xi,I and yi,I are the interface compositions; and xi,b and yi,b
denote the bulk phase compositions.
108

Figure 6.2 Composition profiles near the interface during inter-phase mass transfer, after Taylor and Krishna
(1993).
Therefore, the gas-oil mass transfer for across the interface is governed by two sets of equations:
Continuity of molar fluxes at the interface
𝑁 𝑖 𝑥 =𝑁 𝑖 𝑦 , 𝑖 =1,…,𝑛 𝑐 , (6.7)
and thermodynamic equilibrium at the interface
𝜇 𝑖 ,𝐼 𝑥 =𝜇 𝑖 ,𝐼 𝑦 , 𝑖 =1,…,𝑛 𝑐 , (6.8)
where, by assuming a positive direction from left to right, 𝑁 𝑖 𝑥 denotes the total molar flux of
component i from the “x” phase into the interface, 𝑁 𝑖 𝑦 denotes the total molar flux of component
i from interface into the “y” phase, and 𝜇 𝑖 ,𝐼 𝑥 (𝜇 𝑖 ,𝐼 𝑦 ) is the chemical potential of component i in the
“x” (“y”) phase at the interface.
The overall molar flux, which consists of advective and diffusive fluxes, is always continuous. It
should be noted that individual continuity of the diffusive and advective fluxes is generally not a
requirement, as a discontinuity of the advective flux can be compensated by a discontinuity in the
diffusive flux; The continuity of the overall molar flux can be formulated as follows
109

𝐽 𝑖 ,𝑚 +𝑥 𝑖 𝑁 𝑡 =𝐽 𝑖 ,𝑓 +𝑦 𝑖 𝑁 𝑡
, 𝑖 =1,…,𝑛 𝑐 −1 . (6.9)
In the above equation, 𝐽 𝑖 ,𝑚 and 𝐽 𝑖 ,𝑓 are diffusive fluxes from matrix to interface, and from fracture
to interface respectively. 𝑁 𝑡 is the total molar flux and 𝑥 𝑖 𝑁 𝑡 and 𝑦 𝑖 𝑁 𝑡 are the advective fluxes of
each component on either side of the interface. Note that the above equation can be written for nc-
1 components. The total molar flux is defined as
𝑁 𝑡 =𝜌 𝑥 𝑢 𝑥 =𝜌 𝑦 𝑢 𝑦 , (6.10)
where 𝑢 𝑥 (𝑢 𝑦 ) is the average velocity of components in the oil (gas) phase. The diffusive flux on
either side of the interface can be represented in weak form as
𝐽 𝑖 ,𝑚 =−𝜌 𝑥 ∑ 𝐷 𝑖𝑘
𝑥 𝑥 𝑘 ,𝐼 −𝑥 𝑘 ,𝑏 𝐿 𝑥 𝑛 𝑐 −1
𝑘 =1
, 𝑖 =1,…,𝑛 𝑐 −1 ,
(6.11)
𝐽 𝑖 ,𝑓 =−𝜌 𝑦 ∑ 𝐷 𝑖𝑘
𝑦 𝑦 𝑘 ,𝑏 −𝑦 𝑘 ,𝐼 𝐿 𝑦 𝑛 𝑐 −1
𝑘 =1
, 𝑖 =1,…,𝑛 𝑐 −1 ,
(6.12)
where Lx (Ly) is the distance between the interface and the center of the gridblock containing the
“x” (“y”) phase. To calculate the molar gas-oil diffusive fluxes, Eqs. (6.8) and (6.9) must be solved
simultaneously. The interface equilibrium calculations (Eq. 6.8) can be done using an EOS at
interpolated temperature and pressure. Multicomponent diffusion coefficients can be calculated
using the following equation, which is obtained by comparing the generalized Fick’s law and the
MS model
𝐷 =𝐵 −1
γ . (6.13)
Here, the matrix B is a function of the inverse of the MS coefficients and γ represents the non-ideal
behavior of the mixture. The elements of the matrices B and γ are given by the following equations
110

𝐵 𝑖𝑘
={
𝑥 𝑖 Ð
𝑖 𝑛 𝑐 +∑
𝑥 𝑙 Ð
𝑖𝑙
𝑛 𝑐 𝑙 =1
𝑙 ≠𝑖 𝑓𝑜𝑟 𝑖 =𝑘 −𝑥 𝑖 (
1
Ð
𝑖𝑘
−
1
Ð
𝑖 𝑛 𝑐 ) 𝑓𝑜𝑟 𝑖 ≠𝑘 (6.14)
and
𝛾 𝑖𝑘
=𝛿 𝑖𝑘
+𝑥 𝑖 𝜕𝑙𝑛 𝜑 𝑖 𝜕 𝑥 𝑘 |
𝑇 ,𝑝 ,𝑥 𝑙 ,𝑙 ≠𝑘 =1…𝑛 𝑐 −1
,
(6.15)
where Đik represent the MS coefficients; δik is the Kronecher delta and φi is the fugacity coefficient
of component i. The MS coefficients matrix is symmetric and its diagonal elements do not exist.
Therefore, there are nc(nc-1)/2 MS coefficients which can be obtained from (Wesselingh and
Krishna, 1990; Kooijman and Taylor, 1991)
Ð
𝑖𝑘
=(Ð
𝑖𝑘
𝑜 )
𝑥 𝑘 (Ð
𝑘𝑖
𝑜 )
𝑥 𝑖 ∏ (Ð
𝑖𝑙
𝑜 Ð
𝑘𝑙
𝑜 )
𝑥 𝑙 /2
𝑛 𝑐 𝑙 =1
𝑙 ≠𝑖 ,𝑘 , 𝑖 ,𝑘 =1,…,𝑛 𝑐 𝑖 ≠𝑘 .
(6.16)
Đoik is the diffusion coefficient of component i when infinitely diluted in component k. We use the
correlation by Leahy-Dios and Firoozabadi (2007) to calculate the infinite dilution diffusion
coefficients for both vapor and liquid phases.
An iterative approach is required to find interface compositions and to calculate the diffusive fluxes
on either side of the interface. We define the vector of independent variables as
V=[ln𝑦 𝐼 ,1
,…,ln𝑦 𝐼 ,𝑛𝑐
,ln𝑥 𝐼 ,1
,…,ln𝑥 𝐼 ,𝑛𝑐
,𝑢 𝑥 ,𝑢 𝑦 ]
𝑇 . (6.17)
There are 2nc+2 independent variables. The system of nonlinear equations can be written as
follows:
nc equilibrium equations at interface
111

𝑓 1,..,𝑛𝑐
(𝑉 )= 𝑙𝑛𝑓𝑢𝑔 𝑦 𝑖 ,𝐼 −𝑙𝑛𝑓𝑢𝑔 𝑥 𝑖 ,𝐼 , 𝑖 =1,…,𝑛 𝑐 , (6.18a)
nc-1 overall flux equations across interface
𝑓 𝑛𝑐 +1,..,2𝑛𝑐 −1
(𝑉 )= 𝐽 𝑖 ,𝑚 +𝑥 𝑖 ,𝐼 𝜌 𝑥 𝑢 𝑥 −𝐽 𝑖 ,𝑓 +𝑦 𝑖 ,𝐼 𝜌 𝑦 𝑢 𝑦 , 𝑖 =1,…,𝑛 𝑐 −1 , (6.18b)
One total molar flux equation
𝑓 2𝑛𝑐
(𝑉 )=𝜌 𝑥 𝑢 𝑥 −𝜌 𝑦 𝑢 𝑦 , (6.18c)
and two constraints on mole fractions at the interface (sum to unity)
𝑓 2𝑛𝑐 +1
(𝑉 )=∑𝑥 𝑖 ,𝐼 𝑛𝑐
𝑖 =1
−1 (6.18d)
𝑓 2𝑛𝑐 +2
(𝑉 )=∑𝑦 𝑖 ,𝐼 𝑛𝑐
𝑖 =1
−1 (6.18e)
The above system of equations can be solved using a Newton’s approach.
6.3.2. Gas-Gas and Oil-Oil Diffusion
Diffusion between pairs of the same fluid phase across the matrix-fracture interface can be
expressed as follows
𝐽 𝑖 ,𝑙 =−𝜌 𝑥 ∑ 𝐷 𝑖𝑘
𝑥 𝑥 𝑘 ,𝑚 −𝑥 𝑘 ,𝑏 𝐿 𝑛 𝑐 −1
𝑘 =1
, 𝑖 =1,…,𝑛 𝑐 −1 ,
(6.19a)
𝐽 𝑖 ,𝑔 =−𝜌 𝑦 ∑ 𝐷 𝑖𝑘
𝑦 𝑦 𝑘 ,𝑚 −𝑦 𝑘 ,𝑏 𝐿 𝑛 𝑐 −1
𝑘 =1
, 𝑖 =1,…,𝑛 𝑐 −1 ,
(6.19b)
where Dik,j are the multicomponent diffusion coefficients in phase j. L is the distance between the
center of the matrix block and the center of the fracture block. In the above formulation, the sum
112

of diffusive fluxes of all components in phase j must be zero (∑ 𝐽 𝑖 =0
𝑛𝑐
𝑖 =1
); and the diffusive flux
for the last component (Jnc) is obtained accordingly (Taylor and Krishna, 1993).
6.4. Transfer Functions
Since the introduction of dual-porosity models in the 1960s, various transfer functions with
different levels of sophistication have been suggested for M-F interactions. Warren and Root
(1963) proposed a pseudo-steady state (PSS) transfer function for single-phase flow between
matrix and fracture. The application of a shape factor in reservoir simulation was first introduced
by Kazemi et al. (1976) for three-phase flow subject to a PSS assumption. Coats (1989) extended
the dual-porosity model for compositional simulation and included molecular diffusion in the
transfer function.
The transfer rate of a component i from matrix to fracture is given by
Γ
i
=Γ
i,e
+Γ
i,d
+Γ
i,gd
, 𝑖 =1,…,𝑛 𝑐 , (6.20)
where Γi,e, Γi,d and Γi,gd are transfer rates due to expansion (e), diffusion (d) and gravity drainage
(gd), respectively. Since matrix permeability is low, we do not consider the gravity drainage term
in our analysis. The expansion term is given by
Γ
i,e
=𝜎 1
𝐾 𝑚 ∑ 𝜆 𝑗𝑚
𝑥 𝑖𝑗 ,𝑚 𝜌 𝑗𝑚
(𝑝 𝑚 −𝑝 𝑓 )
𝑛𝑝
𝑗 =1
, 𝑖 =1,…,𝑛 𝑐 ,
(6.21)
where σ1 is the shape factor, Km is the (equivalent isotropic) matrix permeability, λjm is the mobility
of phase j in the matrix, xij,m is the mole fraction of component i in phase j in the matrix, ρjm is the
molar density of phase j in the matrix, and pm and pf are the pressure in matrix and fracture,
respectively. The shape factor proposed by Lim and Aziz (1995) for systems with anisotropic
rectangular matrix blocks is used in our model.
113

𝜎 1
=
𝜋 2
𝐾 𝑚 (
𝐾 𝑚 ,𝑥 𝐿 𝑥 2
+
𝐾 𝑚 ,𝑦 𝐿 𝑦 2
+
𝐾 𝑚 ,𝑧 𝐿 𝑧 2
) , (6.22)
where Km,x, Km,y and Km,z are the matrix permeabilities in x-, y-, and z-direction, respectively, and
Lx, Ly, Lz are the dimensions of matrix block in x-, y-, and z-direction, respectively. The equivalent
isotropic matrix permeability is defined as
𝐾 𝑚 =(𝐾 𝑥 𝐾 𝑦 𝐾 𝑧 )
1/3
.
(6.23)
To calculate the transfer rate of each component from matrix to fracture due to molecular diffusion,
we sum over the diffusion rates of that component. For a two-phase gas-oil system, the transfer
rate of a component i is calculated from
Γ
i,d
=Γ
i,d,g
+Γ
i,d,o
+Γ
i,d,c
, 𝑖 =1,…,𝑛 𝑐 , (6.24)
where Γi,d,g, Γi,d,o and Γi,d,c are transfer rates due to gas-gas, oil-oil and gas-oil molecular diffusion,
respectively. The transfer rate of a component i due to gas-gas diffusion is given by
Γ
i,d,g
=𝜎 1
𝜙 𝑚 𝜌 𝑔 𝑆 𝑔 ∑ 𝐷 𝑖𝑘 ,𝑔 𝑛 𝑐 −1
𝑘 =1
(𝑦 𝑘𝑚
−𝑦 𝑘𝑓
), 𝑖 =1,…,𝑛 𝑐 −1 ,
(6.25)
where ϕm is the matrix porosity and Sg=min(Sgm,Sgf) is the fraction of the total M-F surface area
that is available for gas-gas diffusion. The transfer rate of a component i due to oil-oil diffusion is
calculated similarly to gas-gas diffusion. The transfer rate of each component from matrix to
fracture due to gas-oil diffusion is calculated based on film theory, as mentioned in section 6.3.1.
After calculating the diffusive fluxes on both sides of the interface, we use an upwinding scheme
for calculation of the transfer rates. The transfer rate of a component due to gas-oil diffusion can
be either of the following, based on the upwind direction.
114

Γ
i,d,c,m
=−𝜎 1
𝜙 𝑚 𝜌 𝑚 𝑆 𝑔𝑜
∑ 𝐷 𝑖𝑘
𝑥 (𝑥 𝑘 ,𝐼 −𝑥 𝑘 ,𝑏 )
𝑛 𝑐 −1
𝑘 =1
, 𝑖 =1,…,𝑛 𝑐 −1 ,
(6.26a)
Γ
i,d,c,f
=−𝜎 2
𝜌 𝑓 𝑆 𝑔𝑜
∑ 𝐷 𝑖𝑘
𝑦 (𝑦 𝑘 ,𝑏 −𝑦 𝑘 ,𝐼 )
𝑛 𝑐 −1
𝑘 =1
, 𝑖 =1,…,𝑛 𝑐 −1 ,
(6.26b)
where Sgo=1- min(Sgm,Sgf) - min(Som,Sof) is the fraction of the total M-F surface area that is available
for gas-oil diffusion; σ1 is obtained from an isotropic version of Eq. (6.22); and σ2 is given by
𝜎 2
=𝜋 2
(
1
𝐿 𝑥 𝑏 𝑥 +
1
𝐿 𝑦 𝑏 𝑦 +
1
𝐿 𝑧 𝑏 𝑧 ) , (6.27)
where bx, by, and bz are the fracture openings in x-, y-, and z-direction, respectively. In order to
ensure the numerical stability of the system of equations presented in Eq. (6.18a-e), we have
assumed that 𝜎 2
=10𝜎 1
.
6.5. Calculation Examples
Two numerical examples are presented in this section to demonstrate the significance of accurate
modeling of molecular diffusion and M-F interactions during gas injection in fractured reservoirs.
In both cases, the size of reservoir in x-, y-, and z-direction is 250 m, 250 m, and 10 m respectively.
25×25×5 gridblocks are used for both flowing and stagnant domains. The flowing domain has an
effective permeability of 500 md in all three directions and a porosity of 1%. In the stagnant
domain, the permeabilities in horizontal directions (Kx and Ky) are 1 md while the vertical
permeability is 0.001 md. The stagnant domain has a porosity of 20%. The low vertical
permeability of the matrix blocks was selected to prevent a significant contribution from gravity
drainage and to focus on the representation of diffusive mass transfer. The fracture spacing and
opening in both x- and y-direction are 10 m and 1 mm respectively.
115

In both examples, one injection and one production well are located at opposite corners of the
reservoir. The injection well is completed at the top of formation while the production well is
completed at the bottom of formation. The bottomhole pressure is kept constant for the production
well while the injection well is operated at constant rate.
Linear relative permeabilities are used for the flowing domain (fractures), while Corey-type
relative permeability curves, with parameters given in Table 6.1, are used for the stagnant domain
(matrix).
Table 6.1 Parameters for Corey-type relative permeabilities
S
gc
S
or
n
g
n
o
k
rg,e
k
ro,e

Flowing Domain 0.00 0.00 1.0 1.0 1.0 1.0
Stagnant Domain 0.00 0.25 3.0 4.0 0.75 1.0

Table 6.2 provides a description of the fluid components used in this study, listing their critical
properties, acentric factors, molecular weights, volume shift parameters and parachors (used for
evaluation of the interfacial tension).
Table 6.2 Fluid properties
Comp Tc [K] Pc [atm]
Acentric
factor
MW
[g/mole]
Volume
shift
CO2 304.21 72.86 0.2236 44.01 -0.0817
CH4 190.56 45.39 0.0115 16.043 -0.154
nC4 425.12 37.46 0.202 58.123 -0.06413
116

nC12 658 17.96 0.5764 170.338 0.1150

In the following examples, we compare the conventional approach for representing diffusive mass
transfer to the proposed approach discussed above. The key differences between the two
approaches are summarized below:
Conventional approach:
Molecular diffusion is modeled using the classical Fick’s law in which off-diagonal diffusion
coefficients are neglected. Moreover, constant diffusion coefficients are used and hence the model
ignores the dependency of diffusion coefficients on pressure, temperature and composition. The
diagonal elements of the diffusion coefficient matrix at initial (oil) and injection (gas) conditions
are used for the entire simulation for oil and gas phases, respectively.
The gas-oil diffusion is modeled by assuming an average composition at the interface. By
assuming an average interface composition in Eq. (6.9), one can obtain expressions similar to Eqs.
(6.26a) and (26.26b) for gas-oil diffusion.
Proposed approach:
Molecular diffusion is modeled using the generalized Fick’s law. Diffusion coefficients (full
matrix) are obtained from the MS model (Eq. 6.13) and are pressure-, temperature- and
composition-dependent.
The gas-oil diffusion is modeled based on film theory. According to film theory, when gas and oil
phases are placed in contact, a two-phase interface is formed between them in which gas and oil
117

are in thermodynamic equilibrium. In addition, the molar flux of each component across the
interface must be continuous.
6.5.1. Example 1
In the first example, the initial reservoir pressure at the top of the formation is 157 bar. The
bottomhole pressure for the production well is fixed at 157 bar and gas is injected at rate of 165
Rm3/day. The injected gas is composed of 90% CO2 and 10% CH4 by mole. The reservoir is
initially saturated with oil (30% CH4, 40% nC4H10 and 30% nC12H26 by mole). The reservoir
temperature is 373.15 K at which the oil has a bubble point pressure of 91.6 bar and the minimum
miscibility pressure (MMP) for the oil/gas pair is 150.5 bar. Accordingly, a multicontact miscible
displacement may develop during the displacement process (see more discussion later).
Fig. 6.3 compares the simulation results for example 1 in terms of oil recovery and producing gas-
oil ratio (GOR) for the cases when diffusion is not considered (red), when diffusion is modeled
using our approach (blue), and when diffusion is modeled using the conventional approach (black).
The left tab in Fig. 6.3 shows the liquid recovery (%) while the right tab shows the producing GOR
versus time for ten years of injection. There are minor oscillations in the gas-oil ratio in early times.
Our IMPES simulator uses a volume balance approach (Watts 1986, Aziz and Wong 1988-1989).
As a result, a slight mismatch between the cell volume and the calculated fluid volume, cause a
volume error in the simulator. The volume error results in minor fluctuations in phase velocities
which in turn cause the GOR to oscillate.
In the absence of diffusion, the injected gas mainly sweeps the oil residing in the flowing domain
(fractures), which in turn leads to a very low recovery and a very high producing GOR. However,
when molecular diffusion is included, diffusion of gas components from fractures into the matrix
118

extracts oil components from matrix (swelling and diffusion) and leads to a substantially higher
recovery and a lower producing GOR. Fig. 6.3 also shows that the simulation results obtained from
the conventional approach are significantly different from the results based on the proposed
approach: Significantly higher recoveries and lower gas-oil ratios are obtained using the
conventional approach. We defer further analysis of the differences between predictions by these
two approaches until the discussion section.

Figure 6.3 Oil recovery (left) and producing GOR (right) versus time for example 1 when diffusion is not
considered (red), when diffusion is modeled using the conventional approach (black) and when diffusion is modeled
using our approach (blue).

119

Figure 6.4 Oil recovery (left) and producing GOR (right) versus time for example 1 when diffusion is not
considered (solid red), when our approach is used (solid blue) and when our approach is used and the off-diagonal
diffusion coefficients are neglected (dashed blue).

Next, we investigate how much the results will change for the proposed approach if the off-
diagonal diffusion coefficients are neglected (i.e. the classical Fick’s law). Fig. 6.4 compares the
simulation results for example 1 when diffusion is not considered (solid red), when our approach
is used (solid blue) and when our approach is used and the off-diagonal diffusion coefficients are
neglected (dashed blue).
From Fig. 6.4 we observe that neglecting the off-diagonal diffusion coefficients introduce a ~10%
difference (2.8% out of 26%) in the predicted recovery over the time period investigated. That is
a notable difference, but a much smaller one than the ~ 38% difference (10% out of 26%) observed
between the cumulative recoveries calculated using our approach (which uses film theory for gas-
oil diffusion) and the conventional approach (which assumes an average composition at gas-oil
interface) (see Fig. 6.3). Therefore, a key component in accurate modeling of gas injection in
fractured reservoirs, when molecular diffusion is the main recovery mechanism, is using an
appropriate physical model for gas-oil diffusion (i.e. film theory), and not the contribution from
off-diagonal diffusion coefficients.
From Fig. 6.4 we also observed that when off-diagonal diffusion coefficients are taken into
account, an increase in recovery is predicted due to the contributions from positive (in this case)
off-diagonal coefficients.
The CO2 composition (mole %) in the stagnant domain (matrix) is reported in Fig. 6.5 after ten
years of injection for the case with diffusion based on calculations with the proposed approach
120

(left) and the difference (absolute values) between the proposed approach and the conventional
approach (right). We observed that the differences in the amount of CO2 that entered the stagnant
domain between the two approaches are considerable. It should be noted that the black lines in
Fig. 6.5 represent the boundaries of simulation cells, not the actual fracture network.

Figure 6.5 The CO
2
composition (mole %) in the stagnant domain after 10 years from calculations with diffusion
based on the proposed approach (left) and absolute difference (mole %) between the proposed approach and the
conventional approach (right).
6.5.2. Example 2
In the second example, the initial reservoir pressure is set at 150 bar (top of formation), the
bottomhole pressure of the production well is kept constant at 150 bar and gas is injected at rate
of 180 Rm3/day. The injected gas is composed of 50% CO2 and 50% CH4 by mole. The reservoir
is initially saturated with oil (30% CH4, 40% nC4H10 and 30% nC12H26 by mole) at a temperature
of 336.15K (bubble point of oil = 78.6 bar). The MMP of the oil/gas pair is 179 bar and
accordingly, this displacement cannot develop multicontact miscibility.
Figure 6.6 compares the calculation results for example 2 when diffusion is not considered (red),
when diffusion is modeled the proposed approach (blue), and when diffusion is modeled using the
Producer
Injector
Producer
Injector
121

conventional approach (black). The left tab shows the oil recovery (%) while the right tab presents
producing GOR versus time for ten years of injection.
Similar to the previous example, a low recovery and a high producing GOR is observed in the
absence of molecular diffusion. However, when molecular diffusion is included, substantially
higher oil recoveries and lower producing GOR’s are obtained. Again, we observe that a
significantly higher liquid recovery and a lower GOR are predicted by the conventional approach
relative to the proposed approach.

Figure 6.6 Oil recovery (left) and producing GOR (right) versus time for example 2 when diffusion is not
considered (red), when diffusion is modeled using the conventional approach (black) and when diffusion is modeled
using the proposed approach (blue).
The impact of the off-diagonal diffusion coefficients on the calculation results for example 2 are
demonstrated in Fig. 6.7. The same color codes and line styles of Fig. 6.4 are used to represent
each physical model in Fig. 6.7.
122

Figure 6.7 Oil recovery (left) and producing GOR (right) versus time for example 2 when diffusion is not
considered (solid red), when our approach is used (solid blue) and when our approach is used but the off-diagonal
diffusion coefficients are neglected (dashed blue).
Similar to the previous calculation example, neglecting the off-diagonal diffusion coefficients
introduce notable differences among the results relative to the proposed approach. However, the
differences are, again, much smaller than what we observe between the proposed approach and the
conventional approach (see Fig. 6.6): Lower liquid recoveries are obtained when off-diagonal
diffusion coefficients are taken into account due to the contributions of the off-diagonal
coefficients.
Compared to example 1, we observe a significant reduction in the oil recovery (~16% relative to
~26%) when the injection gas composition is changed and the BHP is reduced (157 bar to 150
bar). The change in the injection gas composition is the main reason for the observed reduction
in recovery. This is primarily due to the higher solubility of the CO2-rich stream (example 1) that
promotes swelling of the oil. The decrease in recovery observed between example 1 and 2 is
hence directly related to the decrease in solubility of the injectant (lower solubility and less
swelling of the oil in example 2).
123

6.6. Discussion
To further understand the dynamics of oil recovery by molecular diffusion during gas injection in
fractured reservoirs, we look at the saturation and composition paths/profiles of a matrix block and
its corresponding fracture over time in Fig. 6.8.

Figure 6.8 Saturation variation in a matrix block (dotted line) and its corresponding fracture (solid line) near the
injection well for example 1 as calculated from the proposed approach (top) and the conventional approach
(bottom).
The top panel of Fig. 6.8 shows the saturation profile observed during the first ten years of injection
in a matrix block (dotted line) near the injection well and its corresponding fracture (solid line) for
example 1 as calculated using our approach. The fracture becomes fully saturated with gas in a
short period of time, while the gas saturation in the matrix remains zero for almost two years, after
which it gradually increases.
124

The bottom panel of Fig. 6.8 shows the saturation profile observed during the first ten years of
production at the same matrix block (dotted line) and its corresponding fracture (solid line) as
calculated from the conventional approach. Again, the fracture becomes fully saturated with gas
in a short period of time. However, the gas saturation in the matrix remains zero for almost seven
years, after which it gradually increases. We note that a gaseous phase was formed in the matrix
after ~ two years when the proposed approach was used in the calculation.
To further investigate this difference between the conventional approach and the proposed
approach, we consider the changes in phase compositions with time that are generated from the
different modeling approaches (Figure 6.9).
The four panels of Fig. 6.9 report the compositions (mass fraction) versus time for each of the four
components of a matrix block near the injection well for example 1. In each panel, the solid line
and the dashed line represent observed compositions as calculated from the conventional approach
and the proposed approach, respectively.
Fig. 6.9 demonstrates why a gaseous phase does not form in the matrix until after ~ seven years
for calculations based on the conventional approach. During the first seven years, the CO2
composition largely increases with small changes in the (relative) compositions of other
components. This acts to keep the fluid inside the matrix in a single-phase state for nearly seven
years. The associated swelling of the oil results in larger liquid recoveries relative to the proposed
approach.
125

Figure 6.9 Composition paths (mass fraction) of a matrix block near the injection well for example 1 during the first
10 years as calculated using the conventional approach (solid line) and the proposed approach (dashed line).
The composition profile in the matrix, as calculated by the conventional approach, differs
significantly from what is calculated using the proposed approach. This difference between
calculated composition paths may introduce significant differences in the predicted recovery
dynamics when gravity drainage contributes to the displacement process. In this work, we have
neglected capillary forces that depend on the interfacial tension. As the interfacial tension depends
on the two-phase compositions (and reservoir pressure and temperature) that prevail in a given
126

system, we expect that the different mass transfer models may result in different gravity drainage
dynamics as well. Additional work is currently in progress to explore the potential impact of mass
transfer modeling on the gravity drainage process.
6.7. A Different Approach to Find Interface Compositions
The formulations presented in this chapter for calculating the interface composition (Eq. 6.17,
6.18a-e) are more rigorous, but computationally expensive. The calculations based on the proposed
model requires significantly more CPU time than the conventional model used in commercial
reservoir simulators, since an iterative approach is used to find the interface compositions. One
way to overcome this concern is to find an approximation of the interface composition that, while
preserving the accuracy of calculations, eliminates the need for iterative computations.
One such approximation can be introduced by assuming that the mass transfer limitation on the
fracture side is marginal. This is a reasonable assumption since the gas phase in fractures is flowing
at relatively high velocities. Based on this assumption, the interface composition would be a
mixture of the gas phase in the fracture and a liquid phase in the matrix, that will correspond to an
equilibrium state. We can define a mixing ratio (f) which represents the fraction of the liquid phase
in the interface mixture. The overall composition at the interface is then given by
𝑧 i
=f𝑥 𝑖 +(1−f)𝑦 𝑖 , (6.28)
The mixture at interface is then flashed to obtain the liquid and vapor composition at the interface.
The mixing ratio (f) may need to be adjusted to ensure that the overall mixture at interface
corresponds to a two-phase system at equilibrium. The fracture aperture is much smaller than the
matrix block length. As a result, the composition at the interface is more similar to the fracture
fluid composition than the matrix fluid composition. For practical purposes, f may be set to a value
127

of 0.975. If this mixture of bulk phase fluids doesn’t reach a two-phase equilibrated mixture, the
value of f must be adjusted (shifted to lower gas fractions).
In order to verify the validity of this approximation, we repeated the calculations for example 1
based on Generalized Fick’s law using the interface flash approach (presented in this section).
Figure 6.10 displays the calculation results when diffusion is not accounted for in calculations (red
line), when diffusion is accounted for, and the interface composition is obtained using the approach
suggested in section 6.3.1 (blue line), and when diffusion in accounted for, and the interface
composition is obtained using a flash of bulk fluid mixtures (dashed blue line).

Figure 6.10 Oil recovery (left) and producing GOR (right) versus time for example 1 when diffusion is not
considered (solid red), when our approach presented in sec. 6.3.1 is used (solid blue) and when our approach
presented in sec. 6.7 is used (dashed blue)
The calculation results presented in Figure 6.10 show that the differences between the iterative
approach and the interface flash approach are negligible (~2% difference in recovery out of ~26%).
As a result, it is suggested to use the simplified approximation to find the interface composition
when the available resources for simulation pose a time limitation on calculations.
128

6.8. Conclusions
Base on the examples and discussion presented in the previous sections, we arrive at the following
conclusions:
1. Molecular diffusion can substantially enhance the oil recovery by gas injection in fractured
reservoirs. Diffusion of gas components from fractures into the matrix extracts oil components
from matrix (via diffusion and associated swelling) and leads to substantially higher liquid
recoveries and lower producing GOR’s.
2. The proposed approach for modeling mass transfer, which is based on a rigorous physical
model for gas-gas, oil-oil and gas-oil diffusion (film theory), provides significantly different
results in terms of oil recovery and GOR relative to the conventional approach, which assumes
an average composition at the gas-oil interface.
3. The dragging effect (impact of off-diagonal diffusion coefficients) has a moderate impact on
the oil recovery dynamics during gas injection in fractured reservoirs.
4. The composition path that is dictated by the mass transfer model controls the variation in IFT
during gas injection in fractures systems and may impact the dynamics of a gravity drainage
process.
5. A satisfactory approximation of the interface composition can be obtained by using a mixture
of fracture and matrix fluids at the interface.

129

Chapter 7 Coarse Scale Modeling of Multicomponent
Diffusive Mass Transfer in Dual-porosity Models
7.1. Introduction
Naturally fractured reservoirs contain a substantial amount of hydrocarbons worldwide. However,
recovery from such reservoirs has been rather low. The low oil recovery from such reservoirs
highlight the need for accurate characterization, modeling and simulation of naturally fractured
reservoirs to determine the flow patterns within these reservoirs and the location of bypassed oil.
The main mechanisms of recovery in naturally fractured reservoirs include molecular diffusion,
viscous displacement and gravity drainage. The relative significance of these mechanisms depends
on several factors including matrix permeability, fracture intensity, fluid properties, injection rate
and reservoir pressure and temperature. (Shojaei and Jessen, 2015). Molecular diffusion can play
an important role in fractured reservoirs because of the large surface area available for diffusion.
During gas injection, diffusion of gas components from a fracture into the matrix extracts oil
components from matrix and delays, to some extent, the gas breakthrough via high-permeability
fractures. As a result, both displacement and sweep efficiencies are increased. Different
investigators have demonstrated the role of molecular diffusion in fractured reservoirs (e.g. da
Silva and Belery, 1989; Hu and Whitson, 1991; Darvish et al., 2006; Hoteit and Firoozabadi, 2009;
Vega et al., 2010; Jamili et al., 2011).
Dual-porosity models are widely used in simulation of naturally-fractured reservoirs. These
models assume the existence of two overlapping and continuous media; The fractures (flowing
domain) provide the flow path and the matrix (stagnant domain) provides the storage region for
the fluid. In dual-porosity models, it is also assumed that no viscous flow occurs between
130

neighboring matrix blocks. The flow between the flowing domain and the stagnant domain is
defined by appropriate transfer rates. The main advantage of dual-porosity models is the
computational efficiency compared to discrete fracture models, while providing acceptable
accuracy in modelling results.
The concept of dual-porosity models was first introduced by Barenblatt et al (1960). It was later
extended into petroleum engineering applications by Warren and Root (1963), primarily for
pressure transient analysis. Warren and Root formulated the transfer function between the flowing
and stagnant domain based on the assumption of pseudo-steady state flow between the two
domains. Kazemi et al. (1976) extended the Warren and Root model to multiphase flow by
assuming a Darcy-like flux between fracture and matrix. the application of a shape factor in
reservoir simulation was also first introduced by Kazemi et al. (1976) for multiphase flow subject
to a pseudo-steady state assumption.
Since the rise in the popularity of dual-porosity models, many researchers have aimed their work
on improving the accuracy of such models by comparing their results to fine scale simulation, and
deriving relevant formulations for transfer rates and shape factors for a range of fluid compositions,
matrix geometries and flow settings. Sonier et al. (1988) and Litvak (1985) provided a dynamic
approach to improve the modeling of the interaction of gravity and capillary forces in
matrix/fracture systems, while Gilman and Kazemi (1988) presented a formulation to account for
the viscous displacement in the matrix blocks in dual-porosity models. Saidi (1983) developed a
fully implicit reservoir simulator with discretization in matrix blocks and was able to match
production data from Haft Kel field in Iran. Kazemi and Gilman (1993) performed an analysis of
the gravity drainage effects on the transfer rate between matrix and fracture in dual porosity
131

simulation and concluded that the gravity terms should be a function of both the height of the
matrix block and the contrast between fluid dynamic forces between the fracture and matrix.
By solving the single-phase pressure equation and approximating the exact solution with
simplified but reasonably accurate functions, Zimmerman et al. (1993) and Lim and Aziz (1995)
derived matrix-fracture transfer functions without assuming a PSS condition. Based on the work
of Zimmerman et al. (1993), Lu et al. (2008) proposed a general transfer function that accounts for
both early- and late-time behavior of multiphase matrix-fracture interactions. The general transfer
function formulation contains a correction factor that is relatively high in early-time and
approaches unity at later times, when a pseudo-steady state assumption holds. The formulation for
the general transfer function is developed for pressure diffusion, and adapted for molecular
diffusion on the grounds of similarity.
While the approach suggested by Lu et al. (2008) works well for a single component system, it
does not provide a basis for adoption of the same formulation in multicomponent settings.
Diffusive fluxes of components across the interface are highly coupled, and an extension of the
single-phase pressure equation to these settings is non-trivial. In fact, despite the rich literature on
the subject of transfer functions in dual-porosity models, there are no works in the literature, to the
best of our knowledge, that have focused on generalized transfer functions for multicomponent
molecular diffusion. In this chapter, we aim to present time-dependent transfer functions for
multicomponent molecular diffusion that accurately represent both early-time and late-time mass
transfer. The proposed transfer functions are obtained using a semi-analytical approach in solving
the partial differential equations that govern matrix-fracture transfer.
132

7.2. Methodology
In this section, we discuss the mass balance equations for both the stagnant domain and the flowing
domain. We will characterize the transfer function between the two domains and provide
estimations to the analytical solution, and define a “correction factor” to the Warren and Root
formulation for the transfer function. The approach used in this work is adopted from previous
work by Zimmerman et al. (1993).
The mass conservation equation for a component i in the flowing domain (fracture) containing np
phases and nc components can be written as (Jessen et al., 2008)
𝜙 𝜕 𝐶 𝑖 (𝑥 𝑓 ,𝑡 )
𝜕𝑡
+∇∙𝐹 𝑖 +∇∙𝐻 𝑖 −Γ
i
=𝑞 𝑖 , 𝑖 =1,…,𝑛 𝑐 ,
(7.1)
while the mass conservation equation for a component i in the stagnant domain can be written as
𝜙 𝜕 𝐶 𝑖 (𝑥 𝑚 ,𝑡 ;𝑥 𝑓 )
𝜕𝑡
+Γ
i
=0, 𝑖 =1,…,𝑛 𝑐 .
(7.2)
With the overall molar density Ci given by
𝐶 𝑖 =∑𝑥 𝑖𝑗
𝜌 𝑗 𝑆 𝑗 𝑛 𝑝 𝑗 =1
, 𝑖 =1,…,𝑛 𝑐 . (7.3)
In the above equations, ϕ is the porosity; xij is the mole fraction of component i in phase j; ρj and
Sj represent the molar density and saturation of phase j, respectively; vj is the velocity vector of
phase j; qi is the source term of component i; and Γi represents the transfer of component i from
matrix to fracture. The transfer rate of a component i from matrix to fracture is given by
Γ
i
=Γ
i,e
+Γ
i,d
+Γ
i,gd
, 𝑖 =1,…,𝑛 𝑐 , (7.4)
133

where Γi,e, Γi,d and Γi,gd are transfer rates due to expansion (e), diffusion (d) and gravity drainage
(gd) respectively. Transfer rates due to each of these mechanisms is treated separately in
compositional simulation. We will aim our focus on mass transfer due to molecular diffusion.
Using a Fickian approach, the mass balance in the matrix block can be described as
𝜕𝐶 (𝑥 𝑚 ,𝑡 ;𝑥 𝑓 )
𝜕𝑡
=D∇
2
𝐶 (𝑥 𝑚 ,𝑡 ;𝑥 𝑓 ) .
(7.5)
In the Warren and Root model, no attempt is made to solve the diffusion equation within each
block, and the distribution within the block is governed by an ordinary, rather than partial
differential equation. Zimmerman et al. (1993) explain the derivation process for the Warren and
Root model in detail. Assuming an average concentration within the matrix block, their model
results in the following ordinary differential equation
𝜕 𝐶 ̅ 𝑚 (𝑡 ;𝑥 𝑓 )
𝜕𝑡
=𝜎𝐷 (𝐶 𝑓 −𝐶 ̅ 𝑚 ) .
(7.6)
The mass transfer rate due to diffusion is then defined as
Γ
d
=𝜑𝜎𝐷 (𝐶 𝑓 −𝐶 ̅ 𝑚 ) . (7.7)
In the above equation, 𝜎 is the “shape factor” that depends on block shape and geometry, and has
dimensions of area
-1
. This representation of the transfer rate is commonly known as the “quasi-
steady state” approximation. In the next two sections, we will discuss the exact solution to mass
transfer within each matrix block for a binary system followed by a more general multicomponent
fluid system. We will provide semi-analytical approximations of the exact solution in each of these
settings, and propose a correction factor to the Warren and Root model for calculation of the
transfer rate in a dual-porosity setting.
7.2.1. Binary Fluid System
134

In a binary fluid component system, the diffusion calculations only need to be performed for one
of the species. Since the sum of the diffusive fluxes within the system is zero, the diffusive flux of
the second component can be calculated accordingly. Since there are only two components in the
system, we can rewrite Eq. 7.5 in terms of molar concentration of the first component in the system
as follows
𝜕𝑐 (𝑥 𝑚 ,𝑡 ;𝑥 𝑓 )
𝜕𝑡
=D∇
2
𝑐 (𝑥 𝑚 ,𝑡 ;𝑥 𝑓 ) .
(7. 8)
As mentioned in the previous section, the Warren and Root formulation assumes pseudo-steady
state conditions within the matrix block. This assumption is valid for late-time behavior, but does
not capture the early- and intermediate-time behavior. The goal of this section is to arrive at an
approximation of the analytical solution that is valid for all times. For most of the following
sections, we assume that the matrix block is a sphere of radius am. We will discuss the extension
of the proposed formulation to other geometries later in this section. Under the assumption of a
spherical matrix block, the exact solution to the density distribution within the matrix block can
be found by solving Eq. 7.8, subject to the following initial and boundary conditions
𝑐 (𝑥 𝑚 ,𝑡 =0)=𝑐 𝑖𝑛𝑖𝑡 , (7. 9)
𝑐 (|𝑥 𝑚 |=𝑎 𝑚 ,𝑡 )=𝑐 𝑓 , (7. 10)
where cinit is the initial concentration of the first component within the matrix block and cf is the
concentration at the block’s outer boundary. The analytical solution to the above set of equations
can be found in many textbooks (e.g. Crank, 1975) and can be written as follows
𝑐 ̅−𝑐 𝑖𝑛𝑖𝑡 𝑐 𝑓 −𝑐 𝑖𝑛𝑖𝑡 =1−
6
𝜋 2
∑
1
𝑛 2
exp(−
𝐷𝑛
2
𝜋 2
𝑡 𝑎 𝑚 2
)
∞
𝑛 =1
. (7.11)
135

For sufficiently large times, all terms in the above summation beyond the first term are negligible.
Dropping all the other terms, integrating and eliminating t results in the Warren and Root
approximation of the transfer rate.
𝜕 𝑐 ̅ 𝜕𝑡
=
𝜋 2
𝐷 𝑎 𝑚 2
(𝑐 𝑓 −𝑐 ̅) . (7.12)
Comparing Eq. 7.12 and 7.6 shows that the shape factor in Warren and Root’s formulation for a
spherical matrix block is 𝜎 =
𝜋 2
𝑎 𝑚 2
.
As mentioned previously, this approximation is only valid for late-time behavior and does not
mimic the early time behavior within the matrix block. A better approximation of the analytical
solution over all time scales was suggested by Vermeulen (1953) and presented by Zimmerman et
al. (1993) as follows
𝑐 ̅−𝑐 𝑖𝑛𝑖𝑡 𝑐 𝑓 −𝑐 𝑖𝑛𝑖𝑡 =[1− exp(−
𝐷 𝜋 2
𝑡 𝑎 𝑚 2
)]
1/2
.
(7.13)
Differentiating Eq. 7.13 with respect to time and eliminating t from the result leads to
𝜕 𝑐 ̅ 𝜕𝑡
=
𝜋 2
𝐷 2𝑎 𝑚 2
(𝑐 𝑓 −𝑐 𝑖𝑛𝑖𝑡 )
2
−(𝑐 ̅−𝑐 𝑖𝑛𝑖𝑡 )
2
(𝑐 ̅−𝑐 𝑖𝑛𝑖𝑡 )
. (7.14)
Rearranging the above equation results in
𝜕 𝑐 ̅ 𝜕𝑡
=
𝜋 2
𝐷 𝑎 𝑚 2
(𝑐 𝑓 −𝑐 ̅)
(|𝑐 𝑓 −𝑐 𝑖𝑛𝑖𝑡 |)+(|𝑐 ̅−𝑐 𝑖𝑛𝑖𝑡 |)
2(|𝑐 ̅−𝑐 𝑖𝑛𝑖𝑡 |)
.
(7.15)
By comparing the above formulation with the Warren and Root formulation (Eq. 7.6), we can see
that the approximation to the analytical solution has resulted in a “correction factor” based on
initial, fracture, and matrix concentrations, to the Warren and Root formulation. This correction
factor is close to unity in late time, when the concentration within the matrix reaches that of the
136

surrounding fracture the transfer function then reverts to the traditional Warren and Root model.
In the early time, however, the concentration within the matrix is close to its initial value and the
correction factor is large. We can write Eq. 7.15, including the correction factor (B) as
𝜕 𝑐 ̅ 𝜕𝑡
=
𝜋 2
𝐷 𝑎 𝑚 2
𝐵 (𝑐 𝑓 −𝑐 ̅) , (7.16)
𝐵 =
(|𝑐 𝑓 −𝑐 𝑖𝑛𝑖𝑡 |)+(|𝑐 ̅−𝑐 𝑖𝑛𝑖𝑡 |)
2𝑀𝑎𝑥 (𝘀 𝑐 ̅,|𝑐 ̅−𝑐 𝑖𝑛𝑖𝑡 |)
,
(7.17)
and the mass transfer rate between matrix and fracture due to diffusion as
Γ
d
=𝜑𝜎𝐵𝐷 (𝑐 𝑓 −𝑐 ̅) , (7.18)
where 𝘀 is a small parameter (0.01-0.1) that is introduced to ensure the stability of numerical codes
at very early times, otherwise the correction factor diverges at initial conditions, causing
convergence problems (Lu et al., 2008).
7.2.2. General Formulation for Multicomponent Fluid System
In this section we will expand the above formulation for the binary case to a more general
multicomponent formulation. In a multicomponent fluid system, the mass conservation within the
matrix block is governed by
𝜕 𝐶 𝑖 𝜕𝑡
=∇.(D∇𝐶 𝑖 ) , 𝑖 =1,..,𝑛 𝑐 −1 ,
(7.19)
where D is a matrix of diffusion coefficients. The elements of D can be calculated from Maxwell-
Stefan diffusion coefficients as follows
𝐷 = 𝐵 −1
𝛾 , (7.20)
137

where matrix B is a function of the inverse of the MS coefficients and γ represents the non-ideal
behavior of the mixture. Calculation of the multicomponent diffusion coefficients is explained in
length in the previous chapter.
Eq. 7.19 represents a set of nc-1 coupled partial differential equations. As a result, one cannot easily
obtain an exact analytical solution to this set of equations, or define correction factors to the
Warren and Root model to account for early-time behavior. In order to decouple this set of
equations and obtain analytical solutions to mass transfer in the multicomponent case, we use the
linearized theory of Toor et al. (1965). The basis of this theory is that Ct (overall molar density in
matrix) and D can be considered constant. Based on this assumption, and using matrix notation,
Eq. 7.19 reduces to
𝐶 𝑡 𝜕 𝒙 𝜕𝑡
=𝐶 𝑡 D∇
2
𝒙 ,
(7.21)
or equivalently
𝜕 𝒙 𝜕𝑡
=D∇
2
𝒙 .
(7.22)
In the above equations, x is the vector of mole fractions. The linearized theory proposed by Toor
exploits the properties of the modal matrix P whose columns are eigenvectors of D (Taylor and
Krishna, 2003). The matrix product
𝑃 −1
𝐷𝑃 =𝐷̂
(7.23)
is a diagonal matrix with elements that are the eigenvalues of D. Multiplying both sides of Eq.
7.22 by P
-1
results in
𝜕 𝒙̂
𝜕𝑡
=D
̂
∇
2
𝒙̂ .
(7.24)
138

where we have defined “pseudocompositions” 𝒙̂ by
𝒙̂ = 𝑃 −1
𝒙 (7.25)
The above transformation reduces the original set of nc-1 coupled partial differential equations to
nc-1 uncoupled partial differential equations in pseudocompositions. We can now write the
solution to the uncoupled multicomponent mass transfer problem in terms of pseudomole fractions.
Specifically, if a binary diffusion problem has the solution
𝑥 1
−𝑥 1.𝑖𝑛𝑖𝑡 𝑥 1,𝑓 −𝑥 1,𝑖𝑛𝑖𝑡 =𝑓 (𝑡 ,𝐷 12
)
(7.26)
then the corresponding multicomponent problem has the solution
𝑥̂
𝑖 −𝑥̂
𝑖 .𝑖𝑛𝑖𝑡 𝑥̂
𝑖 ,𝑓 −𝑥̂
1,𝑖𝑛𝑖𝑡 =𝑓 (𝑡 ,𝐷̂
𝑖 ) (7.27)
where 𝑥̂
𝑖 .𝑖𝑛𝑖𝑡 and 𝑥̂
𝑖 ,𝑓 are suitably transformed initial and boundary conditions. In order to recover
the solution to the original problem in terms of real mole fractions and diffusive fluxes, we can
apply the inverse transformations
𝒙 = 𝑃 𝒙̂ (7.28)
and
𝑱 = 𝑃 𝑱̂
(7.29)
Once the system of equations is decoupled given the above transformations, we can define the
correction factor to the Warren and Root formulation for each component similarly to the binary
system. The correction factor, based on pseudocompositions, can be calculated as
139

𝐵 𝑖 =
(|𝑥̂
𝑖 ,𝑓 −𝑥̂
𝑖 ,𝑖𝑛𝑖𝑡 |)+(|𝑥̂
𝑖 −𝑥̂
𝑖 .𝑖𝑛𝑖𝑡 |)
2𝑀𝑎𝑥 (𝘀 𝑥̂
𝑖 ,|𝑥̂
𝑖 −𝑥̂
𝑖 .𝑖𝑛𝑖𝑡 |)
, 𝑖 =1,…,𝑛 𝑐 −1 .
(7.30)
The pseudo-diffusive flux for each component can then be calculated as
𝐽̂
𝑖 =𝐵 𝑖 𝐷̂
𝑖 ,𝑖 (𝑥̂
𝑖 ,𝑓 −𝑥̂
𝑖 ) , 𝑖 =1,…,𝑛 𝑐 −1 (7.31)
And the diffusive flux can be calculated by transforming the pseudo flux using Eq. 7.29.
7.2.3. Extension to Other Geometries
The equations presented in the previous section were derived for a spherical geometry. Generally,
the spherical geometry presents a good estimate of a 3-D matrix block, while reducing the
complexity of modeling a 3-D problem. The presented approximation of the analytical solution to
the mass transfer problem based on Vermeulen’s formulation does not easily generalize to other
geometries. However, Zimmerman et al. (1993) show that Vermeulen’s approximation can be
extended to other block geometries with a high accuracy if the block has a known shape and
suitable “shape factor” values are used with each shape. Warren and Root left σ as an adjustable
parameter and finding the correct σ has been the subject of numerous research studies (e.g. Kazemi
et al. 1976, Thomas et al 1983, Chang 1993, Lim and Aziz 1995, Hassanzadeh et al. 2009). We
will further discuss the application of the presented approach to a linear 1-D block geometry in the
results section below.
7.3. Results and Discussion
In this section, we report calculation examples that highlight the importance of the proposed
correction factor and demonstrate the accuracy of the calculations based on the correction factor
for both binary and multicomponent fluid systems.
140

7.3.1. Example 1 – Binary Fluid System in Spherical Geometry
In the first example, we consider a case with one spherical matrix block with a radius of am=1 cm,
surrounded by fractures. We assume that the block is initially filled with CO2 (100% initial
concentration of CO2) and that the surrounding fracture has a zero concentration of CO2. We
further assume that the diffusion coefficient of the CO2 in this system is D = 10
-5
cm
2
/sec. The
exact analytical solution for the concentration within the matrix block follows Eq. 7.11 with the
following initial and boundary conditions
𝑐 𝑖𝑛𝑖𝑡 =1 (7.32)
𝑐 𝑓 =0 (7.33)
At constant pressure and temperature and no external forces, the main mechanism for mass transfer
between the matrix block and the surrounding fracture is molecular diffusion. Next, we perform
calculations for average concentration within matrix block (
𝑐 ̅−𝑐 𝑖𝑛𝑖𝑡 𝑐 𝑓 −𝑐 𝑖𝑛𝑖𝑡 ) and diffusive flux of CO2
versus time using various calculation methods. The calculation methods include
1. Analytical solution: results for the analytical solution are obtained based on Eq. 7.11, with
initial and boundary conditions listed in Eq. 7.32 and 7.33
2. Fine grid model: In this model, the matrix block is divided into 100 gridblocks and an
explicit, finite difference approach is used to solve the partial differential equation in Eq.
7.8 numerically.
3. Warren and Root model: Results from this model are calculated based on the Warren and
Root model as explained in Eq. 7.12. The approximation of the analytical solution is based
on a pseudo-steady state assumption.
141

4. Warren and Root model with correction factor: Calculations are made with a correction
factor multiplied by the Warren and Root model, as introduced in Eq. 7.16-7.18. At each
timestep, the diffusive flux and mass transfer rate to the fracture is calculated.
Concentrations are then updated accordingly for the next timestep.
The left panel of Figures 7.1 reports the average concentration within the matrix block (
𝑐 ̅−𝑐 𝑖𝑛𝑖𝑡 𝑐 𝑓 −𝑐 𝑖𝑛𝑖𝑡
)
versus time for each of the calculation methods mentioned above. The solid black line represents
the results of the exact analytical solution. The dashed black line represents the result of the finite-
difference model, and the blue and red lines represent the average concentration within the matrix
block based on the WR models, and WR model using the correction factor, respectively. As can
be seen in the left panel of Fig. 7.1, the Warren and Root model fails to capture the early time
behavior of molecular diffusion in this system. Adding the correction factor to the calculation
improves the results of the dual-porosity model significantly. We also observe that the fine grid
model (a single porosity model with further discretization inside the matrix block) also results in
an accurate approximation of the analytical solution. However, such calculations are very costly
in a large-scale reservoir simulator.

142

Figure 7.1 Average concentration vs time (left) and correction factor vs time (right) for binary fluid system in
spherical geometry
The right panel of Fig. 7.1. displays the correction factor applied to transfer rate calculations versus
time. Initially, the concentration of CO2 within the matrix block is different from the fracture
concentration and the correction factor assumes a large value. An 𝘀 value of 0.1 is used in these
calculations for numerical stability at early times. As the matrix concentration approaches the
fracture concentration, the correction factor is reduced, and ultimately reaches unity when pseudo-
steady state conditions are reached.
7.3.2. Example 2 – Binary Fluid System in Linear Geometry
As mentioned in the previous section, the same set of calculations obtained for a spherical
geometry can be used in a linear geometry, with suitable shape factors in the Warren and Root
equations to incorporate block shape and geometry. Zimmerman et al. suggest a shape factor of
3𝜋 2
𝐿 2
⁄ for a 1-D linear slab of thickness L. This example depicts the accuracy of the suggested
correction factor based on Vermeulen’s formulation where the matrix block is a 1-D linear system
143

of length L. Initial and boundary conditions of the system are the same as the spherical system,
and calculations in example 1 are repeated for a linear geometry.
The left panel of Figure 7.2 shows the average concentration in the matrix block versus time as
obtained from the analytical solution (solid black line), fine grid numerical simulation with n=100
(dashed black line), Warren and Root formulation (solid blue line), and Warren and Root
formulation incorporating the proposed correction factor (solid red line), while the right panel of
Fig. 7.2 depicts the change in the proposed correction factor with time. The correction factor
provides a near-perfect match with the analytical solution in this case as well, demonstrating the
accuracy of the approximation in linear geometry, provided that a suitable shape factor is used.

Figure 7.2 Average concentration vs time (left) and correction factor vs time (right) for binary fluid system in
spherical geometry

144

7.3.3. Example 3 – Experimental Data from a CO
2
Diffusion Experiments
In this example, we will use the proposed correction to the Warren and Root model to simulate
experimental results of CO2 diffusion in a water saturated porous medium. The experimental data
used in this example was published by Shi et al. (2017). The experimental setup consisted of a
reference cell and a sample cell, both kept in an air bath to ensure a constant temperature of 50 °C.
After loading a known mass of water and quartz in the sample cell the reference cell was charged
with CO2 and allowed to establish thermal equilibrium. The diffusion experiment was then
initiated by opening a valve between the reference and the sample cells to allow CO2 to enter the
sample cell. When the pressure in the sample cell reached a certain value, the valve between the
two cells was closed. Past this point, the pressure in the sample cell decreases due to CO2 diffusion
into the liquid phase. A detailed explanation of the interpretation of the pressure and temperature
data for calculation of the mass transfer rates of CO2 into the liquid phase is provided in their
paper.
The diffusion process in the sample cell can be represented as a 1-D linear matrix block as
illustrated in Figure 7.3. The mass transfer of CO2 from the vapor phase into the liquid phase is
described by Eq. 7.8, subject to the following initial and boundary conditions
𝑐 (𝑧 ,𝑡 =0)=0 (7.34)
𝑐 (0,𝑡 )=𝑐 𝑒𝑞
(𝑝 𝑠 ,𝑇 𝑠 ) (7.35)
𝜕𝑐
𝜕𝑧
(𝐿 ,𝑡 )=0
(7.36)
Here, c denotes the concentration of CO2 in the liquid phase (mol/cm
3
), ceq is the solubility
(mol/cm
3
) of CO2 in the liquid phase at the pressure (ps) and temperature (Ts) of the sample cell,
D is the effective CO2 diffusivity (cm
2
/s - diffusion coefficient divided by the tortuosity) in the
145

liquid. The quartz used in this example to represent a porous media has a particle size of 125-
150µm and a porosity of 0.45. The length of the porous medium is L = 5.83 cm and the effective
CO2 diffusivity is 3.37×10
-5
cm
2
/s.
Figure 7.4 illustrates the amount of CO2 in the vapor phase versus time from experimental data
(circles), fine grid simulation results with 100 cells in the matrix block (dashed black line), Warren
and Root model (blue line), and Warren and Root model with the proposed correction factor (red
line). It should be noted that it is not possible to obtain an exact analytical solution to this problem
in the form of Eq. 7.11, as the boundary condition is not constant and is changing with time, based
on the solubility of CO2 at the interface.

Figure 7.3 Schematic of the assumed 1-D diffusion process
146

Figure 7.4 Amount of CO
2
in vapor phase vs time
As can be seen, the Warren and Root model does not capture the behavior of the diffusion process
correctly. However, adding the correction factor substantially increase the accuracy of calculations
and provides for an acceptable agreement with the experimental data (without discretization of the
“matrix block”).
7.3.4. Example 4 – Multicomponent Fluid System in Spherical Geometry
Next, we demonstrate the applicability of the correction factor in a multicomponent fluid system.
The matrix block and surrounding fracture are defined in the same way as example 1. The fluid
used for this study is a 4-component fluid consisting of CO2, CH4, nC4, and nC12. The fluid
properties, along with the composition of the fluid in the matrix block and at the block interface
are listed in Table 7.1.

147

Table 7.1 Fluid description
Comp Tc [K] Pc [atm]
Acentric
factor
M
W

[g/mole]
Volume
shift
Matrix block
composition
Interface
composition
CO
2
304.21 72.86 0.2236 44.01 -0.0817 0.0 0.57
CH
4
190.56 45.39 0.0115 16.043 -0.154 0.3 0.14
nC
4
425.12 37.46 0.202 58.123 -0.06413 0.4 0.18
nC
12
658 17.96 0.5764 170.338 0.1150 0.3 0.11

The system is assumed to have a pressure of 150 bar and a temperature of 373.15 K. The matrix
of diffusion coefficients is calculated at the given pressure, temperature and composition based on
the formulae presented in Chapter 6 (Eq. 6.13-6.16). Since the matrix of diffusion coefficients are
assumed to be constant (basis for linearized theory of Toor), the elements of the matrix only need
to be calculated once. Eigenvalue decomposition of the matrix is performed and transformation of
variables are performed following Eq. 7.23 – 7.25.
The simulation is continued for 2×10
7
seconds and at each timestep, the diffusive flux and transfer
rate are calculated based on Eq. 7.29-7.30. Compositions are updated after each step and average
compositions in the block in each timestep is calculated as
𝑥 𝑖 −𝑥 𝑖 .𝑖𝑛𝑖𝑡 𝑥 𝑖 ,𝑓 −𝑥 𝑖 ,𝑖𝑛𝑖𝑡
(7.37)
It should be noted that the flux and transfer rates can only be calculated for nc-1 components.
Accordingly, the analytical and numerical results presented within this section are for the first nc-
1 components. The flux and transfer rate of the last component is calculated based on the fact
that there is no net diffusive flux within the system (∑ 𝐽 𝑖 =0
𝑛𝑐
𝑖 =1
). Figure 7.5 reports the average
composition within the matrix block for each of the nc-1 components versus simulation time.
Similar to the previous examples, the calculations with the proposed correction factors results in
148

higher accuracy with the analytical solution, compared to the original Warren and Root model.

Figure 7.5 Average composition of CO2 (top), nC4 (middle) and nC12 (bottom) in matrix block vs time
The correction factors used in the multicomponent example are defined independently for each
component. However, this does not pose any inconsistencies in the system, since the flux of the
149

last component is calculated in a way to counteract the diffusive fluxes of the other components,
such that there is no net diffusive flux in the system.
Based on the observations from the numerous examples presented in this section, we conclude that
the proposed correction factor be incorporated in molecular diffusion calculation in
multicomponent fluid systems in reservoir simulation to improve the accuracy of the predictions
in simulation of naturally fractured reservoirs.
7.3.5. Implementation in Reservoir Simulation
The suggested correction factor can be implemented in a dual-porosity compositional reservoir
simulator as a multiplier to the mass transfer rate between each set of matrix and fracture blocks
at every timestep. In this section, we will mention some considerations to be taken into account
for such an implementation.
The correction factor is dependent on the initial matrix composition, as well as matrix boundary
composition. Knowledge of initial compositions needs to be stored for each block. In the examples
presented in this chapter we have assumed that the boundary condition (composition at matrix-
fracture interface) is known. However, in a reservoir simulator, one needs to solve for the interface
composition between matrix and fracture at every timestep. This can be achieved by using film
theory (as explained in section 6.3.1) or by assuming an equilibrated mixture of matrix and fracture
fluids at the interface (as explained in section 6.7)
In film theory calculations, the correction factor appears in both sides of the equation (Eq. 6.18b)
as a multiplier to the diffusive flux. In addition, since the correction factor is dependent on the
interface composition, its value should be updated after each iteration on the interface composition.
However, this may cause numerical instabilities and singularity of the Hessian matrix. In contrast,
150

when a mixture of fracture and matrix fluids is used to represent the interface composition, the
correction factor does not appear in interface composition calculations. Once the interface
composition is found (using Eq. 6.28 and flashing the composition at average temperature and
pressure), the correction factor acts as a multiplier to the diffusive flux (e.g. Eq. 7.31)
Furthermore, in multicomponent fluid settings, the boost factor is defined based on
pseudocompositions. As a result, one needs to implement eigenvalue decomposition of the
diffusion coefficients between each set of matrix-fracture blocks at each timestep. Suitable
transformations of variables need to be performed to find pseudocompositions (Eq. 7.25).
7.4. Conclusions
Based on the formulations, examples and discussion presented in the previous sections, we arrive
at the following conclusions
• Molecular diffusion is an important mechanism of oil recovery in fractured reservoirs.
• Dual-porosity modeling of fractured reservoirs based on Warren and Root model does not
correctly represent the early-time behavior of molecular diffusion
• The proposed approach for calculating a correction factor to the Warren and Root model
provides significant improvements in capturing all-time behavior of the transfer rates between
fracture and matrix
• The proposed approach was tested against analytical solutions, fine scale numerical models,
and observed experimental data
• In a multicomponent system, eigenvalue decomposition of the diffusion coefficients matrix is
needed for calculation of correction factors. The decomposition is, however, less costly than
fine scale simulation.
151

Chapter 8 Summary and Future Research Directions
8.1. Summary
This research project focuses on improving the accuracy of compositional simulation in many
different aspects. In Chapters 3 and 4, we investigated the role of counter-current flow in
simultaneous water and gas injection (SWAG) and water alternating gas injection (WAG) EOR
processes. We showed that counter-current flow plays an important role in the flow mechanisms
during these processes and that the commercial simulation tools often overlook the importance of
counter-current flow in their calculations.
We used two different sets of 2-phase relative permeability curves for co-current flow and counter-
current flow and discussed how we transition from co-current to counter-current relative
permeability curves based on the direction of fluid velocities. We presented two models for
representing three-phase relative permeability calculations and accounted for hysteresis and IFT
scaling in our relative permeability models.
From several calculation examples, both in a 2-D reservoir model and a heterogeneous 3-D
reservoir models, we show that accounting for counter-current flow in our calculations greatly
impacts the simulation results in terms of saturation distribution, recovery and GOR. We then
performed a detailed sensitivity analysis on important design factors for an EOR project and
displayed the impact of accounting for counter-current flow for each setting.
We did not account for capillary pressure in our calculations in these two chapters. In order to
demonstrate the relative importance of capillary pressure, we performed calculations using Eclipse
300 for our 3-D reservoir model and showed that including capillarity does not significantly affect
152

the results. We concluded that the reduction in relative phase mobilities caused by the transition
between co-current and counter-current flow settings plays an important role in the dynamics of
SWAG and WAG injection processes, especially if the slug size is relatively small, and the vertical
permeability is relatively high.
In Chapter 5, we focused on improving phase behavior calculation tools in modeling of
unconventional reservoirs. Phase behavior of fluids in unconventional reservoirs is greatly
impacted by the presence of high capillary pressures caused by the interfacial tension exhibited in
micro- and meso-pores. Although the impact of capillarity on phase behavior is well-documented
by experimental measurements, available simulation tools still lack the appropriate
tools/algorithms to account for this effect.
We presented algorithms for stability testing for fluid systems affected by large values of capillary
pressure. We then demonstrated the efficiency and accuracy of the proposed algorithms for a range
of temperature and pressure conditions in oil and condensate systems. The analysis presented in
this chapter assumes a single average pore size to represent unconventional resources. Calculations
will be more complicated in practice, as a wide range of pore sizes are observed in reservoir rocks.
In Chapter 6, we presented a dual-porosity model in which the generalized Fick’s law is used for
molecular diffusion; and gas-oil diffusion at the fracture-matrix interface is modeled based on film
theory. A novel matrix-fracture transfer function was introduced for gas-oil diffusion based on
film theory.
We used field-scale examples to show that our approach, which is based on a sophisticated physical
model for gas-oil diffusion (film theory), gives significantly different results from the conventional
approach, in which classical Fick’s law is used to represent molecular diffusion and an average
153

interface composition is assumed. It was also demonstrated that the off-diagonal diffusion
coefficients can moderately impact the oil recovery during gas injection in fractured reservoirs.
Capillary forces were neglected in our formulations by assuming that IFT effects are minimal for
gas injection processes in oil reservoirs. When capillary forces are present, they act as a barrier to
gas-oil gravity drainage which is driven by density difference between oil in the matrix and gas in
fracture. Neglecting capillary effects will lead to overestimating oil recovery due to gravity
drainage. To overcome this, we have used small matrix permeabilities in the vertical direction to
minimize the role of gravity drainage in our example calculations.
In Chapter 7 we demonstrated the need for time-dependent transfer function for molecular
diffusion in a dual-porosity model. Dual-porosity models based on Warren and Root formulation
assume a pseudo-steady state flow condition between fracture and matrix and make no attempt in
solving the composition distribution within the matrix block. This assumption, while reducing the
complexity of calculations, results in erroneous outcomes that do not correctly predict the mass
transfer rates.
We presented the exact analytical solution to the concentration distribution within the matrix block
in a binary fluid system, provided an estimation to the analytical solution, and derived a “correction
factor” to the Warren and Root formulation for mass transfer rates between matrix and fracture.
We then exhibited the calculations for a multicomponent system, by using linearized theory of
Toor and eigenvalue decomposition of the diffusion coefficients.
We showed numerous examples in 1-D spherical and 1-D linear geometries and compared the
results obtained from the analytical solution to that of the proposed correction to the Warren and
Root model and concluded that the results of the proposed formulation lead to better
154

approximations of the analytical solution. We further validated the accuracy of the proposed model
by comparison of the results from our model with CO2 diffusion experimental data.
8.2. Future Research
In Chapters 3 and 4, accurate evaluation of counter-current relative permeability curves is essential
in improving the physical models used in WAG and SWAG displacements calculations. The
results of laboratory experiments for counter-current flow are, however, greatly impacted by the
boundary conditions imposed by the setup. Further investigation in lab experiments is
recommended in order to ensure the accuracy of the counter-current relative permeability
measurements. As an alternative, we suggest the use of pore-scale modeling to calculate counter-
current relative permeabilities, where the impact of the imposed boundary conditions can be
evaluated.
For future research in Chapter 5, we suggest using a J-function representation for capillary pressure
to account for the pore size distribution in a reservoir. A detailed study on the coupled impact of
phase behavior and mass transfer in unconventional reservoirs can also be performed.
In Chapter 6, capillary effects can be included in our dual-porosity formulation in future research
activities, where matrix-fracture transfer function needs to be re-defined carefully. The
competition among molecular diffusion, capillary forces and gravity drainage in terms of oil
recovery from matrix during gas injection in fractured reservoirs can be studied using such a dual-
porosity model.
In section 7.3.5 of Chapter 7, we discuss the implementation of the proposed formulation in a
compositional reservoir simulation. Future research in this area can be focused on implementing
155

the correction factor in a dual-porosity reservoir simulator and investigating the computational
efficiency and accuracy of the proposed model.

156

Nomenclature
am Radius of spherical block
A Cross sectional area
B Correction factor to Warren and Root formulation
c Concentration
ceq Solubility of CO2 in water
ct Total fluid compressibility
Ct Total fluid molar density
𝐶 𝑔 Land’s trapping coefficient for gas
𝐶 ℎ
Land’s trapping coefficient for hydrocarbon phase
𝐶 𝑜 Land’s trapping coefficient for oil
D Diffusion coefficient
Đ MS coefficient
Đ
o
Infinite dilution diffusion coefficient
𝐸 General exponent for interfacial tension model
f Objective function as defined by Eq. 6.18
𝐹 Interfacial tension scaling factor
𝑔 𝑖 Gradient of tm function
𝐻 𝑖𝑗
Hessian of tm function
𝐽 Leverett-J function
𝐽 𝑖𝑗
Jacobian of objective function, as described by Eq. 5.23
𝐽 𝑚 ,𝑓 Molar diffusive flux
𝐾 𝑖 Equilibrium constant
𝑘 𝑟 Relative permeability
𝑘 𝑟𝑔
Relative permeability to gas
𝑘 𝑟𝑔 ,𝑒 End-point relative permeability of gas
𝑘 𝑟𝑔
𝑛𝑜 _𝑖𝑓𝑡 Relative permeability to gas when interfacial tension between gas and oil is zero
𝑘 𝑟𝑔
𝑟𝑒𝑓 _𝑖𝑓𝑡 Relative permeability to gas at the reference state interfacial tension
𝑘 𝑟𝑜
Relative permeability to oil
𝑘 𝑟𝑜 ,𝑒 End-point relative permeability of oil
𝑘 𝑟𝑜𝑔 Relative permeability to oil in two-phase oil gas setting
𝑘 𝑟𝑜𝑤 Relative permeability to oil in two-phase water oil setting
𝑘 𝑟𝑜𝑤𝑚𝑎𝑥 Maximum relative permeability to oil in two-phase water oil setting
𝑘 𝑟𝑤
Relative permeability to water
n Number of moles
𝑛 𝑐 Number of components
ng Corey-exponent for gas
no Corey-exponent for oil
157

Nt Total molar flux
𝑃 Pressure
𝑃 𝑐 Capillary pressure
𝑃 𝑐 ,𝑖 Critical pressure
𝑝 𝑥 Liquid phase pressure
𝑃 𝑦 Vapor phase pressure
𝑃 𝑤 Trial phase pressure
𝑃 𝑧 Reference pressure
q Source term in the mass conservation equation
𝑅 Radius of curvature of interface, nm
𝑟 𝑐 Pore radius, nm
S Saturation
𝑆 𝑔𝑐
Critical gas saturation
𝑆 𝑔𝑓
Flowing gas saturation
𝑆 𝑔 𝑚𝑎𝑥
Maximum gas saturation obtained in a displacement
𝑆 𝑔𝑛𝑜𝑟 Normalized gas saturation
𝑆 𝑔𝑟
Residual gas saturation
𝑆 ℎ𝑓 Flowing hydrocarbon saturation
𝑆 ℎ
𝑚𝑎𝑥
Maximum hydrocarbon saturation
𝑆 ℎ𝑟 Residual hydrocarbon saturation
𝑆 𝑜𝑚
Residual oil saturation in a 3-phase system
𝑆 𝑜𝑛𝑜𝑟 Normalized oil saturation
𝑆 𝑜𝑟
Residual oil saturation
𝑆 𝑤𝑐
Connate water saturation
𝑆 𝑤𝑛𝑜𝑟 Normalized water saturation
t Time
𝑇 Temperature
𝑇 𝑐 ,𝑖 Critical temperature
𝑡𝑚 Reduced modified tangent plane distance
tpd Reduced tangent plane distance
TPD Tangent Plane Distance as defined by Eq. 5.5
𝑢 𝑥 Overall velocity in liquid phase
𝑢 𝑦 Overall velocity in vapor phase
𝑉 𝑡 Total fluid volume
𝑤 𝑖 Mole fraction of component i in the trial phase
𝑊 𝑖 Mole number of component i in the trial phase
𝑊 𝑡 Total mole numbers of the trial phase
𝑥 𝑖 Mole fraction of component i in the liquid phase
𝑦 𝑖 Mole fraction of component i in the gas phase
𝑧 𝑖 Mole fraction of component i in reference phase
158

𝛽 Scaling component for IFT scaling
Γd M-F transfer rate due to diffusion
Γe M-F transfer rate due to expansion
Γgd M-F transfer rate due to gravity drainage
δ Kronecher delta function
𝜃 Contact angle for wetting phase
𝜇 Chemical potential
ρ Density
𝜎 Interfacial tension (dyne/cm)
𝜎 0
Reference interfacial tension (dyne/cm)
σ
1,2
Shape factor
𝜑̂ Fugacity coefficient
ϕ Porosity
𝜒 Parachor coefficient

159

References
Alkandari H. A. 2002: Numerical Simulation of Gas-Oil Gravity Drainage for Centrifuge
Experiments and Scaled Reservoir Matrix Blocks (Doctoral dissertation, Colorado School of
Mines. Arthur Lakes Library).
Al-Kobaisi M., Kazemi H., Brown J. S. 2013: Compositional Phase Trapping in CO2 WAG
Simulation. Paper SPE 165983 presented at SPE Reservoir Characterization and Simulation
Conference and Exhibition, Abu Dhabi, UAE.
Alharthy N. S., Nguyen T.; Kazemi, H., Teklu, T., Graves, R. 2013: Multiphase Compositional
Modeling in Small-Scale Pores of Unconventional Shale Reservoirs. SPE Annual Technical
Conference and Exhibition. No. SPE-166306-MS
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, Paper SPE
63146 presented at SPE Annual Technical Conference and Exhibition, Dallas, Texas.
Aziz K., Settari A. 1979: Petroleum Reservoir Simulation. Applied Science Publishers, London.
Aziz K. and Wong T.W., 1988: Considerations in the Development of Multipurpose Reservoir
Simulation Models, proceedings of the 1st and 2nd International Forum on Reservoir
Simulation, Alpbach, Austria.
Baker L.E. 1998: Three-Phase Relative Permeability Correlations, Paper SPE 17369 presented at
the SPE/DOE Enhanced Oil Recovery Symposium, Tulsa, Oklahoma, 17–20 April.
Barenblatt G.I., Zheltov 1u.P. and Kochina, I.N. 1960: Basic concepts on the theory of seepage of
homogeneous liquids in fissured rocks. Prikladnaya Matematikai Mekhanika, Akad. Nauk,
SSSR 24(5): 852-864.
Bensten R., Manai A. 1993: On the Use of Conventional Co-current 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
Bhambri P., Mohanty K. K. 2005: Streamline Simulation of Four-Phase WAG Processes. Paper
SPE 96940 presented at Annual Technical Conference and Exhibition, Dallas, Texas
Blunt M. J. 2000: An empirical model for three-phase relative permeability. SPE Journal 5 (04),
435-445.
Bourbiaux B., Kalaydjian F., 1990: Experimental Study of Cocurrent and Countercurrent Flows in
Natural Porous Media, SPE Reservoir Engineering, 5 (3), 361-368
Brooks R.H. and Corey A.T. 1964: Hydraulic Properties of Porous Media. Colorado State
University, Hydrology Papers (5).
Carman, P. C. 1997: Fluid Flow through Granular Beds. Chemical Engineering Research and
Design, 75, S32-S48.
Caudle B. H., Dyes A. B. 1958: Improving miscible displacement by gas-water injection. Paper
presented at 32nd Annual Fall Meeting of SPE, Dallas, Texas.
Chang M.M. 1993: Deriving the Shape Factor of a Fractured Rock Matrix. Technical Report
NIPER-696 (DE93000170), NIPER, Bartlesville, Oklahoma.
Christensen J. R., Stenby E. H., Skauge A. 2001: Review of WAG field experience. SPE Reservoir
Evaluation & Engineering 4 (02), 97-106.
Coats K.H., 1980: An Equation of State Compositional Model. SPE Journal 20 (05), 363-376
160

Coats, K.H. 1989. Implicit Compositional Simulation of Single-Porosity and Dual-Porosity
Reservoirs. Paper SPE 18427 presented at the SPE Symposium on Reservoir Simulation,
Houston, TX, 6-8 February.
Crank J. 1975: The Mathematics of Diffusion: 2d Ed. Clarendon Press.
Cruichshank J., McDougall S. R., Sorbie K. S. 2002: Anchoring Methodologies for Pore-Scale
Network Models: Application to Relative Permeability and Capillary Pressure
Prediction. Petrophysics 43 (04), 365-375.
da Silva F.V ., and Belery P. 1989: Molecular Diffusion in Naturally Fractured Reservoirs: A
Decisive Recovery Mechanism. Paper SPE 19672 presented at the SPE Annual Technical
Conference and Exhibition, San Antonio, TX, 8-11 October.
Danesh A., Dandekar A. Y ., Todd A. C., Sarkar 1991: R. A Modified Scaling Law and Parachor
Method Approach for Improved Prediction of Interfacial Tension of Gas-Condensate Systems.
SPE Annual Technical Conference and Exhibition, No. SPE-22710-MS
Danesh, A. 1998: PVT and phase behavior of petroleum fluids. Elsevier, Amsterdam
Darvish G., Lindeberg E., Holt T., Kleppe,J., and Arild Utne S. 2006: Reservoir Conditions
Laboratory Experiments of CO2 injection into Fractured Cores. Paper SPE 99650 presented at
the SPE Europec/EAGE Annual Conference and Exhibition, Vienna, Austria, 12-15 June.
De la Cruze V ., Spanos T., 1991: Steady-state Countercurrent Flow in One Dimension. Transport
in Porous Media, 6 (2), 173-182.
Delshad M., Pope G.A. 1989: Comparison of the Three-Phase Oil Relative Permeability Models
Transport in Porous Media, 4 (1), 59-83.
Du L., Chu L. 2012: Understanding Anomalous Phase Behavior in Unconventional Oil Reservoirs.
In SPE Canadian Unconventional Resources Conference. No. SPE-161830-MS
Duncan J.B., and Toor H.L. 1962: An Experimental Study of Three Component Gas
Diffusion. AIChE Journal 8(1): 38-41.
Eastwood J., Spanos T., 1991: Steady-state Countercurrent Flow in One Dimension. Transport in
Porous Media, 6 (2), 173-182
Evans R., Marconi U. M. B., Tarazona P. 1986: Fluids in Narrow Pores: Adsorption, Capillary
Condensation, and Critical Points. The Journal of chemical physics, 84(4), 2376-2399.
Fai-Yengo V . A. 2014: Impact of Light Component Stripping During CO2 Injection in Bakken
Formation. Paper SPE 1922932 presented in Unconventional Resources Technology
Conference (URTEC), Denver, Colorado.
Firincioglu T., Ozkan E., Ozgen C. 2012: Thermodynamics of Multiphase Flow in Unconventional
Liquids-Rich Reservoirs. SPE Annual Technical Conference and Exhibition, No. SPE-159869-
MS
Firoozabadi A. 1999: Thermodynamics of Hydrocarbon Reservoirs, 1st Edition, McGraw-Hill.
Gilman J. R., & Kazemi H. 1988: Improved Calculations for Viscous and Gravity Displacement
in Matrix Blocks in Dual-Porosity Simulators. Journal of Petroleum Technology, 40(01), 60-
70.
Giovambattista N., Rossky P. J., Debenedetti P. G. 2006: Effect of Pressure on The Phase Behavior
and Structure of Water Confined Between Nanoscale Hydrophobic and Hydrophilic
Plates. Physical Review E 73(4), 041604.
Guo P., Wang Z., Shen P., and Du J. 2009: Molecular Diffusion Coefficients of the Multicomponent
Gas–Crude Oil Systems under High Temperature and Pressure. Ind. Eng. Chem. Res. 48(19):
9023–9027.
161

Hamouda A. A., Tabrizy V . A. 2013: The Effect of Light Gas on Miscible CO2 Flooding to Enhance
Oil Recovery from Sandstone and Chalk Reservoirs, Journal of Petroleum Science and
Engineering, 108, 259-266.
Hassanzadeh H., Pooladi-Darvish M., & Atabay S. 2009: Shape Factor in the Drawdown Solution
for Well Testing of Dual-Porosity Systems. Advances in Water Resources, 32(11), 1652-1663.
Haugen K., Firoozabadi A. 2009: Composition at the Interface between Multicomponent Non-
Equilibrium Phases. J. Chem. Phys. 130(6): 64707.
Hebden M. D. 1973: An Algorithm for Minimization Using Exact Second Derivatives. UKAEA
Theoretical Physics Division Harwell.
Hoteit H. 2013: Modeling Diffusion and Gas-Oil Mass Transfer in Fractured Reservoirs. Journal
of Petroleum Science and Engineering 105: 1-17.
Hoteit H., and Firoozabadi A. 2009: Numerical Modeling of Diffusion in Fractured Media for Gas-
Injection and –Recycling Schemes. SPE Journal 14: 323-337.
Hough E. W., Stegemeier G. L. 1961: Correlation of Surface and Interfacial Tension of Light
Hydrocarbons in the Critical Region. Society of Petroleum Engineers Journal, 1(04), 259-263.
Hu H., Whitson C.H., and Qi Y . 1991: A Study of Recovery Mechanisms in a Nitrogen Diffusion
Experiment. Paper SPE 21893 presented at the 6th European Symposium on Improved Oil
Recovery.
Irani M., Bozorgmehry R., Pishvaie M.R., Tavasoli A. 2009: A Study on the Effects of
Thermodynamic Non-Ideality and Mass Transfer on Multi-Phase Hydrodynamics Using CFD
Methods. World Academy of Science, Engineering and Technology 34: 627-632.
Jamili A., Willhite G.P., and Green D. 2011: Modeling Gas-Phase Mass Transfer between Fracture
and Matrix in Naturally Fractured Reservoirs. SPE Journal 16(4): 795-811.
Javaheri M., Jessen K., 2010: Integration of Counter-current Relative Permeability in the
Simulation of CO2 Injection into Saline Aquifers. International Journal of Greenhouse Gas
Control, 5 (5), 1272-1283.
Jenkins M. K. 1984: An analytical model for water/gas miscible displacements. SPE Enhanced Oil
Recovery Symposium. Society of Petroleum Engineers.
Jessen K., Gerritsen M., Mallison, B. 2008: High-Resolution Prediction of Enhanced Condensate
Recovery Processes. SPE Journal 13(2): 257-266.
Jin Z., Firoozabadi A. 2016: Thermodynamic Modeling of Phase Behavior in Shale Media. SPE
Journal, 21(1), 190-207.
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), pp.537-552
Kazemi H., & Gilman J. R. 1993: Multiphase Flow in Fractured Petroleum Reservoirs. In Flow
and Contaminant Transport in Fractured Rock (pp. 267-323).
Kazemi H., Merrill L.S., Porterfield K.L., and Zeman P.R. 1976: Numerical Simulation of Water-
Oil Flow in Naturally Fractured Reservoirs. SPE Journal 16(6): 317-326.
Killough J. E. 1976: Reservoir Simulation with History-dependent Saturation Functions, Trans.
AIME 261: 37-48.
Kooijman H.A., and Taylor R. 1991: Estimation of Diffusion Coefficients in Multicomponent
Liquid Systems. Industrial & Engineering Chemistry Research 30(6): 1217-1222.
Kozeny J. 1927: Über kapillare Leitung des Wassers im Boden: (Aufstieg, Versickerung und
Anwendung auf die Bewässerung). Hölder-Pichler-Tempsky.
Krishna R., and Standart G. 1976: Determination of Mass and Energy Fluxes for Multicomponent
Vapour–Liquid Systems. Letts Heat Mass Transfer 3: 173–182.
162

Lake L.W. 1989: Enhanced Oil Recovery. Preentice Hall, New Jersey.
Land C.S., 1968: Calculation of Imbibition Relative Permeability for Two- and Three-phase Flow
from Rock Properties, Society of Petroleum Engineers Journal, 8(02), 149-156.
Langaas K, Papatzacos P., 2001: Numerical Investigations of the Steady State Relative
Permeability of a Simplified Porous Medium. Transport in Porous Media, 45 (2), 241-266.
Larsen J. A., Skauge A. 1999: Simulation of the Immiscible WAG Process Using Cycle-Dependent
Three-Phase Relative Permeabilities. Paper SPE 56475 presented at SPE Annual Technical
Conference and Exhibition, Houston, Texas.
Leahy-Dios A., and Firoozabadi A. 2007: Unified Model for Nonideal Multicomponent Molecular
Diffusion Coefficients. AIChE Journal 53(11): 2932-2939.
Lee S. T., Chien M. C. H. 1984: A New Multicomponent Surface Tension Correlation Based on
Scaling Theory. SPE Enhanced Oil Recovery Symposium. No. SPE-12643-MS
Lelievre R.F. 1966: Etude D’ecoulement Diphasiques a Contre-courant en Milieu Poreux,
Comparaison avec les Ecoulements de Mems Sens. University of Toulouse, France, PhD
Dissertation
Li D., Kumar K., Mohanty K. K. 2005: Compositional Simulation of WAG Processes for a Viscous
Oil. Paper SPE 84074 presented at SPE Annual Technical Conference and Exhibition, Denver,
Colorado.
Li G., Karpyn Z. T., Halleck P. N., Grader A. S., 2005: Numerical Simulation of a CT Scanned
Counter-current Flow Experiment. Transport in Porous Media, 60 (2), 225-240.
Lim K.T., and Aziz K. 1995: Matrix-Fracture Transfer Shape Factors for Dual-Porosity
Simulators. Journal of Petroleum Science and Engineering 13(3): 169-178.
Litvak, B. L. 1986: Simulation and Characterization of Naturally Fractured Reservoirs.
In Reservoir Characterization (pp. 561-584).
Lu H., Di Donato G., and Blunt M. 2008: General Transfer Functions for Multiphase Flow in
Fractured Reservoirs. SPE Journal 13(3): 289-297.
Macleod D.B. 1923: On a Relation between Surface Tension and Density, Transactions of the
Faraday Society, 19, 38-41
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), pp.72–77
Mani V ., Mohanty K. K. 1998: Pore-Level Network Modeling of Three-Phase Capillary Pressure
and Relative Permeability Curves. SPE Journal, 3 (03), 238-248.
Michelsen M. L. 1982: The Isothermal Flash Problem. Part I. Stability. Fluid Phase
Equilibria, 9(1), 1-19.
Moortgat J., and Firoozabadi A. 2013: Fickian Diffusion in Discrete-Fractured Media from
Chemical Potential Gradients and Comparison to Experiment. Energy & Fuels 27: 5793-5805.
Nojabaei B., Johns R. T., Chu L. 2013: Effect of Capillary Pressure on Fluid Density and Phase
Behavior in Tight Rocks and Shales. SPE Reservoir Evaluation and Engineering,16(03), 281-
289.
Nojabaei B., Siripatrachai N., Johns R. T., Ertekin T. 2014: Effect of Saturation Dependent
Capillary Pressure on Production in Tight Rocks and Shales: A Compositionally-Extended
Black Oil Formulation. SPE Eastern Regional Meeting. No. SPE-171028-MS
Odeh A., 1959: Effect of viscosity ratio on relative permeability. Trans. AIME, 216, pp.346–352
163

Orbach H., Crowe C.M. 1971: Convergence Promotion in the Simulation of Chemical Processes
with Recycle – The Dominant Eigenvalue Method. The Canadian Journal of Chemical
Engineering, 49, 509-513.
Orr F.M. 2007: Analytical Theory of Gas Injection Processes. Tie Line Publications, Denmark.
Pedersen K.S., Christensen P.L. 2007: Phase Behavior of Petroleum Reservoir Fluids, CRC Press,
Taylor & Francis Group.
Pejic D., Maini B. B. 2003: Three-phase Relative Permeability of Petroleum Reservoirs. Paper
SPE 81021 presented at SPE Latin American and Caribbean Petroleum Engineering
Conference, Port-of-Spain, Trinidad and Tobago.
Perkins T.K., and Johnston O.C. 1963: A Review of Diffusion and Dispersion in Porous Media.
SPE Journal 3: 70-84.
Polanyi M. 1932: Theories of the Adsorption of Gases: A General Survey and Some Additional
Remarks. Transactions of the Faraday Society 28: 316–333.
Ramirez B. A., Kazemi H., Al-Kobaisi M., Ozkan E., & Atan S. 2007: A Critical Review for Proper
Use of Water-Oil-Gas Transfer Functions in Dual-Porosity Naturally Fractured Reservoirs-Part
I. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers.
Rasmussen C.P., Krejbjerg K., Michelsen M.L., Bjurstrom K.E. 2006: Increasing the
Computational Speed of Flash Calculations with Applications for Compositional, Transient
Simulations. SPE Reservoir Evaluation and Engineering, 9, 2-38
Restagno F., Bocquet L., Biben T. 2000: Metastability and Nucleation in Capillary
Condensation. Physical Review Letters, 84(11), 2433.
Rossen W. R., Van Duijn C. J., Nguyen Q. P., Shen C., Vikingstad A. K. 2010: Injection Strategies
to Overcome Gravity Segregation in Simultaneous Gas And Water Injection into
Homogeneous Reservoirs. SPE Journal, 15 (01), 76-90.
Saidi A. M. 1983: Simulation of Naturally Fractured Reservoirs. In SPE Reservoir Simulation
Symposium. Society of Petroleum Engineers.
Sandoval Lemus D.R., Yan W., Michelsen M. L., Stenby E. H. 2016: The phase envelope of
multicomponent mixtures in the presence of a capillary pressure difference. Industrial &
Engineering Chemistry Research, 55: 6530-6538
Schechter D. S., Guo B. 1998: Parachors Based on Modern Physics and Their Uses in IFT
Prediction of Reservoir Fluids. SPE Reservoir Evaluation & Engineering, 1(03), 207-217.
Shapiro A. A., Potsch K., Kristensen J. G., Stenby E. H. 2000: Effect of Low Permeable Porous
Media on Behavior of Gas Condensates. SPE European Petroleum Conference. No. SPE-
65182-MS
Shapiro A.A., and Stenby E.H. 1997: Kelvin Equation for a Non-ideal Multicomponent Mixture.
Fluid Phase Equilibria 134 (1-2) pp 87-101
Shapiro A.A., and Stenby, E.H. 1998. Potential Theory of Multicomponent Adsorption. Journal of
Colloid and Interface Science 201: 146–157.
Sherafati M., Javaheri M., Jessen K. 2014: The Role of Counter-Current Flow in Simultaneous
Water and Gas Injection Processes. Paper SPE 169553 presented at SPE Western North
American and Rocky Mountain Joint Meeting, Denver, Colorado.
Shi Z., Wen, B.,Hesse, M., Tsotsis T. T., & Jessen K. 2017: Measurement and Modeling of CO2
Mass Transfer in Brine at Reservoir Conditions. Advances in Water Resources.
Shojaei H., & Jessen K. 2015: Diffusion and Matrix-Fracture Interactions During Gas Injection in
Fractured Reservoirs. In IOR 2015-18th European Symposium on Improved Oil Recovery.
164

Shojaei H., Rastegar R., and Jessen K. 2012: Mixing and Mass Transfer in Multicontact Miscible
Displacements. Transport in Porous Media 94(3): 837-857.
Soave G. 1972: Equilibrium Constants from a Modified Redlich-Kwong Equation of State.
Chemical Engineering Science, 27 (6), 1197-1203.
Song Z., Li Z., Wei M., Lai F., Bai B. 2014: Sensitivity Analysis of Water-alternating-CO2
Flooding for Enhanced Oil Recovery in High Water Cut Reservoirs, Computers and Fluids 99:
93-103
Sonier F., Souillard P., & Blaskovich F. T. 1988: Numerical Simulation of Naturally Fractured
Reservoirs. SPE Reservoir Engineering, 3(04), 1-114.
Stone H. L. 1982: Vertical Conformance in an Alternating Water-Miscible Gas Flood. Paper SPE
11130 presented at SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana.
Sugden S. 1924: The Variation of Surface Tension with Temperature and Some Related Functions,
Journal of Chemical Society, Transactions, 125, 32-41
Taylor R., and Krishna R. 1993: Multicomponent Mass Transfer. John Wiley & Sons, New York.
Thomas L. K., Dixon T. N., & Pierson R. G. 1983: Fractured Reservoir Simulation. SPEJ, Soc.
Pet. Eng. J.;(United States), 23(1).
Toor H. L., Seshadri C. V ., & Arnold K. R. 1965: Diffusion and Mass Transfer in Multicomponent
Mixtures of Ideal Gases. AIChE Journal, 11(4), 746-755.
Vega B., O’Brien W.J. Kovscek A.R. 2010: Experimental Investigation of Oil Recovery from
Siliceous Shale by Miscible CO2 Injection. Paper SPE 135627 presented at the SPE Annual
Technical Conference and Exhibition, Florence, Italy, 19-22 September.
Vermeulen T. 1953: Theory for irreversible and constant-pattern solid diffusion. Industrial &
Engineering Chemistry, 45(8), 1664-1670.
Wang L., Parsa E., Gao Y ., Ok J. T., Neeves K., Yin X., Ozkan E. 2014: Experimental Study and
Modeling of the Effect of Nanoconfinement on Hydrocarbon Phase Behavior in
Unconventional Reservoirs. SPE Western North American and Rocky Mountain Joint Meeting.
No. SPE-169581-MS
Wang Y ., Yan B., Killough J. 2013: Compositional Modeling of Tight Oil Using Dynamic
Nanopore Properties. SPE Annual Technical Conference and Exhibition, No. SPE-166267-MS
Warren J.E., and Root P.J. 1963: The Behavior of Naturally Fractured Reservoirs. SPE
Journal 3(3): 245-255.
Watts J.W. 1986: A Compositional Formulation of the Pressure and Saturation Equations. Paper
SPE 12244, SPE Reservoir Engineering, 1(3): 243-252
Wei Yan, et al. 2004: Three-phase Compositional Streamline Simulation and Its Application to
WAG. SPE/DOE Symposium on Improved Oil Recovery.
Weinaug C. F., Katz D. L. 1943: Surface Tensions of Methane-Propane Mixtures. Industrial &
Engineering Chemistry, 35(2), 239-246.
Wesselingh J.A., and Krishna R. 1990: Mass Transfer. Ellis Horwood, Chichester, England.
Whitaker S., 1986: Flow in porous media II: The governing equations for immiscible, two-phase
flow. Transport in Porous Media, 1(2), pp.105-125.
Yan W., Michelsen M. L., Stenby E. H., Berenblyum R. A., Shapiro A. A. 2004: Three-phase
Compositional Streamline Simulation and Its Application to WAG. Paper SPE 89440 presented
at SPE/DOE Symposium on Improved Oil Recovery, Tulsa, Oklahoma.
Yan B., Alfi M., Wang Y., Killough J. E. 2013: A New Approach for the Simulation of Fluid Flow
in Unconventional Reservoirs through Multiple Permeability Modeling. SPE Annual Technical
Conference and Exhibition. No. SPE-116173-MS
165

Zhou D., Fayers F. J., Orr Jr F. M. 1997: Scaling of Multiphase Flow in Simple Heterogeneous
Porous Media. SPE Reservoir Engineering, 12 (03), 173-178.
Zimmerman R.W., Chen G., Hadgu T., and Bodvarsson G.S. 1993: A Numerical Dual-Porosity
Model with Semianalytical Treatment of Fracture/Matrix Flow. Water resources
research 29(7): 2127-2137.

Abstract (if available)

Abstract Compositional reservoir simulation is a powerful tool commonly used to estimate the potential hydrocarbon recovery from Enhanced Oil Recovery methods. Successful design and implementation of EOR methods are highly dependent on the accuracy by which simulation tools can represent the physics that govern the displacement behavior in a reservoir. In this work, we investigate, and aim to improve the accuracy of some physical models that are frequently used to describe and interpret multiphase flows via compositional reservoir simulation. ❧ In the simulation of gas injection EOR methods, it is common to use relative permeability values obtained from co-current experiments and it is assumed that the viscous coupling between flowing phases is negligible. In this work we show that Water Alternating Gas (WAG) and Simultaneous Water and Gas (SWAG) injection processes are dominated by counter-current flow in the vertical direction. Previous experimental and simulation studies on counter-current flow have shown that counter-current relative permeability is always less than the co-current relative permeability. Most of the previous studies on the subject of counter-current flow have focused on two-phase flow. In this work we focus on counter-current flow in a three-phase flow system such as WAG and SWAG injection processes. ❧ We use both co-current and counter-current relative permeability functions in the simulation of these gas injection EOR processes (co-current flow in the horizontal direction and counter-current flow in the vertical direction) and show that the saturation profiles, oil recoveries and gas oil ratios are significantly influenced by the transition of flow from co-current to counter-current mode. Accounting for counter-current relative permeability results in higher oil recovery and lower GOR in all cases. The extent to which the conventional flow model and the flow model that accounts for counter-current relative permeability differ, relies on the properties of the reservoir (i.e. vertical to horizontal permeability ratio) and the history of production and flooding of the reservoir (e.g. initial fluid saturation in the reservoir at the beginning of the WAG/SWAG flood). ❧ Next, we focus on the phase behavior of fluids in unconventional reservoirs. Despite the large potential of unconventional resources, many unknowns still exist regarding the physics controlling the extraction processes in these settings. These include accurate representation of phase equilibrium in tight formations and effective implementation of relevant models in simulation tools. When a fluid is confined in pore spaces of nanometer size, significant interfacial curvatures may occur that can result in large capillary pressures between the liquid and vapor phases. The pressure difference between the two phases will likely affect the vapor-liquid equilibrium state. Previous efforts have shown that this effect is negligible for conventional reservoirs (with pores in the micron range) and current commercial reservoir simulators commonly ignore the effect of capillary pressure in the VLE calculations. However, experimental and modeling efforts have shown that ignoring capillary pressure in the VLE calculations will not be a valid approximation for unconventional (tight) reservoirs. We analyze the numerical aspects of including capillarity phenomena in VLE calculations, in an effort, to arrive at robust and efficient algorithms for stability analysis that can be used in compositional modeling/simulation of unconventional reservoirs. ❧ Next, we aim our focus on molecular diffusion and mass transfer during gas injection in fractured reservoirs. Molecular diffusion can play a significant role in oil recovery during gas injection in fractured reservoirs, especially when matrix permeability is low and fracture intensity is high. Diffusion of gas components from a fracture into the matrix extracts oil components from matrix and delays, to some extent, the gas breakthrough. This in turn increases both sweep and displacement efficiencies. ❧ In current simulation models, molecular diffusion is commonly modeled using a classical Fick’s law approach with constant diffusion coefficients. In the classical Fick’s law approach, the dragging effects (off-diagonal diffusion coefficients) are neglected. In addition, the gas-oil diffusion at the fracture-matrix interface is normally modeled by assuming an average composition at the interface which does not have a sound physical basis. ❧ We present a dual-porosity model in which the generalized Fick’s law is used to describe molecular diffusion to account for dragging effects. Gas-oil diffusion across the fracture-matrix interface is modeled based on film theory in which the gas in fracture and oil in the matrix are assumed to be at equilibrium at the interface. A novel and consistent matrix-fracture (M-F) transfer function is introduced for gas-oil diffusion based on film theory. Diffusion coefficients are calculated using the Maxwell-Stefan approach and are pressure, temperature and composition dependent. ❧ Finally, we will discuss a common problem in coarse scale modeling of molecular diffusion in naturally fractured reservoirs. In commercial reservoir simulation tools, transfer rates between the flowing and the stagnant domain in a dual-porosity model are modeled using the Warren and Roots model, in which no attempt is made to solve diffusion equations within each block and quasi-steady state transfer is assumed between the two domains. Previous researchers have shown that such a simplification does not correctly represent the mass transfer between the domains in a dynamic displacement. ❧ A more accurate approach to modeling transfer rates in such reservoirs would be to further discretize the matrix blocks but such an approach lacks computational efficiency for the purposes of large scale reservoir simulation, because orders of magnitude more computational cells are required for this approach. ❧ In the final chapter, we present a semi-analytical approach to model multicomponent molecular diffusion within each matrix block in a coarse-scale simulation model and develop equations for time-dependent transfer functions between the flowing and stagnant domains. The time-dependency of the transfer functions are characterized based on the composition of domains. Generalized Fick’s law is used to describe the diffusive fluxes to account for dragging effects between components. Analytical and semi-analytical solutions to the multicomponent mass transfer equations are obtained using the linearized theory of Toor and eigenvalue decomposition of the diffusion coefficient matrix. ❧ To demonstrate the accuracy of the proposed semi-analytical approach, various examples are presented. In these examples the results of the suggested approach are compared to the Warren-root model, analytical solution and high-resolution fine-scale models for 1-D linear and spherical geometries for a variety of fluid compositions. We observe that our approach (coarse-scale model using time-dependent transfer functions) accurately represents the analytical solution, with far lower CPU time requirement compared to high-resolution fine-scale models. Furthermore, we investigate experimental measurements of mass transfer rates in binary mixtures and further validate the precision of the proposed model by comparison of the results obtained from our model with the experimental data. ❧ The novel and consistent matrix-fracture transfer function for coarse-scale models, based on semi-analytical solutions to multicomponent mass transfer, overcome the inaccuracies of conventional treatment of diffusion in dual porosity models, while eliminating the need for further discretization of the matrix blocks. This allows for accurate, yet efficient simulation of mass transfer in fractured reservoirs.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Mass transfer during enhanced hydrocarbon recovery by gas injection processes

PDF

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

PDF

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

PDF

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

$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

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

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

PDF

Reactivation of multiple faults in oilfields with injection and production

PDF

Investigation of gas transport and sorption in shales

$An extended finite element method based modeling of hydraulic fracturing$

PDF

An extended finite element method based modeling of hydraulic fracturing

PDF

Application of data-driven modeling in basin-wide analysis of unconventional resources, including domain expertise

PDF

Real-time reservoir characterization and optimization during immiscible displacement processes

$A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow$

PDF

A coupled poromechanics-damage mechanics approach to model fracturing in multiphase flow

PDF

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

$History matching for fractured reservoirs with Hough-transform-based parameterization$

PDF

History matching for fractured reservoirs with Hough-transform-based parameterization

PDF

Flood front tracking and continuous recording of time lag in immiscible displacement

$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

Asset Metadata

Creator Sherafati, Marjan (author)

Core Title Modeling and simulation of complex recovery processes

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Publication Date 07/29/2018

Defense Date 03/23/2018

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

Tag coarse-scale modeling,compositional reservoir simulation,molecular diffusion,OAI-PMH Harvest,phase behavior,unconventional reservoirs,viscous coupling,water-alternating-gas

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

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

Creator Email marjan.sherafati@gmail.com,sherafat@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-39191

Unique identifier UC11671368

Identifier etd-SherafatiM-6451.pdf (filename),usctheses-c89-39191 (legacy record id)

Legacy Identifier etd-SherafatiM-6451.pdf

Dmrecord 39191

Document Type Dissertation

Format application/pdf (imt)

Rights Sherafati, Marjan

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