University of Southern California Dissertations and Theses

Mass transfer during enhanced hydrocarbon recovery by gas injection processes

(USC Thesis Other)

Mass transfer during enhanced hydrocarbon recovery by gas injection processes

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

MASS TRANSFER DURING ENHANCED HYDROCARBON RECOVERY BY
GAS INJECTION PROCESSES
by
Hasan Shojaei

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)

May 2014

Copyright 2014 Hasan Shojaei
i

“Yesterday I was clever, so I wanted to change the world. Today I am wise, so I am changing
myself.” - Rumi

ii

Dedication

To my parents
For their unconditional love, support and encouragement

iii

Acknowledgements

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 Roger Ghanem, as my dissertation committee members. Their valuable
comments and recommendations have significantly improved the quality of my doctoral
dissertation.
Professor Jessen and Professor Ershaghi have been more than great mentors to me. Every
time I was stuck in my research, or disappointed by hardships that an international student with
Iranian passport and single-entry visa (!) may encounter in the U.S., I just needed to talk to either
of them to regain my full enthusiasm and energy. Winning the first place in SPE student paper
contest (PhD division) in Western North America, would have not been possible without their
continuous encouragement and critical comments during rehearsal sessions.
iv

I would like to extend my gratitude to the Graduate School at USC for supporting my
doctoral studies through the Provost’s PhD Fellowship, one of the most generous scholarships
one can find in a top U.S. university. I also owe my gratitude to the Graduate School, 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 Occidental
Petroleum for the internship opportunity.
I am grateful to my amazing friends in the greater Los Angeles area with whom I have shared
wonderful memories. My special thanks go to Shahram Farhadi, Arman Khodabakhsh, Reza
Rastegar, Cyrus Ashayeri, Siavash Hakim-Elahi, Mohammad Javaheri, Mohammad Evazi and
Mehran Rahmani who have always been there to support and encourage me. I have also enjoyed
the company of Hamidreza Jahangiri, Hamed Haddadzadegan, Ehsan Tajer,
Dalad Nattwongasem, Mehran Hosseini, Mehrdad Mahdavi, Bahar Hosseini, Sahra
Homayounian, Asal Rahimi, Sanaz Norouzi, Tayeb Ayatollahi, Mahshad Samnejad, Marjan
Sherafati, Nadia Kadkhodayan, Gunel Murad and Lawrence Bustos. I am grateful to all of them,
and to many more friends not mentioned here, for their constant support and encouragement.
Having been born and raised in a far-off village in south of Iran, where even some
fundamental educational services were not available at the time, I could have not managed to
achieve the highest academic degree without the love, encouragement and support of my family.
My deepest gratitude goes to my lovely parents, Hossein and Madineh, my brilliant brother,
Mohammad Reza, and my beautiful sisters, Sakineh and Salimeh.

v

Abstract

In order to estimate the potential incremental hydrocarbon recovery by gas injection,
compositional reservoir simulators are commonly used in the industry. Successful design and
implementation of gas injection processes (mainly CO
2
) rely in part on the accuracy by which
the available simulation tools can represent the physics that govern the displacement behavior in
a reservoir.
In the first part of this thesis, we investigate the accuracy of some physical models that are
frequently used to describe and interpret dispersive mixing and mass transfer in compositional
reservoir simulation. The calculations from compositional simulation are compared with the
experimental observations.
A quaternary analog fluid system (alcohol-water-hydrocarbon) that mimics the phase
behavior of CO
2
-hydrocarbon mixtures at high pressure and temperature has been designed in
our research group. A porous medium was designed using Teflon materials to ensure that the
analog oil acts as the wetting phase, and the properties of the porous medium were characterized
in terms of porosity and permeability. Relative permeability and interfacial tension
measurements were also performed to delineate interactions between the fluid system and the
porous medium. Displacement experiments at First-contact miscible (FCM), near-miscible and
multicontact miscible (MCM) conditions were consequently performed (Rastegar, 2010).
The effluent concentrations from two-component FCM displacement experiments exhibit a
tailing behavior that is attributed to imperfect sweep of the porous medium: A feature that is not
vi

captured by normal dispersion models. To represent this behavior in displacement calculations,
we use dual-porosity (DP) models including mass transfer between flowing and stagnant
porosities. The 4-component two-phase displacement experiments (near-miscible and MCM) are
consequently simulated using the DP models constructed based on observations from FCM
displacements.
We demonstrate that the accuracy of our displacement calculations relative to the
experimental observations is sensitive to the selected models for dispersive mixing, mass transfer
between flowing and stagnant porosities, and IFT scaling of relative permeability functions. We
also demonstrate that numerical calculations substantially agree with the experimental
observations for some physical models with limited need for model parameter adjustment.
The second part of this thesis is devoted to diffusion and matrix-fracture interactions 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.
In this work, we present a dual-porosity model in which the generalized Fick’s law is used
for molecular diffusion to account for the dragging effects; and gas-oil diffusion at the fracture-
vii

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. A novel matrix-fracture transfer function is introduced for gas-
oil diffusion based on film theory. Diffusion coefficients are calculated using the Maxwell-
Stefan model and are pressure, temperature and composition dependent. A time-dependent
transfer function is used for matrix-fracture exchange in which the shape factor is adjusted using
a boost factor to differentiate between the transfer rate at early and late times.
Field-scale examples are used 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. It is also demonstrated that the dragging effects (off-diagonal diffusion
coefficients) and time-dependency of matrix-fracture transfer function can moderately impact the
oil recovery during gas injection in fractured reservoirs. We also show that miscibility is not
developed in the matrix blocks even at pressures above minimum miscibility pressure (MMP)
when molecular diffusion is the main recovery mechanism during gas injection in fractured
reservoirs.
In recent years, coalbed methane has become an important source of energy in the United
States. Since primary production techniques typically recover less than half of the methane in a
coalbed, enhanced coalbed methane (ECBM) recovery processes are needed in which CO
2

and/or N
2
are injected into the coalbed to recover more CH
4
.
One of the main mechanisms that govern the dynamics of ECBM recovery processes is the
sorption of gases onto the coal surfaces. Despite the well-documented complexity of
multicomponent sorption phenomena, adsorption and desorption of CH
4
/CO
2
/N
2
mixtures in the
porous coal is commonly modeled with the extended Langmuir model. The extended Langmuir
viii

model has been proven unable to accurately describe the multicomponent sorption that is central
to ECBM recovery processes and, therefore, more sophisticated sorption models are needed.
In the third part of this thesis we apply potential theory to describe the multicomponent
sorption of relevance to ECBM processes. In this approach for modeling multicomponent
sorption, each component is assumed to be affected by a characteristic potential field emitted by
the coal surface. We discuss the implementation of potential theory with emphasis on the
simulation of ECBM processes where computational efficiency and accuracy must be balanced.
The model must be solved by an iterative scheme, and is hence more computationally expensive
than the extended Langmuir approach.
The results and analysis presented in this work demonstrate that the application of potential
theory of sorption to modeling of ECBM recovery processes can improve the accuracy of
calculations. However, the additional complexity of the model and the associated increase in the
computational efforts may not balance the gain in accuracy sufficiently to warrant application in
general purpose reservoir simulation.

ix

Table of Contents
Abstract ........................................................................................................................................... v
List of Tables ............................................................................................................................... xiii
List of Figures .............................................................................................................................. xiv
Chapter 1 Introduction .................................................................................................................... 1
1.1 Background ........................................................................................................................... 1
1.2 Motivation ............................................................................................................................. 4
1.3 Objectives ............................................................................................................................. 5
1.4 Manuscript Organization ...................................................................................................... 7
Chapter 2 Literature Review ........................................................................................................... 9
2.1 Miscibility ............................................................................................................................. 9
2.2 Mixing Mechanisms............................................................................................................ 10
2.3 Effect of Mixing on Miscibility .......................................................................................... 12
2.4 Mass Transfer in Dual-Porosity (DP) Systems ................................................................... 13
2.5 Sorption in Coalbeds ........................................................................................................... 15
2.6 Mathematical Modeling ...................................................................................................... 16
Chapter 3 Experimental and Modeling Study of Multicontact Miscible Displacements ............. 19
3.1 Introduction ......................................................................................................................... 19
3.2 Fluid System ....................................................................................................................... 20
3.2.1 Phase Behavior............................................................................................................. 22
x

3.2.2 Mixture Densities ......................................................................................................... 22
3.2.3 Mixture Viscosities ...................................................................................................... 24
3.2.4 Interfacial Tension ....................................................................................................... 25
3.3 Packed Column ................................................................................................................... 27
3.4 Solid-Fluid Interactions ...................................................................................................... 27
3.4.1 Wettability .................................................................................................................... 27
3.4.2 Relative Permeability ................................................................................................... 28
3.5 Displacement Experiments ................................................................................................. 29
3.5.1 FCM Displacements ..................................................................................................... 30
3.5.2 Near-miscible and MCM Displacements ..................................................................... 30
3.6 Mathematical Modeling ................................................................................................. 36
3.7 Displacement Calculations.................................................................................................. 38
3.7.1 FCM Displacements ..................................................................................................... 39
3.7.2 Near-miscible and MCM Displacements ..................................................................... 40
3.8 Effect of Mixing on Near-miscible and MCM Displacements ........................................... 41
3.9 Importance of Selected Physical Models ............................................................................ 45
3.10 Discussion ......................................................................................................................... 49
3.11 Conclusions ....................................................................................................................... 54
Chapter 4 Diffusion and Matrix-Fracture Interactions during Gas Injection in Fractured
Reservoirs ..................................................................................................................................... 56
xi

4.1 Introduction ......................................................................................................................... 56
4.2 Mathematical Model ........................................................................................................... 60
4.3 Molecular Diffusion ............................................................................................................ 63
4.3.1 Intra-phase Diffusion ................................................................................................... 64
4.3.2 Cross-phase Diffusion .................................................................................................. 65
4.4 Transfer Function ................................................................................................................ 68
4.5 Results ................................................................................................................................. 72
4.5.1 Example 1 .................................................................................................................... 74
4.5.2 Example 2 .................................................................................................................... 79
4.6 Discussion ........................................................................................................................... 82
4.7 Conclusions ......................................................................................................................... 87
Chapter 5 Application of Potential Theory of Adsorption to Modeling of ECBM Recovery ...... 89
5.1 Introduction ......................................................................................................................... 89
5.2 Potential Theory .................................................................................................................. 91
5.3 Numerical Approach ........................................................................................................... 93
5.4 Application of MPTA to Coal............................................................................................. 95
5.5 Discussion and Conclusions ............................................................................................. 103
Chapter 6 Summary and Future Research Directions ................................................................. 105
Nomenclature .............................................................................................................................. 108
References ................................................................................................................................... 111
xii

Appendix A: UNIQUAC Model for Phase Equilibria ............................................................ 124
Appendix B: UNIQUAC Viscosity Model ............................................................................. 126
Appendix C: Numerical Dispersion ........................................................................................ 127
Appendix D: Diffusion Coefficients ....................................................................................... 128
D.1 Example (Binary Mixtures) ............................................................................................. 130
D.2 Example (Ternary Mixture) ............................................................................................. 131
Appendix E: Shape Factors for Gas-Oil Diffusion ................................................................. 133
Appendix F: Transfer Rate due to Gravity Drainage ............................................................. 137

xiii

List of Tables
Table 3.1 Pure component property data for the analog fluid system .......................................... 22
Table 3.2 Parameters for UNIQUAC viscosity model (interactions and structural parameters) .. 24
Table 3.3 Summary of displacement experiments ........................................................................ 31
Table 3.4 Model parameters obtained from vertical dispersivity experiment .............................. 39
Table 4.1 Initial and injection compositions and the MMP of the gas/oil pair used in example 1.
....................................................................................................................................................... 74
Table 4.2 Initial and injection compositions and the MMP of the gas/oil pair used in example 2.
....................................................................................................................................................... 79
Table 5.1 MPTA parameters from pure component isotherms. .................................................... 97
Table 5.2 MPTA parameters from pure component isotherms using composite potential function.
..................................................................................................................................................... 100
Table 5.3 Model parameters for Extended Langmuir Model (ELM) ......................................... 102
Table 5.4 Average absolute relative error (%) in prediction of binary and ternary excess sorption
..................................................................................................................................................... 104
Table A.1 Parameters for UNIQUAC model (interactions and structural parameters) .............. 125
Table D.1 Calculated diffusion coefficient compared with experimental data from Kett and
Anderson (1969) for a ternary liquid mixture at 298.15 K and 0.1 MPa. ................................... 132

xiv

List of Figures
Figure 1.1 Evolution of CO
2
injection projects and oil prices in the U.S. Data from Oil and Gas
Journal EOR surveys 1980-2010 and U.S. EIA 2010 (Alvarado and Manrique, 2010). ................ 2
Figure 1.2 Percentage of current EOR projects in the U.S. by EOR methods. Data from Oil and
Gas Journal (Koottungal, 2008). ..................................................................................................... 3
Figure 3.1 Quaternary phase diagrams (mass fractions). Top: CO2-CH4-nC4-C12 at 2280 psi
and 212°F as calculated from the PR EOS. Bottom: Water-MeOH-IPA-iC8 at 68°F and 14.7 psi
as calculated from the UNIQUAC model. .................................................................................... 21
Figure 3.2 Comparison of mixture densities (25°C and 1 atm) from Otero et al. (2000) with ideal
mixing calculations as a function of Water mass fraction along the Water-IPA-iC8 binodal curve.
....................................................................................................................................................... 23
Figure 3.3 Experimental and calculated viscosity for the binary mixtures (25°C and 1 atm). Data
from Tanaka et al., 1987 (MeOH-Water and Water-IPA), Ku, 2008 (IPA-iC8) and Soliman and
Marschall, 1990 (MeOH-IPA). ..................................................................................................... 25
Figure 3.4 Simplistic IFT model based on tie-line lengths (mass fractions), IFT in mN/m. Data
from Morrow et al., 1988 (Water-IPA-iC8) and Garcia-Flores et al., 2007 (Water-MeOH-iC8). 26
Figure 3.5 Interface between iC8 and MeOH in a PTFE capillary tube. While both pure
substances spread on PTFE, iC8 preferentially wets the surface in this binary setting. ............... 28
Figure 3.6 iC8 (wetting phase) - Water drainage relative permeability: Steady state observations
and Corey-type model. .................................................................................................................. 29
Figure 3.7 Effluent concentrations (in mass fractions) from FCM displacement experiments with
IPA and iC8: Experimental data (circles), single-porosity model (solid line), DP I (broken line),
and DP II (dotted line). Top: Vertical displacement. Bottom: Horizontal displacement. ............ 31
xv

Figure 3.8 Natural logarithm of equilibrium K-values and saturation of analog gas along
displacement length (z
D
=z/L) for experiment III at the dispersion-free limit after 0.7 PVI. ........ 32
Figure 3.9 Effluent concentrations (in mass fractions) for displacement experiment III from dual-
porosity model I (solid line), dual-porosity model II (dotted line), and experimental observations
(circles). ........................................................................................................................................ 33
Figure 3.10 Natural logarithm of equilibrium K-values for the 4 components and saturation of
analog gas along the displacement length (z
D
=z/L) for experiment IV at the dispersion-free limit
after 0.7 PVI. ................................................................................................................................. 34
Figure 3.11 Effluent concentrations (in mass fractions) for displacement experiment IV from
dual-porosity model I (solid line), dual-porosity model II (dotted line), and experimental
observations (circles). ................................................................................................................... 35
Figure 3.12 Natural logarithm of equilibrium K-values for the 4 components and saturation of
analog gas along the displacement length (zD=z/L) for experiment III with dispersion and mass
transfer after 0.7 PVI. .................................................................................................................... 42
Figure 3.13 Natural logarithm of equilibrium K-values for the 4 components and saturation of
analog gas along the displacement length (zD=z/L) for experiment IV with dispersion and mass
transfer after 0.7 PVI. .................................................................................................................... 43
Figure 3.14 Composition path (in mass fractions) for Exp. III from dispersion-free calculations
(broken line) and calculations with dispersive mixing and mass transfer (dotted line). ............... 44
Figure 3.15 Composition path (in mass fractions) for Exp. IV from dispersion-free calculations
(broken line) and calculations with dispersive mixing and mass transfer (dotted line). ............... 44
xvi

Figure 3.16 Component recovery for Exp. III from dispersion-free calculations (IPA: solid line,
iC8: broken line) and calculations with dispersive mixing and mass transfer (IPA: broken dotted
line, iC8: dotted line). ................................................................................................................... 45
Figure 3.17 Component recovery for Exp. IV from dispersion-free calculations (IPA: solid line,
iC8: broken line) and calculations with dispersive mixing and mass transfer (IPA: broken dotted
line, iC8: dotted line). ................................................................................................................... 46
Figure 3.18 Effluent concentrations (in mass fractions) for displacement experiment III from DP
I with IFT scaling of relative permeabilities (solid line), DP I without IFT scaling of relative
permeabilities (broken dotted line), DP I with linear relative permeabilities (dotted line) and
experimental observations (circles). ............................................................................................. 48
Figure 3.19 Effluent concentrations (in mass fractions) for displacement experiment IV from DP
I with IFT scaling of relative permeabilities (solid line), DP I without IFT scaling of relative
permeabilities (broken dotted line), DP I with linear relative permeabilities (dotted line) and
experimental observations (circles). ............................................................................................. 49
Figure 3.20 Relative impacts of longitudinal dispersion and transverse mass transfer for FCM
displacements. Top: Vertical displacement. Bottom: Horizontal displacement. .......................... 52
Figure 4.1 Schematic of two neighboring gridblocks containing oil and gas phases. .................. 63
Figure 4.2 Composition profiles near the interface during inter-phase mass transfer, after Taylor
and Krishna (1993). ...................................................................................................................... 66
Figure 4.3 Schematic of the reservoir and wellbore locations. White lines represent the fracture
system. .......................................................................................................................................... 73
Figure 4.4 Relative permeability curves for flowing (solid line) and stagnant (dashed line)
domains. Gas and oil relative permeabilities are shown with red and green colors respectively. 73
xvii

Figure 4.5 Liquid recovery (left) and producing GOR (right) versus time for example 1 when
diffusion is not considered (red), when diffusion is modeled using our approach (blue) and when
diffusion is modeled using the conventional approach (black). ................................................... 76
Figure 4.6 Liquid 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), when our
approach is used but PSS condition is assumed for M-F transfer function (dashed blue) and when
our approach is used but PSS condition is assumed for M-F transfer function and off-diagonal
diffusion coefficients are neglected (dotted blue). ........................................................................ 77
Figure 4.7 The CO
2
composition (%) in the stagnant domain after 10 years (example 1) for the
case without diffusion (left) and the case with diffusion based on calculations using our approach
(right). ........................................................................................................................................... 79
Figure 4.8 Liquid recovery (left) and producing GOR (right) versus time for example 2 when
diffusion is not considered (red), when diffusion is modeled using our approach (blue) and when
diffusion is modeled using the conventional approach (black). ................................................... 80
Figure 4.9 Liquid 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), when our
approach is used but PSS condition is assumed for M-F transfer function (dashed blue) and when
our approach is used but PSS condition is assumed for M-F transfer function and off-diagonal
diffusion coefficients are neglected (dotted blue). ........................................................................ 81
Figure 4.10 The CO
2
composition (%) in the stagnant domain after 10 years (example 2) for the
case without diffusion (left) and the case with diffusion based on calculations using our approach
(right). ........................................................................................................................................... 82
xviii

Figure 4.11 The saturation profile of a matrix block (dotted line) near the injection well and its
corresponding fracture (solid line) for example 1 as calculated using our approach. .................. 83
Figure 4.12 Composition profiles (mass fraction) of a matrix block (green) near the injection
well and its corresponding fracture (red) for example 1 during the first 10 years as calculated
using our approach. The blue curves show the two-phase boundaries in each ternary plane at 157
bar and 373.15 °K. The dashed brown line represents the dilution line. ...................................... 84
Figure 4.13 The saturation profile of a matrix block (dotted line) near the injection well and its
corresponding fracture (solid line) for example 1 as calculated using the conventional approach.
....................................................................................................................................................... 85
Figure 4.14 Composition profiles (mass fraction) of a matrix block (green) near the injection
well and its corresponding fracture (red) for example 1 during the first 10 years as calculated
using the conventional approach. The blue curves show the two-phase boundaries in each ternary
plane at 157 bar and 373.15 °K. The dashed brown line represents the dilution line. .................. 86
Figure 4.15 The CO
2
composition (%) in the stagnant domain after 10 years (example 1) as
calculated using the conventional approach (left) and calculations using our approach (right). .. 86
Figure 5.1 Pure component isotherms at 45 °C. Data from Ottiger et al. (2008b). Solid lines:
MPTA after regression. ................................................................................................................. 97
Figure 5.2 Experimental observations and MPTA model predictions for a ternary mixture of
CH4, CO2 and N2 at 45 °C and pressures from 5-188 bars. Excess sorption data from Ottiger et
al. (2008b). .................................................................................................................................... 98
Figure 5.3 Experimental observations and MPTA model predictions for a ternary mixture of
CH4, CO2 and N2 at 45 °C and pressures from 5-188 bars. Excess sorption data from Ottiger et
al. (2008b). .................................................................................................................................. 101
xix

Figure 5.4 Experimental observations and ELM predictions for a ternary mixture of CH4, CO2
and N2 at 45 °C and pressures from 5-188 bars. Excess sorption data from Ottiger et al. (2008b).
..................................................................................................................................................... 102
Figure D.1 Calculated diffusion coefficient (solid lines) compared with experimental data
(circles) from Sigmund (1976) for methane-propane mixtures at different values of pressure,
temperature and composition. ..................................................................................................... 130
Figure D.2 Calculated diffusion coefficient (solid lines) compared with experimental data
(circles) from Lo (1974) for binary liquid mixtures of alkanes at 298.15 K and 0.1 MPa with
varying compositions. ................................................................................................................. 131
Figure E.1 Schematic of a matrix block surrounded by fractures ............................................... 133

1

Chapter 1 : Introduction

1.1 Background
With the increase in energy demands, there is a need for increased oil/gas production in the
upcoming decades while alternative energy sources (e.g. solar, wind, geothermal, nuclear fusion,
biofuels) are developed and implemented. Since the oil and gas reserves are limited, improved
oil/gas recovery techniques are expected to play a crucial role in meeting the energy needs in the
years to come by increasing recovery from existing reserves.
Improved oil/gas recovery techniques include intelligent well monitoring and control,
secondary recovery methods and tertiary recovery processes. In secondary recovery methods
water (or gas) is injected into aquifer (or gas cap) to maintain the reservoir pressure, while
tertiary recovery processes encompass thermal methods (steam flooding and in-situ combustion),
gas injection (hydrocarbon, N
2
and CO
2
) and chemical flooding (polymer, alkaline and
surfactant).
Injection of CO
2
into oil/gas reservoirs, which is the main focus of this research project, has
gained tremendous industrial application over the last three decades. Figure 1.1 depicts evolution
of CO
2
injection projects in the U.S. over the last three decades based on the data from Oil and
Gas Journal EOR surveys 1980-2010 and U.S. EIA 2010 (Alvarado and Manrique, 2010). Even
though the oil price dropped during the 1980s and 2008-2010, the number of CO
2
injection
projects increased during these time intervals. This is a very interesting observation and shows
the economic feasibility of CO
2
injection projects even at low oil prices.
2

One of the main reasons for the increase in CO
2
injection projects in the U.S. is the cheap
CO
2
available from natural sources (Moritis, 2001). In addition, the miscible conditions for CO
2

injection processes are achieved at significantly lower pressures compared to natural gas or
Nitrogen. The solubility of CO
2
in oil causes the oil to swell and flow at lower viscosities. The
reduction in oil viscosity caused by CO
2
offsets, to some extent, the unfavorable mobility ratios
which result from low viscosity of CO
2
at reservoir conditions (Stalkup, 1983). All of these
factors have made CO
2
injection the most favorable EOR process for medium and light oil
reservoirs with API gravity greater than 25 (Taber, 1997; Alvarado and Manrique, 2010).

Figure ‎ 1.1 Evolution of CO
2
injection projects and oil prices in the U.S. Data from Oil and Gas
Journal EOR surveys 1980-2010 and U.S. EIA 2010 (Alvarado and Manrique, 2010).

In 2008, CO
2
and other gas injection processes accounted for more than 60% of the EOR
projects in the U.S. while approximately 30% of the EOR projects were based on thermal
methods (Koottungal, 2008). Most of these gas injection projects are performed on non-fractured
3

reservoirs. But there are huge fractured reservoirs in different parts of the world which may have
a potential for gas injection. Gas injection in fractured reservoirs is one of the focus areas of this
research project.

Figure ‎ 1.2 Percentage of current EOR projects in the U.S. by EOR methods. Data from Oil and
Gas Journal (Koottungal, 2008).

CO
2
injection is not limited to oil reservoirs. In recent years, coalbed methane has become an
important source of energy in the United States. Since primary production techniques typically
recover less than half of the methane in a coalbed, enhanced coalbed methane (ECBM) recovery
processes are used in which CO
2
and/or N
2
are injected into the coalbed to recover more CH
4
.
Climate change due to greenhouse gas emissions has become an important concern over the
last decade. With CO
2
being one of the main greenhouse gases in the earth’s atmosphere, CO
2

injection into depleted oil reservoirs and coalbeds provides additional benefit of sequestering
carbon in the subsurface.
4

1.2 Motivation
The initiation of an EOR project depends on the economic evaluation of the project which, in
turn, depends on the potential incremental hydrocarbon recovery from the project. In order to
estimate the potential incremental hydrocarbon recovery by gas injection, compositional
reservoir simulators are commonly used in the industry. In addition, successful design and
implementation of these processes rely, in part, on the accuracy by which the available
simulation tools can represent the physics that govern the displacement behavior in a reservoir.
Despite tremendous efforts and progress by different investigators over the past decades (see
Chapter 2), the interaction of flow and phase behavior in porous media and its proper modeling
in the context of CO
2
injection (and more generally miscible displacements) is not completely
understood due to the complex nature of these multiphase, multicomponent displacement
processes. This has motivated us to investigate the accuracy of some physical models that are
frequently used in compositional simulation of miscible displacements by comparing numerical
calculations with observations from miscible and near-miscible displacement experiments
(Chapter 3).
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 4).
5

The sorption of gases onto the coal surfaces is one of the main mechanisms that govern the
dynamics of ECBM recovery processes. Despite the well-documented complexity of
multicomponent sorption phenomena, adsorption and desorption of CH
4
/CO
2
/N
2
mixtures in the
porous coal is commonly modeled with the extended Langmuir model. The extended Langmuir
model has been proven unable to accurately describe the multicomponent sorption that is central
to ECBM recovery processes and, therefore, more sophisticated sorption models are needed.
This has motivated us to apply the potential theory to describe the multicomponent sorption of
relevance to ECBM processes (Chapter 5).
1.3 Objectives
In the first part of this research, the objective is to investigate the accuracy of key physical
models that are frequently used to describe and interpret dispersive mixing and mass transfer in
compositional simulation of miscible displacements with emphasis on CO
2
injection processes.
A quaternary analog fluid system (alcohol-water-hydrocarbon) that mimics the phase
behavior of CO
2
-hydrocarbon mixtures at high pressure and temperature has been designed in
our research group (Rastegar 2010). A porous medium was designed using
PolyTetraFlouroEthylene (PTFE) materials to ensure that the analog oil acts as the wetting
phase, and the properties of the porous medium were characterized in terms of porosity and
permeability. Relative permeability and interfacial tension (IFT) measurements were also
performed to delineate interactions between the fluid system and the porous medium.
Displacement experiments at First-contact miscible (FCM), near-miscible and multicontact
miscible (MCM) conditions were consequently performed (Rastegar 2010).
The effluent concentrations from two-component FCM displacement experiments exhibit a
tailing behavior that is attributed to imperfect sweep of the porous medium: A feature that is not
6

captured by normal dispersion models. To represent this behavior in displacement calculations,
we use dual-porosity (DP) models including mass transfer between flowing and stagnant
porosities. The 4-component two-phase displacement experiments (near-miscible and MCM) are
consequently simulated using the DP models constructed based on observations from FCM
displacements.
We demonstrate that the accuracy of our displacement calculations relative to the
experimental observations is sensitive to the selected models for dispersive mixing, mass transfer
between flowing and stagnant porosities, and IFT scaling of relative permeability functions. We
also demonstrate that numerical calculations substantially agree with the experimental
observations for some physical models with limited need for model parameter adjustment.
In the second part of this research, the objective is to formulate a dual-porosity model in
which the generalized Fick’s law is used for molecular diffusion to account for the dragging
effects (component interactions); and gas-oil diffusion at 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. A novel matrix-fracture transfer function is introduced for gas-oil diffusion based
on film theory. Diffusion coefficients are calculated using the Maxwell-Stefan model and are
pressure, temperature and composition dependent. A time-dependent transfer function is used for
matrix-fracture exchange in which the shape factor is adjusted using a boost factor to
differentiate between the transfer rate at early and late times.
Field-scale calculation examples are used to demonstrate 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. It is also demonstrated that the dragging effects (off-
diagonal diffusion coefficients) and time-dependency of matrix-fracture transfer function can
7

moderately impact the oil recovery during gas injection in fractured reservoirs. We also show
that miscibility is not developed in the matrix blocks even at pressures above minimum
miscibility pressure (MMP) when molecular diffusion is the main recovery mechanism during
gas injection in fractured reservoirs.
The objective, in the third part of this research, is to apply potential theory to describe the
multicomponent sorption of relevance to ECBM processes. In this approach for modeling
multicomponent sorption, each component is assumed to be affected by a characteristic potential
field emitted by the coal surface. We discuss the implementation of potential theory with
emphasis on the simulation of ECBM processes where computational efficiency and accuracy
must be balanced. The model must be solved by an iterative scheme, and is hence more
computationally expensive than the extended Langmuir approach.
The results and analysis presented in this work demonstrate that the application of potential
theory of sorption to modeling of ECBM recovery processes can improve the accuracy of
calculations. However, the additional complexity of the model and the associated increase in the
computational efforts may not balance the gain in accuracy sufficiently to warrant application in
general purpose reservoir simulation.
1.4 Manuscript Organization
The remainder of this manuscript is organized as follows. In Chapter 2, a brief literature review
on miscibility, mixing mechanisms, the effect of mixing on miscible displacements, mass
transfer in dual-porosity systems, and the mathematical modeling of these processes is presented.
Chapter 3 includes the experimental procedure and data together with mathematical modeling,
displacement calculations and interpretations of the experiments performed in our research group
to study mixing and mass transfer in multicontact miscible displacements.
8

In chapter 4, 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. A time-dependent transfer
function is presented that accounts for both early- and late-time behavior of matrix-fracture
interactions. 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.
In Chapter 5, we apply the potential theory to describe the multicomponent sorption of
relevance to ECBM processes. We discuss the implementation of potential theory with emphasis
on the simulation of ECBM processes where computational efficiency and accuracy must be
balanced. The calculations from potential theory are compared with experimental sorption data
of CH
4
, CO
2
, N
2
and their binary and ternary mixtures on real coal samples. Finally, the future
direction of this research project is presented in Chapter 6.

9

Chapter 2 : Literature Review

2.1 Miscibility
During gas/oil displacements in porous media, since the injected gas is not initially at
equilibrium with the oil, the contact(s) between the two phases results in mass transfer; and
therefore changes the properties of both oil and gas phases. If the pressure is high enough, the
displacing gas and displaced oil become similar; and the displacement process becomes highly
efficient. This condition is generally referred to as development of miscibility.
First-contact miscibility is a condition where two fluids form a single phase at their first
contact when they are mixed at any proportion at a given pressure and temperature. If two fluids
form a single phase after a few contacts, they are called to be multicontact miscible at the given
pressure and temperature. Multicontact miscibility can be achieved through three mechanisms:
vaporizing, condensing and combined vaporizing/condensing drives. In the vaporizing drive
mechanism, light and intermediate components transfer from the oil into gas phase (lean gas like
CH
4
, N
2
and CO
2
). In contrast, intermediate components transfer from the enriched gas phase
into the oil when condensing drive mechanism takes place (Stalkup, 1983). Combined drives
exhibit both features.
Ternary diagrams have been used traditionally to explain development of miscibility. In
ternary systems, a process is MCM by a vaporizing (or condensing) drive only if the initial (or
injection) composition lies outside the region of tie-line extensions. However, it has been shown
that a combination of vaporizing and condensing derives is responsible for development of
multicontact miscibility for reservoir fluids where more than three components are present (Zick,
10

1986 and Stalkup, 1987). Johns et al. (1993) showed that at least four components are needed to
describe combined condensing/vaporizing drive (quaternary diagrams).
The minimum miscibility pressure (MMP) is the pressure above which the injected gas and
initial oil become miscible by a multicontact mechanism. The MMP can be determined in the lab
by slimtube experiments in which gas is injected into a long coiled slimtube that has been packed
with sand and saturated with oil. The slimtube experiment is performed at increasing pressures;
and the resulting oil recoveries at 1.2 pore volumes injected (PVI) are plotted versus
displacement pressure. The break-point on the plot is considered as the MMP.
Different approaches have been used to estimate the MMP including single mixing cell,
multiple mixing cells, 1D slimtube simulation and the method of characteristics (MOC). A
description of these methods and comparisons between calculations with different methods and
experimental data can be found in Jessen (2000) and Abedini et al. (2014).
2.2 Mixing Mechanisms
Diffusion and dispersion are often used to represent mixing mechanisms during displacement
processes in porous media at scales where permeability heterogeneity is not represented
explicitly. Perkins and Johnston (1963) presented an extensive review of studies of diffusion and
dispersion in porous media. Diffusion takes place due to random motion of molecules when two
miscible fluids are placed in contact. While diffusion is driven by gradients in concentrations (or
more generally by gradients in the chemical potential), additional apparent mixing is caused by
local variations in the fluid velocity in a porous medium. The variations in fluid velocity and the
associated mixing are generally referred to as dispersion (Perkins and Johnston, 1963).
Diffusion will oftentimes not play a significant role in regions of a porous medium where
fluids velocities are high because the characteristic time for diffusion is relatively large.
11

However, in regions of the porous media where flow velocities are low (e.g. stagnant or
bypassed pore space) diffusion can dominate mass transfer over viscous flow. For simplicity, the
contribution of diffusion and dispersion are commonly lumped into an effective dispersion
coefficient (or dispersivity) as discussed later. Several parameters are known to influence the
effective mixing in a porous medium including: flow rate, heterogeneity, turbulence, viscosity
and density differences, grain size distribution, and grain shapes (Perkins and Johnston, 1963).
Coats and Smith (1964) observed asymmetric effluent concentrations from tracer
experiments on consolidated cores in contrast to nearly symmetrical effluent concentrations for
unconsolidated porous materials. They used a differential capacitance model to match their
experimental observations. This model accounts for dead-end pore space (stagnant porosity) in a
porous medium and includes mass transfer between flowing and stagnant pore space. Their
approach was subsequently modified by Baker (1977) in a study of dispersion and mass transfer
in FCM displacement experiments.
Bretz and Orr (1987) used the Coats-Smith model and a porous-sphere model to interpret
their experimental observations for FCM displacements in carbonate cores. The porous-sphere
model assumes an assemblage of porous spheres with internal diffusive mass transfer limitations.
They found that a porous-sphere model provided for a better match of their experimental
observations and explained this by the actual structure of the cores being more similar to a
porous-sphere model (from thin-section analysis).
Arya et al. (1988) presented a numerical study of the effect of aspect ratio, heterogeneity and
diffusion coefficients on the effective dispersive mixing at larger scales. They demonstrated that
dispersive mixing in heterogeneous systems with large aspect ratio behaves like equivalent
layered systems.
12

Peters et al. (1996) presented a method for measuring longitudinal dispersion coefficients in
porous media through computed tomography imaging while Schulze-Makuch (2005) reviewed
109 dispersion related studies and compiled values of dispersivity for different formations/scales.
Bijeljic and Blunt (2006) introduced a measure of transit times for particles between
neighboring pores to overcome the limitations of conventional time independent dispersion
coefficient to reproduce experimental observations at small scale.
Garmeh et al. (2009) used pore-scale simulation to study dispersion in porous media. The
results of their simulations align well with the classical work of Perkins and Johnston that dictate
a relationship between the pore-scale Peclet number and the observed effective longitudinal
dispersivity.
Despite a broad industrial applicability, less effort has been devoted to the study of
dispersion in multiphase systems. Sahimi et al. (1982) investigated dispersion in two-phase
immiscible displacements using Monte Carlo simulation combined with percolation theory to
determine phase distributions and mixing in artificial pore structures. They demonstrated that
dispersive mixing depends strongly on phase saturation and saturation history.
Delshad et al. (1985) performed multiphase experiments in sandpacks and Brea sandstone to
measure relative permeabilities and dispersive mixing. They generalized single phase dispersion
theory and applied it to interpret their experimental observations and arrived at dispersivities that
differed from the equivalent single phase dispersivities.
2.3 Effect of Mixing on Miscibility
It has long been known that numerical dispersion affects the efficiency in MCM displacements
(Gardner, 1981 and Stalkup, 1988). Walsh and Orr (1990) studied the effect of dispersion
(physical and numerical) on the composition path for one-dimensional flow of ternary, two-
13

phase mixtures. They demonstrated that the composition path passes through the two-phase
region as a result of numerical/physical dispersion and that oil recovery is reduced as the Peclet
number is decreased. They also demonstrated that the effect of dispersion on composition path,
and therefore recovery, depends on the shape and size of two-phase region on ternary diagram.
Johns et al. (1994) demonstrated that numerical dispersion causes some two-phase flow in a
4-component displacement that would be MCM in the absence of dispersion. They observed that
as numerical dispersion is increased, the oil recovery is reduced, especially near the minimum
miscibility enrichment (MME). Other investigators have also shown that the efficiency of
enhanced oil recovery by multicontact miscible gas injection depends, in part, on the level of
mixing that occurs between the injected gas and the oil in the reservoir (e.g. Solano et al., 2001;
Johns et al., 2002; Jessen et al., 2004).
2.4 Mass Transfer in Dual-Porosity (DP) 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 and demonstrated that diffusion can have a significant effect on oil
recovery away from the miscibility pressure.
14

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 CO
2
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 CO
2
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 CO
2
injection. They observed that diffusion played a more
significant role than gravity drainage in recovering oil from a long vertical tight matrix block
when CO
2
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
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.
15

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 CO
2
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.
2.5 Sorption in Coalbeds
In coalbeds, gas is present in two phases: a bulk gaseous phase that occupies the pore space and a
liquid-like adsorbed phase on the coal surfaces/pores. In primary recovery, the coalbed is
dewatered to reduce the overall reservoir pressure which causes CH
4
to be desorbed from coal
surfaces. Since primary production typically recovers less than half of the methane in a coalbed
(Stevens et al., 1998), enhanced coalbed methane (ECBM) recovery processes are needed in
which CO
2
and/or N
2
are injected into the coalbed to recover more CH
4
. Injection of CO
2
into
coalbeds also provides additional benefit of sequestering carbon in the subsurface.
16

Gas injection in ECBM recovery provides a method to maintain the overall coalbed pressure.
In addition, injecting a second gas, or a mixture of gases, decreases the partial pressure of CH
4
in
the free gas. As a result, desorption of CH
4
from coal surface is enhanced. The convective flow
of injected gas sweeps desorbed CH
4
through the coalbed. Therefore the sorption of gases onto
the coal surfaces is one of the main mechanisms that govern the dynamics of ECBM recovery
processes.
Despite the well-documented complexity of multicomponent sorption phenomena (Stevenson
et al., 1991; DeGance et al., 1993; Chaback et al., 1996), adsorption and desorption of
CH
4
/CO
2
/N
2
mixtures in ECBM recovery is usually modeled with the extended Langmuir model
because of its simplicity and associated low computational cost (Guo, 2003; Zhu et al., 2003;
Smith et al., 2005; Seto et al., 2009). The extended Langmuir has been proven unable to
accurately describe the multicomponent sorption that is relevant to ECBM recovery processes
(Clarkson, 2003; Wei et al., 2005). Jessen et al. (2008) demonstrated that extended Langmuir
was able to model the sorption process in binary displacements; but failed to describe the
behavior of ternary displacements.
2.6 Mathematical Modeling
Description of flow in a porous medium is typically started by writing the mass conservation
equation (also referred to as continuity equation). When convection, diffusion and dispersion are
the main physical mechanisms that govern the flow dynamics, the multicomponent multiphase
continuity equation for a single-porosity system can be written as

(2.1)
17

where x
ij
is the mole fraction of component i in phase j, is the porosity, S
j
is the saturation of
phase j, ρ
j
is the molar density of phase j,

is the velocity of phase j,

represents the effective
dispersion tensor (combined impact of dispersion and diffusion) for component i in phase j, t
represents time, n
p
represents the number of phases, and n
c
is the number of components.
It should be noted that the effects of chemical reaction, adsorption and temperature variation
have been neglected in derivation of Eq. (2.1) (Lake, 1989; Orr, 2007). The longitudinal (along
the main direction of flow) and transverse (perpendicular to the main direction of flow)
dispersion coefficients can be obtained from Eq. (2.2) and (2.3) respectively (Bear, 1972; Orr,
2007)

(2.2)

(2.3)
where α denotes a material constant known as dispersivity, subscripts l and t refer to longitudinal
and transverse directions, and D
ij
is the effective diffusion coefficient in the porous medium. We
note that the molecular diffusion coefficient is multiplied by

to estimate the effective
diffusion coefficient in the porous medium, where FF is the formation factor (Perkins and
Johnston, 1963).
The phase velocities in Eq. (2.1) are calculated from the extended Darcy’s law for multiphase
flow in porous media

,
(2.4)
18

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, and k represents the absolute permeability. The
relationship between phase pressures are represented by capillary pressure functions

(2.5)
The phases are assumed to be in thermodynamic equilibrium which can be represented by
equality of chemical potentials

(2.6)
It can be shown that the equality of chemical potentials leads to the equality of fugacities, and the
fugacities are calculated using the equations of state (EOS) for reservoir fluids (Danesh, 1998;
Orr, 2007). The phase saturations and mole fractions in each phase must sum to one

(2.7)

(2.8)
The equations presented in this section provide enough information to obtain the solution to a
flow and transport problem that models the effects of convection, diffusion, dispersion and phase
equilibrium (Orr, 2007). Note that the governing equations presented in this section are valid for
single-porosity systems. The governing equations for dual-porosity systems will be presented in
the subsequent chapters.

19

Chapter 3 : Experimental and Modeling Study of Multicontact Miscible
Displacements
1

3.1 Introduction
In this chapter, we investigate the accuracy of some physical models that are commonly used to
describe dispersive mixing and mass transfer in compositional reservoir simulation. We combine
experimental and numerical efforts to investigate the effect of mixing and mass transfer in near-
miscible and MCM displacements.
We use an analog fluid system that exhibit liquid-liquid phase behavior at ambient conditions
that is comparable to the vapor-liquid equilibrium behavior of mixtures of CO
2
and hydrocarbons
at high pressure and temperature. Two-component FCM, and 4-component near-miscible and
MCM displacement experiments were performed in a synthetic porous medium in our research
group to form the basis of our modeling study with emphasis on physical models for representing
mixing and mass-transfer in these relevant displacement processes (Rastegar, 2010).
We start by presenting the analog fluid system and the design/characterization of our porous
medium and solid-fluid interactions. We then present our observations from FCM, near-miscible
and MCM displacement experiments and provide an interpretation of the observations based on
numerical calculations with different physical models. A discussion and analysis of the presented
results and observations concludes the chapter.

1
Most of the results presented in this chapter have been published in Transport in Porous Media:
Shojaei, H., Rastegar, R., & Jessen, K.: Mixing and mass transfer in multicontact miscible displacements. Transport
in Porous Media 94, 837-857 (2012)
20

3.2 Fluid System
To select a suitable analog system for studying mixing and mass transfer in the context of CO
2
injection processes, the quaternary system of CO
2
, Methane (CH
4
), Butane (nC
4
) and Dodecane
(C
12
) at high pressure was initially examined. Figure 3.1 (top) shows the quaternary phase
diagram for this system at 2280 psi and 212°F as calculated from the Peng-Robinson (PR)
equation of state (EOS) (Peng and Robinson, 1976).
CO
2
represents the injected gas composition and a representative oil composition is shown
with a circle in the base of the quaternary phase diagram. At the given conditions, the
displacement of oil by CO
2
will be MCM as compositions that are formed during the
displacement process (assuming one-dimensional disperion-free flow) will pass through the
critical locus close to the front ternary of the quaternary diagram (Orr, 2007).
To facilitate experimental work at ambient conditions, analog solvents that mimic the phase
behavior of high pressure CO
2
/hydrocarbon systems were used. The quaternary system of Water,
Methanol (MeOH), Isopropanol (IPA), and Isooctane (iC
8
) was selected due to similarities with
the high-pressure CO
2
/hydrocarbon system (Rastegar, 2010). Pure component properties of the
analog fluids are summarized in Table 3.1.
Other authors have used water-alcohol-hydrocarbon systems to study MCM displacements in
porous media at ambient conditions (e.g. Batycky, 1994; Al-Wahaibi et al., 2007; and AlHamdan
et al., 2011). However, all previous studies have utilized simplified representations of the
associated phase behavior and transport properties that complicate the interpretation of the
displacement experiments.
21

Figure ‎ 3.1 Quaternary phase diagrams (mass fractions). Top: CO2-CH4-nC4-C12 at 2280 psi
and 212°F as calculated from the PR EOS. Bottom: Water-MeOH-IPA-iC8 at 68°F and 14.7 psi
as calculated from the UNIQUAC model.
22

Table ‎ 3.1 Pure component property data for the analog fluid system
Component M
w
(g/mol) ρ (g/cm
3
)* Viscosity (cP)* Boiling Point (
o
C)**
Water 18.01 0.997 0.890 99.9
MeOH 32.04 0.787 0.544 64.7
IPA 60.10 0.781 1.960 82.3
iC
8
114.23 0.687 0.473 99.3
*) density and viscosity at 25°C, **) Normal boiling point

Figure 3.1 (bottom) shows the phase diagram of the analog system at ambient conditions as
calculated by the Universal Quasichemical Activity Coefficient (UNIQUAC) model (see
Appendix A for more details). For the analog system, mixtures of Water and MeOH represent an
injection gas composition (CO
2
or CO
2
-rich gas) while compositions in the bottom ternary of the
phase diagram represent an initial oil.
3.2.1 Phase Behavior
To characterize and model the phase behavior of the analog fluid system, a range of ternary and
quaternary equilibrium mixtures was created and analyzed by gas chromatography (GC). The
UNIQUAC model, introduced by Abrams and Prausnitz (1975), was subsequently used to model
the phase behavior of the analog fluid system and a set of binary interaction coefficients was
extracted from the experimental observations. Additional details related to the phase behavior of
the analog fluid system are presented in Rastegar and Jessen (2011).
3.2.2 Mixture Densities
Mixture densities are commonly provided by the EOS model that is used in the equilibrium
calculations for hydrocarbon systems. As the UNIQUAC model does not provide us with
23

mixture density information, a separate density model is required. We assume that the excess
volume is negligible and calculate mixture densities based on the assumption of ideal mixing.
This assumption is supported by the experimental density data reported by Otero et al. (2000) for
mixtures of Water-IPA-iC
8
at 25°C.

Figure ‎ 3.2 Comparison of mixture densities (25°C and 1 atm) from Otero et al. (2000) with ideal
mixing calculations as a function of Water mass fraction along the Water-IPA-iC8 binodal curve.

Figure 3.2 compares the data from Otero et al. (2000) with ideal-mixing calculations along the
binodal curve of the Water-IPA-iC
8
ternary. The maximum error between calculated and
observed mixture densities is less than 2% and is observed for water-rich equilibrium mixtures.
Accordingly, we assume that the assumption of ideal mixing is sufficiently accurate for our
modeling work.
24

3.2.3 Mixture Viscosities
The viscosity of mixtures of the analog fluids is a highly nonlinear function of the composition
due to the strong nonideal behavior of water-alcohol systems. Simple viscosity correlations are
hence not suitable for this system. The UNIQUAC viscosity model, as proposed by Martin et al.
(2001), was selected and used throughout our modeling work (see Appendix B for more details).
The binary interaction coefficients of the UNIQUAC viscosity model were obtained from
viscosity data for binary mixtures of MeOH-Water (Tanaka et al. 1987), Water-IPA (Tanaka et
al. 1987), IPA-iC
8
(Ku 2008), and MeOH-IPA (Soliman and Marschall 1990). Experimental
observations are compared to model calculations for the relevant binary systems in Fig. 3.3 and
the model parameters are summarized in Table 3.2.

Table ‎ 3.2 Parameters for UNIQUAC viscosity model (interactions and structural parameters)
Comp Water MeOH IPA iC
8

Water 0 -276.03 -536.44 0.00 0.9945 0.28486
MeOH 423.54 0 398.30 0.00 1.4320 1.4311
IPA 999.08 -210.98 0 -155.96 2.2571 3.3915
iC
8
0.00 0.00 -8.32 0 5.0080 5.8463

The viscosity model was subsequently tested with viscosity data for 11 mixtures of MeOH-IPA-
Water (Soliman and Marschall, 1990) and a satisfactory agreement was observed with an
average error of less than 3%. We note that the interaction coefficient between MeOH and iC
8
is
set to zero due to lack of experimental data for this binary. This results in a monotonic and
almost linear compositional dependence of the viscosity for the associated binary system. We
25

believe that this is a reasonable representation for this polar-nonpolar binary as cross-association
will not contribute to a nonmonotonic compositional dependence as seen for example for the
IPA-Water binary.

Figure ‎ 3.3 Experimental and calculated viscosity for the binary mixtures (25°C and 1 atm). Data
from Tanaka et al., 1987 (MeOH-Water and Water-IPA), Ku, 2008 (IPA-iC8) and Soliman and
Marschall, 1990 (MeOH-IPA).

3.2.4 Interfacial Tension
Pendant drop measurements were performed for equilibrated ternary mixtures of Water, IPA and
iC
8
to validate the experimental observations of Morrow et al. (1988) and our observations were
in good agreement with reported values (Rastegar, 2010). Initial attempts to correlate the data
from Morrow et al. (1998) including the Parachor method and the more complex model proposed
by Bahramian and Danesh (2004) failed in predicting the IFT with reasonable accuracy. To
26

include IFT effects in our modeling efforts, we use a simple correlation derived from the data of
Morrow et al. (1988). The proposed correlation has the form
b, (3.1)
where x denotes the Euclidian norm of a tie line based on mass fractions. A comparison of the
observed and calculated IFT (along with the coefficients of Eq. 3.1) is provided in Fig. 3.4. The
correlation was subsequently validated for a mixture of Water-MeOH-iC
8
(circle in Fig. 3.4)
from Garcia-Flores et al. (2007), and was found to be in reasonable agreement. We note that the
correlation predicts a value of IFT equal to 4e-4 mN/m for critical mixtures which is sufficiently
accurate for our modeling purposes. However, additional work is warranted to develop more
rigorous IFT models for these complex mixtures.

Figure ‎ 3.4 Simplistic IFT model based on tie-line lengths (mass fractions), IFT in mN/m. Data
from Morrow et al., 1988 (Water-IPA-iC8) and Garcia-Flores et al., 2007 (Water-MeOH-iC8).

27

3.3 Packed Column
A porous material was selected to allow for the hydrocarbon-rich phase (oil compositions) to act
as the wetting phase while water-rich mixtures (injection gas compositions) act as the non-
wetting phase. To achieve this, PTFE materials were used to design the packed column. PTFE
powder from Sigma-Aldrich with a 100 micron average particle diameter was used as the
packing material. The column was designed with an inner diameter of 0.375 inches and length of
12 inches. This diameter is sufficiently large to minimize wall effects in the packed column.
The porosity of the column was measured gravimetrically to 47% by saturating the column
with iC
8
under moderate vacuum. The permeability of the packed column was calculated to 0.39
Darcy from Darcy’s law using the stabilized pressure drop across the column when iC
8
was
injected at a constant rate (Rastegar, 2010).
3.4 Solid-Fluid Interactions
3.4.1 Wettability
Additional considerations were made regarding the wettability of the individual components of
the analog fluid system. From the pure components, Water is the only component that does not
spread on a PTFE surface. Mixtures of Water and MeOH at sufficiently high MeOH
concentrations (above 65% by mass) also spread on a PTFE surface and hence, the use of Water-
MeOH mixtures to mimic injection gas compositions was questionable.
With iC
8
as the primary component in the analog oil phase, an additional experiment with a
PTFE capillary tube was performed to identify the preferential wetting component of the MeOH-
iC
8
pair. Figure 3.5 shows the interface between MeOH and iC
8
at equilibrium in a PTFE
capillary tube that confirms that iC
8
will preferentially wet the surface and that MeOH will act as
28

a nonwetting phase in displacements where iC
8
is present in the porous medium (Rastegar,
2010).
3.4.2 Relative Permeability
Steady state relative permeability experiments were performed for the immiscible pair of iC
8
and
Water to establish the maximum residual oil saturation in the porous medium at the maximum
attainable IFT of 39 mN/m (longest tie-line in Fig. 3.1). Figure 3.6 reports the relative
permeabilities of the porous medium to Water and iC
8
phases as a function of the iC
8
saturation.
From the experiments, a residual oil saturation of 29% is observed (Rastegar 2010).
To model the observed relative permeability data, we use Corey-type relative permeability
functions with saturation exponents for wetting and non-wetting phases of 1.5 and 3,
respectively.

Figure ‎ 3.5 Interface between iC8 and MeOH in a PTFE capillary tube. While both pure
substances spread on PTFE, iC8 preferentially wets the surface in this binary setting.
29

Figure ‎ 3.6 iC8 (wetting phase) - Water drainage relative permeability: Steady state observations
and Corey-type model.

3.5 Displacement Experiments
A series of displacement experiments was performed with the analog fluid system and the PTFE
column. In all experiments, 3 droplets of effluent (less than 1% of the pore volume) were
diverted to a clean sample vial, every 5-10 minutes, diluted with Ethanol and analyzed for
composition by gas chromatography. This approach was used to minimize the additional
apparent mixing during sampling.
Displacement experiments were, unless otherwise noted, performed vertically to promote
gravity stable fronts and to avoid phase segregation (override/underride) due to gravity: A more-
dense and more-viscous phase was injected to displace a less-dense and less-viscous phase from
the bottom to the top of the column. According to Lake (1989), the upward vertical
displacements performed in this work are unconditionally stable since the end-point mobility
ratios are less than one and the injection phases are denser than the initial phases.
30

3.5.1 FCM Displacements
Two FCM displacement tests (experiments I and II) were performed to study the effective
mixing in the porous material. In the first displacement experiment, iC
8
was displaced

by IPA in
the vertical direction (bottom to top) at 23°C at a volumetric flow rate of 0.05 ml/min. In a
second FCM experiment, IPA was used to displace iC
8
in the horizontal direction also at a rate of
0.05 ml/min (Rastegar, 2010). Based on the analysis of the scaling groups suggested by Zhou et
al. (1997), the horizontal FCM test (exp. II) is not dominated by gravity forces. A 2D simulation
with the conditions in experiment II verifies that gravity underride does not happen in this
experiment.
A flow rate of 0.05 ml/min corresponds to an average linear velocity of 3.3 ft/day. The
effluent concentrations from these FCM displacements are reported in Fig. 3.7 (top and bottom
figures are for vertical and horizontal displacements, respectively) as a function of pore volumes
injected (PVI). We note that both FCM displacement experiments show a distinct tailing
behavior and defer additional discussion/interpretation to the modeling section below.
3.5.2 Near-miscible and MCM Displacements
Two 4-component near-miscible and MCM displacements were subsequently performed in the
vertical direction (bottom to top). For these displacement experiments, the initial and injection
compositions were selected to provide for a near-miscible and a MCM displacement in the
absence of dispersion. The initial and injected compositions are reported in Table 3.3 (Rastegar,
2010).

31

Figure ‎ 3.7 Effluent concentrations (in mass fractions) from FCM displacement experiments with
IPA and iC8: Experimental data (circles), single-porosity model (solid line), DP I (broken line),
and DP II (dotted line). Top: Vertical displacement. Bottom: Horizontal displacement.

Table ‎ 3.3 Summary of displacement experiments
Initial composition by mole Injection composition by mole
Run T (°C) q (ml/min) Water MeOH IPA iC8 Water MeOH IPA iC8
I , II 23 0.05 0.00 0.00 0.00 1.00 0.00 0.00 1.00 0.00
III 19 0.05 0.00 0.00 0.25 0.75 0.25 0.75 0.00 0.00
IV 21 0.05 0.08 0.00 0.29 0.63 0.08 0.90 0.02 0.00

32

In the first 4-component displacement experiment (Exp. III), an initial oil composition of 75%
iC
8
and 25% IPA (by mole) was displaced by 75% MeOH and 25% water (by mole). This
corresponds to a near-miscible displacement as demonstrated by a numerical calculation in 1D
with 10,000 cells in the absence of physical dispersion/diffusion (hyperbolic form). Figure 3.8
reports the variation in equilibrium K-values (K
i
= y
i
/x
i
) and analog gas saturation along the
displacement length after 0.7 PVI. We observe that the trailing vaporization shock propagates at
a specific velocity that is less than unity and that near-miscible conditions exist upstream of the
leading edge of the displacement.

Figure ‎ 3.8 Natural logarithm of equilibrium K-values and saturation of analog gas along
displacement length (z
D
=z/L) for experiment III at the dispersion-free limit after 0.7 PVI.

Figure 3.9 reports the effluent concentration of the four components in experiment III as a
function of PVI. At the dispersion-free limit for MCM displacements, all compositions would
propagate at a unit characteristic velocity through the porous medium. However, the impacts of
mixing and mass transfer in the experiment at near-miscible conditions cause the characteristic
33

velocities of the various compositions to depart from unity. In addition, we observe that a bank
of IPA is formed at the leading edge of the displacement process. This is, in part, due to the
curvature of the two-phase boundary in the quaternary diagram that forces the compositional
path of the displacement process to pass through a region of higher IPA concentrations relative
to both initial and injection compositions (Orr, 2007).

Figure ‎ 3.9 Effluent concentrations (in mass fractions) for displacement experiment III from
dual-porosity model I (solid line), dual-porosity model II (dotted line), and experimental
observations (circles).

In a second 4-component displacement experiment (Exp. IV), a small amount of Water was
added to the oil composition to place it on the Water–IPA–iC
8
ternary while the injected
PVI
34

composition was enriched with IPA (see Table 3.3). The enrichment of the injected composition
was performed to further promote the development of multicontact miscibility. To confirm the
development of miscibility, a numerical calculation in 1D with 10,000 cells that represents the
dispersion-free limit was performed. Figure 3.10 reports the variation in the equilibrium K-
values and saturation along the displacement length after 0.7 PVI. The development of
miscibility by a vaporizing/condensing mechanism is evident from the hour-glass shape observed
in the top panel of Fig. 3.10.

Figure ‎ 3.10 Natural logarithm of equilibrium K-values for the 4 components and saturation of
analog gas along the displacement length (z
D
=z/L) for experiment IV at the dispersion-free limit
after 0.7 PVI.

The observed effluent concentrations for the four components in the MCM displacement
experiment (IV) are reported in Fig. 3.11. As for experiment III, we observe a departure from
unit characteristic velocity for all components. Again, this is attributed to the mixing and mass
transfer that influences the displacement characteristics relative to a hypothetical dispersion-free
limit. In this experiment, we observe that banks of both Water and IPA are formed at the leading
35

edge of the displacement: Water and IPA are both more soluble in MeOH than in iC
8
and are
hence gradually stripped from the initial oil in the system to form these banks.

Figure ‎ 3.11 Effluent concentrations (in mass fractions) for displacement experiment IV from
dual-porosity model I (solid line), dual-porosity model II (dotted line), and experimental
observations (circles).

36

3.6 Mathematical Modeling
We observe a tailing effluent concentration in both FCM displacement experiments (see Fig. 3.7)
that signifies an imperfect sweep of the column. This, in turn, suggests that a 1D single-porosity
(SP) model will be insufficient to capture the dynamics of the displacement experiments. A 1D
dual-porosity (DP) model with mass transfer between flowing and stagnant porosities (Coats and
Smith, 1964), or a two (or three)-dimensional single-porosity model with explicit representation
of heterogeneity can be used to capture the tailing behavior observed in dispersivity experiments.
Since explicit information regarding the heterogeneity of the porous material is not available,
we use the Coats-Smith model to interpret our displacement experiments. Advection, diffusion,
dispersion, and gravity effects are the main physical mechanisms that control the displacement
dynamics in the presented experiments. To represent the relevant physics, we write the 1D
continuity equation for multicomponent two-phase flow in the packed column as

(3.2)
where x
ij
is the mole fraction of component i in phase j, is the porosity, S
j
is the saturation of
phase j, ρ
j
is the molar density of phase j,

is the velocity of phase j,

represents the effective
dispersion (combined impact of dispersion and diffusion) for component i in phase j, f denotes
the flowing fraction, t represents time, denotes the direction along the packed column, n
p

represents the number of phases, and n
c
is the number of components. The mass transfer between
the flowing and the stagnant porosities (q
i
) is represented by a simple linear mass-transfer model

(3.3)
37

where θ
i
is the mass transfer coefficient of component i. Superscript * denotes the quantities of
mixtures in the stagnant porosity. Phase velocities are calculated from Darcy’s law and we
assume that capillary pressure is negligible for our displacement experiments.
The effective dispersion coefficients in Eq. 3.2 are assumed to be identical for all
components within a phase (by assuming that molecular diffusion is negligible compared to
velocity dependent dispersion) and are calculated from the phase saturations and velocities as

(3.4)
where is the effective longitudinal dispersivity of the packed column. This simplification is
introduced to reduce the number of model parameters.
The phases are assumed to be in thermodynamic equilibrium and, as discussed in section
3.2.1, we use the UNIQUAC model to predict the phase behavior of the quaternary system used
in our experiments. The UNIQUAC viscosity model (see section 3.2.3) is used for viscosity
calculations while density calculations are based on ideal mixing (see section 3.2.2). To calculate
IFT, we apply the simple correlation presented in section 3.2.4 while Corey-type relative
permeability functions (see section 3.4.2) with IFT scaling are used in our displacement
calculations to account for the effects of IFT on relative permeabilities (see Coats, 1980;
Amaefule and Handy, 1982). A scaling factor, F, in the range of 0 to 1 that indicates how close a
mixture is to a critical point (miscible condition) is calculated from

(3.5)
where σ is the mixture IFT, σ
0
is a reference IFT, and β represent an adjustable scaling exponent.
The scaling factor interpolates between straight lines and reference (σ
0
) relative permeability.
According to Coats (1980), the scaling exponent β is generally in the range of 0.1 to 0.25. Since
38

our relative permeability experiments were performed at the maximum attainable IFT of 39
mN/m, we use this value as the reference IFT in Eq. 3.5. Residual oil saturations and critical gas
saturations are then modified by the scaling factor F as (Coats, 1980)

(3.6)
Normalized oil and gas saturations are subsequently calculated from

(3.7)
and Oil and gas relative permeabilities are finally calculated from

(3.8)
where n
o
and n
g
are Corey exponents (see section 3.4.2) for oil and gas respectively. Throughout
our displacement calculations, we assume that S
gc
is zero.
An explicit finite volume formulation was used to simulate our displacement experiments. All
simulations were performed in one dimension (as discussed above) with 1,000 grid blocks to
minimize the effects of numerical diffusion.
3.7 Displacement Calculations
Two representations of a dual-porosity porous medium were considered in this work. In the first
model (DP I), the stagnant porosity is evenly distributed along displacement length, while in the
second model (DP II), the stagnant (bypassed) porosity is located at the inlet and outlet of the
packed column. This representation attempts to mimic imperfect distribution of the fluids at the
inlet and outlet manifold.

39

3.7.1 FCM Displacements
Three model parameters were obtained by matching the numerical calculations with DP I to the
effluent concentrations from the vertical dispersivity experiment (Exp. I). The model parameters,
reported in Table 3.4, include the flowing fraction (f = 0.965), the effective longitudinal
dispersivity (α = 0.0015 m) and the mass transfer coefficient (θ = 0.0003 sec
-1
). The agreement
between experimental observations and numerical calculations is reported in the top panel of Fig.
3.7. We note that a dispersivity of 0.0015 m corresponds to a Peclet number of 203 that indicates
a moderate level of dispersive mixing in the porous medium.

Table ‎ 3.4 Model parameters obtained from vertical dispersivity experiment
Flowing fraction: f Dispersivity: α (m) Mass transfer: θ (1/sec)
DP I 0.965 0.0015 0.0003
DP II 0.650 0.0015 0.0050

In DP II, we allocate 5% of the column length at the inlet and outlet of the system (10% of the
total column length) to represent the stagnant porosity. This representation attempts to mimic
imperfect distribution of the fluids at the inlet and outlet manifolds. The flowing fraction at the
inlet and outlet of the packed column is calculated to 0.65 by assuming that the overall flowing
porosity is the same for DP I and DP II.
We apply the value of dispersivity obtained from DP I in our displacement calculations with
DP II and match DP II model to the experimental data from vertical dispersivity test by adjusting
the mass transfer coefficient (θ). Since the characteristic length for transverse mass transfer is
larger for DP II (as compared to DP I), we expect a larger value of mass transfer coefficient (θ).
40

By matching DP II to the effluent concentrations from the vertical dispersivity experiment, a
value of θ=0.005 sec
-1
is obtained that is consistent with our expectations (see Table 3.4).
Figure 3.7 (top) compares the calculated effluent concentrations from SP and DP models (I
and II) with experimental effluent data from the vertical dispersivity experiment. From Fig. 3.7
(top) we observe that the tailing behavior is captured accurately by the DP models for vertical
FCM experiment. In order to validate the matches obtained for DP models using the vertical
dispersivity test, we predict the effluent concentrations for the horizontal dispersivity test and
compare the predictions with the experimental data in Fig. 3.7 (bottom). From Fig. 3.7 (bottom)
we observe that both DP models are able to accurately predict the effluent concentrations for the
horizontal dispersivity test.
We note also that there is no significant difference between the calculated effluent
compositions from DP I and DP II. For comparison, Fig. 3.7 reports the calculations with a SP
model and we observe, as expected, that the tailing behavior of the dispersivity experiments
cannot be reproduced. We note that the value of dispersivity as obtained from DP I (α = 0.0015
m) was used with the SP model.
3.7.2 Near-miscible and MCM Displacements
Next, we consider the near-miscible (Exp. III) and MCM (Exp. IV) displacement experiments.
The model parameters from the dispersivity experiments (see Table 3.4) are used to simulate the
four-component displacement experiments. The value of θ, as estimated from the dispersivity
experiments, was used as the effective mass transfer coefficient for all components in the four-
component displacement calculations (see Eq. 3.3).
From our modeling efforts, a value of β equal to 0.15 (see Eq. 3.5) was found to provide the
best agreement between experimental and calculated effluent concentrations for experiments III
41

and IV. We note that the IFT scaling exponent (β=0.15) is the only parameter that was adjusted
to match the four-component displacement experiments. The additional three model parameters
obtained from the dispersivity experiments were used directly as input parameters to the
displacement calculations without further adjustments.
Figure 3.9 compares calculations with DP I and DP II to the experimental effluent
concentrations from experiment III. Despite the very limited parameter adjustment, we find a
reasonable agreement between calculated and measured effluent concentrations for all four
components. A minor mismatch is observed for the Water effluent concentrations at later times
and we attribute this to experimental errors introduced in the sampling and analysis of the
effluent concentrations. We also observe that the numerical calculations with DP I and DP II are
close to identical.
Numerical calculations are compared to the experimental effluent concentrations for
experiment IV in Fig. 3.11. The main features of the effluent concentrations are captured well by
the calculations. The moderate mismatch observed for Water and IPA effluent concentrations is
again attributed to experimental errors associated, in part, with sampling and subsequent
compositional analysis. We observe, as for experiment III, that DP I and DP II provide close to
identical effluent concentrations.
3.8 Effect of Mixing on Near-miscible and MCM Displacements
For MCM displacements it has been shown that dispersive mixing (or alternatively mass transfer
between flowing and stagnant regions) causes the composition paths to intersect the two-phase
region (e.g. Walsh and Orr, 1990; Johns et al., 1994).
Our simulations confirm that two-phase flow occur in both experiments III and IV. Figure
3.12 (3.13) reports the K-values and gas saturation along the displacement length (in flowing
42

porosity) when dispersive mixing and mass transfer are included in DP I calculations using 1,000
grid blocks for experiment III (IV). By comparing Fig. 3.12 (3.13) with Fig. 3.8 (3.10) we
observe that dispersive mixing and mass transfer between flowing and stagnant porosities cause
the displacement to depart from MCM conditions and promote two-phase flow inside the porous
medium.

Figure ‎ 3.12 Natural logarithm of equilibrium K-values for the 4 components and saturation of
analog gas along the displacement length (zD=z/L) for experiment III with dispersion and mass
transfer after 0.7 PVI.

43

Figure ‎ 3.13 Natural logarithm of equilibrium K-values for the 4 components and saturation of
analog gas along the displacement length (zD=z/L) for experiment IV with dispersion and mass
transfer after 0.7 PVI.

Figure 3.14 (3.15) compares the composition path for the dispersion-free calculations (SP
simulation without physical dispersion using 10,000 cells) with the calculations with dispersive
mixing and mass transfer (simulation with DP I using 1,000 cells) for experiment III (IV). As
observed from Figures 3.14 and 3.15, dispersive mixing and mass transfer between flowing and
stagnant porosities cause the composition path to move inside the two-phase region for both
experiments III and IV.
The compoent recoveries for the dispersion-free calculations and the calculations with
dispersive mixing and mass transfer are compared for experiments III and IV in Figures 3.16 and
3.17 respectively. Since dispersion and mass transfer cause the displacements to depart from the
MCM conditions, the analog oil recovery is reduced when these mechanisms are present.
44

Figure ‎ 3.14 Composition path (in mass fractions) for Exp. III from dispersion-free calculations
(broken line) and calculations with dispersive mixing and mass transfer (dotted line).

Figure ‎ 3.15 Composition path (in mass fractions) for Exp. IV from dispersion-free calculations
(broken line) and calculations with dispersive mixing and mass transfer (dotted line).
45

Figure ‎ 3.16 Component recovery for Exp. III from dispersion-free calculations (IPA: solid line,
iC8: broken line) and calculations with dispersive mixing and mass transfer (IPA: broken dotted
line, iC8: dotted line).

3.9 Importance of Selected Physical Models
To investigate the impact of the selected dispersion model as well as the impact of IFT scaling of
relative permeabilities on the calculated effluent concentrations, additional displacement
calculations were performed, and the results are presented in this section.
First we investigate the formulation of the dispersive mixing and replace the phase-velocity
dependent model (see Eqs. 3.2 and 3.4) with a model that is based on the total fluid velocity. The
effective dispersion coefficients for the latter model are calculated from the phase saturations and
total velocity as
46

Figure ‎ 3.17 Component recovery for Exp. IV from dispersion-free calculations (IPA: solid line,
iC8: broken line) and calculations with dispersive mixing and mass transfer (IPA: broken dotted
line, iC8: dotted line).

(3.9)
where

is the total linear velocity (summation of phase velocities). We note that when the
dispersion model is based on the total fluid velocity (Eq. 3.9), the continuity equation (Eq. 3.2) is
rewritten as

(3.10)
47

We found that the use of total velocity in the dispersion model does not have a significant effect
on the simulation results or agreement with the experimental observations for both experiments
III and IV.
Experiment III is a near-miscible displacement and hence some two-phase flow is expected
inside the column. In addition, as discussed in the previous section, dispersive mixing and mass
transfer between flowing and stagnant porosities promote two-phase flow in both experiments III
and IV. Hence the observation that the use of total or phase velocities in the dispersion model
lead to substantially identical results is attributed to the moderate level of mixing observed in
these experiments (α=0.0015 m). This argument was subsequently validated with additional
simulations of experiments III and IV using larger values of dispersivity (smaller Pe numbers)
exhibiting a more notable effect of the selected dispersion model.
Next, we probe the sensitivity of our displacement calculations to the IFT scaling of relative
permeability functions for experiment III. Figure 3.18 compares the experimental observations
with 3 different numerical calculations (DP I) with values of F equal to 0.0, 1.0 and as calculated
from Eq. 3.5 with β equal to 0.15. F = 0 corresponds to calculations with straight-line relative
permeability functions (valid for fully miscible condition) while F = 1 defaults the relative
permeability to the steady-state measurements reported in Fig. 3.6 (valid for immiscible
condition). From Fig. 3.18, we observe that the calculated effluent is very sensitive to the scaling
of the relative permeabilities and that the IFT scaling of relative permeabilities is required to
match the experimental observations. The observed sensitivity can be explained by the variations
in the interfacial tension due to composition changes that arise from the interaction of flow and
phase behavior inside the column.
48

Similar observations regarding the sensitivity of the IFT scaling were obtained for
displacement experiment IV as shown in Fig. 3.19. We note that additional steady-state relative
permeability measurements can be performed with pre-equilibrated phase compositions at an IFT
different from 39 mN/m to determine the scaling exponent (β) directly. This would allow for
predictive calculations of the multicomponent displacements without any need for parameter
adjustments.

Figure ‎ 3.18 Effluent concentrations (in mass fractions) for displacement experiment III from DP
I with IFT scaling of relative permeabilities (solid line), DP I without IFT scaling of relative
permeabilities (broken dotted line), DP I with linear relative permeabilities (dotted line) and
experimental observations (circles).
49

Figure ‎ 3.19 Effluent concentrations (in mass fractions) for displacement experiment IV from DP
I with IFT scaling of relative permeabilities (solid line), DP I without IFT scaling of relative
permeabilities (broken dotted line), DP I with linear relative permeabilities (dotted line) and
experimental observations (circles).
3.10 Discussion
In the previous sections, we presented the experimental and modeling aspects of our work on
mixing and mass transfer in near-miscible and MCM displacement processes. Both dispersivity
(FCM) experiments display a tailing behavior in the effluent concentrations that is attributed to
imperfect sweep in the column. A 1D SP model was not able to capture this behavior and,
50

therefore, a 1D DP model with mass transfer between flowing and stagnant porosities (Coats-
Smith model) was applied to represent the tailing behavior observed in dispersivity experiments.
From our interpretation of the dispersivity experiments, we find that dispersive mixing is
moderate in the porous medium with a Peclet number of approximately 200. For comparison, the
value of longitudinal dispersivity can be estimated using the correlation proposed by Perkins and
Johnston (1963). They reported the following correlation to estimate the dispersivity in
unconsolidated porous media

(3.11)
where

,

, FF, , , and

denote longitudinal dispersion coefficient, molecular diffusion,
formation factor, porosity, interstitial velocity, and particle diameter, respectively. Using a
typical value for molecular diffusion in liquids (D
o
=1e-5 cm
2
/sec), a formation factor for
unconsolidated porous media (FF=2.5), and additional parameters as described in previous
sections, we arrive at a value of α = 2.1e-4 m from the correlation. This is slightly lower than the
value obtained from matching our experimental observations (α = 1.5e-3 m): According to the
summary of experimental observations presented by Lake (1989), the value of dispersivity from
laboratory-scale experiments are typically in the range of 1e-4 to 1e-2 m. Thus the agreement
between the values of dispersivity obtained from our experiments and the correlation suggested
by Perkins and Johnston (1963) is reasonable.
We believe that the mass transfer between flowing and stagnant porosities in our experiments
is diffusive. To verify this, we compare the characteristic time for diffusion in a porous medium
to the characteristic time for mass transfer between flowing and stagnant porosities estimated
from our experiments. The characteristic time for diffusion can be estimated as
51

(3.12)
where L is the characteristic length and D is the diffusion coefficient in the porous medium.
Assuming that the characteristic length for mass transfer is 10-20% of the column diameter
combined with typical values of D
o
=1e-5 cm
2
/sec and FF=2.5, we estimate a characteristic time
of 959-3837 sec for diffusion in the porous medium. The mass transfer coefficient obtained by
matching DP I to the experimental data from dispersivity tests is θ = 3e-4 sec
-1
; thus the
characteristic time for mass transfer between flowing and stagnant porosities is approximately

which is in the same range as the estimated characteristic time for diffusion in
the porous medium. It is therefore reasonable to assume that the mass transfer between flowing
and stagnant porosities in DP I model is diffusive.
To investigate the relative impacts of longitudinal dispersion (along streamlines) and
transverse mass transfer (between stagnant and flowing regions), we compare simulations with
the individual mechanisms turned on and compare the results with the experimental data for the
FCM displacements in Fig. 3.20.
From Fig. 3.20, we observe that the transverse mass transfer has a more significant impact on
the agreement between calculated and experimental effluent concentrations for the vertical FCM
displacement. However, for the horizontal FCM displacement, an explicit representation of both
the longitudinal dispersion and the transverse mass transfer is important. The observation that
longitudinal dispersion has a smaller impact on the effluent concentrations for the vertical
experiment is attributed to the gravity stabilization that reduces the local variations in fluid
velocity.
52

Figure ‎ 3.20 Relative impacts of longitudinal dispersion and transverse mass transfer for FCM
displacements. Top: Vertical displacement. Bottom: Horizontal displacement.

To further establish that the level of numerical diffusion is small compared to the physical
dispersion, we approximate numerical Peclet numbers for the relevant calculations (see
Appendix C for details on estimating numerical diffusion in single- and two-phase flow). The
numerical Peclet number is estimated in the range of 2200 to 3000 for experiments I through IV
for simulations with 1,000 grid blocks. Accordingly, numerical Peclet numbers are one order of
magnitude larger that the physical Peclet number of approximately 200 estimated from our
experiments, and we can assume that the impact of numerical artifacts is marginal.
From the displacement experiments, we observe that DP I and DP II provide for almost
identical results. We note that in order to match the experimental data, separate parameter
53

adjustments were performed for DP I and DP II using the effluent concentrations from the
vertical dispersivity experiment.
The Coats-Smith model used in our interpretations is a 1D model that attempts to accounts
for some unresolved 2D/3D effects by means of a mass transfer term. The two representations of
the dual-porosity system (DP I and II) used in this work attempts to include additional spatial
information regarding the location of stagnant porosity in our application of the Coats-Smith
model. However, since both DP I and DP II are able to reproduce the experimental observations,
we cannot argue which model (DP I or DP II) provides for a better representation of the
unresolved heterogeneity. This, in turn, points out the need for experiments where the flowing
and stagnant regions are well defined by e.g. CT imaging or by careful design of heterogeneity in
bead packs.
We observe that the calculated effluent concentrations are insensitive to the velocity
dependence of the effective dispersion coefficient (phase velocities or total velocity). This is
despite the fact that two phases are present inside the column due to near-miscible conditions
(for experiment III) and dispersive mixing and mass transfer between flowing and stagnant
porosities (for experiments III and IV). We attribute this insensitivity to the moderate level of
mixing in our experiments. At larger scale, where the dispersivity is expected to increase, the
impact of the dispersion model is likely to be more significant.
IFT scaling of relative permeabilities was demonstrated to be of great importance in our
modeling efforts. This importance is caused by the variations in the IFT that arise from
interactions between flow and phase behavior during the displacements. This is in contrast to the
theory of gas injection processes at the dispersion-free limit (Orr, 2007) where development of
miscibility is independent of relative permeability. For MCM displacements in the absence of
54

dispersion, two-phase flow will not occur. However, the effect of dispersive mixing and mass
transfer between flowing and stagnant porosities in our experiments results in some two-phase
flow during both 4-component displacements.
We show that if IFT scaling is not included, the numerical calculations fail in reproducing the
experimental observations: Breakthrough time is underestimated and the cumulative effluent of
Water and MeOH is overestimated (compared to the case where IFT scaling is performed) at the
cost of iC
8
. In addition, the IPA is predicted to reside preferentially in the Water-rich phase
rather than in the hydrocarbon-rich phase as observed in the experiments.
Some degree of mismatch was observed between calculated and experimental effluent
concentrations of Water and IPA. We attribute this to experimental errors introduced in effluent
sampling and GC analysis. Additional error can be introduced via the phase behavior model used
in the displacement calculations. However, we believe that the latter source of error is marginal
based on the experiments and analysis presented by Rastegar and Jessen (2011).
Finally, the use of a simple linear mass-transfer model appears to be sufficient to represent
the component fluxes in and out of the stagnant porosities. While more sophisticated models
exist (e.g. Maxwell-Stephan diffusion as discussed by Taylor and Krishna (1993)) and may be
warranted in some displacement processes, the additional work required to determine component
specific coefficients, does not appear to be crucial in the context of the displacement experiments
reported in this work.
3.11 Conclusions
Based on the results and analysis presented in this chapter, we conclude that:
 The use of phase velocities or total velocity in dispersion modeling of near-miscible and
MCM displacement processes with a moderate level of mixing does not have a significant
55

impact on the calculation results. This is despite the fact that two phases are flowing
simultaneously inside the porous medium as a result of dispersive mixing and mass transfer
between the stagnant and flowing porosities.
 The two representations of a DP system used in our calculations provide for almost identical
results. Therefore, it is not clear which model (DP I or DP II) provides for a better
representation of the unresolved heterogeneity inside the column using the Coats-Smith
model (A 1D model that accounts for some 2D/3D effects via a mass transfer term). This
observation illustrates the need for 2D/3D displacement experiments where heterogeneity is
carefully characterized to further test/validate the physical models that are used in
compositional simulation.
 IFT scaling of relative permeabilities is required in our calculations to match the
experimental observations from displacement experiments where dispersive mixing and mass
transfer between flowing and stagnant porosities promote some two-phase flow inside the
porous medium. This observation contrasts the theory of gas injection processes that studies
the dispersion-free limit of 1D flows where development of miscibility is independent of
relative permeability.

a)
56

Chapter 4 : Diffusion and Matrix-Fracture Interactions during Gas
Injection in Fractured Reservoirs
1

4.1 Introduction
Gravity drainage, molecular diffusion and viscous displacement are known as the main recovery
mechanisms 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 oil recovery because the injected gas channels
through high-permeability fractures which comprise only a few percentage 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. Different investigators have demonstrated the efficiency 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).

1
Most of the results in this chapter have been presented in the SPE Improved Oil Recovery Symposium:
Shojaei, H., and Jessen, K.: Diffusion and matrix-fracture interactions during gas injection in fractured reservoirs.
Paper SPE 169152 presented at the SPE Improved Oil Recovery Symposium, Tulsa, OK, 12-16 April (2014)
57

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. Gravity drainage has a similar
effect on displacement and sweep efficiencies. In this chapter we focus on molecular diffusion
and matrix-fracture exchange while more details on other drive mechanisms in fractured
reservoirs can be found elsewhere (e.g. Lemonnier and Bourbiaux, 2010; Chordia and Trivedi,
2010; Rezaveisi et al., 2012).
Diffusive mass transfer in porous media can be attributed to three main mechanisms:
Knudsen (free-molecule) diffusion, molecular (continuum) diffusion and surface diffusion.
Knudsen diffusion occurs when the pore size is smaller than the mean free path of the molecules:
Knudsen diffusion is dominant when the pore size is in the range of nanometers. Molecular
diffusion takes place when the pore sizes are relatively large compared to the mean free path of
the molecules. Surface diffusion describes the transport of matter in an adsorbed layer and
involves interactions between surface and molecules (Mason and Malinauskas, 1983). In this
work Knudsen and surface diffusions are neglected because it is assumed that the pores are in the
range of micrometers and that there is no adsorbed layer on the walls of the pores.
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 media 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).
58

Single-phase multicomponent molecular diffusion can be modeled using different
approaches: The classical Fick’s law, the generalized Fick’s law and the Maxwell-Stefan (MS)
model. In the classical Fick’s law approach, the diffusion flux of each component is only a
function of its own concentration gradient. In other words, the interactions among different
species (dragging effects) are neglected. The generalized Fick’s law approach, in contrary, takes
into account the component interactions; i.e. the diffusion 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., 2011), multicomponent molecular diffusion is modeled
using a classical Fick’s law approach with effective diffusion coefficients. This means the
diffusion 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 (reverse diffusion) (Duncan and Toor, 1962).
Therefore ignoring the dragging effects (off-diagonal diffusion coefficients) may lead to
considerable errors in simulation results as will be shown later.
Effective diffusion coefficients for multicomponent mixtures (classical Fick’s law) are
usually calculated using the approach of da Silva and Belery (1989) which is an extension of
Sigmund’s correlation for binary mixtures (Sigmund, 1976). These coefficients correspond to
diagonal elements of the diffusion coefficient matrix. The full-matrix diffusion coefficients
59

(generalized Fick’s law) can be calculated using the MS model (e.g. Taylor and Krishna, 1993;
Ghorayeb and Firoozabadi, 2000; Leahy-Dios and Firoozabadi, 2007) as will be explained later.
The gas-oil diffusion at the fracture-matrix interface (cross-phase diffusion) 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 gas-oil
diffusion coefficient calculations do not have a sound physical basis.
In chemical engineering literature, the cross-phase 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 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 to model gas-oil
diffusion in lab 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 to represent the molecular
diffusion; and gas-oil diffusion at the fracture-matrix interface is modeled based on film theory.
A novel matrix-fracture transfer function is introduced for gas-oil diffusion based on film theory.
Diffusion coefficients are calculated using the MS model and are pressure, temperature and
composition dependent. A time-dependent transfer function is used for the matrix-fracture
interactions in which the shape factor is adjusted using a boost factor to differentiate between the
transfer rate at early and late times.
Field-scale examples are used 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
60

conventional approach. It is also demonstrated that the dragging effects (off-diagonal diffusion
coefficients) and time-dependency of matrix-fracture transfer function can moderately impact the
oil recovery during gas injection in fractured reservoirs. We also show that miscibility is not
developed in the matrix blocks even at pressures above minimum miscibility pressure (MMP)
when molecular diffusion is the main recovery mechanism during gas injection in fractured
reservoirs.
The remainder of this chapter is organized as follows: The governing equations for transport
of fluid in a dual-porosity reservoir are first presented. Intra-phase and cross-phase diffusion are
then discussed based on generalized Fick’s law and film theory. A generalized matrix-fracture
transfer function is subsequently presented. Different physical models for molecular diffusion
and matrix-fracture interactions are then compared using field-scale examples. A discussion and
analysis of the presented results concludes the chapter.
4.2 Mathematical Model
We use a dual-porosity approach to model gas injection in fractured reservoirs. Dual-porosity
models assume there are two communicating domains in a fractured reservoir: a flowing domain
(fracture) and a stagnant domain (matrix). Mass (or volume) balance equations are solved
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 dual-porosity formulation because it provides a practical representation of
fractured reservoirs with reasonable accuracy and computational efficiency. Other approaches
such as fine-grid single-porosity and discrete fracture models are computationally expensive and
need extensive data which are not generally available (e.g. fracture distribution in the reservoir).
61

In addition, a dual-porosity model can be constructed from detailed discrete fracture
characterization using proper upscaling techniques (e.g. Karimi-Fard et al., 2006; Gong et al.,
2008).
The mass conservation equation for a component i in the flowing domain (fracture)
containing n
p
phases and n
c
components, where advection, diffusion and gravity are the main
physical mechanisms, can be written as (Jessen et al., 2008)

(4.1)
while the mass conservation equation for a component i in the stagnant domain (matrix) is
written as

(4.2)
where the overall molar density C
i
is given by

(4.3)
the overall molar advective flux F
i
is given by

(4.4)
and the overall molar diffusive flux H
i
is given by

(4.5)
where H
i,j
denotes the molar diffusive flux of component i within phase j (intra-phase diffusion)
and H
i,jk
represents the molar diffusive flux of component i at the interface between phases j and
62

k (cross-phase diffusion). More details on calculation of molar diffusive fluxes due to intra-phase
and cross-phase diffusion are given in the subsequent section.
In Equations (4.1) through (4.4), ϕ is the porosity; x
ij
is the mole fraction of component i in
phase j; ρ
j
and S
j
represent the molar density and saturation of phase j respectively; v
j
is the
velocity vector of phase j; q
i
is the source term of component i; and Γ
i
represents the transfer of
component i from matrix to fracture.
The phase velocities in Eq. (4.1) are evaluated from the pressure field which is obtained by
solving the following volume-balance equation

(4.6)
where c
t
is the total fluid compressibility; p is the pressure; V
t
is the total fluid volume; V
cell
is the
pore volume; and n
i
is the number of moles of component i. Capillarity is not included because
we focus on gas injection processes in which IFT effects are less significant.
We use a finite volume IMPES (implicit pressure explicit saturation/composition)
formulation with a Cartesian grid in this work. This means at each new time level, Eq. (4.6) is
solved with coefficients fixed at the old time level. The discrepancy between cell volume and
fluid volume is minimized by carrying errors forward in time (Trangenstein and Bell, 1989;
Jessen et al., 2008).
Once the pressure field is obtained at the new time level using Eq. (4.6), phase velocities at
the gridblock interfaces are obtained using Darcy’s law

(4.7)
63

where K is the permeability tensor; λ
j
and ρ
mj
are the mobility and mass density of phase j; and g
is the gravitational vector. The total number of moles of each component in a gridblock l is then
updated using

(4.8)
where A
lm
denotes the area connecting gridblocks l and m; (F+H)
i,lm
is the total molar flux
(advective+diffusive) of component i out of gridblock l at the interface lm, and q
i,l
represent the
source term. The phase-equilibrium calculations are performed using Peng-Robinson (PR)
equation of state (EOS). The advective flux is calculated using the standard single point upwind
(SPU) scheme.
4.3 Molecular Diffusion
Molecular diffusion may occur within a single phase (intra-phase) or between two different
phases (cross-phase). For two neighboring gridblocks (l and m) that contain oil and gas phases
(Fig. 4.1), the cross sectional areas that are available for gas-gas, oil-oil and gas-oil diffusion are
A
g
=A
lm
×min(S
g,l
,S
g,m
), A
o
=A
lm
×min(S
o,l
,S
o,m
) and A
c
=A
lm
-A
g
-A
o
respectively.

Figure ‎ 4.1 Schematic of two neighboring gridblocks containing oil and gas phases.

64

In the following subsections we explain how to calculate intra-phase and cross-phase molecular
diffusion based on generalized Fick’s law and film theory.
4.3.1 Intra-phase Diffusion
The molar diffusive flux of a component i in phase j (generalized Fick’s law) is given by

(4.9)
where D
ik,j
are the multicomponent diffusion coefficients in phase j. The sum of diffusive fluxes
of all components in phase j must be zero (

); and the diffusive flux for the last
component (H
nc,j
) is obtained accordingly (Taylor and Krishna, 1993). The diagonal elements in
the diffusion coefficients matrix D are called the main diffusion coefficients; while the off-
diagonal elements are called the cross-diffusion coefficients which are generally non-zero and
non-symmetric (i.e. D
ik
≠ D
ki
for i≠k). There are (n
c
-1)
2
Fickian diffusion coefficients which are
dependent on the numbering of components.
Multicomponent diffusion coefficients can be calculated using the following equation which
is obtained by comparing the generalized Fick’s law and MS model

(4.10)
where 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

(4.11)
65

and

(4.12)
where Đ
ik
represent the MS coefficients; δ
ik
is the Kronecher delta function 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 n
c
(n
c
-1)/2 MS coefficients which can be obtained using
(Wesselingh and Krishna, 1990; Kooijman and Taylor, 1991)

(4.13)
where Đ
o
ik
is the diffusion coefficient of component i infinitely diluted in component k. In this
work we use the correlation by Leahy-Dios and Firoozabadi (2007) to calculate the infinite
dilution diffusion coefficients for both vapor and liquid phases. See Appendix D for relevant
equations and also for comparisons between the calculated and experimental diffusion
coefficients.
4.3.2 Cross-phase Diffusion
Let us consider two different phases which have been placed in contact as shown in Fig. 4.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. 4.2, x
i,I
and y
i,I
are the interface
compositions; and x
i,b
and y
i,b
denote the bulk phase compositions (Taylor and Krishna, 1993).
66

Figure ‎ 4.2 Composition profiles near the interface during inter-phase mass transfer, after Taylor
and Krishna (1993).

Therefore the cross-phase mass transfer for a stationary interface is governed by two sets of
equations: Continuity of molar fluxes at the interface

(4.14)
and thermodynamic equilibrium

(4.15)
where, assuming the positive direction is from left to right,

denotes the total molar flux of
component i from “x” phase into the interface;

denotes the total molar flux of component i
from interface into “y” phase; and

(

) is the chemical potential of component i in “x” (“y”)
phase at the interface.
The molar flux consists of advective and diffusive fluxes. Since the advective flux is
computed using an upwind scheme, the continuity of advective flux will always be honored.
Therefore the continuity of the molar flux (Eq. 4.14) reduces to the continuity of diffusive flux
which, based on the generalized Fick’s law, can be written in the following finite difference form
67

(4.16)
where L
x
(L
y
) is the distance between the interface and the center of the gridblock containing the
“x” (“y”) phase. To calculate the molar cross-phase diffusive fluxes, Equations (4.15) and (4.16)
need to be solved simultaneously. The interface equilibrium calculations (Eq. (4.15)) can be done
using an EOS at interpolated temperature and pressure. Multicomponent diffusion coefficients
are obtained using Eq. (4.10).
The solution method described in the above paragraph, is an iterative and hence time-
consuming approach. Hoteit (2013) proposed a non-iterative solution based on gradients of the
chemical potential (instead of concentration gradients) which does not need explicit knowledge
of interface compositions. In this approach, which is used in our work, the continuity of diffusive
fluxes is given by

(4.17)
where

; B is the matrix given by Eq. (4.11); X is a diagonal matrix with diagonal
elements [1/x
k
]
k=1,…,nc-1
; ψ
k
=ln(f
k
); and f
k
is the fugacity of component k. We note that because of
the chemical equilibrium condition at the interface (Eq. (4.15)), we have

=

. Eq.
(4.17) can be rearranged to obtain the following expression in matrix form

(4.18)
where

and

. Once

is calculated, it can be substituted in either left-
hand side or right-hand side of Eq. (4.17) to calculate the molar diffusive flux of each component
across the interface. The sum of diffusive fluxes of all components across the interface must be
zero; and the diffusive flux of the last component is obtained accordingly.
68

4.4 Transfer Function
Since the introduction of dual-porosity models in the 1960s, various transfer functions with
different levels of sophistication have been suggested for matrix-fracture interactions. Warren
and Root (1963) proposed a pseudo-steady state (PSS) transfer function for single-phase flow
between matrix and fracture. In their model, the transfer rate per unit bulk volume is obtained by
multiplying the average pressure difference between matrix and fracture by mobility and a
parameter which is commonly referred to as the shape factor. The shape factor has a dimension
of reciprocal area and is a characteristic of the fractured rock.
The application of the shape factor in reservoir simulation was first introduced by Kazemi et
al. (1976) for three-phase flow subject to a pseudo-steady state assumption. A gravity term was
later included in the transfer function by Gilman and Kazemi (1983), Gilman (1986) and Sonier
et al. (1988). However, it was added for flow between matrix and fracture in all three directions,
which could lead to an overestimation of the transfer rate. Quandallet and Sabathier (1989)
resolved this issue by treating flow in horizontal and vertical directions separately.
Coats (1989) extended the dual-porosity model for compositional simulation and included
molecular diffusion in the transfer function. He suggested the use of shape factors that are
exactly twice those of Kazemi et al. (1976). By solving the single-phase pressure diffusion
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 pseudo-steady state condition. While pseudo-steady state transfer functions
are only good for late times, the transfer functions developed by Zimmerman et al. (1993) and
Lim and Aziz (1995) are able to provide reasonable accuracy for both early and late times.
69

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 multi-phase matrix-fracture
interactions. Their approach, which is used in our work with some modifications, is explained in
the following paragraphs. The transfer rate of a component i from matrix to fracture is given by

(4.19)
where Γ
i,e
, Γ
i,d
and Γ
i,gd
are transfer rates due to expansion, diffusion and gravity drainage
respectively. We note again that capillary effects are assumed to be negligible for gas injection
processes. The expansion term is given by

(4.20)
where σ
1
is the shape factor; B
e
is the boost or correction factor for the expansion term; K
m
is the
(equivalent isotropic) matrix permeability; λ
jm
is the mobility of phase j in the matrix; x
ij,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 p
m
and p
f
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 this work

(4.21)
where K
m,x
, K
m,y
and K
m,z
are the matrix permeabilities in x-, y-, and z-direction respectively; and
L
x
, L
y
, L
z
are the dimensions of matrix block in x-, y-, and z-direction respectively. The
equivalent isotropic matrix permeability is defined as

(4.22)
70

The boost factor for expansion is given by the following equation (Zimmerman et al., 1993; Lu
et al., 2008)

(4.23)
where p
m,init
stands for the initial pressure in the matrix. At late times, when matrix and fracture
pressures are similar, the boost factor approaches unity and Eq. (4.20) reduces to Kazemi’s
pseudo-steady state transfer function. However, at early times when the matrix pressure is close
to its initial value, the correction factor is larger than one.
To calculate the transfer rate of each component from matrix to fracture due to molecular
diffusion, we sum the intra-phase and cross-phase diffusion rates of that component. For a two-
phase gas-oil system, the transfer rate of a component i is calculated from

(4.24)
where Γ
i,d,g
, Γ
i,d,o
and Γ
i,d,c
are transfer rates due to gas-gas (intra-phase), oil-oil (intra-phase) and
gas-oil (cross-phase) molecular diffusion respectively. The transfer rate of a component i due to
gas-gas diffusion is given by

(4.25)
where B
d,g
is the boost factor for gas-gas diffusion; ϕ
m
is the matrix porosity and S
g
=min(S
gm
,S
gf
)
is the fraction of the total matrix-fracture surface area that is available for gas-gas diffusion.
Since the molecular diffusion equation is analogous to the pressure diffusion equation, the
following expression can be used to calculate the boost factor for gas-gas diffusion

(4.26)
71

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 the film theory as explained earlier. For the matrix-fracture system, the
continuity of diffusive flux across the interface can be written as

(4.27)
where S
go
=1- min(S
gm
,S
gf
) - min(S
om
,S
of
) is the fraction of the total matrix-fracture surface area
that is available for gas-oil diffusion; σ
1
is obtained from an isotropic version of Eq. (4.21); and
σ
2
is given by

(4.28)
where b
x
, b
y
, and b
z
are the fracture openings in x-, y-, and z-direction respectively. The
derivations of Eqs. (4.27) and (28) are given in Appendix E.
By rearranging Eq. (4.27), we arrive at an equation similar to Eq. (4.18) in matrix form

(4.29)
where

and

. Once

is calculated, it can be substituted in
either left-hand side or right-hand side of Eq. (4.27) to calculate the cross-phase molar diffusive
flux of each component from matrix to fracture. The boost factor for gas-oil diffusion can be
obtained from

(4.30)
72

The transfer rate of each component from matrix to fracture due to gravity drainage is calculated
using the equations given in Appendix F.
4.5 Results
Two numerical examples are presented in this section to demonstrate the significance of proper
modeling of molecular diffusion and matrix-fracture 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 a permeability of 500 md in all three directions; and a porosity of 1%. In the stagnant
domain, the permeabilities in horizontal directions (K
x
and K
y
) are 1 md while the vertical
permeability is 0.001 md. The stagnant domain has a porosity of 20%. The fracture spacing and
opening in both x- and y-direction are 10 m and 1 mm respectively.
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 both
injection and production wells, with values that will be specified for each example. A schematic
of the reservoir and wellbore locations is shown in Fig. 4.3.
73

Figure ‎ 4.3 Schematic of the reservoir and wellbore locations. White lines represent the fracture
system.

Figure ‎ 4.4 Relative permeability curves for flowing (solid line) and stagnant (dashed line)
domains. Gas and oil relative permeabilities are shown with red and green colors respectively.

Linear relative permeabilities are used for the flowing domain (fractures); while Corey-type
relative permeability curves are used for the stagnant domain (matrix). The saturation exponents
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
rel perm
Sg
kro_m
krg_m
kro_f
krg_f
74

for gas and oil phases are 3.0 and 4.0 respectively; and the critical gas saturation and residual oil
saturation are 0.0 and 0.25 respectively for the stagnant domain. Relative permeabilities for both
flowing and stagnant domains are shown in Fig. 4.4.
4.5.1 Example 1
In this example, the initial reservoir pressure at the top of formation is 157 bar. The bottomhole
pressure for the production well is kept constant at 157 bar and gas is injected at an average rate
of 165 Rm
3
/day. The injected gas is composed of 90% CO
2
and 10% CH
4
by mole. The reservoir
is initially saturated with oil of 30% CH
4
, 40% nC
4
H
10
and 30% nC
12
H
26
by mole. The reservoir
temperature is 373.15 K at which the minimum miscibility pressure (MMP) of the oil/gas pair is
150.5 bar. The oil has a bubble point pressure of 91.6 bar at the reservoir temperature. The initial
and injection compositions together with the MMP of the gas/oil pair used in example 1 are
summarized in Table 4.1.
Table ‎ 4.1 Initial and injection compositions and the MMP of the gas/oil pair used in example 1.

Composition (mole %) MMP (bar)

CO
2

CH
4
nC
4
H
10
nC
12
H
26

Initial 0 30 40 30
} 150.5
Injection 90 10 0 0

Fig. 4.5 compares the simulation results for example 1 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 differences between the conventional
approach and our proposed approach are summarized below:
In the conventional approach,
75

1- Molecular diffusion is modeled using the classical Fick’s law in which the off-diagonal
diffusion coefficients are neglected. Moreover, constant diffusion coefficients are used by
ignoring the dependency of diffusion coefficients on pressure, temperature and
composition. The diagonal elements of the diffusion coefficient matrix obtained from Eq.
(4.10) at initial and injection conditions are used for the entire simulation (constant
diffusion coefficients) for oil and gas phases respectively.
2- The gas-oil diffusion is modeled by assuming an average composition at the interface. By
assuming an average interface composition in Eq. (4.16), one can obtain expressions
similar to Eqs. (4.27) and (4.29) for gas-oil diffusion.
3- Pseudo-steady state condition is assumed for the matrix-fracture transfer function. This is
equivalent to setting the boost factor equal to one. Pseudo-steady state assumption is only
valid for late times.
In our proposed approach,
1- Molecular diffusion is modeled using the generalized Fick’s law. Diffusion coefficients
(full matrix) are obtained based on the MS model (Eq. (4.10)) and are pressure-,
temperature- and composition-dependent.
2- 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 are at thermodynamic equilibrium. In addition, the molar flux of each
component across the interface must be continuous.
3- A time-dependent transfer function is used for matrix-fracture interactions. The boost
factor that is included in the transfer function differentiates between the early- and late-
time behaviors of matrix-fracture interactions.
76

The left tab in Fig. 4.5 shows liquid recovery (%) while the right tab shows producing gas-oil
ratio (GOR) versus time for ten years. In the absence of diffusion, the injected gas mainly sweeps
the oil residing in the flowing domain (fracture) which in turn leads to very low liquid recovery
and very high producing GOR. However, when molecular diffusion is present, diffusion of gas
components from fracture into the matrix extracts oil components from matrix and leads to
substantially higher liquid recoveries and lower producing GOR.

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

Fig. 4.5 also shows that the simulation results obtained using the conventional approach are
significantly different from the predictions based on our approach: Significantly higher liquid
recoveries and lower gas-oil ratios are obtained using the conventional approach. We defer
further analysis of the differences between predictions by these two different approaches until
the discussion section.
77

The next step is to see how much the results would change if pseudo-steady state (PSS)
condition was assumed for the matrix-fracture (M-F) interactions (i.e. setting the boost factor
equal to one) and off-diagonal diffusion coefficients were neglected (i.e. the classical Fick’s law)
in our approach. Fig. 4.6 compares the simulation results for example 1 when diffusion is not
considered (solid red), when our approach is used (solid blue), when our approach is used but
PSS condition is assumed for M-F transfer function (dashed blue) and when our approach is used
but PSS condition is assumed for M-F transfer function and off-diagonal diffusion coefficients
are neglected (dotted blue).

Figure ‎ 4.6 Liquid 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), when our
approach is used but PSS condition is assumed for M-F transfer function (dashed blue) and when
our approach is used but PSS condition is assumed for M-F transfer function and off-diagonal
diffusion coefficients are neglected (dotted blue).

From Fig. 4.6 we observe that assuming PSS condition in M-F transfer function and neglecting
the off-diagonal diffusion coefficients cause ~ %10 (%3 out of %28) difference in the predicted
78

liquid recoveries using different variants of our approach over the time period investigated. This
is a notable difference, but much smaller than the ~ %60 difference (%17 out of %28) we
observed between the cumulative liquid 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. 4.5). Therefore the key factor in proper modeling of
gas injection in fractured reservoirs when molecular diffusion is the main recovery mechanism,
is using the appropriate physical model for gas-oil diffusion (i.e. film theory), not the time-
dependency of M-F transfer function or the off-diagonal diffusion coefficients.
From Fig. 4.6 it is also observed that when off-diagonal diffusion coefficients are taken into
account, higher liquid recoveries are obtained because of the larger contributions by the positive
off-diagonal coefficients in this case. The effect of including boost factors in the M-F transfer
function is more evident during the earlier times according to Fig. 4.6.
The CO
2
composition (mole %) in the stagnant domain (matrix) after ten years (example 1) is
shown in Fig. 4.7 for the case with no diffusion (left) and the case with diffusion based on
calculations using our approach (right). It is observed that, in the absence of molecular diffusion,
very small amounts of CO
2
have entered the stagnant domain due to viscous displacement and
gravity drainage. However when molecular diffusion is present, considerable amounts of CO
2

enter the stagnant domain and extract the oil components. It should be noted that the black lines
in Fig. 4.7 represent the boundaries of simulation cells, not the actual natural fractures.

79

Figure ‎ 4.7 The CO
2
composition (%) in the stagnant domain after 10 years (example 1) for the
case without diffusion (left) and the case with diffusion based on calculations using our approach
(right).
4.5.2 Example 2
In this example, the initial reservoir pressure is 150 bar, the bottomhole pressure of the
production well is kept constant at 150 bar and gas is injected at an average rate of 180 Rm
3
/day.
The injected gas is composed of 50% CO
2
and 50% CH
4
by mole. The reservoir is initially
saturated with oil of the same composition as in example 1 (30% CH
4
, 40% nC
4
H
10
and 30%
nC
12
H
26
by mole) at a temperature of 336.15 °K. The MMP of the oil/gas pair is 178.8 bar; and
the oil has a bubble point pressure of 78.6 bar. The initial and injection compositions together
with the MMP of the gas/oil pair used in example 2 are summarized in Table 4.2.
Table ‎ 4.2 Initial and injection compositions and the MMP of the gas/oil pair used in example 2.

Composition (mole %) MMP (bar)

CO
2

CH
4
nC
4
H
10
nC
12
H
26

Initial 0 30 40 30
} 178.8
Injection 50 50 0 0

Injector
Producer
Producer
Injector
80

Fig. 4.8 compares the simulation results for example 2 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 shows liquid recovery (%) while
the right tab presents producing GOR versus time for ten years.
Similar to the previous example, low liquid recovery and high producing GOR is observed in
the absence of molecular diffusion. However, when molecular diffusion is present, substantially
higher liquid recoveries and lower producing GOR are obtained. Significantly higher liquid
recoveries and lower gas-oil ratios are calculated using the conventional approach as compared
to our approach.

Figure ‎ 4.8 Liquid recovery (left) and producing GOR (right) versus time for example 2 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 effect of PSS assumption in M-F transfer function and the off-diagonal diffusion coefficients
on simulation results for example 2 is investigated in Fig. 4.9. The same color codes and line
styles as in Fig. 4.6 are used to represent each physical model.
81

Similar to the previous example, assuming PSS condition in M-F transfer function and
neglecting the off-diagonal diffusion coefficients cause notable differences among the results
obtained using different variants of our approach, but much smaller than the difference we
observed between the results of our approach and the conventional approach (see Fig. 4.8).
Similar to the previous example, higher liquid recoveries are obtained when off-diagonal
diffusion coefficients are taken into account, because of the larger contributions by the positive
off-diagonal coefficients. According to Fig. 4.9, the effect of including boost factors in M-F
transfer function is more evident during the earlier times.

Figure ‎ 4.9 Liquid 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), when our
approach is used but PSS condition is assumed for M-F transfer function (dashed blue) and when
our approach is used but PSS condition is assumed for M-F transfer function and off-diagonal
diffusion coefficients are neglected (dotted blue).

Fig. 4.10 shows the CO
2
composition (mole %) in the stagnant domain (matrix) after ten years
(example 2) for the case without diffusion (left) and the case with diffusion as calculated using
82

our approach (right). Very small amounts of CO
2
have entered the stagnant domain (due to
viscous displacement and gravity drainage) when molecular diffusion is absent. However,
considerable amounts of CO
2
enter the stagnant domain and extract the oil components when
molecular diffusion is present.

Figure ‎ 4.10 The CO
2
composition (%) in the stagnant domain after 10 years (example 2) for the
case without diffusion (left) and the case with diffusion based on calculations using our approach
(right).

4.6 Discussion
To understand how the oil is recovered by molecular diffusion during gas injection in fractured
reservoirs, we analyze the saturation and composition profiles of a matrix block and its
corresponding fracture. Fig. 4.11 shows the saturation profile observed during the first ten years
of production 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 one year, after which it gradually increases.
Injector Injector
Producer
Producer
83

Figure ‎ 4.11 The saturation profile of a matrix block (dotted line) near the injection well and its
corresponding fracture (solid line) for example 1 as calculated using our approach.

Fig. 4.12 shows the composition routes (mass fraction) of a matrix block (green) near the
injection well and its corresponding fracture (red) for example 1 during the first ten years of
production as calculated using our approach. The blue curves show the two-phase boundaries in
each ternary plane at 157 bar and 373.15 °K; while the dashed brown line represents the dilution
line.
Although the pressures in example 1 are above the MMP, miscibility is not developed in the
matrix block. As observed from Fig. 4.12, molecular diffusion generates a composition profile
deep into the two-phase region. Initially, the CH
4
composition increases with a slight increase in
CO
2
composition; while the nC
4
H
10
(nC
4
) composition decreases and nC
12
H
26
(C
12
) composition
remains almost unchanged. Later, mainly CO
2
composition increases and C
12
composition
decreases until the composition in the matrix approaches injection composition. It is also
observed that the composition profile in the fracture departs from the multi-contact miscibility
(MCM) condition due to exchange with the matrix (and the limited number of cells used in the
representation of the domain).
84

Figure ‎ 4.12 Composition profiles (mass fraction) of a matrix block (green) near the injection
well and its corresponding fracture (red) for example 1 during the first 10 years as calculated
using our approach. The blue curves show the two-phase boundaries in each ternary plane at 157
bar and 373.15 °K. The dashed brown line represents the dilution line.

Fig. 4.13 shows the saturation profile observed during the first ten years of production at the
same matrix block (dotted line) as in Fig. 4.11 and its corresponding fracture (solid line) as
calculated using the conventional approach. As before, the fracture becomes fully saturated with
gas in a short period of time; while 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 ~ one year when our approach was used (see Fig. 4.11).

85

Figure ‎ 4.13 The saturation profile of a matrix block (dotted line) near the injection well and its
corresponding fracture (solid line) for example 1 as calculated using the conventional approach.

Fig. 4.14 can explain why a gaseous phase does not form in the matrix until after ~ seven years
as calculated using the conventional approach. During the first seven years, the CO
2
composition
largely increases with small changes in compositions of other components. This keeps the fluid
inside the matrix in the single phase region for nearly seven years. The resulting swelling of oil
cause large liquid recoveries. The composition profile in the matrix, as calculated by the
conventional approach, differs significantly from what is calculated using our approach. This
explains why the recovery calculations using these two different approaches differ significantly.
The CO
2
composition (mole %) in the stagnant domain (matrix) after ten years (example 1)
as calculated by the conventional approach (left) is compared with calculations using our
approach (right) in Fig. 4.15. It is observed that the two different approaches predict significantly
different distributions of CO
2
composition in the stagnant domain after 10 years.

86

Figure ‎ 4.14 Composition profiles (mass fraction) of a matrix block (green) near the injection
well and its corresponding fracture (red) for example 1 during the first 10 years as calculated
using the conventional approach. The blue curves show the two-phase boundaries in each ternary
plane at 157 bar and 373.15 °K. The dashed brown line represents the dilution line.

Figure ‎ 4.15 The CO
2
composition (%) in the stagnant domain after 10 years (example 1) as
calculated using the conventional approach (left) and calculations using our approach (right).

Producer
Injector
Producer
Injector
87

Capillary forces have been 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.
4.7 Conclusions
In this chapter, we presented a dual-porosity model in which the generalized Fick’s law is used
to represent molecular diffusion; and the 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. Diffusion coefficients are calculated based on MS model and
are pressure, temperature and composition dependent. In the matrix-fracture transfer function, a
boost factor is used to adjust the shape factor in order to account for the differences between the
transfer rate at early and late times. The following conclusions can be drawn based on the results
presented in this chapter:
 Molecular diffusion can substantially enhance the oil recovery by gas injection in fractured
reservoirs. Diffusion of gas components from fracture into the matrix extracts oil components
from matrix and leads to substantially higher liquid recoveries and lower producing GOR.
 Our approach, which is based on a sophisticated physical model for gas-oil diffusion (film
theory), gives significantly different results from the conventional approach which assumes
an average composition at the gas-oil interface.
 The dragging effects (off-diagonal diffusion coefficients) and time-dependency of matrix-
fracture transfer function can moderately impact the oil recovery during gas injection in
88

fractured reservoirs.
 Miscibility is not developed in the matrix block even at pressures above MMP when
molecular diffusion is the main recovery mechanism during gas injection in fractured
reservoirs. In this case, molecular diffusion pushes the composition profile deep into the two-
phase region.

89

Chapter 5 : Application of Potential Theory of Adsorption to Modeling of
ECBM Recovery
1

5.1 Introduction
In coalbeds, gas is present in two phases: a bulk gaseous phase that occupies the pore space and a
liquid-like adsorbed phase on the coal surfaces/pores. In primary recovery, the coalbed is
dewatered to reduce the overall reservoir pressure which causes CH
4
to be desorbed from coal
surfaces. Since primary production typically recovers less than half of the methane in a coalbed
(Stevens et al., 1998), enhanced coalbed methane (ECBM) recovery processes are needed in
which CO
2
and/or N
2
are injected into the coalbed to recover more CH
4
. Injection of CO
2
into
coalbeds also provides additional benefit of sequestering carbon in the subsurface.
Gas injection in ECBM recovery provides a method to maintain the overall coalbed pressure.
In addition, injecting a second gas, or a mixture of gases, decreases the partial pressure of CH
4
in
the free gas. As a result, desorption of CH
4
from coal surface is enhanced. The convective flow
of injected gas sweeps desorbed CH
4
through the coalbed. Therefore the sorption of gases onto
the coal surfaces is one of the main mechanisms that govern the dynamics of ECBM recovery
processes.
Despite the well-documented complexity of multicomponent sorption phenomena (Stevenson
et al., 1991; DeGance et al., 1993; Chaback et al., 1996), adsorption and desorption of
CH
4
/CO
2
/N
2
mixtures in ECBM recovery is usually modeled with the extended Langmuir model

1
Most of the results in this chapter have been presented in the SPE Western Regional Meeting:
Shojaei, H. and Jessen, K.: Application of potential theory to modeling of ECBM recovery. Paper SPE 144612
presented at the SPE Western Regional Meeting, Anchorage, Alaska, 7-11 May (2011)
90

because of its simplicity and associated low computational cost (Guo, 2003; Zhu et al., 2003;
Smith et al., 2005; Seto et al., 2009). The extended Langmuir has been proven unable to
accurately describe the multicomponent sorption that is relevant to ECBM recovery processes
(Clarkson, 2003; Wei et al., 2005). Jessen et al. (2008) demonstrated that extended Langmuir
was able to model the sorption process in binary displacements; but failed to describe the
behavior of ternary displacements.
The ideal adsorbate solution (IAS) theory has been used to improve the agreement between
calculated and measured values in ternary sorption processes. Although improvements have been
attained over the extended Langmuir, still considerable mismatch between the calculated and
measured values remained (Manik, 1999; Jessen et al., 2007). In addition, Manik et al. (2002)
showed that at high pressures, predictions from IAS theory deviated from measured adsorption
capacities due to non-ideality in the adsorbed phase. The real adsorbate solution (RAS) theory
may be used to account for the non-ideality of adsorbate; but due to the complexity of required
calculations and the large number of adjustable parameters involved in RAS modeling, it may
not be suitable for ECBM recovery calculations.
In recent years, a new model for the sorption behavior of gases and liquids in porous materials
has been proposed (Shapiro and Stenby, 1998). The so called “multicomponent potential theory
of adsorption” (MPTA) is based on the potential concept originally suggested by Polanyi (1932).
The successful application of MPTA for sorption calculations in microporous materials (Shapiro
and Stenby, 1998; Monsalvo and Shapiro, 2007; Monsalvo and Shapiro, 2009), and the lack of
accuracy in previous ECBM calculations (Jessen et al., 2007; Jessen et al., 2008), motivated us to
test this model in the context of ECBM recovery calculations. To the best of our knowledge, this
is the first time that MPTA is used for sorption calculations in this context.
91

The remainder of this chapter is organized as follows: First, we describe the multicomponent
potential theory of adsorption (MPTA) and its implementation with emphasis on ECBM
recovery calculations. Subsequently, we apply the MPTA approach to multicomponent gas
sorption on real coal and compare predictions to experimental observations. Several example
calculations illustrate the improved predictive capability of the MPTA model. A discussion and
conclusions section completes the manuscript.
5.2 Potential Theory
Shapiro and Stenby (1998) developed the MPTA model based on the potential concept originally
suggested by Polanyi (1932). In this theory, the mixture is considered to be a segregated phase in
the potential field emitted by the adsorbent. Each component in the mixture is affected by the
adsorption potential

. The isothermal equilibrium state between the bulk and adsorbed phases
in the external potential field is described by

(5.1)
where is interpreted as porous volume, is the pressure in the sorbed phase while and x are
mole fractions of the bulk phase and adsorbed phase, respectively. The above equation can be
written in a more convenient form using fugacities

that can be estimated by an appropriate
equation of state:

(5.2)
where R and T represent the gas constant and the temperature respectively. Given pressure

and composition in the bulk phase, the distribution of pressures and mole fractions
are uniquely determined by Eq. (5.2). To solve the system of equations (5.2), we need a
92

representation of the adsorption potentials

. The generalized Dubinin dependence for
adsorption potentials in porous media has the following form:

(5.3)
where

is the total porous volume,

is the characteristic potential for component , and is
the Dubinin exponent. A value of corresponds to the standard Dubinin-Radushkevich
(DR) potential and is usually used for the activated carbon. ,

and

are commonly used as
adjustable parameters and are obtained by matching the MPTA model to relevant single-
component sorption isotherms.
The surface excess is defined as the difference between the actual amount adsorbed on the
rock surface and the amount present in the bulk phase of the pores pace if there was no
adsorption potential. The value of surface excess Γ

for each component is determined using:
Γ

(5.4)
where

is total pore volume, and is the molar density (obtained from an equation of state).
Since the integrand in Eq. (5.4) cannot be written explicitly in terms of porous volume z, the
above integral must be evaluated by numerical integration. After calculating the surface excess
for each component, the total molar excess is calculated from
Γ Γ

, (5.5)
and the average mole fraction of each component in the adsorbed phase is calculated from:

Γ

Γ

. (5.6)

93

5.3 Numerical Approach
Application of the MPTA model involves the solution of Eqs. (5.2) to (5.6) by a suitable
numerical scheme. The evaluation of Eq. (5.6) requires numerical integration and hence repeated
solution of the equilibrium problem stated by Eq. (5.2) at the nodes of any selected integration
rule. The equilibrium problem stated by Eq. (5.2) is essentially similar to a dew point calculation
where the solution is given by an incipient phase (here the adsorbed phase) and the associated
equilibrium pressure (in the adsorbed phase). Given the similarity to dew point calculations, we
adapt the approach by Michelsen (1985) for calculation of saturation points.
At any given node (z) of the integration rule, we know the temperature (T), the bulk phase
composition (y) and the pressure of the bulk phase (p
y
). From these, we wish to calculate the
composition (x) and pressure (p) of the adsorbate. We start by rewriting Eq. (5.2) in terms of the
fugacity coefficients, φ

. (5.7)
Next, we introduce the equilibrium K-value

γ

, (5.8)
with
γ

. (5.9)
The requirement for equilibrium can now be stated as

(5.10)
94

and we can proceed by solving Eq. (5.10) for the equilibrium pressure (p) and the incipient phase
composition (x) following the ideal solution based method proposed by Michelsen (1985). The
iterative procedure is as follows: Given an initial estimate of p and x, we can calculate the K-
values from Eq. (5.8) at iteration level k

γ

. (5.11)
Next, we evaluate the trial function (Eq. 5.10) and the derivative wrt p

,
(5.12)

.
(5.13)
The pressure is then updated to the next level by a Newton correction

, (5.14)
and the mole fractions of the adsorbed phase are finally updated by direct substitution

.
(5.15)
Equations (5.11) to (5.15) are then repeated until convergence. This approach for solving the
vapor adsorbate equilibrium problem is simple and relatively inexpensive as the derivatives of
the fugacity coefficients wrt composition are not needed. The solution strategy, does not
guarantee that the equilibrium adsorbate phase is stable as pointed out by Shapiro and Stenby
(1998). Stability testing (Michelsen, 1982) can be applied to the adsorbate phase but we have
found this to be unnecessary for temperatures above the critical temperature of CO
2
(most
coalbeds).
The overall algorithm for calculation of excess sorption is as follows:
95

1- Select a number of interior points, z, within the integration interval

as dictated by
the selected integration rule. In this work we use Simpson’s rule for simplicity.
2- Start at the upper integration limit and obtain initial estimates for and in the
adsorbate. For

, p
y
is a good initial estimate for p, and the extended Langmuir
approach provides initial estimates of x.
3- For subsequent (smaller) values of z, the converged values of p and x from previous
integration node provide good initial estimates. The initial estimate of p can be refined by
using the derivative of Eq. (5.10) wrt z from the previously converged level
. (5.16)
4- Evaluate the surface excess for each component using Eq. (5.4) and integration rule.
5- Evaluate the mole fraction of each component in the adsorbate phase using Eq. (5.6).
Initial tests of our implementation were performed for experimental observations of pure
component (CH
4
, CO
2
and N
2
) and mixture sorption on well characterized activated carbon. We
used the experimental data of Dreisbach et al. (1999) to compare our calculations with the work
of Monsalvo and Shapiro (2007) and found a good agreement. In the next section, we apply the
MPTA model to simulate the sorption behavior of CH
4
, CO
2
, N
2
and their mixtures on real coal
samples.
5.4 Application of MPTA to Coal
Ottiger et al. (2008a, 2008b) presented a comprehensive study of the sorption characteristics of
CH
4
, CO
2
, N
2
and their binary and ternary mixtures on dried coal samples from the Sulcis Coal
Province in Italy. They applied lattice density functional theory (DFT) and arrived at a good
agreement between modeling and experimental observations.
96

Some of the key challenges in modeling sorption on real coal samples are related to the multi-
porosity and highly heterogeneous nature of coals. From careful characterization of the coal
samples including low pressure CO
2
sorption and Hg porosimetry, Ottiger et al. (2008b), report
significant contributions from both micro- and meso- to macro-pores to the overall porosity of
the coal with the micropores making up approximately 0.07 cc/gram of the 0.202 cc/gram of
pore space. In their lattice DFT modeling efforts, they select three discrete pore sizes: Two in the
microporous region (1.2 and 1.6 nm) and one to represent the entire meso- to macro-porous
region (20nm). In addition, they include the assumption that CH
4
, CO
2
and N
2
have access to
different portions of the pore volume in their modeling efforts.
Both the multi-modal pore size distribution and the need for restricting species for some
ranges of porosity suggests that MPTA may not be able to fully capture the complex sorption
characteristics of coal given the inherent continuum scale equilibrium assumption. In addition, a
simple potential function (e.g. Eq. (5.3)) may not be sufficient to represent the interactions
between coal surfaces and gases over the entire range of the pore size distribution.
With this in mind, we apply the MPTA model in its original form (Shapiro and Stenby, 1998)
by fitting the model parameters (characteristic energies and pore volume) to the pure component
isotherms. We use the Peng-Robinson equation of state with volume shift throughout this work.
All binary interaction coefficients are set to zero to avoid additional parameter estimation.
Figure 5.1 shows the pure component isotherms for the three gases at 45 °C in terms of the
excess sorption at pressures of 5-188 bars. The associated model parameters are reported in
Table 5.1.

97

Figure ‎ 5.1 Pure component isotherms at 45 °C. Data from Ottiger et al. (2008b). Solid lines:
MPTA after regression.

Table ‎ 5.1 MPTA parameters from pure component isotherms.
CH
4
CO
2
N
2

ε
0i
/R

(K) 425.1 602.8 269.3
beta (-) 0.59
z
0
(cc/gram) 0.122

We note that the pore volume calculated from simultaneous regression of the model parameters
to the pure component isotherms is less than what is reported by Ottiger et al. (2008b). This can
be interpreted in a multitude of ways including: a) the entire pore volume does not contribute to
the excess sorption for this sample or b) the potential function that we use in this work does not
adequately represent the interactions between solid and the gases and c) the equation of state
does not capture the non-ideal behavior of the sorbed phase density accurately.
While b) and c) are most likely both true, the fact that the pore volume needed in the potential
function to match the experimental observations differs significantly from the measured pore
0
0.5
1
1.5
2
2.5
0 50 100 150 200
Excess sorption (mmol/gram)
Pressure (bar)
CH4
CO2
N2
98

volume suggests that the MPTA model, in the original form, may not be suitable for modeling
sorption on complex coal materials. In addition, we observe that the value of beta is somewhat
low relative to what has been used for more homogeneous microporous materials such as
activated carbon. We return to this after looking at the performance of the MPTA model for
prediction of sorption for ternary mixtures of CH
4
, CO
2
and N
2
.
Based on the parameters listed in Table 5.1, we have applied the MPTA model, in predictive
mode, to represent the sorption behavior of binary and ternary mixtures of CH
4
, CO
2
and N
2
at
45 °C. Figure 5.2 compares the experimental observations with the model predictions for a
ternary mixture of equimolar bulk phase composition.

Figure ‎ 5.2 Experimental observations and MPTA model predictions for a ternary mixture of
CH4, CO2 and N2 at 45 °C and pressures from 5-188 bars. Excess sorption data from Ottiger et
al. (2008b).

We observe that the overall excess sorption is predicted accurately by the MPTA model but the
excess sorption of CO
2
is overestimated at the cost of both CH
4
and N
2
. In this calculation
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
0 50 100 150 200
Excess Sorption (mmol/gram)
Pressure (bar)
Total
CO2
CH4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
Calculated excess sorption
(mmol/gram)
Experimental excess sorption
(mmol/gram)
Total
CH4
N2
CO2
99

example, the excess sorption of N
2
is predicted to be negative at moderate to high pressures
causing a relatively large underestimation of the N
2
concentration in the sorbed phase. Similar
behavior is observed for binary mixtures (CO
2
-CH
4
and CO
2
-N
2
) where the total excess sorption
is reasonably accurate while the excess sorption of CO
2
is overestimated.
To improve the predictive performance of the MPTA model and to address the above
concerns regarding the pore volume used in the potential function, we attempted to re-estimate
the pure component parameters by fixing z
0
to the reported value of 0.202cc/gram. However, it
was not possible to arrive at a reasonably accurate description of the pure component isotherms
by this approach.
In the interest of integrating some of the pore volume information from the coal
characterization into the MPTA model via the potential function, we propose to divide the
potential function into a microporous section covering the pore volume from 0 to 0.07 cc/gram
and a meso- to macro-porous section covering the pore volume from 0.07 to 0.202 cc/gram. In
this initial attempt, we keep the functional form of the potential in both sections but allow the use
of different characteristic parameters in each section. This can be done by introducing

(5.17)

,
(5.18)
where b
0
represents the pore volume of the microporous region. Eq. (5.17) ensures continuity in
the chemical potential in the transition from the microporous to the mesoporous region. The new
composite potential function is used with fixed values of b
0
and z
0
set to 0.07 and 0.202 cc/gram,
respectively. In this form, a total of 8 parameters must be estimated simultaneously from the 3
pure-component isotherms relative to 5 parameters in the original version of the MPTA model.
100

Table 5.2 reports the estimated parameters for the composite potential function and Figure 5.3
compares model predictions with experimental observations for the ternary mixtures. The
composite potential function improves the overall prediction for the ternary system slightly.
Primarily, we see a better agreement between model and experiments for CO
2
at lower pressures
and for CH
4
over the entire pressure range. However, the modified potential function does not
improve the prediction of the excess sorption for N
2
to any significant extent.

Table ‎ 5.2 MPTA parameters from pure component isotherms using composite potential
function.
CH
4
CO
2
N
2

(ε
0i
I)/R (K) 600.3 722.4 485.6
(ε
0i
II)/R (K) 162.4 247.2 52.5
beta I (-) 0.688
beta II (-) 0.122
z
0
(cc/gram) 0.202
b
0
(cc/gram) 0.070

The modest improvement in the accuracy of the model predictions between the original MPTA
model and composite version, suggests that additional segments may be needed. This would
require additional information regarding the distribution of the pore volume in the meso- to
macro-porous region. Such information is, however, not available and we choose instead to
compare the MPTA predictions to those of the extended Langmuir model (ELM).

101

Figure ‎ 5.3 Experimental observations and MPTA model predictions for a ternary mixture of
CH4, CO2 and N2 at 45 °C and pressures from 5-188 bars. Excess sorption data from Ottiger et
al. (2008b).

An initial attempt to apply the ELM to predict the pure component isotherms suffered the same
problem as our initial MPTA attempt: a significant adjustment of the pore volume was required
to match the excess sorption of the individual gases. To overcome this problem, we apply the
following composite version of ELM and force the pore volumes of the micro- and meso- to
macro-porous regions to be honored.

. (5.19)
The parameters in Eq. (5.19) were obtained by regression of the pure component isotherms and
are reported in Table 5.3.

0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
0 50 100 150 200
Excess sorption (mmol/gram)
Pressure (bar)
Total
CO2
CH4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
Calculated excess sorption
(mmol/gram)
Experimental excess sorption
(mmol/gram)
Total
CH4
CO2
N2
102

Table ‎ 5.3 Model parameters for Extended Langmuir Model (ELM)
CH
4
CO
2
N
2

V
i,I
(mmol/cc) 1.21E-02 1.89E-02 5.34E-03
b
i,I
(1/atm) 0.246389 0.443316 0.600095
V
i,II
(mmol/cc) 4.57E-02 6.00E-02 3.88E-02
b
i,II
(1/atm) 2.69E-03 4.19E-03 2.94E-03
b
0
(cc/gram) 0.070
z
0
(cc/gram) 0.202

Figure 5.4 compares the experimental observations with ELM predictions for the ternary mixture
of CH
4
/ CO
2
/ N
2
.

Figure ‎ 5.4 Experimental observations and ELM predictions for a ternary mixture of CH4, CO2
and N2 at 45 °C and pressures from 5-188 bars. Excess sorption data from Ottiger et al. (2008b).

From Figure 5.4, we observe a significant difference in the predictive capabilities of ELM
relative to MPTA. The prediction of the excess sorption of CH
4
is improved at the cost of
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
0 50 100 150 200
Excess sorption (mmol/gram)
Pressure (bar)
Total
CO2
CH4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
Calculated excess sorption
(mmol/gram)
Experimental excess sorption
(mmol/gram)
Total
CH4
CO2
N2
103

underestimation of the total excess and the CO
2
excess while the N
2
excess is now slightly
overestimated. The applied version of the ELM requires 12 parameters to be estimated from pure
component data as compared to 5 parameters for the original MPTA and 8 parameters for the
composite; suggesting that a 3 stage MPTA model (with 12 adjustable parameters) might be
more valuable in terms of gains in accuracy.
5.5 Discussion and Conclusions
In the previous sections, we have discussed the numerical and application aspects of the MPTA
model and presented a comparison between calculations and experimental adsorption data from a
real coal sample. In this work, we have used the MPTA in a purely predictive mode based on
pure component isotherms alone. The model is able to accurately predict the total excess sorption
over a range of binary and ternary mixtures but is found to overestimate the CO
2
excess sorption
and underestimate the CH
4
and N
2
excess sorptions. Table 5.4 summarizes the accuracy of the
presented models for all experimental observations presented by Ottiger et al. (2008a, 2008b).
The summary reports the errors in prediction of the component excess sorption as average
absolute relative errors. The number of data points for each system and the number of adjustable
parameters in the model are listed in parenthesis.
We believe that the predictive capability of the MPTA model can be improved by improving
the interactions between species in the adsorbed phase. Using more sophisticated model (relative
to a standard cubic EOS) for density calculations in the high-pressure sorbed phase can also
enhance the agreement between the model and experimental observations.
The potential function(s) used in the MPTA model also deserve(s) additional attention as and
more flexible functional form may need to be used. In addition, a careful consideration of the
relevant pore size distribution seems essential in this context. From a comparison of the results
104

from the MPTA and the ELM models, we observe that the use of the MPTA model improves the
agreement with experimental data and trends significantly. At the same time, however, the
computational cost is substantially increased. Therefore the use of MPTA or ELM for sorption
calculations is a matter of balancing the accuracy versus CPU time. We believe that at the
current state, the MPTA is not suitable for large scale general purpose simulation of ECBM
recovery processes; while it’s suggested for standalone sorption calculations and/or simulation of
ECBM experiments at the laboratory scale.

Table ‎ 5.4 Average absolute relative error (%) in prediction of binary and ternary excess sorption

Mixture

CH
4
/N
2
(30) CH
4
/CO
2
(41) CO
2
/N
2
(40) CH
4
/CO
2
/N
2
(10)
MPTA (5 par) % % % %
Ntot 4.1 3.9 4.3 2.4
nCH4 54.3 70.8 0.0 55.0
nCO2 0.0 20.9 28.6 23.1
nN2 70.4 0.0 94.4 97.5
MPTA (8 par)
Ntot 2.6 3.8 6.5 2.9
nCH4 42.9 58.6 0.0 34.4
nCO2 0.0 15.2 24.0 20.4
nN2 62.0 0.0 91.7 96.6
ELM (12 par)
Ntot 8.4 9.0 17.8 12.7
nCH4 31.3 50.9 0.0 15.6
nCO2 0.0 15.8 38.9 30.5
nN2 31.1 0.0 81.2 78.2

105

Chapter 6 : Summary and Future Research Directions

This research project focuses on the role of mass transfer mechanisms and their impact on the
displacement dynamics in enhanced oil/gas operations. In Chapter 3, we investigated the
accuracy of some of the physical models that are frequently used to represent dispersive mixing
and mass transfer in a 1D dual-porosity system by comparing the results from compositional
simulations with experimental observations. We showed that the accuracy of our displacement
calculations relative to the experimental observations was sensitive to the selected models for
dispersive mixing, mass transfer between flowing and stagnant porosities, and IFT scaling of
relative permeability functions.
The two representations of a dual-porosity model used in our calculations lead to almost
identical results. Therefore, we were not able to conclude which model (DP I or DP II) provides
for a better representation of the unresolved heterogeneity inside the column using the Coats-
Smith model. Based on this observation, we recommend performing 2D/3D displacement
experiments by carefully defining the heterogeneity in the porous medium. The porous medium
should consist of two layers with high permeability contrast to represent a dual-porosity system.
Such experiments will allow for further testing/validating the physical models that are used in
compositional simulation of miscible displacements in dual-porosity systems, and coming up
with the simplest accurate physical models to be used in this context.
For the new set of experiments, it is suggested that steady-state relative permeability
measurements with pre-equilibrated phase compositions at two different values of IFT be
performed. This will allow one to determine the scaling exponent (β) and hence test/validate the
106

IFT scaling method suggested by Coats (1980) by comparing the experimental data and
displacement calculations in a predictive mode (i.e. without adjusting β).
As was discussed in section 2.4, several attempts by different investigators have failed in
accurately reproducing the experimental observations for miscible gas injection into dual-
porosity systems using commercial compositional simulators. Therefore combined experimental
and modeling study of multicontact miscible displacements in dual-porosity systems is needed
for investigating the accuracy of the physical models that are commonly used in compositional
simulations for dual-porosity systems.
In Chapter 4, 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. A time-dependent transfer function is used for matrix-fracture exchange in which
the shape factor is adjusted using a boost factor to differentiate between the transfer rate at early
and late times.
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 gas-oil diffusion is modeled by simply assuming an average
interface composition. It was also demonstrated that the off-diagonal diffusion coefficients and
time-dependency of matrix-fracture transfer function 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
107

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.
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 Chapter 5, we discussed the numerical and application aspects of the MPTA model and
presented a comparison between calculations and experimental adsorption data from a real coal
sample. We used the MPTA in a purely predictive mode based on pure component isotherms
alone. The model was able to accurately predict the total excess sorption over a range of binary
and ternary mixtures but was found to overestimate the CO
2
excess sorption and underestimate
the CH
4
and N
2
excess sorptions.
The predictive capability of the MPTA model can be improved by improving the interactions
between species in the adsorbed phase in future research activities. Using more sophisticated
model (relative to a standard cubic EOS) for density calculations in the high-pressure sorbed
phase can also enhance the agreement between the model and experimental observations.
The potential function(s) used in the MPTA model also deserve(s) additional attention as and
more flexible functional form may need to be used. In addition, a careful consideration of the
relevant pore size distribution seems essential in this context.

108

Nomenclature
A Cross sectional area
b Langmuir parameter
B
d
Boost factor for diffusion term
B
e
Boost factor for expansion term
B
gd
Boost factor for gravity drainage term
c
t
Total fluid compressibility
C Overall molar density
d
p
Particle diameter
D Diffusion coefficient
Đ

MS coefficient
Đ
o

Infinite dilution diffusion coefficient
f Flowing fraction, fugacity
F IFT scaling factor, molar advective flux
FF Formation factor
g Gravitational vector
H Molar diffusive flux
k
r
Relative permeability
k
rg,e
End-point relative permeability of gas
k
ro,e
End-point relative permeability of oil
K
i
K-value of component i

Effective dispersion coefficient of component i in phase j
K
l
Longitudinal dispersion coefficient
109

n

Number of moles
n
g
Corey-exponent for gas
n
o
Corey-exponent for oil
N Total molar flux
p Pressure
q Source term in the mass conservation equation
R Universal gas constant
S

Saturation
S
gc
Critical gas saturation
S
or
Residual oil saturation
t Time
T Temperature
V

Langmuir parameter
V
t
Total fluid volume
x Length of tie line (mass fraction)
x
ij
Mole fraction of component i in phase j
z Direction along packed column, porous volume
α Longitudinal dispersivity
β IFT scaling exponent, Dubinin exponent
Γ Matrix-fracture transfer rate, surface excess
Γ
d
Matrix-fracture transfer rate due to diffusion
Γ
e
Matrix-fracture transfer rate due to expansion
Γ
gd
Matrix-fracture transfer rate due to gravity drainage
110

δ

Kronecher delta function
ε

Sorption potential
θ

Mass transfer coefficient
λ

Mobility
μ Viscosity, chemical potential
v Velocity
ρ Density
σ Interfacial tension, shape factor
σ
0
Reference interfacial tension
φ Fugacity coefficient
ϕ Porosity
ψ Natural log of fugacity

111

References
Abedini, A., Mosavat, N., and Torabi, F.: Determination of minimum miscibility pressure of
crude oil–CO
2
system by oil swelling/extraction test. Energy Technology. 2, 431-439 (2014)
Abrams, D.S., Prausnitz, J.M.: Statistical thermodynamics of liquid mixtures: a new expressionn
for the excess Gibbs energy of partly or completely miscible systems. AIChE Journal. 21,
116-128 (1975)
AlHamdan, M.R., Cinar, Y., Suicmez, V.S., Dindoruk, B.: Experimental and numerical study of
compositional two-phase displacements in layered porous media. Paper SPE 147967
presented at the SPE reservoir characterization and simulation conference and exhibition,
Abu Dhabi, UAE, 9-11 October (2011)
Alvarado, V., Manrique, E.: Enhanced oil recovery: an update review. Energies. 3, 1529-1575
(2010)
Al-Wahaibi, Y.M., Muggeridge, A.H., Grattoni C.A.: Experimental and numerical studies of
gas/oil multicontant miscible displacements in homogeneous and crossbedded porous media.
SPE Journal. 12, 62-76 (2007)
Amaefule, J.O., Handy, L.L.: The effect of interfacial tensions on relative oil/water
permeabilities of consolidated porous media. SPE Journal. 22, 371-381 (1982)
Arya, A., Hewett, T.A., Larson, R.G., Lake, L.W.: Dispersion and reservoir heterogeneity. SPE
Reservoir Engineering. 3, 139-148 (1988)
Asghari, K., Torabi, F.: Effect of miscible and immiscible CO
2
flooding on gravity drainage:
Experimental and simulation results. Paper SPE 110587 presented at the SPE/DOE improved
oil recovery symposium, Tulsa, Oklahoma, USA, 19-23 April (2008)
112

Bahramian, A., Danesh, A.: Prediction of liquid-liquid interfacial tension in multi-component
systems. Fluid Phase Equilibria. 221, 197-205 (2004)
Baker, L.E.: Effects of dispersion and dead-end pore volume in miscible flooding. SPE Journal.
17, 219-227 (1977)
Barenblatt, G.I., Zheltov, 1u.P. and Kochina, I.N.: Basic concepts on the theory of seepage of
homogeneous liquids in fissured rocks. Prikladnaya Matematikai Mekhanika, Akad. Nauk,
SSSR 24, 852-864 (1960)
Batycky, R.P.: Experimental verification of MOC theory for three and four component systems.
MS thesis, Stanford University (1994)
Bear, J.: Dynamics of fluids in porous media. Dover Publications, New York (1972)
Bijeljic, B., Blunt, M.J.: A physically based description of dispersion in porous media. Paper
SPE 102869 presented at the SPE annual technical conference and exhibition, San Antonio,
Texas, USA, 24-27 September (2006)
Bretz, R.E., Orr, F.M.: Interpretation of miscible displacements in laboratory cores. SPE
Reservoir Engineering. 2, 492-500 (1987)
Burger, J.E., Mohanty, K.K.: Mass transfer from bypassed zones during gas injection. SPE
Reservoir Engineering. 12, 124-130 (1997)
Chaback, J.J., Morgan, W.D., and Yee, D.: Sorption of nitrogen, methane, carbon dioxide and
their mixtures on bituminous coals at in-situ conditions. Fluid Phase Equilibria 117, 289–296
(1996)
Chordia, M., and Trivedi, J.: Diffusion in naturally fractured reservoirs-a-review. Paper SPE
77339 presented at the SPE Asia Pacific Oil and Gas Conference and Exhibition, Brisbane,
Australia, 18-20 October (2010)
113

Clarkson, C.R.: Application of a new multicomponent gas adsorption model to coal gas
adsorption systems. SPEJ 8, 236–251 (2003)
Coats, K.H.: An equation of state compositional model. SPE Journal. 20, 363-376 (1980)
Coats, K.H.: 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 (1989)
Coats, K.H., Smith, B.D.: Dead-end pore volume and dispersion in porous media. SPE Journal.
4, 73-84 (1964)
da Silva, F.V., Belery, P.: Molecular diffusion in naturally fractured reservoirs: A decisive
recovery mechanism. Paper SPE 19672 presented at the SPE annual technical conference and
exhibition, San Antonio, Texas, USA, 8-11 October (1989)
Danesh, A.: PVT and phase behavior of petroleum fluids. Elsevier, Amsterdam (1998)
Darvish, G.R., Lindeberg, E., Holt, T., Utne, S.A., Kleppe, J.: Reservoir conditions laboratory
experiments of CO
2
injection into fractured cores. Paper SPE 99650 presented at the
SPE/EAGE annual conference and exhibition, Vienna, Austria, 12-15 June (2006)
DeGance, A.E., Morgan, W.D., and Yee, D.: High pressure adsorption of methane, nitrogen and
carbon dioxide on coal substrates. Fluid Phase Equilibria 82, 215–224 (1993)
Delshad, M., Macalister, D.J., Pope, G.A., Rouse, B.A.: Multiphase dispersion and relative
permeability experiments. SPE Journal. 25, 524-534 (1985)
Di Donato, G., Tavassoli, Z., and Blunt, M.J.: Analytical and numerical analysis of oil recovery
by gravity drainage. Journal of Petroleum Science and Engineering 54, 55-69 (2006)
114

Dreisbach, F., Staudt, R., and Keller, J.U.: High pressure adsorption data of methane, nitrogen,
carbon dioxide and their binary and ternary mixtures on activated carbon. Adsorption 5, 215–
227 (1999)
Duncan, J.B., and Toor, H.L: An experimental study of three component gas diffusion. AIChE
Journal 8, 38-41 (1962)
Fuller, E.N., Schettler, P.D., and Giddings, J.C.: A new method for prediction of binary gas-
phase diffusion coefficients. Ind. Eng. Chem. 58, 18–27 (1966)
Fuller, E.N., Ensley, K., and Giddings, J.C.: Diffusion of halogenated hydrocarbons in helium:
The effect of structure on collision cross sections. J. Phys. Chem. 73, 3678–3685 (1969).
Garcia-Flores, B.E., Trejo, A., Aguila-Hernandez, J.: Liquid–liquid phase behaviour, liquid–
liquid density, and interfacial tensions of multicomponent systems at 298. Fluid Phase
Equilibria. 255, 147-159 (2007)
Gardner, J.W., Orr, F.M., Patel, P.D.: The effect of phase behavior on CO
2
flood displacement
efficiency. SPE Journal. 21, 2067-2081 (1981)
Garmeh, G., Johns, R.T., Lake, L.W.: Pore-scale simulation of dispersion in porous media. SPE
Journal. 14, 559-567 (2009)
Ghorayeb, K., and Firoozabadi, A.: Molecular, pressure, and thermal diffusion in nonideal
multicomponent mixtures. AIChE Journal 46, 883-891 (2000)
Gilman, J.R.: An efficient finite-difference method for simulating phase segregation in the
matrix blocks in double-porosity reservoirs. SPE Reservoir Engineering 1, 403-413 (1986)
Gilman, J., and Kazemi, H.: Improvements in simulation of naturally fractured reservoirs. SPE
Journal 23, 695-707 (1983)
115

Gong, B., Karimi-Fard, M., and Durlofsky, L.: Upscaling discrete fracture characterizations to
dual-porosity, dual-permeability models for efficient simulation of flow with strong
gravitational effects. SPE Journal 13, 58-67 (2008)
Guo, X., Du, Zh., and Li, Sh.: Computer modeling and simulation of coalbed methane reservoirs.
Paper SPE 84815 presented at the SPE Eastern Regional/AAPG Eastern Section Joint
Meeting, Pittsburgh, Pennsylvania, 6–10 September (2003)
Guo, P., Wang, Z., Shen, P., and Du, J.: Molecular diffusion coefficients of the multicomponent
gas–crude oil systems under high temperature and pressure. Ind. Eng. Chem. Res. 48, 9023–
9027 (2009)
Haugen, K., Firoozabadi, A.: Composition at the interface between multicomponent non-
equilibrium phases. J. Chem. Phys. 130: 64707 (2009)
Hoteit, H.: Modeling diffusion and gas-oil mass transfer in fractured reservoirs. Journal of
Petroleum Science and Engineering 105, 1-17 (2013)
Hoteit, H., Firoozabadi, A.: Numerical modeling of diffusion in fractured media for gas-injection
and –recycling schemes. SPE Journal. 14, 323-337 (2009)
Hu, H., Whitson, C.H., and Qi, Y.: A study of recovery mechanisms in a nitrogen diffusion
experiment. Paper SPE 21893 presented at the 6th European Symposium on Improved Oil
Recovery (1991)
Irani, M., Bozorgmehry R., Pishvaie, M.R., and Tavasoli, A.: 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 (2009)
Jamili, A., Willhite, G.P., and Green, D.: Modeling gas-phase mass transfer between fracture and
matrix in naturally fractured reservoirs. SPE Journal 16, 795-811 (2011)
116

Jessen, K.: Effective algorithms for the study of miscible gas injection processes. PhD thesis,
Technical University of Denmark, Lyngby, Denmark (2000)
Jessen, K., Gerritsen, M., and Mallison, B.: High-resolution prediction of enhanced condensate
recovery processes. SPE Journal 13, 257-266 (2008)
Jessen, K., Lin, W., and Kovscek, A.R.: Multicomponent sorption modeling in ecbm
displacement calculations. Paper SPE 110258 presented at the SPE Annual Technical
Conference and Exhibition, Anaheim, California, 11–14 November (2007)
Jessen, K., Stenby, E.H., Orr, F.M.: Interplay of phase behavior and numerical dispersion in
finite-difference compositional simulation. SPE Journal. 9, 193-201 (2004)
Jessen, K., Tang, G-Q., and Kovscek, A.R.: Laboratory and simulation investigation of enhanced
coalbed methane recovery by gas injection. Transp. Porous Med. 73, 141–159 (2008)
Johns, R.T., Dindoruk, B., Orr, F.M.: Analytical theory of combined condensing/vaporizing gas
drives. SPE ATS. 1, 7-16 (1993)
Johns, R.T., Fayers, F.J., Orr, F.M.: Effect of gas enrichment and dispersion on nearly misible
displacements in condensing/vaporizing drives. SPE ATS. 2, 26-34 (1994)
Johns, R.T., Sah, P., Solano, R.: Effect of dispersion on local displacement efficinecy for
multicomponent enriched-gas floods above the MME. SPE Reservoir Evaluation and
Engineering. 5, 4-10 (2002)
Karimi-Fard, M., Gong, B., and Durlofsky, L.J.: Generation of coarse-scale continuum flow
models from detailed fracture characterizations. Water Resources Research 42, W10423
(2006)
Kazemi, H., Merrill, L.S., Porterfield, K.L., and Zeman, P.R.: Numerical simulation of water-oil
flow in naturally fractured reservoirs. SPE Journal 16, 317-326 (1976)
117

Kett, T.K. and Anderson, D.K.: Ternary isothermal diffusion and the validity of the onsager
reciprocal relations in nonassociating systems. J. Phys. Chem. 73, 1268–1274 (1969)
Kooijman, H.A., and Taylor, R.: Estimation of diffusion coefficients in multicomponent liquid
systems. Industrial & Engineering Chemistry Research 30, 1217-1222 (1991)
Koottungal, L.: Special report: 2008 worldwide EOR survey. April 21, 47-59 (2008)
Krishna, R., and Standart, G.: Determination of mass and energy fluxes for multicomponent
vapour–liquid systems. Letts Heat Mass Transfer 3, 173–182 (1976)
Ku, H.C.: Densities, viscosities, refractive indexes, and surface tensions for binary and ternary
mixtures of Tetrahydofuran, 2-propanol, and 2,2,4-Trimethylpentane. J. Chem. Eng. Data.
53, 566-573 (2008)
Lake, L.W.: Enhanced Oil Recovery. Preentice Hall, New Jersey (1989)
Lantz, R.B.: Quantitative evaluation of numerical dispersion (truncation error). SPE Journal, 11,
315-320 (1971)
Leahy-Dios, A., and Firoozabadi, A.: Unified model for nonideal multicomponent molecular
diffusion coefficients. AIChE Journal 53, 2932-2939 (2007)
Lemonnier, P., and Bourbiaux, B.: Simulation of naturally fractured reservoirs. State of the art-
Part 1–Physical mechanisms and simulator formulation. Oil & Gas Science and Technology–
Revue de l’Institut Français du Pétrole, 65, 239-262 (2010)
Li, H., Putra, E., Schechter, D.S., Grigg, R.B.: Experimental investigation of CO2 gravity
drainage in a fractured system. Paper SPE 64510 presented at the SPE Asia Pacific oil and
gas conference and exhibition, Brisbane, Australia, 16-18 October (2000)
Lim, K.T., and Aziz, K.: Matrix-fracture transfer shape factors for dual-porosity
simulators. Journal of Petroleum Science and Engineering 13, 169-178 (1995)
118

Lo, H.Y.: Diffusion coefficients in binary liquid n-alkane mixtures. J. Chem. Eng. Data. 19, 236–
241 (1974)
Lohrenz, J., Bray, B.G., and Clark, C.R.: Calculating viscosities of reservoir fluids from their
compositions. Journal of Petroleum Technology 16, 1-171 (1964)
Lu, H., Di Donato, G., and Blunt, M.: General transfer functions for multiphase flow in fractured
reservoirs. SPE Journal 13, 289-297 (2008)
Manik, J.: Compositional modeling of enhanced coalbed methane recovery. PhD Dissertation,
Pennsylvania State University, University Park, Pennsylvania, United States (1999)
Manik, J., Ertekin, T., and Kohler, T.E.: Development and validation of a compositional coalbed
simulator. JCPT 41, 39–45 (2002)
Martins, R.J., Cardoso, J.E., Barcia, O.E.: Calculation of viscosity of ternary and quaternary
liquid mixtures. Ind. Eng. Chem. Res. 40, 1271-1275 (2001)
Mason, E.A., and Malinauskas, A.P.: Gas Transport in Porous Media: The Dusty-Gas Model.
Elsevier, Amsterdam (1983)
Michelsen, L.M.: The isothermal flash problem. Part I. Stability. Fluid Phase Equilibria 9, 1–19
(1982)
Michelsen, L.M.: Saturation point calculations. Fluid Phase Equilibria 23, 181–192 (1985)
Monsalvo, M.A., and Shapiro, A.A.: Prediction of adsorption from liquid mixtures in
microporous media by the potential theory. Fluid Phase Equilibria 261, 292–299 (2007)
Monsalvo, M.A., and Shapiro, A.A.: Study of high-pressure adsorption from supercritical fluids
by the potential theory. Fluid Phase Equilibria 283, 56–64 (2009)
Moortgat, J., and Firoozabadi, A.: Fickian diffusion in discrete-fractured media from chemical
potential gradients and comparison to experiment. Energy & Fuels 27, 5793-5805 (2013)
119

Moritis, G.: Future of EOR & IOR: New companies, infrastructure, projects reshape landscape
for CO2 EOR in the US. Oil and Gas J. 99, 68–69, 72–73 (2001)
Morrow, N.R., Chatzis, I., Taber, J.J.: Entrapment and mobilization of residual oil in bead packs.
SPE Reservoir Engineering. 3, 927-934 (1988)
Orr, F.M.: Analytical Theory of Gas Injection Processes. Tie Line Publications, Denmark (2007)
Otero, J.J., Comesana, J.F., Correa, J.M., Correa, A.: Liquid−liquid equilibria of the system
Water+ 2-Propanol + 2,2,4-Trimethylpentane at 25°C. J. Chem. Eng. Data. 45, 898-901
(2000)
Ottiger, S., Pini, R., and Storti, G.: Competitive adsorption equilibria of CO
2
and CH
4
on a dry
coal. Adsorption 14, 539–556 (2008a)
Ottiger, S., Pini, R., Storti, G., and Mazzotti, M.: Measuring and modeling the competitive
adsorption of CO2, CH4, and N2 on a dry coal. Langmuir 24, 9531-9540 (2008b)
Peng, D.Y. and Robinson, D.B.: A new two-constant equation of state. Ing. Eng. Chem. Fundam.
15, 59–64 (1976)
Perkins, T.K., and Johnston, O.C.: A review of diffusion and dispersion in porous media. SPE
Journal 3, 70-84 (1963)
Peters, E.K., Gharbi, R., Afzal, N.: A look at dispersion in porous media through computed
tomography imaging. Journal of Petroleum Science and Engineering. 15, 23-31 (1996)
Polanyi, M.: Theories of the adsorption of gases: a general survey and some additional remarks.
Transactions of the Faraday Society 28, 316–333 (1932)
Poling, B.E., Prausnitz, J.M., O’Connell, J.P.: The properties of gases and liquids (5
th
edition).
McGRAW-Hill, New York (2001)
120

Quandalle, P., and Sabathier, J.C.: Typical features of a multipurpose reservoir simulator. SPE
Reservoir Engineering 4, 475-480 (1989)
Rastegar, R.: A study of dispersive mixing and flow based lumping/delumping in gas injection
processes. PhD Thesis, University of Southern California (2010)
Rastegar, R., Jessen, K.: Measurement and modeling of liquid-liquid equilibrium for ternary and
quaternary mixtures of Water, Methanol, 2-Propanol, and 2,2,4-Trimethylpentane at 293.2 K.
J. Chem. Eng. Data. 56, 278-281 (2011)
Rezaveisi, M., Ayatollahi, S., and Rostami, B.: Experimental investigation of matrix wettability
effects on water imbibition in fractured artificial porous media. Journal of Petroleum Science
and Engineering 86, 165-171 (2012)
Sahimi, M., Heiba, A.A., Hughes, B.D., Davis, H.T., Scriven, L.E.: Dispersion in flow through
porous media. Paper SPE 10969 presented at the SPE annual technical conference and
exhibition, New Orleans, Louisiana, 26-29 September (1982)
Schulze-Makuch, D.: Longitudinal dispersivity data and implications for scaling behaviour.
Ground Water. 43, 443-456 (2005)
Seto, C.J., Jessen, K., and Orr, F.M. Jr.: A multicomponent, two-phase-flow model for CO2
storage and enhanced coalbed-methane recovery. SPEJ 14, 30–40 (2009)
Shapiro, A.A., and Stenby, E.H.: Potential theory of multicomponent adsorption. Journal of
Colloid and Interface Science 201, 146–157 (1998)
Shojaei, H., Rastegar, R., and Jessen, K.: Mixing and mass transfer in multicontact miscible
displacements. Transport in Porous Media 94, 837-857 (2012)
121

Sigmund, P.: Prediction of molecular diffusion at reservoir conditions. Part 1-Measurement and
prediction of binary dense gas diffusion coefficients. Journal of Canadian Petroleum
Technology 15, 48-57 (1976)
Smith, D.H., Bromhal, G., Sams, W.N., Jikich, S., and Ertekin, T.: Simulating carbon dioxide
sequestration/ECBM production in coal seams: Effects of permeability anisotropies and the
diffusion-time constant. SPEREE 8, 156–163 (2005)
Soave, G.: Equilibrium constants from a modified Redlich-Kwong equation of state. Chemical
Engineering Science 27, 1197-1203 (1972)
Solano, R., Johns, R.T., Lake, L.W.: Impact of reservoir mixing on recovery in enriched-gas
drives above the minimum miscibility enrichment. SPE Reservoir Engineering and
Evaluation. 4, 358-365 (2001)
Soliman, K., Marschall, E.: Viscosity of selected binary, ternary, and quaternary liquid mixtures.
J. Chem. Eng. Data. 35, 375-381 (1990)
Sonier, F., Souillard, P., and Blaskovich, F.T.: Numerical simulation of naturally fractured
reservoirs. SPE Reservoir Engineering 3, 1114-1122 (1988)
Stalkup, F.I.: Miscible displacement. SPE Monograph 8, New York (1983)
Stevenson, M.D., Pinczewski, W.V., Somers, M.L., and Bagio, S.E.: Adsorption/desorption of
multicomponent gas mixtures at in-seam conditions. Paper SPE 23026 presented at the SPE
Asia-Pacific Conference, Perth, Western Australia, 4–7 November (1991)
Stevens, S.H., Spector, D., and Riemer, P.: Enhanced coalbed methane recovery using CO
2

injection: Worldwide resources and CO
2
sequestration potential. Paper SPE 48881 presented
at the SPE International Conference and Exhibition, Beijing, China, 2–6 November (1998)
122

Stiel, L.I. and Thodos, G.: The viscosity of nonpolar gases at normal pressures. AICHE J. 7,
611–615 (1961)
Stulkup, F.I.: Displacement behavior of the condensing/vaporizing gas drive processes. Paper
SPE 16715 presented at the SPE annual technical conference and exhibition, Dallas, Texas,
27-30 September (1987)
Stulkup, F.I.: Effect of gas enrichment and numerical dispersion on compositional simulator
predictions of oil recovery in reservoir condensing and condensing/vaporizing gas drives.
Paper SPE 18060 presented at the Fall meeting of SPE of AIME, Houston, Texas, 2-5
October (1988)
Taber, J.J., Martin, F.D., Seright, R.S.: EOR screening criteria revisited—Part II: Applications
and impact of oil price. SPE Reservoir Engineering. 12, 199–205 (1997)
Tanaka, Y., Matsuda, Y., Fujiwara, H., Kubota, H., Makita, T.: Viscosity of (Water+Alcohol)
mixtures under high pressure. International Journal of Thermophysics. 8, 147-163 (1987)
Taylor, R., and Krishna, R.: Multicomponent Mass Transfer. John Wiley & Sons, New York
(1993)
Trangenstein, J.A., and Bell, J.B.: Mathematical structure of compositional reservoir
simulation. SIAM Journal on Scientific and Statistical Computing 10, 817-845 (1989)
Vega, B., O’Brien, W.J., Kovscek, A.R.: Experimental investigation of oil recovery from
siliceous shale by miscible CO
2
injection. Paper SPE 135627 presented at the SPE annual
technical conference and exhibition, Florence, Italy, 19-22 September (2010)
Walsh, B.W., Orr, F.M.: Prediction of miscible flood performance: the effect of dispersion on
composition paths in ternary systems. InSitu. 14, 19-47 (1990)
123

Warren, J.E., and Root, P.J.: The behavior of naturally fractured reservoirs. SPE Journal 3, 245-
255 (1963)
Wei, X.R., Wang, G.X., and Massarotto, P.: A review on recent advances in the numerical
simulation for coalbed methane recovery process. Paper SPE 93101 presented at the Asia-
Pacific Oil and Gas Conference and Exhibition, Jakarta, Indonesia, 5–7 April (2005)
Wesselingh, J.A., and Krishna, R.: Mass Transfer. Ellis Horwood, Chichester, England (1990)
Zhou, D., Fayers, F.J., Orr, F.M.: Scaling of multiphase flow in simple heterogeneous porous
media. SPE Reservoir Engineering 12, 173-178 (1997)
Zhu, J., Jessen, K., Kovscek, A.R., and Orr, F.M. Jr.: Analytical theory of coalbed methane
recovery by gas injection. SPEJ 8, 371–379 (2003)
Zick, A.A.: A combined condensing/vaporizing mechanism in the displacement of oil by
enriched gas. Paper SPE 15493 presented at the SPE annual technical conference and
exhibition, New Orleans, Louisiana, 5-8 October (1986)
Zimmerman, R.W., Chen, G., Hadgu, T., and Bodvarsson, G.S.: A numerical dual-porosity
model with semianalytical treatment of fracture/matrix flow. Water resources research 29,
2127-2137 (1993)

124

Appendix A: UNIQUAC Model for Phase Equilibria
The Universal Quasichemical Activity Coefficient (UNIQUAC) model is an activity coefficient
model used in description of phase equilibria. In this model the activity coefficient of component
i in a multicomponent mixture is described by a combinatorial and a residual contribution
(Abrams and Prausnitz, 1975; Poling et al., 2001)

(A.1)
where

is an entropic term which quantifies the deviation from ideal solubility due to
differences in the shape of molecule and

is an enthalpic correction resulting from the change
in interacting forces between different molecules upon mixing. The combinatorial term is
calculated using the relative Van der Waals volumes r
i
and surface areas q
i
of the pure
components

(A.2)
where

is the mole fraction, n
c
is the number of components, z is the coordination number,

is
the area fraction,

is the segment fraction, and

is defined as

(A.3)
A value of 10 is frequently used for the coordination number (number of close interacting
molecules around a central molecule). The values of

and

are obtained from the following
equations
125

(A.4)

(A.5)
The residual term is calculated using Eq. (A.6)

(A.6)
Where

is defined as

(A.7)
The binary interaction parameters,

, are obtained using the experimental phase behavior data.
The binary interaction parameters together with r
i
and q
i
for the analog fluid system used in this
work are given in Table A.1 (Rastegar and Jessen 2011).
Table ‎ A.1 Parameters for UNIQUAC model (interactions and structural parameters)
Comp Water MeOH IPA iC
8

Water 0 -760.7 301.2 429.0 1.4000 0.9200
MeOH 2175.8 0.0 146.3 11.5 1.4320 1.4311
IPA -241.1 -143.2 0.0 -83.2 3.1240 3.2491
iC
8
2512.6 598.5 169.0 0.0 5.0080 5.8463

126

Appendix B: UNIQUAC Viscosity Model
The UNIQUAC viscosity model has been developed based on the Eyring’s theory of viscous
flow and the UNIQUAC phase behavior model (Martins et al., 2001). In this model, the mixture
viscosity is calculated from

(B.1)
where is the mixture viscosity,

is the viscosity of component i,

is the mole fraction of
component i,

is the molar volume of component i,

is the molar volume of the mixture, z is
the coordination number, q
i
is the surface area,

is the segment fraction (obtained from Eq.
(A.4)),

is the area fraction (obtained from Eq. (A.5)), and

is the binary interaction
parameter between component k and i. The model parameters for the analog fluid system used in
this work are given in Table 3.2.

127

Appendix C: Numerical Dispersion
The truncation errors in finite-difference numerical calculations produce numerical diffusion.
Lantz (1971) showed that the numerical Peclet number for single-phase displacements can be
estimated from

(C.1)
where and represent dimensionless distance and time, respectively. For simulations with
1000 grid locks, numerical Peclet numbers of 2636 and 2196 are estimated for experiments I and
II respectively.
The numerical Peclet number for two-phase displacements can be estimated from (Walsh and
Orr, 1990; Orr, 2007)

(C.2)
where df
1
/dS
1
represents the slope of fractional flow curve. We use 3.33 and 3.13 for df
1
/dS
1
(the
maximum slope estimated from immiscible fractional flow curve) for experiments III and IV,
respectively. For simulations with 1000 grid blocks, numerical Peclet numbers of 3057 and 2611
are estimated for experiments III and IV, respectively.
The estimated values of numerical Peclet number for simulation of all experiments are one
order of magnitude larger that the physical Peclet number of approximately 200 observed in our
experiments. Hence, numerical artifacts should be marginal.
128

Appendix D: Diffusion Coefficients
Fickian diffusion coefficients can be obtained from MS coefficients which are less sensitive to
composition changes and hence easier to measure in the lab (Taylor and Krishna 1993). The
relationship between Fickian diffusion coefficients and MS coefficients are given in section 4.3.1
of the manuscript (Eqs. 4.10 through 4.12).
Multicomponent MS coefficients are commonly calculated using binary infinite dilution
diffusion coefficients (see Eq. 3.13). In this work we use the correlation by Leahy-Dios and
Firoozabadi (2007) for infinite dilution diffusion coefficients for both vapor and liquid phases as
explained in the following paragraphs.
By performing non-linear least-squares minimization on various relationships, Leahy-Dios
and Firoozabadi (2007) found the following expression to best describe the experimental data for
infinite dilution diffusion coefficients:

(D.1)
where Đ
o
21
is the diffusion coefficient of component 2 infinitely diluted in component 1; ρ is the
molar density (mol/m
3
) of component 1; μ is the viscosity (Pa.s) of component 1; (ρĐ)
0
and μ
0
are
the dilute gas density-diffusion coefficient product (mol/m.s) and viscosity (Pa.s) respectively; T
r

and p
r
are the reduced temperature and pressure respectively; and ω is the acentric factor.
The parameters A
0
to A
1
are given by

(D.2)
129

where a
1
= ─ 0.0472; a
2
= 0.0103; a
3
= ─ 0.0147; a
4
= ─ 0.0053; a
5
= ─ 0.337; a
6
= ─ 0.1852; a
7
= ─ 0.1914.
The approach of Fuller et al. (1966, 1969) is used to calculate the dilute gas density-diffusion
coefficient product

(D.3)
where M is the molecular weight (g/mol); T is the absolute temperature (K); and Συ is the so-
called “diffusion volume increments” which is calculated by summing the atomic diffusion
volumes (see Poling et al., 2001).
The low-pressure viscosity for each component is calculated using the correlation by Stiel
and Thodos (1961)

(D.4)
where

(D.5)
The dilute gas viscosity of the mixture is obtained using a weighted average of dilute gas
viscosity of the components

(D.6)
To validate our diffusion coefficient calculations (Eqs. 4.10 through 4.13, and D.3 through D.6),
we compare our calculations with experimental data for binary and ternary mixtures. We note
that SRK equation of state (Soave, 1972) with volume translation is used for density
130

calculations; and viscosities are obtained using LBC correlation (Lohrenz et al., 1964). We also
note that similar comparisons can be found in Leahy-Dios and Firoozabadi (2007).
D.1 Example (Binary Mixtures)
Figure D.1 compares the calculated Fickian diffusion coefficients with experimental data from
Sigmund (1976) for methane-propane binary gas mixtures at different temperatures and
pressures. Very good agreements between the calculated and experimental values are observed.

Figure ‎ D.1 Calculated diffusion coefficient (solid lines) compared with experimental data
(circles) from Sigmund (1976) for methane-propane mixtures at different values of pressure,
temperature and composition.

131

Calculated Fickian diffusion coefficients are compared with experimental data from Lo (1974)
for binary liquid mixtures of alkanes (n-C
7
, n-C
10
, n-C
12
, nC
8
and nC
14
) at 298.15 K and 0.1 MPa
in Fig. D.2. The agreements between calculated and experimental values are substantial.

Figure ‎ D.2 Calculated diffusion coefficient (solid lines) compared with experimental data
(circles) from Lo (1974) for binary liquid mixtures of alkanes at 298.15 K and 0.1 MPa with
varying compositions.

D.2 Example (Ternary Mixture)
Table D.1 compares the calculated Fickian diffusion coefficients with experimental data from
Kett and Anderson (1969) for a ternary liquid mixture of nC
6
(33.3 mole %), nC
12
(35 mole %)
132

and nC
16
(31.7 mole %) at 298.15 K and 0.1 MPa. Reasonable agreements between the
calculated and experimental values are observed.

Table ‎ D.1 Calculated diffusion coefficient compared with experimental data from Kett and
Anderson (1969) for a ternary liquid mixture at 298.15 K and 0.1 MPa.
D
ij
(10
-9
m
2
/sec)
nC
12
(1), nC
16
(2), nC
6
(3)
Subscripts i j Exp. Model ∆ (%)
1 1 0.968 1.05 8
1 2 0.266 0.143 46
2 1 0.225 0.108 52
2 2 1.031 0.988 4

133

Appendix E: Shape Factors for Gas-Oil Diffusion
Consider a matrix block that is surrounded by fractures (Fig. E.1). Assume matrix and fractures
are filled with different phases (e.g. oil in the matrix and gas in the fractures); and there is gas-oil
equilibrium at the fracture-matrix interface at all six faces of the matrix block. Also assume the
fluids in all fractures surrounding the matrix block have the same pressure and composition. The
transfer rate of each component due to molecular diffusion from the center of matrix block to the
matrix-fracture interface at each face is given by

Figure ‎ E.1 Schematic of a matrix block surrounded by fractures

(E.1)
134

The total transfer rate of each component from the matrix block to the matrix-fracture interface is
obtained by summing the transfer rates in Eq. (E.1)

(E.2)
Dividing both sides of Eq. (E.2) by

will give us the transfer rate per unit volume of
matrix block

(E.3)
where σ
1
is given by

(E.4)
The transfer rate of each component due to molecular diffusion from the center of each fracture
to the matrix-fracture interface is given by
135

(E.5)
The total transfer rate of each component from the fracture to the matrix-fracture interface is
obtained by summing the transfer rates in Eq. (E.5)

(E.6)
Dividing both sides of Eq. (E.6) by

will give us the transfer rate per unit volume of matrix
block

(E.7)
where σ
2
is given by

(E.8)
136

Based on the film theory, the transfer rates given by Eqs. (E.3) and (E.7) must be equal. This
leads to Eq. (27). Since the molecular diffusion equation is analogous to the pressure diffusion
equation, one can replace 4 with π
2
in Eqs. (E.4) and (E.8) and multiply the transfer function by a
boost factor based on the derivations of Zimmerman et al. (1993) and Lim and Aziz (1995).

137

Appendix F: Transfer Rate due to Gravity Drainage
The transfer rate of a component i from matrix to fracture due to gravity drainage is given by

(F.1)
where T
jk
represents the displacement of phase j by phase k due to gravity drainage which acts
only in the vertical direction. Since we are looking at two-phase gas-oil systems, we need to
calculate T
go
and T
og
. Lu et al. (2008) proposed the following equations based on the
observations of Di Donato et al. (2006)

(F.2)
where B
gd
is the correction factor for the gravity drainage term; S
gm,init
is the initial gas saturation
in the matrix block; S
*
gm
is the final gas saturation in the matrix at the end of displacement; b is
the exponent of the oil relative permeability function; and β
go
and F(S
gf
) are parameters that are
defined as follows

(F.3)
where λ
gm
, λ
om
and λ
tm
are the gas, oil and total mobilities in the matrix respectively; ∆ρ
mgo
is the
difference between mass densities of gas and oil; g is the acceleration of gravity; and h denotes
the hight of the matrix block. Since the gas is much more mobile than the oil, λ
gm
λ
om
/λ
tm
can be
approximated by

where

is the maximum oil relative permeability.
The parameter F(S
gf
) is given by

(F.4)
138

where S
gf
is the gas saturation in the fracture; and K
f
is the fracture permeability. Since the
gravity drainage acts in the vertical direction, vertical permeability is used in Eqs. (F.3) and
(F.4).
The correction factor for the gravity drainage term is given by

(F.5)
We note that the correction factor for gravity drainage has a value less than unity and will
approach unity at late times. We also note that T
og
= ─T
go
for incompressible flow settings (Lu et
al., 2008).

Abstract (if available)

Abstract In order to estimate the potential incremental hydrocarbon recovery by gas injection, compositional reservoir simulators are commonly used in the industry. Successful design and implementation of gas injection processes (mainly CO₂) rely in part on the accuracy by which the available simulation tools can represent the physics that govern the displacement behavior in a reservoir. ❧ In the first part of this thesis, we investigate the accuracy of some physical models that are frequently used to describe and interpret dispersive mixing and mass transfer in compositional reservoir simulation. The calculations from compositional simulation are compared with the experimental observations. ❧ A quaternary analog fluid system (alcohol‐water‐hydrocarbon) that mimics the phase behavior of CO₂‐hydrocarbon mixtures at high pressure and temperature has been designed in our research group. A porous medium was designed using Teflon materials to ensure that the analog oil acts as the wetting phase, and the properties of the porous medium were characterized in terms of porosity and permeability. Relative permeability and interfacial tension measurements were also performed to delineate interactions between the fluid system and the porous medium. Displacement experiments at First‐contact miscible (FCM), near‐miscible and multicontact miscible (MCM) conditions were consequently performed (Rastegar, 2010). ❧ The effluent concentrations from two‐component FCM displacement experiments exhibit a tailing behavior that is attributed to imperfect sweep of the porous medium: A feature that is not captured by normal dispersion models. To represent this behavior in displacement calculations, we use dual‐porosity (DP) models including mass transfer between flowing and stagnant porosities. The 4‐component two‐phase displacement experiments (near‐miscible and MCM) are consequently simulated using the DP models constructed based on observations from FCM displacements. ❧ We demonstrate that the accuracy of our displacement calculations relative to the experimental observations is sensitive to the selected models for dispersive mixing, mass transfer between flowing and stagnant porosities, and IFT scaling of relative permeability functions. We also demonstrate that numerical calculations substantially agree with the experimental observations for some physical models with limited need for model parameter adjustment. ❧ The second part of this thesis is devoted to diffusion and matrix‐fracture interactions 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. ❧ In this work, we present a dual‐porosity model in which the generalized Fick's law is used for molecular diffusion to account for the dragging effects

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Modeling and simulation of complex recovery processes

PDF

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

PDF

Asphaltene deposition during co₂ injection in enhanced oil recovery applications

PDF

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

PDF

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

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

PDF

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

PDF

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

PDF

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

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

PDF

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

PDF

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

PDF

Optimization of coupled CO₂ sequestration and enhanced oil recovery

PDF

Reactivation of multiple faults in oilfields with injection and production

PDF

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

PDF

A method for characterizing reservoir heterogeneity using continuous measurement of tracer data

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

Real-time reservoir characterization and optimization during immiscible displacement processes

PDF

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

PDF

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

PDF

Investigation of gas sorption and mass transfer in gas shales

Asset Metadata

Creator Shojaei, Hasan (author)

Core Title Mass transfer during enhanced hydrocarbon recovery by gas injection processes

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Publication Date 07/21/2014

Defense Date 05/30/2014

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

Tag CO₂ injection,dispersion,enhanced oil recovery,fractured reservoirs,molecular diffusion,OAI-PMH Harvest,reservoir simulation

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Jessen, Kristian (committee chair), Ershaghi, Iraj (committee member), Ghanem, Roger G. (committee member)

Creator Email hasan.shojaei@gmail.com,shojaei@usc.edu

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

Unique identifier UC11287704

Identifier etd-ShojaeiHas-2725.pdf (filename),usctheses-c3-446031 (legacy record id)

Legacy Identifier etd-ShojaeiHas-2725.pdf

Dmrecord 446031

Document Type Dissertation

Format application/pdf (imt)

Rights Shojaei, Hasan

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