University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content UNIVERSITY OF SOUTHERN CALIFORNIA

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

By
Mohammad Evazi Yadecuri

PHD DISSERTATION

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

i

Abstract
Oil and gas reservoirs in fluvial and deltaic depositional systems comprise a significant
portion of the so-called conventional resources. In these systems, the inherent reservoir
heterogeneity is often the major challenge of exploitation that complicates the application
and simulation of enhanced oil recovery techniques. Despite the advancements in the
computational technology, fine-scale simulation of stochastically generated geomodels is
not feasible as hundreds of realizations are analyzed to quantify the impact of geological
uncertainty on the reservoir performance. This has motivated the development of
upscaling techniques that, for highly heterogeneous permeability fields, still remains a
challenge.
In the first part of this work, the effect of heterogeneity on displacement efficiency
in coarse scale modeling is studied. Pore space is ranked based on flow contribution to two
levels of porosity and a dual-porosity dual-permeability flow model is adapted for coarse-
scale flow simulation. We use fine-scale streamline information to transform
heterogeneous geomodels into a dual-continuum coarse model that preserves the global
flow pathways adequately.
The proposed technique is applied and tested on two heterogeneous models with
different types of fluid flow modeling (compositional and black oil). In order to evaluate
the performance of the proposed technique, displacement calculations are performed on
the original fine grid and on a uniform coarse grid with single-porosity and dual-porosity
upscaling. We demonstrate that dual-porosity coarse models predict the breakthrough
time accurately and reproduce the post-breakthrough response adequately while single-
porosity models overestimate arrival time and oil recovery. By preserving large scale
ii

heterogeneities, dual-porosity coarse models are demonstrated to be significantly less
sensitive to the level of upscaling, compared to conventional single-porosity upscaling. This
makes the proposed upscaling approach a relevant and suitable technique for upscaling of
highly heterogeneous geomodels.
Modeling of multi-fractured horizontal wells is the subject of the second part of this
dissertation. Recovery from low and ultra-low permeability reservoirs have been unlocked
by the introduction of horizontal well and hydraulic fracturing technology that has
increased the world’s hydrocarbon reserve significantly. In order to maximize the
productivity from these reservoirs, numerical simulation is often used to investigate the
optimized well spacing and stimulation design. Fracture properties vary during the well life
and must be properly modeled to enable realistic production simulation. Production
history matching of a real-life multi-fractured horizontal well in a liquid-rich Monterey
formation reservoir is considered. We infer from numerical simulations that, in these
systems, time-dependent fracture conductivity model is required to match early-time flush
production and later-time low flow rates.
Accurate representation of the near-well dynamics is key to the realistic simulation
of fractured wells. Due to the large size of these reservoirs, efficient numerical simulation
requires use of large gridblocks that may not be adequate to capture the near-well
dynamics. We adapt a multi-phase upscaling workflow for multi-fractured horizontal wells
under natural depletion. Fine-scale simulation of the depletion process on a local well
model is used to calculate a new set of gas-oil pseudo-relative permeability curves for
fractured and non-fractured gridblocks. Simulation results demonstrate improvements in
the accuracy of the coarse scale simulation over existing approaches.
iii

Dedicated to my parents,
Zeinab and Abbasgholi

iv

Acknowledgements

Many people have helped me in the inception and completion of this work. Dr. Kristian
Jessen has been an inspiration and a phenomenal guru that has supported me at desperate
moments, challenged me as a critical thinker and guided me towards becoming an
independent researcher. My gratitude and thanks for him is immense.
I would like to thank Dr. Iraj Ershaghi, director of petroleum engineering program,
for his insightful advices and persistent support throughout my studies at USC. My
appreciation goes also to my PhD dissertation committee, Dr. Iraj Ershaghi and Dr. Roger
Ghanem, as well as my PhD qualifying committee, Dr. Fred Aminzadeh and Dr. Catherine
Shing.
I am truly grateful for receiving Mork Family Department graduate scholarship and
research assistantship that eased me graduate life. I have gained significant industrial
experience through summer internships in Occidental Oil and Gas Corporation (OXY). I am
very thankful of OXY and all my mentors and managers that helped me learn and grow.
Good friends have turned my PhD life into a very delightful experience. I would like
to specially thank Mohsen Taghavifar, Arman Khodabakhshnejad, Hamed Darabi and
Shahram Farhadi for always being there and sharing memorable fun times. I would also
like to appreciate friendship and insightful discussions with Ehsan Tajer, Reza Rastegar,
Nelia Jafroodi, Hamid Reza Jahangiri, Parham Ghods, Dalad Nattwongasem, Mohammad
Javaheri, Tayeb Ayat, Asal Rahimi, Farnaz Eskandari, Hassan Shojaei, Cyrus Ashayeri,
Siavash Hakim Elahi, Mehran Hoseini, Leonardo Gomes Oliviera, Christian Keller, Azarang
Golmohammadi, Marjan Sherafati, Devang Dasani and Hassan Dashtian.
v

Last but not the very least, my endless love and gratefulness goes to my family. I am
indebted to Zeinab and Abbasgholi, my mother and father, for their everlasting love, care
and support. My brothers, Ali, Hamid and Vahid and my sister, Tahereh, deserve matchless
acknowledgement for being source of nonstop encouragement and support.

vi

Contents
Abstract ......................................................................................................................................... i
Acknowledgements ................................................................................................................. iv
List of figures ............................................................................................................................. ix
List of tables ............................................................................................................................ xiii
Chapter 1: Introduction .......................................................................................................... 2
1.1 Upscaling of geological models ..................................................................................................... 3
1.1.1 Problem statement and research approach ........................................................................ 5
1.2 Coarse scale modeling of multi-fractured horizontal wells............................................... 8
1.2.1 Problem statement and research approach ........................................................................ 9
1.3 Organization of the manuscript ................................................................................................. 11
Chapter 2: Upscaling of flow in porous media .............................................................. 14
2.1 Governing equations ...................................................................................................................... 14
2.1.1 Immiscible, incompressible two-phase flow ................................................................... 15
2.2 Upscaling of flow in porous media ........................................................................................... 17
2.2.1 Calculation of fluid properties ............................................................................................... 18
2.2.2 Porosity upscaling ...................................................................................................................... 19
2.2.3 Upscaling of absolute permeability ..................................................................................... 19
2.2.4 Multi-phase flow upscaling .................................................................................................... 25
2.2.5 Upscaling-downscaling ............................................................................................................ 27
2.2.6 Multi-scale techniques ............................................................................................................. 27
2.3 Upscaling errors .............................................................................................................................. 28
2.3.1 Errors in upscaling of the pressure equation .................................................................. 29
2.3.2 Errors in the upscaling of the saturation equation ....................................................... 29
2.3.3 What is the optimum level of upscaling? .......................................................................... 30
Chapter 3: Dual-porosity upscaling of heterogeneous reservoirs .......................... 33
3.1 Dual-porosity upscaling ............................................................................................................... 34
vii

3.2 Modeling of dual porosity systems .......................................................................................... 35
3.2.1 Dual-porosity flow modeling ................................................................................................. 36
3.2.2 Dual-porosity dual-permeability flow modeling ........................................................... 37
3.3 Coarse-scale two-porosity representation of heterogeneous reservoirs ................. 38
3.3.1 Local porosity splitting ............................................................................................................ 41
3.3.2 Global porosity splitting .......................................................................................................... 45
3.4 Calculation of coarse grid properties ...................................................................................... 48
Chapter 4: Results from dual-porosity upscaling ........................................................ 53
4.1 Test Models ....................................................................................................................................... 54
4.1.1 Model I: Channelized Model ................................................................................................... 55
4.1.2 Model II: 10
th
SPE Comparative Solution Project, layer 37
th
..................................... 56
4.2 Comparative study of single-phase upscaling techniques .............................................. 57
4.3 Upscaling of an immiscible displacement ............................................................................. 61
4.3.1 Test model I .................................................................................................................................. 62
4.3.2 Test model II ................................................................................................................................ 67
4.4 Upscaling of a near-miscible displacement .......................................................................... 73
4.4.1 Test model I .................................................................................................................................. 76
4.4.2 Test model II ................................................................................................................................ 81
4.5 Concluding remarks ....................................................................................................................... 85
Chapter 5: Modeling of multi-fractured horizontal wells ......................................... 88
5.1 Unconventional reservoir characterization.......................................................................... 89
5.2 Hydraulic fracture characterization ........................................................................................ 90
5.2.1 Hydraulic fracture geometry ................................................................................................. 91
5.2.2 Hydraulic fracture conductivity ........................................................................................... 91
5.3 Predictive tools for production ................................................................................................. 95
5.3.1 Decline curve analysis .............................................................................................................. 95
5.3.2 Analytic solution for fractured horizontal wells ............................................................ 96
5.3.3 Numerical simulation ............................................................................................................... 97

viii

Chapter 6: Coarse-scale simulation of multi-fractured wells ................................ 102
6.1 Governing equations .................................................................................................................... 104
6.2 Multi-phase flow upscaling under depletion condition ................................................. 106
6.2.1 Workflow for multiphase flow upscaling of multi-fractured well ........................ 108
6.3 Coarse scale simulation of horizontal well with a single fracture ............................. 109
6.3.1 Explicit simulation of hydraulic fracture ........................................................................ 112
6.3.2 Coarse grid with local grid refinement ............................................................................ 112
6.3.3 Coarse grid with transmissibility modification ............................................................ 114
6.3.4 Coarse grid with multi-phase upscaling .......................................................................... 114
6.4 A history matching study ........................................................................................................... 118
6.4.1 Multi-fractured horizontal well history-matching case study ................................ 120
6.5 Concluding remarks ..................................................................................................................... 125
Chapter 7: Future directions ............................................................................................ 128
Bibliography .......................................................................................................................... 131

ix

List of figures

Figure ‎ 2-1: Schematic showing (a) fine and coarse scale grids, (b) fine scale local region
and (c) fine scale extended local region (from Chen et al., 2003) .................................................. 22
Figure ‎ 3-1: Idealization of a fractured porous system (from Warren and Root, 1963) ......... 36
Figure ‎ 3-2:Permeability map and streamline distribution in a five-spot pattern. (a)
Permeability distribution, logarithmic (b) Normalized map of streamline density (c)
Normalized map of inverse of streamline time-of-flight (d) Streamline Index map ............... 40
Figure ‎ 3-3: Local porosity splitting based on span of streamline index (a) Primary (black)
and secondary (white) porosity distribution in fine level with span cut-off=0.3 (d) Primary
porosity ratio in coarse grid .......................................................................................................................... 43
Figure ‎ 3-4: Local porosity splitting based on span of streamline index (a) Primary (black)
and secondary (white) porosity distribution in fine level with span cut-off=0.2 (d) Primary
porosity ratio in coarse grid .......................................................................................................................... 43
Figure ‎ 3-5: Local porosity splitting based on span of streamline index (a) Primary (black)
and secondary (white) porosity distribution in fine level with span cut-off=0.1 (d) Primary
porosity ratio in coarse grid .......................................................................................................................... 44
Figure ‎ 3-6: Water-cut simulated by fine grid, single-porosity (SP) coarse model and dual-
porosity (DP) models generated by local splitting based on span of SI (SC) ............................. 44
Figure ‎ 3-7: Dynamic F-C diagram based on streamline data for the example shown in
Figure 3.2 .............................................................................................................................................................. 46
Figure ‎ 3-8: Global porosity splitting based on dynamic F-C diagram with thresholding of 0.7
(a) Primary (black) and secondary (white) porosity distribution in fine level (d) Primary
porosity ratio in coarse grid .......................................................................................................................... 47
Figure ‎ 3-9: Global porosity splitting based on dynamic F-C diagram with thresholding of 0.5
(a) Primary (black) and secondary (white) porosity distribution in fine level (d) Primary
porosity ratio in coarse grid .......................................................................................................................... 47
Figure ‎ 3-10: Global porosity splitting based on dynamic F-C diagram with thresholding of
0.3 (a) Primary (black) and secondary (white) porosity distribution in fine level (d)
Primary porosity ratio in coarse grid ......................................................................................................... 48
Figure ‎ 3-11: Water-cut simulated by fine grid, single-porosity (SP) coarse model and dual-
porosity (DP) models generated by global splitting based on flow capacity cut-off from
dynamic F-C curve (FC) ................................................................................................................................... 48
Figure ‎ 3-12: Illustration of inter-block and intra-block flow for two coarse gridblocks ....... 50

x

Figure ‎ 4-1: Permeability map for test model 1 (logarithm of permeability in mD) ................ 55
Figure ‎ 4-2: Porosity and logarithmic permeability (mD) distribution of test model 2 .......... 56
Figure ‎ 4-3: Transmissibility calculated by local, global and adaptive local-global upscaling
................................................................................................................................................................................... 57
Figure ‎ 4-4: Oil production rate calculated by reference fine grid and coarse grid with
different transmissibility upscaling techniques. .................................................................................... 58
Figure ‎ 4-5: Oil production rate error calculated by reference fine grid and coarse grid with
different transmissibility upscaling techniques. .................................................................................... 59
Figure ‎ 4-6: Recovery factor and average reservoir pressure calculated by reference fine
grid and coarse grid with different transmissibility upscaling techniques ................................. 60
Figure ‎ 4-7: Saturation map at different injection times predicted by fine and three single
porosity coarse models .................................................................................................................................... 60
Figure ‎ 4-8: Cumulative oil production resulting from water-flood simulation on test model I
(fine and coarse grid) ....................................................................................................................................... 64
Figure ‎ 4-9: Producing water cut resulting from water-flood simulation on test model I ( fine
and coarse grid) .................................................................................................................................................. 64
Figure ‎ 4-10: Well oil production rates simulated by fine and coarse models ............................ 64
Figure ‎ 4-11: Saturation distribution predicted by fine grid and coarse models ....................... 65
Figure ‎ 4-12: : Effect of upscaling level on performance of single-porosity and dual porosity
coarse models: Fractional water production .......................................................................................... 66
Figure ‎ 4-13: Effect of upscaling level on performance of single-porosity and dual porosity
coarse models: Oil production ...................................................................................................................... 67
Figure ‎ 4-14: Permeability map and streamline distribution in a five-spot pattern on test
model 2. (a) Permeability distribution, logarithmic (b) Normalized map of streamline
density (c) Normalized map of inverse of streamline time-of-flight (d) SI map ...................... 68
Figure ‎ 4-15: Global porosity splitting based on F-C diagram with thresholding of 0.4 (a)
Primary (black) and secondary (white) porosity distribution in fine level (d) Primary
porosity ratio in coarse grid .......................................................................................................................... 69
Figure ‎ 4-16: Water fractional production simulated on test model 2 .......................................... 70
Figure ‎ 4-17: Total oil production simulated on test model 2 ........................................................... 70
Figure ‎ 4-18: Average pressure simulated on test model 2 ................................................................ 70
Figure ‎ 4-19: Oil production rate simulated on test model 2 ............................................................. 70
Figure ‎ 4-20: : Saturation distribution predicted by fine grid and coarse grid models ........... 71

xi

Figure ‎ 4-21: Error in prediction of oil saturation by single-porosity and dual-porosity
coarse models ...................................................................................................................................................... 71
Figure ‎ 4-22: : Effect of upscaling level on performance of single-porosity and dual porosity
coarse models models in test case 2: Fractional water production ............................................... 72
Figure ‎ 4-23: : Effect of upscaling level on performance of single-porosity and dual porosity
coarse models in test case 2: Oil production rate .................................................................................. 73
Figure ‎ 4-24: Cumulative oil production simulated by fine and coarse grids of test model I
................................................................................................................................................................................... 77
Figure ‎ 4-25: Production gas-oil ratio simulated by fine and coarse grids of test model I ... 77
Figure ‎ 4-26: Oil production rate simulated by fine and coarse grids of test model I ............. 77
Figure ‎ 4-27: Accuracy of coarse models in the prediction of oil saturation: Root-mean
square error ......................................................................................................................................................... 77
Figure ‎ 4-28: Near-miscible displacement; saturation map at different pore volume injected
simulated with fine and coarse models. .................................................................................................... 78
Figure ‎ 4-29: Effect of upscaling level on the accuracy of single-porosity and dual-porosity
uspcaled models in the simulation of oil production rate .................................................................. 80
Figure ‎ 4-30: Effect of upscaling level on single-porosity and dual-porosity uspcaled models
in the simulation of producing gas-oil ratio ............................................................................................ 80
Figure ‎ 4-31: Production gas-oil ratio simulated by fine and coarse grids of test model II .. 81
Figure ‎ 4-32: Cumulative oil production simulated by fine and coarse grids of test model II
................................................................................................................................................................................... 81
Figure ‎ 4-33: Oil saturation map at different times simulated by fine and coarse grids: ....... 82
Figure ‎ 4-34: Error in prediction of saturation by single-porosity and dual-porosity coarse
models .................................................................................................................................................................... 82
Figure ‎ 4-35: Effect of upscaling level on the accuracy of single-porosity and dual-porosity
uspcaled models in the simulation of oil production rate .................................................................. 83
Figure ‎ 4-36: Effect of upscaling level on the accuracy of single-porosity and dual-porosity
uspcaled models in the simulation of production gas-oil ratio ........................................................ 84
Figure ‎ 6-1: Cross sectional view of a fractured well in a domain with 9 coarse blocks. ...... 106
Figure ‎ 6-2: Gas-oil relative permeability in the test model ............................................................. 110
Figure ‎ 6-3: Water-oil imbibition relative permeability in the test model ................................. 110
Figure ‎ 6-4: Step-wise bottom-hole pressure reduction in the depletion simulation of
homogeneous test model .............................................................................................................................. 110

xii

Figure ‎ 6-5: Oil PVT Properties used in the test model ...................................................................... 111
Figure ‎ 6-6: Gas PVT Properties used in the test model ..................................................................... 111
Figure ‎ 6-7: Oil production rate simulated with explicit and LGR representation of fracture
................................................................................................................................................................................. 113
Figure ‎ 6-8: Production gas-oil ratio simulated with explicit and LGR representation of
fracture ................................................................................................................................................................. 113
Figure ‎ 6-9: Relative permeability in the well gridblock ................................................................... 115
Figure ‎ 6-10: Relative permeability in X direction in all gridblocks.............................................. 115
Figure ‎ 6-11: Relative permeability in Y direction in non-fractured gridblocks ...................... 116
Figure ‎ 6-12: Relative permeability in Y direction in fractured gridblocks ............................... 116
Figure ‎ 6-13: Relative permeability in Z direction in non-fractured gridblocks ...................... 116
Figure ‎ 6-14: Relative permeability in Z direction in fractured gridblocks ................................ 116
Figure ‎ 6-15: Oil production rate simulated with explicit and coarse representation of
fracture with transmissibility alteration (TRAN) and with multi-phase upscaling (MPU) 117
Figure ‎ 6-16: Production gas-oil ratio simulated with explicit and coarse representation of
fracture with transmissibility alteration (TRAN) and with multi-phase upscaling (MPU).
................................................................................................................................................................................. 118
Figure ‎ 6-18: J-function from capillary pressure data of the layers in the case study ........... 121
Figure ‎ 6-19: Oil and gas PVT properties of the field example ........................................................ 122
Figure ‎ 6-20: Stress-dependent rock compressibility in the case study ...................................... 122
Figure ‎ 6-21: Stress-dependent rock porosity and permeability input to the simulator ...... 122
Figure ‎ 6-23: Producing water cut observation and simulated values ........................................ 124
Figure ‎ 6-24: Well bottom-hole pressure from simulation ............................................................... 125
Figure ‎ 6-25: Producing gas-oil ratio simulated and observed on the well ............................... 125

xiii

List of tables

Table ‎ 4-1: Corey-type relative permeability parameters for oil and water ................................ 62
Table ‎ 4-2: Initial fluid composition ............................................................................................................. 76
Table ‎ 4-3: Injection fluid composition ....................................................................................................... 76
Table ‎ 4-4: Physical properties of hydrocarbon mixture components for PR-EOS ................... 76
Table ‎ 6-1: Homogeneous test model properties ................................................................................. 109
Table ‎ 6-2: Petrophysical properties of the case study ...................................................................... 120

1

Chapter 1
Introduction

2

1 Introduction

This dissertation discusses coarse scale multi-phase flow simulation in highly
heterogeneous reservoirs and multi-fractured horizontal wells. Hydrocarbon reservoirs in
fluvial and deltaic depositional systems comprise significant portion of the so-called
conventional resources where reservoir heterogeneity is one of the major challenges of
exploitation. Fine scale multi-phase flow simulation in these reservoirs is computationally
expensive and upscaling is inevitable. However, significant challenges arise when enhanced
oil recovery techniques are simulated on the upscaled models. The first part of this work is
concerned with upscaling of transport processes in highly heterogeneous reservoirs where
a majority of upscaling techniques fail in reproducing the fine scale results.
The second part of the dissertation considers the coarse grid simulation of horizontal
wells with multiple hydraulic fractures. Low permeability and unconventional resources
provide significant increase to the total hydrocarbon reserves of the world that has been
producing commercially by the application of horizontal well and hydraulic fracturing
technology. Several questions regarding the production optimization from these resources
arise, such as optimum well spacing and stimulation design, which can be addressed via
reservoir simulation studies. Use of locally refined grids has been successfully applied to
capture near-fracture dynamics at single well or sector well modeling studies. This
becomes very tedious and computationally inefficient when field-scale studies are

3

addressed. Significant challenges arise as induced fractures are incorporated into the
coarse scale numerical simulation. This topic is investigated as the second part of this work.
1.1 Upscaling of geological models
Reservoir simulation is an invaluable tool in reservoir management for production history
matching, field development and future performance evaluation. Meaningful simulation
study requires extensive multi-disciplinary characterization work by geologists,
petrophysicists, geophysicists and engineers to construct a representative model or set of
equally probable models. The general workflow for reservoir studies is then to use
numerical flow simulation to validate these static models by calibration to available
dynamic data such as production and pressure measurements. Calibrated reservoir models
are utilized subsequently for production forecast under different scenarios such as primary
or enhanced recovery options, well placement potentials, optimization of stimulation
design and much more.
Geological models which are among the major inputs to flow simulators provide a
detailed description of fluid storage and flow properties of the porous system obtained
from different characterization tools, including log data, seismic, core analysis etc. These
data, which describe geological features at varying scales, are integrated into a so-called
geological model through geostatistical techniques. A characterized fluid description is
another input component to the reservoir simulator that can be very detailed, depending
on the nature of fluid flow modeling while physics of flow in the reservoir is represented in
the form of multi-phase relative permeability and capillary pressure curves. Conservation
of mass, momentum and energy are combined to model the flow and transport (Aziz &

4

Settari, 1979). The accuracy of the flow simulation is accordingly dependent on the
appropriate representation of physics of flow and detailed description of geology and fluid
properties. However, flow simulation at high resolution requires extensive computational
efforts that are commonly impractical.
Advances in hardware technology and computational power have been notable and
enables faster simulations [Lu et al. (2007), Fjerstad et al. (2005)]. However, this has not
fully resolved the computation cost due to the parallel improvements in the resolution of
reservoir descriptions. These models may contain a number of gridblocks in the order of
10
8
while computational resources typically handle as many as 10
6
[Chen and Durlofsky,
2006].

Therefore, there have been persistent efforts by researchers for some sort of
upscaling in the representation of fluid flow. Upscaling of geological models to coarser
models has been the focus of interest since the advent of reservoir simulators. In this
process, a detailed heterogeneous model is averaged to construct a coarse model with
larger and fewer simulation cells. The workflow of mapping and averaging of the fine scale
rock and flow description to a coarse scale model is termed upscaling. This procedure
reduces the number of equations to be solved and results in temporal gain in the
simulation process. Due to the averaging process and loss of subgrid information, several
kinds of errors are introduced in the simulation process. When flow is modeled
compositionally, high resolution fluid description may be required for accurate simulation
which introduces extra computation cost to the simulation work load due to equilibrium
(flash) calculations. In this case, lumping of reservoir fluid descriptions into fewer
components may be applied which reduces the computational cost [for example, Hong
(1982), Jessen and Rastegar (2009)].

5

1.1.1 Problem statement and research approach
In (highly) heterogeneous porous media, fluid flow and transport is often channelized and
does not occur equally in all regions. Some regions with connected, high fluid conductivity
contribute more to the general flow while surrounding regions with lower fluid
conductivity have less or no involvement in the viscous flow. This translates into a natural
ranking of porosity to different levels of active storage based on the degree of contribution
to flow. In other words, at the local level, fluids flow preferably in the high permeability
pathways and surrounding lower permeable rocks may contribute minimally to flow and
transport, which results in imperfect sweep and so-called channeling phenomena [Evazi
and Jessen (2013)]. In flow-based upscaling techniques, transmissibility calculation
partially captures these heterogeneity and flow preferences [Durlofsky (2003), Chen et al.
(2003)]. However, transport calculation can be quite erroneous due to the assumption of
thermodynamic (pressure and saturation) equilibration and uniform displacement at the
local level. This problem can get even worse in the case of miscible flow and displacements
with highly undesirable mobility ratio. Thus, defining active porosity (to flow) in the
process of coarse scale modeling and assigning appropriate weights for flow contribution is
expected to improve simulations of multi-phase flow problems.
Non-uniform upscaling have been used to reduce the homogenization error of
heterogeneity by using finer gridblocks at regions of higher importance and coarser
gridblocks at regions of less variation in heterogeneity [Garcia et al. (1992), Li et al. (1995),
Durlofsky et al. (1996), Castellini (2001), Mahani and Muggeridge (2005), Evazi and
Mahani (2010), Evazi and Mahani (2010)]. Coarse grid can be adaptively refined to capture
the impact of heterogeneity on the displacing fluid front and better represent the sweep

6

efficiency [see e.g., Li et al. (1995), Provost et al. (2005)]. Efficient grid generation and
calculation of flow properties for the gridblocks of varying geometry can be challenging
and reduce the attractiveness of non-uniform upscaling. Furthermore, commercial
simulators do not often handle non-uniform grid and use of uniform grid upscaling is
promoted.
Fine scale streamline simulation is commonly used to define important flow
features, such as connectivity and main flow pathways [For example, Saad et al. (1996),
Shook and Mitchel (2009)]. Streamline density and time-of-flight (TOF) provide valuable
information on potential sweepout areas in the displacement processes. Streamline density
highlights areas of high permeability connectivity that contribute more to the general flow
while TOF indicates the time-scale for flow, thus it represents the inverse of the streamline
characteristic velocity. In this work, a Streamline Index (SI) is defined as the ratio of
streamline density to the TOF that provides a combined picture of the flow setting based on
streamline density and velocity. A map of SI is used to partition the pore space into two
levels; primary and secondary porosity. Throughout the manuscript the terms, primary and
active as well as secondary and stagnant are used interchangeably.
With pore space partitioned, flow and transport is represented using a dual-porosity
dual-permeability (DPDP) flow model [Hill and Thomas, 1985]. In this context, flow and
transport is modeled in the primary and secondary porosity as two distinct continua and a
linear transfer function is used to model cross flow between these continua. Global single
phase upscaling is applied to calculate inter-porosity and inter-block transmissibility. This
treatment of heterogeneous reservoir is analogous to DPDP modeling of fractured
reservoirs, where flow and transport is occurring primarily in the fracture networks and

7

matrix-matrix and matrix-fracture flow occurs at a different time scale (slower). DPDP
modeling simulates (1) flow between primary porosity of neighboring gridblocks (similar
to fracture-fracture flow) (2) flow between secondary porosity of neighboring gridblocks
(similar to matrix-matrix flow) and (3) flow between primary and secondary porosity at
each gridblocks (similar to matrix-fracture).
In the workflow proposed as part of this research, a fine scale model is uniformly
coarsened and, based on the streamline index map, two porosities for each coarse
gridblock are defined. Incompressible single phase flow solution at fine scale is used to
calculate coarse grid inter-block and intra-block flow properties (permeability and
transmissibility). This treatment is expected to reduce the impact of the aggressive
assumption of thermodynamic equilibration in the coarse gridblocks and allow for a better
representation of non-uniform sub-grid saturation distribution and displacement
efficiency. Thermodynamic equilibration is assumed to occur in the two pore volumes
within each coarse gridblock and effect of numerical diffusion, smearing displacement front
and incorrect representation of local displacement and sweep efficiency is accordingly
reduced.
Two highly heterogeneous, channelized permeability fields are used to investigate
the performance of the current and proposed upscaling techniques, by comparison with the
fine scale simulation results. Flow and transport is dominated by high permeability
channels with large correlation length between injector and producer(s) that challenge the
performance of the upscaling techniques. Miscible and immiscible displacement
calculations are performed and the process dependency of upscaling is demonstrated.

8

1.2 Coarse scale modeling of multi-fractured horizontal wells
Low permeability and unconventional reservoirs commonly provide uneconomic
production rates unless wells are stimulated. Hydraulic fracturing is a common option that
enhances the reservoir access and increases fluid production to and beyond the economic
level. In order to maximize the productivity of these reservoirs, optimization of well
spacing and stimulation design should be investigated. For this purpose, predictive tools
that provide realistic estimates of production are desirable. Analytic solutions, numerical
simulation and decline curve techniques have been commonly used.
Decline curve techniques forecast future production based on simple decline
functions constrained to the past production [Fetkovich (1980)]. Uses of these techniques
in unconventional reservoirs are limited as it requires considerable production history for
capturing of different flow periods. This technique is mainly used for reserve calculation
where single well analysis is simply generalized to the field. Analytic solutions [Such as
Brown et al. (2009), Guo et al. (2009), Yuan and Zhou (2010)], alternatively, provide fast
estimates of production based on limited data and do not require production history. These
techniques are derived usually based on several simplifying assumptions and should be
used with care. Extension to the multi-well settings and incorporation of different physics
of flow is not straightforward. Numerical simulation, can on the other hand, provide
accurate solutions provided that input data are indeed valid and that relevant physics are
adequately represented. It provides a flexible framework for the integration of real-world
data including reservoir heterogeneity, hydraulic fracture geometry and dynamics
alongside multi-physics flow modeling and is the prime choice for rigorous reservoir study.

9

1.2.1 Problem statement and research approach
Accurate representation of the near-fracture fast-moving pressure transient is essential for
realistic production simulation, which is optimally achieved by grid refinement close to the
fractures and the well. However, it increases simulation run-time significantly and reduces
efficiency of any large-scale simulation studies. Alternatively, use of equivalent skin or well
radius [based on Cinco-Ley and Samaniego (1981)] has been practiced but is proven to fail
in many cases. Therefore, fracture upscaling and development of accurate coarse scale
simulation techniques are desirable.
Accuracy of the predictive tools is controlled by the reliability of input data and
appropriate modeling of the flow in these reservoirs. Therefore, characterization of
reservoir rock and fluid properties and completion/stimulation require significant
attention. Stress dependence of matrix permeability (rock compaction) and fracture
conductivity becomes significant [Clarkson et al. (2013)] and fluid properties may not obey
conventional (continuum scale) rules. As most of (ultra-) low permeability reservoirs
produce at an uneconomic rate unless fracturing technologies are applied, induced fracture
characterization becomes as important (if not more) as the reservoir characterization.
Pressure depletion due to production and fracturing of (older) nearby wells changes
the stress field in the reservoir and thus influences the orientation and dimensions of the
created fracture. This complicates the fracture characterization and may require the
application of direct measurement techniques, such as microseismic monitoring and tilt-
meter mapping. However, due to economic considerations, this is not often feasible and
requires the development of indirect characterization techniques. Furthermore, induced
fractures are dynamic in nature and the flow and storage properties vary with time and

10

pressure changes. Specifically, pressure decline in the fracture causes proppant
crushing/embedment that reduces fracture porosity and conductivity and thereby result in
productivity loss for a well. Realistic productivity modeling, therefore, requires for the
reasonable representation of these dynamics.
This work demonstrates, through several numerical examples, the inadequacy of
different techniques in the accurate coarse-scale simulation of flow dynamics in the
fractured blocks. We show that the use of coarse grid blocks generally result in the over-
prediction of the well productivity unless a special treatment is introduced. In the history-
matching process, it is demonstrated that this becomes problematic as it results in the
over-prediction of the bottom-hole pressure and thereby in unrealistic representation of
fluid production ratios, especially in oil reservoirs.
We adapt a multi-phase upscaling workflow for multi-fractured horizontal wells
under natural depletion. A local well model representative of the near-fracture region and
far-field area of a reservoir is used for fine scale simulation of the depletion process. We
calculate a new set of gas-oil pseudo relative permeability curves for fractured and non-
fractured gridblocks that improve the predictive accuracy of the current coarse scale
simulation techniques considerably.
We show that the dynamic nature of fractures should be correctly modeled in order
to simulate initial flush and late-time low flow rates and thereby propose a time-dependent
conductivity model that approximates in-situ conductivity loss by history matching of two
model parameters. Proposed workflow is applied to match the production history of multi-
fractured horizontal wells in a Monterey formation of California. We demonstrate the

11

inadequacy of conventional coarse scale simulation in matching of the production history
and demonstrate improvement in history matching by application of the proposed
approaches combined with the fracture conductivity model.
1.3 Organization of the manuscript
In the first part of this work, upscaling of flow in porous media, in the context of reservoir
simulation, is addressed. Chapter 2 summarizes flow upscaling techniques and presents the
calculation of effective coarse scale properties. Chapter 3 discusses dual porosity dual
permeability flow modeling as well as partitioning of heterogeneous reservoirs into
equivalent two-porosity systems. In Chapter 4, coarse scale modeling of displacement
processes is considered and the performance of single porosity and dual porosity models
are compared. Two heterogeneous reservoir descriptions are used to investigate the
performance of the proposed techniques.
In the second part of this work, coarse scale dynamic simulation of horizontal wells
with multiple fractures is discussed. Chapter 5 presents a short discussion of
characterization of low permeability and unconventional reservoirs. Characterization and
dynamic modeling of induced fractures and predictive tools for the productivity of multi-
fractured horizontal wells are discussed. Chapter 6 considers the coarse scale simulation of
hydraulically fractured horizontal wells and compares the performance of current
techniques. New approaches are introduced and improvements over existing approaches
are demonstrated.

12

Chapter 7 summarizes the work and presents the accomplishments and findings of the
dissertation. Recommendation for future study and research opportunities are also
discussed.

13

Chapter 2
Upscaling of flow in porous media

14

2 Upscaling of flow in porous media

2.1 Governing equations
Single-phase flow in porous media is modeled via conservation of mass and inclusion of an
empirical form of the momentum equation (Darcy’s law). When more than one phase is
mobile at reservoir conditions, the extended form of Darcy’s law (Muskat, 1949) is used
where relative permeability to each phase defines the mobility of individual phases. In the
case of immiscible multi-phase flow, conservation of mass for phase p is written as:

(

) .(

⃗

)

̅

(2.1)
And the equation of motion for each phase is given below:
u ⃗

(

g ) (2.2)
Here,

u ⃗

g, q ̅

represent porosity, phase mass density, phase
saturation, phase velocity vector, absolute permeability, phase relative permeability, phase
viscosity, the gravity acceleration, volumetric source/sink flow rate per unit bulk volume
and height, respectively. Subscript p refers to the phase. Fluid properties are functions of
pressure and phase relative permeability is related to phase saturation(s). Substituting
equation 2.2 into equation 2.1 yields:

(

) .(

(

g ))

q ̅

(2.3)

15

Equations 2.3 can be written for the existing phases that result in n

equations and 2n

unknowns. To complete the formulation, one needs to include (n

) capillary pressure
equations and one saturation relationship (∑

).
2.1.1 Immiscible, incompressible two-phase flow
In most oil reservoirs, oil and water make up the initial fluids distribution and flow at
immiscible conditions. To model flow in these systems, one writes equation 2.3 for oil and
water, and further assumes incompressible rock and fluids to arrive at the following
equations:

.(

[

g ]) q ̅

=0 (2.4)

.(

[

g ]) q ̅

(2.5)
Subscript o and w refer to oil and water phases and the mobility of the phase p is defined as

. For two phase flow, we need to include capillary pressure definition and
saturation relation to complete the formulation:

(

) (2.6)
and

(2.7)
Substituting equation 2.6 in 2.5 and adding equations 2.4 and 2.5, we arrive at the pressure
equation:
.[

(

)g

] (q ̅

q ̅

) (2.8)

16

The total mobility,

, is defined as the sum of the phase mobilities (

) . Three
terms in the bracket on the left hand side of equation 2.8 specify the total velocity due to
pressure gradients, gravity and capillary forces. Equation 2.8 can be reduced to the shorter
form:
.v

(q ̅

q ̅

) q ̅

, (2.9)
where the total velocity is given by:
v

v

v

[

(

)g

] (2.10)
Solving equation 2.10 for

and substituting it in equation 2.4 results in saturation
equation:
[v

g
(

)

].

.(

)

q ̅

q

(2.11)
In the equation above, fractional flow is defined as f

. Note that the pressure
equation (2.8) is an elliptic partial differential equation and therefore the solution
(pressure) is diffusive and smooth in the defined domain. On the other hand, the saturation
equation is parabolic when the capillary is non-zero (

) and becomes hyperbolic in
the case of no capillarity [Aziz and Settari, 1979]. In the latter case, discontinuities (shocks)
may appear in the solution that necessitates the application of special numerical
treatments.
Typical displacement calculations in reservoir engineering require the solution of
equations of the form 2.8 and 2.11. These equations are solved numerically by, commonly,
dividing the physical domain of interest (hydrocarbon reservoirs) into control volumes and
applying a discretized form of these equations at each control volume. The division of the

17

physical domain is initially provided through a geological model with stochastic
representation of permeability (and porosity) from limited measurements (core, log,
seismic etc.) and knowledge of the geological structure. These models entail high
resolution data and usually cannot be handled efficiently by flow simulation unless they are
upscaled. Therefore, upscaling of the geological model is a common stage of field-scale
simulation studies and is the topic of discussion in the following sections.
2.2 Upscaling of flow in porous media
It is generally assumed that governing equations at coarse scale have the same form as the
original fine scale [for example, Durlofsky (1991, 2003), Chen et al. (2003)]. Upscaling
techniques, however, calculate effective or equivalent properties for the coarse cells with
the aim of minimizing the mismatch of flow results from the coarse and fine models.
Upscaling techniques vary in the way that fluid and flow properties are calculated.
Single-phase upscaling techniques target the pressure equation and aim at
capturing sub-grid permeability heterogeneity. These techniques are commonly adequate
for simple permeability fields that can be represented by two-point geostatistics, e.g.
variogram based, if the upscaling does not destroy the geological structure. As the
permeability field becomes more complex, e.g. fluvial systems, most single-phase upscaling
techniques fail to reproduce the fine scale flow results and other upscaling techniques
(such as multi-phase upscaling) are required. Alternative techniques, such as multi-scale
methods and upscaling/downscaling, have been developed that improve the accuracy of
the coarse scale simulation.

18

In the following sections, calculation of effective (upscaled) rock, fluid and flow
properties are discussed and variations in the calculations by different techniques are
explored. Then, a short summary of alternative methods of coarse scale modeling is
presented. Finally, quantification of upscaling error and its impact on geological
uncertainty analysis is described.
2.2.1 Calculation of fluid properties
The assumption of homogeneity and thermodynamic equilibration in the coarse cells is
equivalent to uniformity in pressure, temperature and composition. This signifies
importance of calculation of fluid properties in upscaling, especially when displacement
problems are considered where sweep and displacement calculations are very sensitive to
this assumption. Density and viscosity are functions of fluid pressure (and composition)
and are affected accordingly by the assumption of thermodynamic equilibration. At the
fine-scale level, displacements are non-uniform (imperfect sweep efficiency) and so are
fluid distributions. Upon upscaling, homogenization assumptions remove non-uniformity
in fluid distribution and introduce extra errors in displacement calculation by inaccurate
representation of fluid properties. We introduce a two-level porosity upscaling that helps
in reducing the effect of homogenization on the calculation of fluid properties by splitting
porosity and the associated fluid distribution into separate but communicating sub-
volumes.
Initialization of fluid saturation and therefore oil in place calculation in upscaled
model deserves special consideration as well. When different facies are present in the
geological setting, different saturation functions (capillary pressure and relative

19

permeability curves) describe the initial saturation and multi-phase flow dynamics. As the
fine scale model is upscaled, use of averaged capillary pressure may result in significant
change in hydrocarbon in place calculation as the capillary pressure is a non-linear
function of saturation. Two-porosity heterogeneity representation at coarse scale clusters
the flowing pore volume and implicitly captures the facies variation to some extent,
thereby reducing the averaging error in initial saturation calculation. Calculation of
porosity in a coarse scale model has not received much attention and simple volume
averaging is assumed to be adequate.
2.2.2 Porosity upscaling
As upscaling techniques commonly assume the original conservation equations combined
with effective properties are valid in coarse scale, a single value for each cell is calculated.
Porosity of the upscaled model has been calculated commonly with volume-weighted
arithmetic averaging of fine scale values such that the pore volume is conserved [Durlofsky,
2003]. Since most of the upscaling techniques assume homogeneous flow and fluid
properties for each cell, this may introduce tremendous error for highly heterogeneous
reservoirs where preferential flow pathways cause the bypass - portions of the coarse cell.
Here we assume two levels of porosity in the upscaling process that helps to capture the
effect of varying contribution of porosities to the viscous flow.
2.2.3 Upscaling of absolute permeability
In contrast to porosity, absolute permeability calculation has received tremendous
attention by various research groups. In the upscaling of (multi-phase) flow in porous
systems it is generally assumed that calculation of effective or equivalent permeability is

20

sufficient for the coarse model to approximate the fine scale response for a given
displacement problem, even though this can be very erroneous in many heterogeneous
reservoirs when the correlation length of the heterogeneity is comparable to the well
spacing [Muggeridge (1991), Durlofsky et al. (1994), Barker and Thibeau (1997),]. In the
simplest form of heterogeneity, arithmetic, geometric or harmonic averaging can be used to
calculate equivalent permeability which generally provide for very low accuracy in
heterogeneous models [for example, Renard and de Marsily(1997), Farmer (2002)]. So-
called (dynamic) flow-based techniques implicitly capture permeability heterogeneity
effects in equivalent permeability or transmissibility calculations which increase the
accuracy of the upscaling. These approaches belong to the group of techniques that target
the calculating of effective coarse grid properties (permeability, transmissibility, well
index, etc.) based on steady-state single-phase flow simulation results. For this reason,
flow-based upscaling of absolute permeability is interchangeably termed as “single -phase
flow upscaling”. When geological models contain higher variability, Chen et al. (2003)
showed that it is preferable to calculate equivalent transmissibility instead of permeability.
In flow-based upscaling, since the pressure equation is used, sub-grid heterogeneity effect
on pressure is accurately captured if the correlation length does not exceed the size of
domain where fine-scale flow is solved. These techniques vary based on the domain of
upscaling and boundary conditions applied, and are discussed in the following subsections.
2.2.3.1 Local upscaling techniques
In these techniques the single flow problem is solved by assuming Dirichlet and/or
Neumann boundary conditions on a local domain including the coarse gridblock. Thus, no
global flow information is considered. The domain where the flow problem is solved is the

21

fine grid underlying the region including the coarse grid block of interest with or without
an extra border consisting of neighboring coarse grid blocks [for example, Durlofsky,
2003]. Different boundary conditions have been applied and tested. In periodic pressure
boundary conditions, when permeability or transmissibility calculation in a specific
direction is desired, unit pressure difference is set in that direction and zero pressure
difference is set in all other directions. No-flow constant pressure boundary conditions,
also referred to as standard boundary conditions, assumes a unit pressure difference in the
direction of interest and no-flow boundaries in all other directions. Depending on the
domain where flow is simulated, the following techniques have been developed:
2.2.3.1.1 Purely Local Upscaling
Traditionally referred to as the “ pressure solver method”, only fine scale grid blocks
underlying the coarse block are considered [e.g. Warren and Price (1961), Durlofsky
(1991)]. Figure ‎ 2-1(a) shows a typical fine grid and a coarse grid with an upscaling ratio of
4 in each direction. A single coarse gridblock with the underlying fine grid is shown in
Figure ‎ 2-1(b). Assuming a boundary condition pattern discussed above, flow is simulated
on the domain and equation 2.12 is used to calculate effective coarse grid permeability.
This technique provides single value of permeability in each direction. Assuming a diagonal
permeability tensor, permeability in x direction is calculated as:

, (2.12)
where

and

are the upscaled permeability, the cross sectional area and the length
of coarse grid block in the x direction, respectively. is the pressure difference across the

22

coarse gridblock in x direction. The transmissibility can be calculated subsequently by
harmonic averaging of coarse grid block permeabilities.

Figure ‎ 2-1: Schematic showing (a) fine and coarse scale grids, (b) fine scale local region
and (c) fine scale extended local region (from Chen et al., 2003)

2.2.3.1.2 Extended Local Upscaling Technique
Flow is affected by local and global heterogeneities and using fine cells limited to the coarse
gridblock cannot capture local inter-block heterogeneity. This issue can be addressed by
assuming a border around the target coarse block that covers a few surrounding grid
blocks (as shown in Figure ‎ 2-1(c)). In these techniques, known as extended local or border
region upscaling [Durlofsky (2003), Chen and Durlofsky (2006)], solution of steady-state
single phase flow with an assumed boundary condition provides the pressure and velocity

23

field. Inter-block transmissibility is then calculated by dividing the total inter-block flow
rate (sum of flow across fine cells at the interface of the target coarse grid block and its
neighbor) by the pressure difference between the coarse gridblocks.

〈 〉

〈 〉

〈 〉

, (2.13)
where 〈q〉

is the total flow at the interface of grid block i and and 〈 〉

is the volume
average pressure of fine cells inside coarse gridblock i. . This technique provides two values
of transmissibility or permeability in each direction.
2.2.3.2 Global upscaling
In highly heterogeneous permeability distributions generated by multi-point geostatistic,
the local upscaling methods lose accuracy and generally fail to reproduce fine scale
behavior adequately (Chen et al., 2003). In these cases the choice of constant boundary
conditions impacts the accuracy of upscaling results. Chen et al. (2003) showed with a
numerical example that the assumption of constant local boundary conditions is violated in
layered or channelized system and upscaling results are affected accordingly. Global
upscaling techniques [e.g. White and Horne (1987), Holden and Nielsen (2000)] use a
steady-state single phase flow solution on the whole domain at the fine scale to calculate
effective permeability and transmissibility. Global upscaling provides accurate single-phase
upscaled properties and effects of (local) boundary conditions are removed. However, the
solution of a single phase flow problem at the fine scale can be computationally demanding.
This has led to the development of Quasi-global techniques.

24

2.2.3.3 Quasi-global Upscaling (Iterative)
Highly detailed three dimensional geocellular models may contain in the order of 10
8
grid
cells and even the solution of a single-phase flow problem at this scale can be
computationally demanding. In these cases, it is possible to increase the accuracy of local
upscaling results by using local boundary conditions interpolated from a coarse scale
solution of single phase flow. This type of techniques, referred to as quasi-global or local-
global [Chen et al. (2003) and Chen and Durlofsky (2006)], attempt to capture the effects of
local and large scale heterogeneity in the transmissibility calculation. Two variations of this
approach are discussed here:
2.2.3.3.1 Coupled Local-Global Upscaling
This technique uses a global steady-state single-phase flow solution on the coarse grid to
represent the local flow problems in a more realistic fashion [Chen et al., 2003]. Pressure
boundary conditions for the local problem are interpolated from generic global flow (side
to side, corner to corner or well-driven flow) and are then used to calculate the coarse grid
block transmissibility. An algorithm for this method can be summarized as:
1. Upscale K or T using a local upscaling technique
2. Solve generic global single phase flow on the coarse model with current upscaled
properties ( from 1 or 3)
3. Upscale permeability or transmissibility by local upscaling under Dirichlet BC
interpolated from global solution
4. Repeat 2-3 until convergence is obtained (for example, until total flow difference or
pressure difference in successive iterations vanishes)

25

The use of generic global flow information makes the upscaled model process-independent
and thus suitable for simulation studies with different well configurations without the need
to update the coarse model.
2.2.3.3.2 Adaptive local-global Upscaling
It is well understood that upscaling is a process-dependent procedure and appropriate
boundary conditions must be applied to calculate accurate transmissibility. Chen and
Durlofsky (2006) showed through numerical tests that accuracy of coarse model obtained
by local-global upscaling can be quite dependent on the imposed flow field. Thus it is
expected to improve the coupled local-global upscaling results if realistic global boundary
conditions relative to the displacement scenario in question are used. These techniques are
referred to Adaptive local-global (ALG) upscaling [Chen and Durlofsky (2006)]. To make
the upscaling algorithm more time-efficient, the authors included a threshold based on flow
rate to update upscaled properties only in high-flow regions. Single phase and two phase
flow simulation results on highly heterogeneous permeability distribution reported by the
authors show significant improvement over local upscaling and coupled local-global
upscaling.
2.2.4 Multi-phase flow upscaling
Single-phase upscaling attempts to capture the gross features of flow by calculating an
equivalent permeability that simulates the same total single-phase flow through the
homogeneous coarse block as will be calculated by the coarse heterogeneous block with
underlying fine cells [Christie and Clifford, 1997]. As the single-phase upscaling only
considers a single-phase pressure equation, parameters involved in saturation equation, i.e.

26

relative permeability and capillary pressure, is not considered. This is more pronounced in
the case of large upscaling ratios where numerical diffusion/dispersion affects the
saturation equation significantly. In two-phase flow upscaling, upscaled (pseudo) relative
permeability (and capillary pressure) is generated that aims at taking into account the
effect of sub-grid heterogeneities on the multi-phase flow (e.g. Christie and Clifford, 1997).
The basic principle is to solve multi-phase flow problem on the fine grid and enforce
equality of phase flow rate calculated by the fine grid solution and coarse scale with
averaged properties:

∑ (

)

, (2.13)
where variables with superscript * indicate coarse grid values and subscript p refers to the
phase. Depending on how coarse block properties are averaged, different techniques have
been developed [for example Jac s et al. ( 973) yte and Berry ( 975) tones’ method
(1991)]]. Barker and Thibeau (1997) concluded that, except for cases where equilibrium
(capillary or gravity) can be assumed at the coarse gridblock scale, pseudos are not reliable
to scale up a fine-grid model to a coarse grid fluid-flow model. Cao and Aziz (1999) also
reported that the use of pseudos, specifically for highly heterogeneous reservoirs, is
unreliable.
Recent progress by the use of effective flux boundary conditions [Wallstrom (2002), Chen
and Durlofsky (2006)] and time-of-flight [Chen et al. (2013)] for pseudo relative
permeability calculation has shown promising results. However, there is not exhaustive
examination of these approaches available to date.

27

Multiphase flow upscaling is not considered in this work. We argue that, when
compared with single porosity upscaling, porosity partitioning and two levels equilibrium
calculations reduce errors due to sub-grid heterogeneity in the coarse scale response.
Furthermore, in displacement processes where IFT changes with composition, the scaling
of relative permeability (and residual saturations) will lose its meaning if pseudo relative
permeability is used in the course model.
2.2.5 Upscaling-downscaling
Use of coarse grid blocks in the simulation of multi-phase flow in porous media results in
the homogenization (loss) of sub-grid heterogeneity and increased numerical dispersion.
These effects smear the saturation front and introduce errors in the recovery calculation.
Upscaling-downscaling techniques [for example, Rame & Killogh (1992), Guerillot and
Verdiere (1995), Gautier et al. (1999), Audigante and Blunt (2004) ], also known as dual-
mesh method, reduce these errors by introducing an extra step of downscaling calculations
to the upscaling procedure that attempts to capture sub-grid heterogeneity effect on the
saturation (transport) equation. In the downscaling step, a fine grid velocity is
reconstructed from the coarse grid velocity and used in the fine-scale transport
calculations. Coarse grid saturation is then calculated from the fine scale solution. This
extra steps improves flow results, however, the added computational cost is considerable
and may not be viable (Babaei, 2013).
2.2.6 Multi-scale techniques
Multi-scale methods have been used to capture sub-grid heterogeneity and reduce the
numerical dispersion arising from truncation error that tends to smear the fluid front.

28

These techniques aim at modeling physical phenomena at coarse scale whilst using scale
operators to preserve sub-grid fine structures [Hesse, 2008]. Flow equations are derived at
the coarse and fine grid level. The pressure equation at fine level is defined with reduced
boundary conditions that result in the calculation of basis functions (for example, Chen and
Hou, 2003). The fine-scale pressure and velocity is then calculated by linear combination of
these basis functions. In the case of multi-phase flow, these functions must be updated to
account for variation in mobility which poses a major computational drawback [Babaei
(2013)].
2.3 Upscaling errors
Geological information is inadequate to build deterministic full-field model and uncertainty
is inherent in the generated models. In order to model the variability of the geological
model and to quantify its effects on reservoir performance, it is common to construct
realizations of the geological model and quantify the uncertainty in flow results by running
flow simulations. These realizations are generally upscaled for simulation and the
introduced upscaling error can render geological uncertainty analysis ambiguous. Bias
introduced to the flow results by upscaling is a function of the level of upscaling and should
be controlled to avoid the annulling of stochastic modeling efforts.
Two major sources of error affect the accuracy of upscaling of geological models:
Coarsening and homogenization. Coarsening error, known as numerical dispersion, is the
result of truncation of aylor’s series expansion of derivatives and triggers smearin g of the
saturation front (discontinuity) while homogenization error refers to the loss of
heterogeneity that tends to smooth the permeability field. Numerical dispersion smears the

29

front and generally causes earlier breakthrough of the displacing fluid. Homogenization, on
the other hand, tends to suppress channeling and makes the saturation front more stable
with lower velocity and delays the breakthrough time. Accordingly, coarsening and
homogenization effects on upscaling of flow (pressure) and transport equations are
different.
2.3.1 Errors in upscaling of the pressure equation
The pressure equation (equation 2.8) is an elliptic partial differential equation and the
solution (pressure distribution) is diffusive and smooth in nature and the discretization
error is generally moderate. Upscaling error due to homogenization is significant,
compared to coarsening, and increases sharply as the upscaling ratio increases to the level
that the geological structure is lost. Total upscaling errors related to the pressure equation
is mainly controlled by the loss of heterogeneity and follows a similar trend as
homogenization errors [Sablok and Aziz, 2006]. Total upscaling errors in the pressure
equation is low compared to the errors in the saturation equation due to the fact that
single-phase upscaling implicitly captures sub-grid heterogeneity in the transmissibility
calculation.
2.3.2 Errors in the upscaling of the saturation equation
The saturation equation is a parabolic partial differential equation and reduces to a
hyperbolic one in the absence of capillary forces [Aziz and Settari (1981)]. Discretization
errors affect the saturation equation considerably and similar to capillary forces smears
any sharp saturation fronts, especially in the absence of capillary forces. Loss of
heterogeneity, also, introduces considerable errors in the saturation solution by stabilizing

30

the saturation front. Sablok and Aziz (2006) observed that the total upscaling error for the
saturation does not show a sharp increase as the upscaling ratio is increased beyond the
loss of geological structure. These authors reported that, similar to the upscaling of
pressure equation, upscaling error for saturation due to loss of heterogeneity increases
sharply as the geological structure is lost in the coarse model. Due to the opposing effect of
numerical dispersion these two sources can negate each other’s effect and result in a total
upscaling error less than the sum of the absolute error of each.
In order to understand the contribution of the various error sources on the total
upscaling error, Sablok and Aziz (2006) performed an extensive study to quantify each
error and presented guidelines to control the level of upscaling and identify the potential
effects of each of these errors on the pressure and saturation equations as well as the well
flow results. Discussions below are based mainly on the work by the abovementioned
authors.
2.3.3 What is the optimum level of upscaling?
The quantification of upscaling errors and identification of the optimum level of upscaling
is a complicated problem that has not received adequate attention. Numerical tests have
shown that homogenization and coarsening errors can work in opposite directions such
that total upscaling error in global flow results (for example, well rates) does not
necessarily increase with upscaling ratios [Sablok and Aziz (2006)]. These authors
observed that for a highly heterogeneous permeability field, homogenization error is
significant even at low level of upscaling and is larger than the total upscaling error. As the
geological structure is lost at larger upscaling ratios, homogenization error increases

31

sharply. They used Q-Q plot of fine and coarse permeability values to identify the optimum
level of upscaling. Q-Q plot is a graphical method to test if two random variables belong to
the same distribution. It is constructed by plotting the values of the two random variables
at same probability values on two axes. If the Q-Q plot is a straight line with unit slope then
the two random variables are from identical distribution and deviation from this line is
indicative of change in statistical properties of the upscaled values from the original
distribution.
A numerical study of two cases by abovementioned authors, a simple isotropic
permeability field with random distribution and a complex heterogeneous permeability
field with a large correlation range, showed that the Q-Q plot is a good indicator for loss of
heterogeneity. Deviation from the unit slope line is indicative of change in statistical
properties of the upscaled values from the original distribution and thus suggests the loss
of permeability structure. They concluded that heterogeneity dominates displacement
performance and should hence be captured adequately in any upscaled model.

32

Chapter 3
Dual-porosity upscaling of heterogeneous
reservoirs

33

3 Dual-porosity upscaling of heterogeneous reservoirs

Flow and transport in highly heterogeneous porous media is non-uniform and does not
occur equally in all regions. Some regions with high fluid conductivity contribute more to
the general flow while neighboring regions with lower fluid conductivity have less or even
no flow involvement. This translates into a ranking of the pore space to different levels of
storage based on the degree of contribution to flow. In other words, at the local level, fluids
flow preferably in the high permeability pathways and surrounding lower permeable rocks
may contribute minimally to the flow and transport. In extreme cases, this setting is
analogous to fractured porous systems where connected fractures act as main flow
pathways and matrix blocks/segments slowly feed into the fracture network. Simulation of
enhanced oil recovery (EOR) techniques in these systems results in non-uniform
displacements and promotes channeling and fingering problems. At core level, Shojaei et al.
(2012) used dual-porosity approach with mass transfer between active and passive
porosities to model imperfect sweep of a heterogeneous core in a miscible displacement
and obtained a good agreement between their simulation and experimental results.
When a detailed heterogeneous model is upscaled, sub-grid heterogeneity is
homogenized into a united pore space with (equilibrium) constant pressure, composition
and saturation. Homogenization eliminates sub-grid effects on the sweepout and
distribution of compositions and, accordingly alters the related phase behavior. Therefore,
sub-scale mass transfer between different levels of porosity is not correctly modeled in
general (single-porosity) up-scaling frameworks. In flow-based upscaling techniques,
transmissibility calculation partially captures these heterogeneity and flow preferences.

34

However, transport calculation can be quite erroneous due to the equilibrium assumption
and uniform displacement at local level. This problem is particularly problematic in
compositional flows and displacements with a highly unfavorable mobility ratio. Defining
effective storage (relative to flow) in the process of coarse-scale modeling and giving
appropriate weight to flow contribution is expected to improve simulations of multi-phase
flow problems. This work introduces a dual-porosity upscaling approach for coarse scale
simulation of highly heterogeneous reservoirs.
In following sections, the workflow for the dual-porosity upscaling of heterogeneous
geomodels is presented. Then, modeling of flow in dual-porosity systems is described. By
making the analogy of highly heterogeneous reservoirs with dual porosity systems, a
transformation of a heterogeneous reservoir into a dual-porosity model based on
streamline information is described. Finally, the calculation of flow properties of the dual-
porosity coarse model is discussed.
3.1 Dual-porosity upscaling
In order to reduce the upscaling errors due to loss of heterogeneity, we propose to upscale
heterogeneous earth models into a dual-porosity coarse model. It is presumed that
averaging of the sub-grid heterogeneity will be less erroneous when performed at two
levels (instead of conventional single continuum approach) by relaxing the over-averaging
of low and high permeability regions. This requires an efficient division of porosity into
primary and secondary sets based on the contribution to the general flow arising from the
imposed flow field. This issue is addressed via streamline simulation that, given the
boundary conditions, highlights the preferential flow pathways. Porosity splitting results in

35

the generation of a dual-porosity coarse model. Flow properties (permeability or
transmissibility) of the coarse models are calculated from the pressure solution provided
by streamline simulation. The workflow proposed here is summarized below:
1. Trace streamlines on the fine grid
2. Uniformly coarsen the fine grid
3. Split pore space into primary and secondary porosity
4. Calculate inter-block and inter-porosity transmissibility for the coarse gridblocks
5. Use dual-porosity dual-permeability approach for flow modeling on coarse grid
The remainder of this chapter describes the individual steps of the proposed workflow.
3.2 Modeling of dual porosity systems
Flow in fractured porous systems occurs primarily in the fracture networks and the matrix
segments act as a source or sink to the fracture network. These systems may be modeled
with a dual porosity model (Warren and Root, 1963) as shown in Figure ‎ 3-1, where it is
assumed that the porous media can be partitioned into two continua with the matrix
system providing the main fluid storage while fluid flow takes place primarily in (high
permeability) fractures. When viscous flow between matrix blocks is significant, it must be
represented by a dual-porosity dual-permeability model. These modeling approaches are
discussed here:

36

Figure ‎ 3-1: Idealization of a fractured porous system (from Warren and Root, 1963)

3.2.1 Dual-porosity flow modeling
In this approach, the flow equations are formulated for the fracture network and linked to
the matrix with matrix-fracture transmissibility (Warren and Root, 1963). For multiphase
flow in fracture system, combining mass conservation equations for each fluid with the
extended version of Darcy’s law following equation is obtained [Kazemi et al. (1976)]:
.[

(

gD)]

(

) q

(

) (3.1)
The second term on the left side of equation 3.1 is the contribution of matrix to the flow in
fracture and can be obtained by writing the continuity equation for the matrix:

(

)

(

) (3.2)
Superscripts f and m refer to the fracture (f) and matrix (m), respectively, while subscripts
represent the phase (oil, gas and water). B and are transmissibility, mobility,

37

formation volume factor, porosity, pressure and saturation, respectively. Phase mobility is
defined as:

(3.3)
Here

is the relative permeability to phase p and

is the phase viscosity. Matrix-
fracture transmissibility is given by:

, (3.4)
where

is the matrix-fracture permeability (generally taken as matrix permeability)
and is a scaling factor (shape factor) that accounts for the matrix/fracture interface area
per unit volume. Kazemi et al. (1976) derived the following expression for matrix blocks of
size

with isotropic permeability:
(

) (3.5)
Note that this formulation assumes disconnected matrix blocks. Recovery mechanisms
from the matrix in these systems, depending on the nature of flow, can be due to
compressibility (oil expansion), imbibition, gravity drainage or diffusion.
3.2.2 Dual-porosity dual-permeability flow modeling
Viscous flow in matrix blocks are neglected in dual-porosity models that can be erroneous
in fractured reservoirs with low fracture intensity and matrix blocks with larger
permeability. Dual porosity, dual permeability (DPDP) models, on the other hand,
considers interconnected matrix segments and therefore allow for flow between matrix
blocks [Hill and Thomas, 1985]. In this approach, the flow equations are formulated for the

38

fracture network and the matrix and linked together with a matrix-fracture
transmissibility:
.[

(

gD)]

(

) q

(

) (3.6)
.[

(

gD)]

(

)

(

) (3.7)
These equations are commonly solved numerically, by spatial and temporal discretization
of the involved derivatives. In this work, with dual-porosity upscaling of heterogeneous
geomodels, flow is modeled with a dual porosity dual permeability approach. Calculation of
coarse gridblock properties are discussed in section 3.4.
3.3 Coarse-scale two-porosity representation of heterogeneous reservoirs
In upscaling of heterogeneous permeability distributions, a multitude of fine cells is
averaged into a single coarse gridblock with constant properties. As the coarse gridblock
involve larger numbers of fine cells with significant permeability variation, preferential
flow pathways emerge. Depending on the level of heterogeneity and existence of different
facies, these preferential pathways behave more or less similar to fractures and an
equilibrium assumption at this scale is quite erroneous. At this scale, splitting of sub-grid
heterogeneity into two-levels of porosity is quite meaningful. Flow modeling is then
analogous to the dual porosity, dual permeability assumption in modeling of fractured
reservoirs, where flow occurs primarily in main flow channels subject to mass transfer
between secondary porosity and primary porosity (analogous to matrix-fracture flow), as
well as mass transfer between secondary porosities (analogous to matrix-matrix flow).

39

Streamline simulation provides invaluable information regarding flow channels and
connectivity in a porous system and highlights the potential sweepout in a given
displacement problem. Under given external and internal (wells) boundary conditions,
streamlines delineate the contributing high flow pathways and heterogeneity in the context
of sweepout of the in-situ fluid. Streamline density (

) is defined here as the number of
streamlines (n

) that passes through a given gridblock and is a representative of flow
capacity:

log(n

) (3.8)
Time-of-flight ( ) , on the other hand, indicates latency of flow and is calculated from the
inverse of the characteristic velocity along a streamline. These properties of streamlines
cover a large range of values and are represented appropriately by a normalized
logarithmic scale. Figure ‎ 3-2 shows an example for the permeability distribution described
in section Figure ‎ 4-1 where streamline simulation is performed for a five spot pattern with
four producing wells on the corners and an injector in the middle of the domain.

40

Figure ‎ 3-2:Permeability map and streamline distribution in a five-spot pattern. (a)
Permeability distribution, logarithmic (b) Normalized map of streamline density (c)
Normalized map of inverse of streamline time-of-flight (d) Streamline Index map

Streamline density highlights areas of high permeability connectivity that contribute
more to the general flow. As shown in Figure ‎ 3-2(b), streamline density is generally larger
in high permeability channels. However, some channels are not highlighted as these are
short-circuits due to flow configuration arising from boundary conditions. Characteristic
velocity (Figure ‎ 3-2(c)), on the other hand, shows how fast flow occurs and assumes larger
values close to the wells and also in the high permeability areas. For a given heterogeneity
and well configuration, streamline TOF reflects saturation/composition front propagation

41

at various times and provides a direct measure of volumetric sweep at a given time [Datta-
Gupta, 2000]. A streamline index (SI) is defined here as the product of these variables:

(3.9)
As shown in Figure ‎ 3-2(d), SI combines the information of streamline density and velocity
to highlight the high flow rate areas. In this wok, SI in the logarithmic normalized form is
used to split the pore space in our presented approach.
The criteria suggested above can be used in different ways to divide the porous
media into a two-porosity system by flagging the fine-scale gridblocks either as primary or
secondary porosity. The level of contribution of porosity to flow can be evaluated based on
local and/or global measures. In the following sections, two threshold criteria at the local
and the global level are examined to identify the most efficient splitting of the pore space.
In order to evaluate the performance of the resulting coarse models, a five-spot
waterflooding scenario described in chapter 4 is simulated on the coarse and fine models.
Water-cut results on the fine grid (with 100x100 cells) and coarse models (with 10x10
cells) are then used to compare the accuracy of resulting coarse models.
3.3.1 Local porosity splitting
Use of local measures promotes the capturing of local heterogeneities. We define porosity
splitting as local when only the SI information of fine cells underlying the coarse gridblock
is used to define the primary and secondary porosity. Fine cells are flagged as primary if
they satisfy a thresholding criterion. Primary pore volume of the coarse gridblock is then
calculated as the summation of the pore volume of primary fine cells. Secondary porosity is
the difference between the total and primary porosities.

42

Large variation in SI is indicative of non-uniformity in flow and thereby justifies
porosity splitting. The span of SI, defined as the difference between the maximum and
minimum SI of the fine cells included in a coarse gridblock, is selected as a direct measure
of porosity splitting: A displacement is locally non-uniform for larger values of the span of
SI. Therefore, by setting a threshold limit, coarse gridblocks are subject to splitting if the
span of SI is larger than the limit. Fine cells with SI larger than the mean of the SI values are
then flagged as primary.
Figure ‎ 3-3, Figure ‎ 3-4 and Figure ‎ 3-5 show the porosity splitting based on the span
of SI for different cut-off values. In each figure, the first plot shows the fine scale primary
and secondary cells while the second plot shows the primary porosity contribution to the
total pore volume of coarse gridblocks. Depending on the values of the cut-off for the span
of SI, different flow systems are constructed. As smaller values of the span of SI is used to
split porous media, more coarse gridblocks are candidate for splitting and thus, the
primary pore volume increases, resulting in decreased fluid movement with larger sweep
efficiency corresponding to a delayed breakthrough of the injected fluids. This suggests
that we may be able to find optimum value of the cut-off that will result in an enhanced
accuracy of the upscaling of flow and transport.
Water-flood displacement calculations on the fine and coarse models are shown in
Figure ‎ 3-6. While conventional single-porosity upscaling (shown by SP in the plots)
overestimates the break-through time, all dual-porosity coarse models (shown by DP in the
plots) predict the break-through time accurately. Comparison of post break-through time
response, however, reveals that porosity splitting with cut-off value of 0.3 outperforms the
other models slightly.

43

Figure ‎ 3-3: Local porosity splitting based on span of streamline index (a) Primary (black)
and secondary (white) porosity distribution in fine level with span cut-off=0.3 (d) Primary
porosity ratio in coarse grid

Figure ‎ 3-4: Local porosity splitting based on span of streamline index (a) Primary (black)
and secondary (white) porosity distribution in fine level with span cut-off=0.2 (d) Primary
porosity ratio in coarse grid

44

Figure ‎ 3-5: Local porosity splitting based on span of streamline index (a) Primary (black)
and secondary (white) porosity distribution in fine level with span cut-off=0.1 (d) Primary
porosity ratio in coarse grid

Figure ‎ 3-6: Water-cut simulated by fine grid, single-porosity (SP) coarse model and dual-
porosity (DP) models generated by local splitting based on span of SI (SC)

45

3.3.2 Global porosity splitting
The use of global measures in porosity splitting is expected to result in a better
representation of larger scale heterogeneities such as long-range channels. We define
porosity splitting as global when the SI information of all fine cells is used to define the
primary and secondary porosity. A flow capacity-storage capacity diagram obtained with SI
as the flow capacity is used for this purpose. The cross-plot of flow capacity versus storage
capacity has been used to model oil recovery from layered reservoirs [For example, Stiles
(1949), Lake (1989)]. Details of the construction of static and dynamic F-C curves are
described in the paper of Shook and Mitchel (2009). We construct an F-C diagram by using
SI as the measure of flow capacity at each gridblock:

∑

∑

∑

∑

(3.10)
and storage capacity is defined as

∑

∑

(3.11)
Here

is the pore volume of cell k and n is total number of fine cells. Note that fine cells
are ordered with decreasing value of SI before the

calculation. Fine cells are flagged as
primary or secondary if the corresponding flow capacity (F) is less than a pre-specified
threshold. Figure ‎ 3-7 shows the dynamic F-C diagram for the example described in
Figure ‎ 3-2. This diagram, for instance, indicates that 50% of the pore space contributes to
75% of the flow. Different values of flow capacity are used to split the pore volume.

46

Figure ‎ 3-8, Figure ‎ 3-9 and Figure ‎ 3-10 show the resulting fine grid with highlighted flow
pathways and the coarse grid with primary porosity ratio. Water-flood displacement
calculations on the fine and coarse models are shown in Figure ‎ 3-11. Dual-porosity coarse
models generated by cut-off values of 0.3 and 0.5 predict break-through accurately while
the accuracy of the coarse model with cut-off of 0.7 is reduced.

Figure ‎ 3-7: Dynamic F-C diagram based on streamline data for the example shown in
Figure 3.2

47

Figure ‎ 3-8: Global porosity splitting based on dynamic F-C diagram with thresholding of 0.7
(a) Primary (black) and secondary (white) porosity distribution in fine level (d) Primary
porosity ratio in coarse grid

Figure ‎ 3-9: Global porosity splitting based on dynamic F-C diagram with thresholding of 0.5
(a) Primary (black) and secondary (white) porosity distribution in fine level (d) Primary
porosity ratio in coarse grid

48

Figure ‎ 3-10: Global porosity splitting based on dynamic F-C diagram with thresholding of
0.3 (a) Primary (black) and secondary (white) porosity distribution in fine level (d)
Primary porosity ratio in coarse grid

Figure ‎ 3-11: Water-cut simulated by fine grid, single-porosity (SP) coarse model and dual-
porosity (DP) models generated by global splitting based on flow capacity cut-off from
dynamic F-C curve (FC)

3.4 Calculation of coarse grid properties
Once the coarse models are constructed, gridblock properties such as porosity,
permeability and transmissibility must be calculated. Splitting of the porous media into

49

primary and secondary porosities was described in previous sections where fine cells were
flagged as primary or secondary. We define a variable α that is set equal to 1 if a cell is
primary and 0 if secondary:
α {
primary porosity
secondary porosity
(3.12)
Steady-state single-phase incompressible flow solution at the fine grid level is used to
calculate coarse gridblock flow properties. Since two porosities are calculated for each
gridblock, transmissibility between primary porosities and secondary porosities of
neighboring grid blocks as well as transmissibility between primary and secondary
porosity at each gridblock must be calculated.
Consider the two coarse grid blocks and the underlying fine cells that are shown in
Figure ‎ 3-12. Further assume that the fine cells are flagged as primary (dark cells) or
secondary (blank cells). The porosity for medium 1 (primary) and 2 (secondary) in a
coarse gridblock I can then be calculated as:

∑ ∑

(3.13)

∑ ∑ (

)

, (3.14)
where

is the porosity of fine cell i and n

and n

are the number of fine cells in x and y
directions inside the coarse gridblock.

50

Figure ‎ 3-12: Illustration of inter-block and intra-block flow for two coarse gridblocks

The inter-block primary transmissibility is calculated as:

(3.15)
Here q

is the sum of the flow between primary fine cells located at the interface of coarse
gridblocks and J.

and

are the pore volume averaged pressure of the primary fine
cells in the coarse blocks and J:

∑ ∑

∑ ∑

(3.16)
and

∑ ∑ (

)

∑ ∑

. (3.17)
Similarly, inter-block transmissibility between secondary porosities is calculated as:

(3.18)
Here q

is the sum of the flow between primary fine cells located at the interface of coarse
gridblocks and J. Flow from primary to secondary fine cells on the interface of coarse

51

gridblocks is equally divided to the inter-block secondary and primary flows. Intra-block
(secondary-primary) transmissibility for the coarse block I is calculated by summing the
flow across the primary and secondary fine cells inside coarse gridblock I:

(3.19)
In this manner, secondary-primary transmissibility is directly calculated and a shape factor
is not needed. In order to use a commercial simulator for dual-porosity dual-permeability
simulation, however, shape factor should be introduced. In this case, one may set shape
factor as 1 and use

.

52

Chapter 4
Results from dual-porosity upscaling

53

4 Results from dual-porosity upscaling

This chapter reports an investigation of the performance of conventional single-porosity
and the proposed dual-porosity upscaling. For this purpose, displacement calculations on
two challenging test models at fine and coarse scale are performed and simulation results
are compared.
Accurate calculation of single-phase flow properties (absolute permeability or
transmissibility) is an important component of upscaling and coarse-scale modeling. Three
categories of flow-based upscaling techniques, i.e. local, global and quasi-global, are used
for transmissibility calculation. Starting with an immiscible displacement problem, we
show that use of global single-phase flow upscaling outperforms the other two techniques.
Since streamline information is used for porosity partitioning, the fine-scale pressure
distribution is available and thus use of a global upscaling approach is promoted.
Coarse-scale modeling of heterogeneous reservoirs as a single-porosity and dual-
porosity model is then discussed with illustrations through simulation results. We argue
that for realistic representation of heterogeneous reservoirs with a coarse-scale dual-
porosity model, mass transfer between primary and secondary porosity should be
correctly modeled. This can be accomplished by dual-permeability assumption, where
viscous flow in the secondary porosity is also considered. As described in Chapter 3, we use
global single-phase flow solution to calculate appropriate transmissibility factors at coarse
scale.

54

Simulation of immiscible displacements is performed with a single-porosity and a
dual-porosity model constructed by the proposed techniques. Two different heterogeneous
domains are considered for upscaling and testing purposes. Simulation of near miscible CO 2
injection is subsequently performed on the fine and coarse models and results are
compared and discussed. The fine scale model results are considered as base solution and
the relevant response of derived coarse models is evaluated.
Selection of optimum upscaling level is a challenging task that has been rarely
discussed in literature. Observation by different authors indicate that single-porosity
upscaling is very sensitive to the upscaling level and the total upscaling error does not
necessarily increase as higher level of upscaling is practiced (Sablok and Aziz, 2006). We
further demonstrate the sensitivity of conventional single-porosity upscaling to the
upscaling level and report significant reduction in the sensitivity of coarse model to the
upscaling level by preserving the heterogeneity via dual-porosity upscaling.
4.1 Test Models
Two permeability models (two dimensional) are considered here. In both test models, flow
is dominated by high permeability channels between injector and producer(s). The
permeability distribution across and along channels is also quite heterogeneous making
these models challenging candidates for evaluating upscaling techniques. Immiscible and
miscible displacement problems in a 5-spot pattern are considered. No flow boundary
conditions are imposed at the left, right, top and bottom of the models. All models are
horizontal and therefore gravity effects are not included.

55

4.1.1 Model I: Channelized Model
This model is one realization of a 2D model of a synthetic channelized system (taken from
Chen et al., 2004) containing 100×100 gridblocks of size 1x1x1 m
3
in a Cartesian
coordinate system. This model was generated using multi-point geostatistics for the
channels and two-point geostatistics for permeability distributions within each facies. This
permeability model, shown in Figure ‎ 4-1, is characterized by high variations from 0.02 to
6.5×10
7
mD, a mean of 7×10
4
mD and isotropic normalized correlation length varying from
0.05-0.5.

Figure ‎ 4-1: Permeability map for test model 1 (logarithm of permeability in mD)

The porosity is assumed to be constant (25%). In this model, the correlated permeabilities
are not aligned with the Cartesian grid lines, but rotated almost 45° with respect to the grid
boundaries. The model represents a case in which the flow solution is highly affected by

56

grid orientation and flow conditions in a finite volume scheme, as the preferred flow paths
are along high permeabilities, and not aligned with the model coordinate system. High
permeability variation, high permeability channels, and complex geometry of
heterogeneities make it a challenging model to upscale.
4.1.2 Model II: 10
th
SPE Comparative Solution Project, layer 37
th

This layer is part of a highly heterogeneous and channelized fluvial system [Christie and
Blunt (2001)], where permeability ranges from 4x10
-3
to 2x10
4
mD and contains 60×220
Cartesian gridblocks (Figure ‎ 4-2). Porosity is variable in this case and varies from 0.001 to
0.4. Because of the high heterogeneity and presence of channels, this model is also very
challenging to upscale, and is a good test for the evaluation of the proposed upscaling
techniques.

Figure ‎ 4-2: Porosity and logarithmic permeability (mD) distribution of test model 2

57

4.2 Comparative study of single-phase upscaling techniques
The performance of coarse scale models in simulation of displacement problems is affected
by the accuracy of equivalent single-phase flow properties. To test this aspect of coarse
models, various upscaling techniques discussed in Chapter 2 were applied to generate
upscaled models. Transmissibility calculated by three different single-phase upscaling
techniques (Local, Global and Quasi-global) are plotted in Figure ‎ 4-3.

Figure ‎ 4-3: Transmissibility calculated by local, global and adaptive local-global upscaling

Transmissibility calculated by local upscaling and ALG are in close agreement, although
ALG calculates larger transmissibility in some points, presumably at regions with highly
varying permeability. The transmissibility calculated by global upscaling, on the other
hand, differs significantly in some points from those calculated by local and ALG upscaling.
An immiscible displacement in the five-spot pattern is initially simulated on the test
model 1. Fine and single-porosity coarse-grid models (with upscaling ratio of 4 in each

58

direction) obtained from single-phase upscaling are used for simulation. Figure ‎ 4-4 reports
the oil production rate in the individual wells as predicted by the fine and coarse models.

Figure ‎ 4-4: Oil production rate calculated by reference fine grid and coarse grid with
different transmissibility upscaling techniques.

Coarse models overestimate the production rate and breakthrough time (observable in
well 4). The error in oil production rates (%), calculated with respect to the fine-scale
production rates, is reported in Figure ‎ 4-5, which highlights the accuracy of these
techniques.
0 0.1 0.2 0.3 0.4 0.5
0
0.1
0.2
0.3
0.4
0.5
Well 1
Pore volume injected
Oil production rate, Sm3/day

Fine grid
Global upscaling
ALG upscaling
Local upscaling
0 0.1 0.2 0.3 0.4 0.5
0
0.2
0.4
0.6
0.8
1
Well 2
0 0.1 0.2 0.3 0.4 0.5
0
5
10
15
20
Well 3
0 0.1 0.2 0.3 0.4 0.5
0
20
40
60
80
100
120
Well 4

59

Figure ‎ 4-5: Oil production rate error calculated by reference fine grid and coarse grid with
different transmissibility upscaling techniques.

Recovery factor and average reservoir pressure calculated by these models are compared
in Figure ‎ 4-6. A coarse model from global upscaling predicts flow and pressure very
accurately and is isolated from the other two techniques. Clearly, global upscaling
outperforms the local and ALG techniques, even though it fails in capturing the exact
breakthrough time. Saturation maps at different injection times are shown in Figure ‎ 4-7. At
fine scale, channeling is evident, while coarse grids tend to suppress these instabilities by
averaging of sub-grid heterogeneity.
0 0.2 0.4 0.6 0.8
-100
-50
0
50
100
Pore volume injected
Oil production rate error, %
Well 1

Global upscaling
ALG upscaling
Local upscaling
0 0.2 0.4 0.6 0.8
-100
-50
0
50
100
Well 2
0 0.2 0.4 0.6 0.8
-100
-50
0
50
100
Well 3
0 0.2 0.4 0.6 0.8
-100
-50
0
50
100
Well 4

60

Figure ‎ 4-6: Recovery factor and average reservoir pressure calculated by reference fine
grid and coarse grid with different transmissibility upscaling techniques

Figure ‎ 4-7: Saturation map at different injection times predicted by fine and three single
porosity coarse models

61

4.3 Upscaling of an immiscible displacement
In coarse-scale simulation of an immiscible displacement, the assumption of homogeneous
coarse-grid porosity combined with the equilibrium assumption in the coarse gridblocks
dampens front instabilities and thus decelerates front velocity, increases sweep and
decreases local displacement efficiency. Numerical dispersion, as discussed in Chapter 2,
acts at opposite direction and promotes earlier breakthrough by smearing of the displacing
fluid front. The net result is over-prediction or under-prediction of the breakthrough time:
The effect of upscaling on prediction of displacement dynamics depends on the
displacement and sweep efficiency, and can result in either over-prediction or under-
prediction of e.g. oil recovery.
The proposed procedure for transforming a highly heterogeneous reservoir to a
coarse model based on dual-porosity consideration (see section ‎ 3.3) is applied to the test
models. When considering immiscible displacements, e.g. water-oil displacement, in a
single-layer horizontal heterogeneous reservoir, mass transfer between primary and
secondary porosity is due to viscous and capillary forces. In order for a dual-porosity
model to mimic the fine scale behavior, these mass transfer mechanisms should be
represented adequately. Effects of capillary forces on displacement are not considered in
this work. Viscous displacement of the fluid in a secondary porosity (matrix) by displacing
fluid in the primary porosity (fracture) as represented in a dual-porosity, dual-
permeability (DPDP) formulation attempt to represent this flow mechanism. In the
following sections, coarse model result with single-porosity upscaling is denoted by SP and
DP represents dual-porosity coarse model responses.

62

The fine-scale model and the coarse-scale dual-porosity model were used to
simulate the displacement process. Immiscible displacement (water-flooding) in a five-spot
pattern was simulated on the test models in a configuration where water is injected from
the well in the middle of domain with 4 producing wells located on the model corners.
Water is injected at a constant rate and production wells operate at a fixed bottom-hole
pressure. Corey type relative permeability model is used to define simultaneous flow of oil
and water:

s̅

(4.1)

s̅

(4.2)
s̅

(4.3)
Parameters of this model, fixed for all grid blocks, are summarized in Table ‎ 4-1.

Table ‎ 4-1: Corey-type relative permeability parameters for oil and water
K roe K rwe S or S wir n o n o
1 0.5 0.2 0.2 3 2

4.3.1 Test model I
The channelized model with a fine grid representation of 100 100 was uniformly upscaled
to a 10×10 coarse grid, resulting in total upscaling ratio of 100. Splitting of total pore space
into primary and secondary porosity was described in section ‎ 3.3. A five-spot water-flood
was simulated on the fine and coarse models, in which water is injected at constant rate of

63

10 Sm
3
/day and production wells produce at constant pressure of 250 barsa. Figure ‎ 4-8
and Figure ‎ 4-9 compare the recovery and producing water-cut versus pore-volume of
water injected and demonstrates the accuracy of the dual-porosity model in prediction of
breakthrough time and subsequent response. We observe that the DP approach performs
very well in reproducing oil recovery observations from the reference model. Simulation
results using conventionally upscaled single-porosity coarse model, however, shows
significant error in prediction of total recovery and the breakthrough time.
Well oil production rates simulated by fine and coarse models are shown in
Figure ‎ 4-10. Note that wells 3 and 4 (at the lower left and right corners) are main
producers that are partially connected to the high permeable channels. Pore space in the
gridblocks involving these wells are subject to splitting and two flow media is defined. This
has led to more realistic representation of near-well heterogeneity in the coarse model and
productivity is more accurately modeled. Wells 1 and 2, on the other hand, are isolated
from main channels and porosity splitting is not promoted on the corresponding
gridblocks. Accordingly, these wells show a similar response in both coarse models.

64

Figure ‎ 4-8: Cumulative oil production
resulting from water-flood simulation on
test model I (fine and coarse grid)

Figure ‎ 4-9: Producing water cut resulting
from water-flood simulation on test model I
( fine and coarse grid)

Figure ‎ 4-10: Well oil production rates simulated by fine and coarse models

65

Saturation distribution simulated on the fine and coarse grids are shown in Figure ‎ 4-11. In
order to make meaningful comparisons, the fine grid oil saturation is averaged over the
corresponding coarse grid. Saturation distributions predicted by the fine and coarse
models illustrate the impact of upscaling on channeling. The coarse model with single-
porosity dampens channeling by homogenizing non-uniform sub-grid saturation
distributions which results in a decreased front velocity, delayed breakthrough and
increased sweep efficiency.

Figure ‎ 4-11: Saturation distribution predicted by fine grid and coarse models

The waterflooding scenario described above was simulated on coarse models with 20x20
and 5x5 grid blocks, corresponding to total upscaling ratios (UR) of 25 and 400. Fractional

66

water production results, shown in Figure ‎ 4-12, indicate the sensitivity of single-porosity
models to the level of upscaling.

Figure ‎ 4-12: : Effect of upscaling level on performance of single-porosity and dual porosity
coarse models: Fractional water production

This complicates the choice of upscaling level (detection of optimal level). In the contrary,
all dual-porosity coarse models predict the breakthrough time accurately and reproduce
the post-breakthrough response more adequately. This can be explained by the splitting of
the heterogeneous porous medium into primary and secondary porosity that captures the
large-scale preferential flow paths effectively. The process of splitting is independent of
upscaling and the output is a global primary and secondary pore volume that can be
represented in coarse models of desirable upscaling ratio. In this case, however, the total
upscaling error is influenced by numerical dispersion.

67

A similar behavior is observed in the calculation of oil recovery where dual-porosity
coarse models with UR of 25 and 100 produce very similar results (see Figure 1-13). At an
UR of 400, early-time response is adequately reproduced and predictions deteriorate only
at late-times. Single-porosity coarse models, however, are extremely sensitive to UR.

Figure ‎ 4-13: Effect of upscaling level on performance of single-porosity and dual porosity
coarse models: Oil production

4.3.2 Test model II
For the second test model, the permeability and porosity distribution are described in
section ‎ 4.1.2. This is a challenging model to evaluate the upscaling techniques due to the
presence of high-permeability channels with varying width. The main flow channels
meander diagonally through the domain and trigger grid orientation effects. Porosity
splitting using a global F-C diagram is summarized in Figure ‎ 4-14 and Figure ‎ 4-15.

68

Figure ‎ 4-14: Permeability map and streamline distribution in a five-spot pattern on test
model 2. (a) Permeability distribution, logarithmic (b) Normalized map of streamline
density (c) Normalized map of inverse of streamline time-of-flight (d) SI map

69

Figure ‎ 4-15: Global porosity splitting based on F-C diagram with thresholding of 0.4 (a)
Primary (black) and secondary (white) porosity distribution in fine level (d) Primary
porosity ratio in coarse grid

The original fine grid model (60x220) is upscaled uniformly to a 6x11 coarse grid model
(total UR of 200). Results from a five-spot water-flood simulation on the fine and coarse
models are reported in Figure ‎ 4-16. In this example water is injected at constant rate of 1
Sm
3
/day and production wells are operated at constant pressure of 250 barsa. Comparison
of producing water-cut, oil rate, cumulative production and average pressure response of
coarse models relative to the fine solution verify consistently the superior performance of
the dual-porosity coarse model over conventional single-porosity upscaling.

70

Figure ‎ 4-16: Water fractional production
simulated on test model 2

Figure ‎ 4-17: Total oil production simulated
on test model 2

Figure ‎ 4-18: Average pressure simulated on
test model 2

Figure ‎ 4-19: Oil production rate simulated
on test model 2

71

Figure ‎ 4-20: : Saturation distribution predicted by fine grid and coarse grid models

Figure ‎ 4-21: Error in prediction of oil saturation by single-porosity and dual-porosity
coarse models

72

In order to evaluate the impact of upscaling ratio for the coarse model accuracy, the
original fine grid is uniformly coarsened to 12x22 and 20x44 grid systems corresponding
to upscaling ratio of 50 and 15. Simulation of waterflooding on these coarse models with
single-porosity and dual-porosity upscaling indicates consistent results to that of test
model I. Single-porosity models are sensitive to the level of upscaling while dual porosity
models are affected less by the upscaling ratio. This is illustrated with producing water cut
and oil recovery in Figure ‎ 4-22 and Figure ‎ 4-23.

Figure ‎ 4-22: : Effect of upscaling level on performance of single-porosity and dual porosity
coarse models models in test case 2: Fractional water production

73

Figure ‎ 4-23: : Effect of upscaling level on performance of single-porosity and dual porosity
coarse models in test case 2: Oil production rate

4.4 Upscaling of a near-miscible displacement
Oil recovery from near-miscible gas flooding is affected by several factors that impact the
sweep and local displacement efficiency. Recovery efficiency (RE) and its components are
defined as below:

(4.4)

(4.5)
Here

are initial oil saturation, average oil saturation and residual oil saturation,
respectively. represents the recovery efficiency,

is local displacement efficiency and

74

is the volumetric sweep efficiency. Several factors influence recovery efficiency,
including reservoir heterogeneity, mobility ratio, mixing of fluids, interplay of phase
behavior and dispersion, gravity segregation, capillary forces etc. (Stalkup, 1983).
The mobility ratio of the displacing fluid and the displaced fluid significantly affects
the volumetric sweepout in miscible displacements. When the mobility ratio is large
(>>1), displacements are unstable and channeling induced by permeability heterogeneity
grow exponentially. Channeling takes place when displacing fluid shortcuts a low
permeability zones by traveling through low-resistant high-permeability pathways and can
be aggravated with adverse mobility ratios. Channeling and viscous fingering reduce sweep
efficiency and therefore has a negative effect on recovery efficiency. Model coarsening
reduces the heterogeneity contrast and dampens the fingering and channeling which
results in overestimation of sweep efficiency. It is expected that splitting of porosity and
representation of high permeability zones with more realistic fluid distribution will help in
providing a better representation of recovery efficiency in coarse scale models.
In miscible gas flooding, fluids can be first contact miscible (FCM) or multiple
contact miscible (MCM) where, theoretically, oil is fully recoverable, or near-miscible
where high, but incomplete, recovery of oil is possible. In real reservoir conditions,
however, mixing of fluids interferes with the development of miscibility and can cause an
otherwise miscible displacement to depart from miscibility (see e.g. Johns et al., 1993). In
coarse-scale modeling of compositional flows, crude homogenization of fluid compositions
and properties introduces tremendous error via phase behavior calculations in the
predicted displacement efficiency: Fluids in all underlying fine cells are assumed to be in

75

equilibrium and of equal (average) composition, thus neglecting the potential bypassing of
oil.
The consequences of assuming local equilibrium can be more severe when miscible
displacement processes are considered, where flow and phase behavior are tightly coupled.
To explore the accuracy of proposed upscaling techniques, near-miscible gas injection is
simulated on fine scale and upscaled models. Injection and initial fluid compositions are
shown in Table ‎ 4-2 and Table ‎ 4-3. The Peng-Robinson equation of state (1976) is used to
represent the phase behavior of the hydrocarbon mixture. Physical properties of the
mixture components and EOS parameters are listed in Table ‎ 4-4. The initial reservoir fluid
has a saturation pressure of 99.3 bars at a reservoir temperature of 344 K. The minimum
miscibility pressure (MMP) for the gas/oil pair at this temperature is 129.7 bars. An initial
reservoir pressure of 130 bars is assumed and the bottom-hole pressure of the production
wells is set to 110 bars while CO 2 is injected at fixed rate of 250 Sm
3
/day and 1000
Sm
3
/day for the first and second test models, respectively.

76

Table ‎ 4-2: Initial fluid composition
Component Molar ratio
CO2 0.01
Methane 0.35
n-butane 0.25
n-dodecane 0.39

Table ‎ 4-3: Injection fluid composition
Component Molar ratio
CO2 1.0
Methane 0.0
n-butane 0.0
n-dodecane 0.0

Table ‎ 4-4: Physical properties of hydrocarbon mixture components for PR-EOS
Component Tc(K) Pc(atm) Omega Mw Zc Shift
CO2 304.21 72.865 0.2236 44.01 0.2742 -0.0817
Methane 190.56 45.389 0.0115 16.043 0.286 -0.1595
n-Butane 425.12 37.464 0.2002 58.123 0.274 -0.0675
nC12 658.00 17.962 0.5764 170.338 0.238 0.2251
4.4.1 Test model I
Fine and coarse models of the permeability distribution described for test model 1 is used
to simulate the near-miscible displacement scenario described above. Porosity splitting of
this model is described in Chapter 3. The fine scale grid was coarsened with an upscaling
ratio of 10 in each direction (total UR of 100). Simulation results shown in Figure ‎ 4-24 and
Figure ‎ 4-25 and demonstrate the improved performance of the dual-porosity model over
the single-porosity upscaling. The dual-porosity model predicts early time production and
the breakthrough time very well. Both coarse models over-predict oil recovery due to over-
estimation of sweep efficiency, however, this is less pronounced in the case of dual-

77

porosity upscaling where porosity splitting helps to reduce the over-averaging of
heterogeneity and to represent the sweep efficiency more accurately.

Figure ‎ 4-24: Cumulative oil production
simulated by fine and coarse grids of test
model I

Figure ‎ 4-25: Production gas-oil ratio
simulated by fine and coarse grids of test
model I

Figure ‎ 4-26: Oil production rate simulated
by fine and coarse grids of test model I

Figure ‎ 4-27: Accuracy of coarse models in
the prediction of oil saturation: Root-mean
square error

78

Saturation maps at different simulation time, shown in Figure ‎ 4-28, reveals interesting
transport related issues. Partitioned porosity helps to reduce front dampening of the
coarse model and thus, the dual-porosity coarse model performs very well in reproducing
the early-time flow instabilities observed from fine scale model: Channeling and finger
growth is better reproduced by the dual porosity model. Local displacement efficiency is
also captured fairly well with the dual porosity model, especially at early times and
outperforms the single porosity response. This is shown in the plot of root-mean-square
deviation of oil saturation in Figure ‎ 4-27.

Figure ‎ 4-28: Near-miscible displacement; saturation map at different pore volume injected
simulated with fine and coarse models.

79

As a fine scale model is upscaled, sub-grid heterogeneity effects vanish and in general, the
accuracy of the coarse scale model in reproducing the fine scale response reduces. In order
to study the impact of the level of coarsening/upscaling on the performance of the resulting
coarse model, we consider the near-miscible displacement discussed in section ‎ 4.4.1. The
original fine model (100x100) is upscaled to single-porosity and dual-porosity coarse
models with 20x20 and 5x5 grid systems, corresponding to total UR of 25 and 400.
Simulation results from the coarse models are compared with the fine model in Figure ‎ 4-29
and Figure ‎ 4-30, and we observe that as the upscaling ratio increases, accuracy of the
single-porosity upscaling decreases dramatically. This is indicated by over-prediction of
breakthrough time (optimistic) and associated over-prediction of the oil recovery. This
example further demonstrates that loss of heterogeneity in single-porosity upscaling plays
a dominant role in decelerating the saturation fronts.
In contrast, the simulation results indicate that dual-porosity coarse models are
significantly less influenced by the level of upscaling as breakthrough time is predicted
very accurately by all U ’s. t can be argued that as the upscaling ratio increases, coarse
grid blocks contain larger number of fine cells where flow connectivity is better
represented by the dual-porosity approach: A larger variation in features related to the
flow makes porosity splitting and transmissibility calculations more meaningful and the
displacement is represented in a more realistic manner.

80

Figure ‎ 4-29: Effect of upscaling level on the accuracy of single-porosity and dual-porosity
uspcaled models in the simulation of oil production rate

Figure ‎ 4-30: Effect of upscaling level on single-porosity and dual-porosity uspcaled models
in the simulation of producing gas-oil ratio

81

4.4.2 Test model II
The near-miscible displacement scenario described in section ‎ 4.4 was next simulated on
Test Model II, in which CO2 is injected at constant rate of 250 Sm
3
/day. Original fine grid
and 12x22 (UR = 50) uniform coarse models with single-porosity and dual-porosity
upscaling were used to simulate the displacement process. Figure ‎ 4-31 and Figure ‎ 4-32
report the calculated recovery and producing GOR and indicate a very good agreement
between the dual-porosity coarse model and the reference fine-grid model. Saturation
maps at different times further verify the performance of the dual-porosity model. The
average saturation error (Figure ‎ 4-34) is less that 10% in the dual-porosity model
compared to 15% in the single-porosity model.

Figure ‎ 4-31: Production gas-oil ratio
simulated by fine and coarse grids of test
model II

Figure ‎ 4-32: Cumulative oil production
simulated by fine and coarse grids of test
model II

82

Figure ‎ 4-33: Oil saturation map at different times simulated by fine and coarse grids:

Figure ‎ 4-34: Error in prediction of saturation by single-porosity and dual-porosity coarse
models

83

As for Test Model I, the effect of upscaling ratio on the accuracy of coarse models was
studied. The reference fine-grid model was upscaled to 20x44 and 6x11 coarse grids and
the near-miscible displacement process was simulated. The oil production rates and
producing GOR, reported in Figure ‎ 4-35 and Figure ‎ 4-36, demonstrate the sensitivity of
single-porosity coarse model performance to the upscaling ratio where oil recovery is
consistently over-predicted.

Figure ‎ 4-35: Effect of upscaling level on the accuracy of single-porosity and dual-porosity
uspcaled models in the simulation of oil production rate

84

Producing gas-oil ratio, as shown in Figure ‎ 4-36, is underestimated by the single-porosity
models as displacement efficiency is overestimated. Dual-porosity models predict gas
breakthrough very accurately while the post-breakthrough accuracy at very high upscaling
ratio is reduced.

Figure ‎ 4-36: Effect of upscaling level on the accuracy of single-porosity and dual-porosity
uspcaled models in the simulation of production gas-oil ratio

85

4.5 Concluding remarks
This work introduces dual-porosity upscaling for coarse-scale modeling of highly
heterogeneous reservoirs. Transformation of fine-scale heterogeneous geomodels into a
dual-continuum model is accomplished by utilizing streamline information to identify key
flow features and define a streamline index. Local and global measures of porosity splitting
are used to construct a dual-continuum representation of the fine-scale model. A dual-
porosity dual-permeability model is adapted to simulate the relevant mass transfer
phenomena. An incompressible single-phase flow solution is used to calculate appropriate
transmissibility factors. Immiscible and miscible displacements were performed on two
heterogeneous models at fine and coarse scales. Single-porosity and dual-porosity coarse
model responses were compared with fine scale solutions and the sensitivity to the
upscaling ratio was studied. From the results and analysis presented above, we arrive at
the following conclusions:
1. Global single-phase flow upscaling outperforms the local and quasi-global
techniques. More accurate transmissibility calculation decrease errors of coarse
models in simulation of displacement processes.
2. Streamline information provides valuable information for estimating active pore
space of highly heterogeneous realizations. Defining the Streamline Index (SI) as the
ratio of streamline density to the time-of-flight helps to highlight important flow
regions and facilitates efficient porosity partitioning.
3. Splitting of the heterogeneous porous medium into primary and secondary porosity
captures the large scale preferential flow pathways effectively. Porosity partitioning
and dual-porosity dual-permeability flow modeling reduces errors introduced due

86

to thermodynamic equilibrium assumption (numerical diffusion) in the coarse grid
model.
4. Dual-porosity coarse models predict the breakthrough time accurately and
reproduce the post-breakthrough response adequately while single-porosity models
overestimate arrival time and oil recovery.
5. By preserving large scale heterogeneities, dual-porosity coarse models are
demonstrated to be significantly less sensitive to the level of upscaling, compared to
conventional single-porosity upscaling, and the total upscaling error is more
influenced by numerical dispersion/diffusion than by loss of heterogeneity.

87

Chapter 5
Modeling of multi-fractured horizontal wells

88

5 Modeling of multi-fractured horizontal wells

Accurate multi-phase flow modeling in low-permeability and unconventional reservoirs
and simulation of multiply fractured horizontal wells is challenging and has recently
received significant attention. The pore structure of these systems is different from the
conventional reservoirs and storage, flow and transport properties can be quite different as
well. In addition, rock compaction can be considerable and dynamic modeling of induced
fractures is vital to a realistic simulation of productivity of these resources. Flow modeling
and simulation of these systems is complicated and has challenged the capabilities of
current models in history matching of the past production as well as forecasting of the
future production.
Induced fractures increase well access to different zones (e.g. layers) of a reservoir
and provide main flow pathway for production. In order to simulate the fluid production in
these reservoirs, therefore, it is vital to estimate an average fracture surface that would
mimic dynamics of near-well pressure drops and production rates. Fracture properties
vary as pressure and other in-situ parameters change during the well life [Palisch et al.
(2007), Weaver et al. (2010)] and must be properly modeled to enable realistic production
simulation. Flow in fracture also requires special treatments as linear multi-phase flow in a
fracture of small pore volume can be numerically challenging.
Field development and production optimization in these systems are usually
challenged by such common questions as optimum well spacing and selection of the
number of fracturing stages. For this purpose, predictive tools that provide accurate
production predictions for multi-fractured horizontal wells are very desirable. Major

89

challenges in the development of these tools are tied with the proper representation of
very thin fractures in a very large reservoir.
This chapter summarizes work towards accurate production modeling from
unconventional reservoirs. In the following sections, characterization of low permeability
and unconventional reservoirs are shortly discussed. Then, characterization and dynamic
modeling of induced fractures is presented. Predictive tools for the productivity of multi-
fractured horizontal wells, including analytic solution and numerical simulation is
discussed subsequently.
5.1 Unconventional reservoir characterization
Conventional techniques for fluid storage and flow behavior characterization of
unconventional reservoirs become very inefficient, if not impractical, due to significant
difference in pore structure. Wettability, transport and fluid PVT properties at nano-pores
is not yet well understood, especially when pore throat sizes are only an order of
magnitude larger than fluid molecules [Elgmati et al. (2011)]. High pressure (up to 60,000
psi) mercury injection capillary measurement (MICP) is commonly used to provide
effective (connected) pore volume and pore-throat distribution as well as air-mercury
drainage capillary pressure curves. Despite the fact that this technique is destructive and
material integrity is under question, it is currently the most common method for
petrophysical analysis of unconventional reservoirs [Elgmati et al., 2011]. Different
techniques have been developed to calculate absolute permeability from the MICP data
[Swanson (1981), Comisky et al. (2007)]. Direct liquid permeability measurements become

90

tiresome and practically impossible for shales with nano-Darcy permeability. Furthermore,
relative permeability measurement is extremely tedious if not impossible.
Inadequacy of laboratory measurement of storage, flow and transport properties
has led to image-based characterization techniques. Among them are high-resolution X-ray
computed tomography (HRXCT) and focused ion beam-scanning electron microscopy (FIB-
SEM) that delineate 2-D and 3-D pore structure and micro fracture at different resolutions.
Computational fluid dynamics techniques such as Lattice-Boltzman Method (LBM) have
been used to calculate absolute permeability and two-phase relative permeability from the
3-D rock images [Grader et al. (2010), Tölke et al. (2010)].
5.2 Hydraulic fracture characterization
Characterization of the induced fractures is a key component in productivity analysis,
stimulation design and well placement studies. Induced fractures extend the well contact
with the reservoir and thus control, in part, the production dynamics. For the convenience
of numerical simulation, there is a general tendency to assume that fractures grow as bi-
wing planes contained in the pay-zone. However, industry experience of fracture diagnosis
indicates that fractures easily grow asymmetrically out of the pay zone into water or gas
zones with irregular patterns, instead of forming a single plane [Barre et al., 2002]. In
addition, the hydraulic fracture is an unconsolidated proppant pack and undergoes
significant changes during production life. In order to simulate initial high production rates
and late-time low flow rates that are commonly observed in these wells, therefore, the
dynamic nature of the fractures should be correctly modeled.

91

5.2.1 Hydraulic fracture geometry
Fracture diagnosis techniques can be categorized as direct measurements or indirect
methods. Direct measurements include near-well or far-field techniques that measure log-
scale (few feet) and large-scale (up to hundreds of feet) features of the fractures,
respectively. Near-well techniques include radioactive tracers, temperature and production
logs that provide detailed information of the fracture height at the wellbore, production
profile at the perforations and proppant concentration in the perforation intervals [Barre
et al., 2002]. Far-field techniques, such as microseismic and tiltmeter mapping, can provide
a gross view of the fractures such as orientation and dimensions.
Due to economic considerations, direct measurements are done only on few select
wells. This has led to the development of indirect fracture diagnosis techniques. Fracture
modeling, well testing and production data analysis are most common indirect techniques.
Indirect diagnostic analysis usually suffers from non-uniqueness and must be verified by
direct measurements. For this reason, it is common to combine direct diagnostic
measurements in few wells to come up with an average picture of the fracture growth and
utilize it to calibrate a fracture model that simulates fracture growth under in-situ
condition. This calibrated model may be used to predict fracture growth in the same
reservoir if operation condition changes.
5.2.2 Hydraulic fracture conductivity
Induced fractures provide the main pathway for fluid flow from reservoir into the well and
are thus desired to be as conductive as possible. Several lab studies have shown that
porosity and conductivity of induced fractures declines considerably during production

92

[Palisch et al. (2007), Weaver et al. (2010)]. It is also well known that conductivities
observed in the field are considerably different from the lab measurements. Different
mechanisms contribute to the conductivity loss, including proppant embedment and
crushing, fracturing fluid damage, multi-phase flow effects, fines invasion and proppant
diagenesis [Weaver et al. (2007), Weaver et al. (2010)]. Accurate flow modeling from
fractured wells, therefore, necessitates dynamic modeling of fracture conductivity and
failure in realistic representation of conductivity can result in inaccurate interpretation of
well-test data and biased numerical simulation. For instance, in the transient analysis of
shale gas wells, Thompson et al. (2010) suggested that productivity degradation due to
(fracture and matrix) permeability loss may be interpreted as the start of boundary
dominated flow, thus resulting in the under-estimation of reserves.
Therefore, optimum fracture and well spacing studies can be influenced by
improper modeling of hydraulic fracture conductivity: An inaccurate simulation model can
result in selection of larger well and fracture spacing if fracture conductivity deterioration
is not properly taken it into consideration
5.2.2.1 Stress dependence of fracture conductivity
Production from fractured wells accompanies pressure decline in the fracture that in turn
increases the effective stress on the fracture. Increased effective stress results in proppant
embedment/crushing and thus in reducing the fracture conductivity and storage
properties that directly affect the well productivity. Stress dependence of propped
hydraulic fractured conductivity is controlled primarily by the type, size and concentration
of the proppants used [Palisch et al. (2007)].

93

In the API adopted conductivity cell, proppants are placed between core slabs and
conductivity is measured at different confining stress by fluid flow across the fracture.
These measurements provide information regarding the permeability decrease due a pore
pressure reduction [Pedroso and Correa, 1997]. Single phase flow is used in the
experiments and proppants used are clean and thus measured values provide the upper
limits of conductivity. The mechanical properties of the rocks are functions of effective
stress [Pedrosa, 1986], and is expressed by the effective stress law:

(5.1)
Here

is the closure stress, is the Biot’s poroelastic constant and

is the pore pressure.
Biot’s constant is unity for completely unconsolidated roc s such as pac ed proppants
inside fractures. Therefore, equation 5.1 simplifies to:

(5.2)
At reservoir conditions, the closure stress is constant except if tectonic activities occur.
Therefore, the pore pressure controls the effective stress and consequently the rock
properties. A permeability modulus ( ) has been used to express stress dependence of
fracture conductivity (Yilmaz and Nur, 1991):

(5.3)
Here K is the fracture permeability. Integrating equation 5.3 and substituting for

we
find

exp [ (

)

(

)

] (5.4)

94

If the closure stress is constant, equation 5.4 reduces to

exp [ (

)] (5.5)
In commercial simulators, therefore, pressure dependence of fracture permeability may be
input in a tabular form [De Souza et al. (2012), Okouma (2011)]. This treatment, however,
assumes conductivity is solely function of effective stress and neglects the effect of other
sources of conductivity loss discussed below.
5.2.2.2 Laboratory vs. in-situ fracture conductivity
Fracture conductivity depends on proppant type, size and concentration and is affected by
effective confining stress, fracturing fluid residue, multiphase flow effects and etc. (Palisch
et al., 2007). The study by Palisch et al. (2007) indicates that in-situ fracture conductivity
can drop to 5% of the laboratory observations. Separate measurements by Davies and
Kuiper (1988) show that fracturing fluid residue may reduce the proppant conductivity by
90% and, in addition, two-phase flow effects may reduce the effective fracture conductivity
significantly (by a factor of about 5). Furthermore, lab measurement by Weaver et al.
(2007, 2010) indicates that geochemical interaction between proppants and the formation
water at high stress and temperature conditions accelerates the proppant dissolution
process and results in porosity and permeability loss, as high as 90%, in fractions of a year.
5.2.2.3 In-situ conductivity modeling
Use of appropriate values of in-situ fracture conductivity is essential for accurate
production simulation. The dependence of fracture conductivity on different parameters
and its dynamic nature complicates the in-situ conductivity modeling. These issues bring
uncertainty to the fracture conductivity representation in flow modeling and therefore

95

should be modeled with the potential uncertainty and considered as a history matching
parameter.
Rate transient analysis has been used to characterize dynamic variability of
hydraulic fracture and matrix permeability [Clarkson (2013), De Souza (2012)]. Clarkson
et al. (2013) used, in the context of rate-transient analysis, modified pseudo-pressure and
pseudo-time to incorporate stress-dependent matrix permeability and approximated the
fracture conductivity changes via time-dependent skin effect. These authors show that due
to loss of conductivity during production, the traditional linear and square-root of time
plots become non-linear unless modified pseudo-variables are used. Otherwise, one may
interpret the start of non-linearity as the indication of pseudo-steady state flow and
therefore underestimates the reserves. To our best knowledge, there is no time-lapse
conductivity model or correlation currently available for productivity modeling.
5.3 Predictive tools for production
Reliable predictive tools are necessary for the optimization of field development, well
placement, stimulation design and investigation of key parameters influencing the
performance of fractured horizontal wells [Lolon, 2007]. Three commonly used tools with
limitations and advantages are discussed below.
5.3.1 Decline curve analysis
Decline curve techniques forecast future production based on simple decline functions
constrained to the past production [Fetkovich, 1980]. Uses of these techniques in
unconventional reservoirs are limited as it requires considerable production history for

96

capturing of different flow periods. This technique is mainly used for reserve calculation
where single well analysis is simply generalized to the field.
5.3.2 Analytic solution for fractured horizontal wells
Analytic and semi-analytic solutions of productivity and pressure transient in multi-
fractured horizontal wells have received considerable attention [Wei and Economides
(2005), Guo et al. (2009), Yuan and Zhou (2010), Medeiros et al. (2008), Brown et al.
(2009), Zhou et al. (2012)]. Among the most recent and general solution, Brown et al.
(2009) derived an analytical model that represents transient tri-linear single-phase flow.
This solution incorporates linear flow from the outer reservoir into the inner reservoir
(where hydraulic fractures were created), linear flow from the inner reservoir to the
fracture and linear flow from the fracture to the wellbore accounting for the near-well
converging radial flow skin. The solution is provided for single-porosity and dual-porosity
reservoirs. Shojaei and Tajer (2013) extended Brown et al.’s solution to the multi -phase
flow. By incorporating solution-gas drive model [Muskat (1949)], they developed a
transient three-phase flow model for a multi-fractured horizontal well.
Analytical solutions provide a computational advantage over numerical simulations
and represent early-time pressure transient accurately; nonetheless these solutions should
be used with care as there are several simplifying assumptions imbedded in the derivation.
For instance, these solutions may lose accuracy in heterogeneous (layered) reservoir, and
may not be applicable when fracture geometries are not well-defined or are asymmetric
(most-likely). Furthermore, these solutions lose applicability as the physics of flow are not
honored or field-scale problems are studied.

97

5.3.3 Numerical simulation
In-depth reservoir studies of fracturing design and field development are only possible via
reservoir simulators and cannot be replaced with analytic solutions [Meyer et al. (2010)].
Governing multi-phase flow equations in porous media (equations 2-8 and 2-11) is solved
numerically to simulate the productivity of fractured horizontal wells. Numerical
simulation provides an accurate solution based on the validity of input data and a response
model that captures the governing physics. It presents a flexible framework for the
integration of real-world data including reservoir heterogeneity, hydraulic fracture
geometry and dynamics and multi-physics flow modeling and is the prime choice for
rigorous reservoir studies.
Simulation of multi-fractured horizontal wells is influenced by the proper modeling of
hydraulic fractures as the main flow pathways. Special numerical treatments are required
to accurately model the fractures. A fracture, in a generic sense, acts as a constant pressure
plane that connects different layers and makes production from different parts of the
reservoir possible. As the physical domain is divided into (rectangular) control volumes,
typically with dimensions in the order of 100 ft, explicit inclusion of fractures with the
typical width of about 0.05 ft is not practical. In addition, dynamics of multi-phase flow in
the fracture is very different (linear relative permeability curves) than the surrounding
rock (nonlinear curves) and further complicates the fracture simulation. Several
techniques have been used for the simulation of fractured horizontal wells that are
summarized below:

98

5.3.3.1 Explicit simulation
In this approach, very fine gridblocks with the size of the well (diameter and length) and
the fractures are used to explicitly represent the high-conductive well and fractures. Large
transmissibility values are then applied to the well and fracture gridblocks. In this manner,
well and block pressure is equal and thereby eliminates the need for use of well index [Wan
et al. (1998)]. This technique is computationally very expensive and is generally used only
to assess the accuracy of other techniques.
5.3.3.2 Local grid refinement
In order to simulate the near-fracture fast moving pressure transient and the initial
production rates accurately, it is common to refine the gridblock size at the vicinity of
fracture in a logarithmic fashion. Use of small gridblock size with large conductivity
adjacent to a larger gridblock size with significantly lower permeability can cause
numerical convergence problems [Sadrpanah, 2006]. In addition, simulation run-time can
increase beyond the limits of computational capacity at larger-scale simulation studies. In
this work, as in many other works, explicit simulation of fractured wells with locally
refined grid is used as the reference solution.
5.3.3.3 Use of negative skin (effective wellbore radius)
The simplest approach to simulate a fractured horizontal well is to use a vertical well with
a negative skin (Bonney et al., 2013). Cinco-Ley and Samaniego-V (1981) derived an
expression for the effective wellbore radius based on the pseudo steady-state solution.
However, this approach does not represent the initial transient period that can be quite
long for low permeability reservoirs. Cvetkovic and Norge (2009) used a standard vertical

99

well model with an effective wellbore radius and single-fracture half-length to model the
productivity of a multi-fractured horizontal shale-gas well. Similarly, this approach does
not capture the initial transient period. This approach produces radial pressure
disturbance around the well while, in reality, a pseudo-elliptical pressure disturbance is
introduced.
Effective wellbore radius approach enables the simulation of multi-fractured horizontal
wells with coarse gridblocks and is computationally very favorable for field scale studies.
However, there are several drawbacks to this technique. It is developed for single phase
flow in homogeneous reservoirs. Use of high negative skin introduces substantial error in
the calculation of production rates due to numerical problems and it does not capture the
physics of flow defined by the presence of hydraulic fractures. In addition, this method is
applicable only if the effective wellbore radius is smaller than the pressure equivalent
radius of the gridblock [Lolon et al., 2007].
5.3.3.4 Permeability/Transmissibility Modification
In uniform coarse grid simulation, permeability or transmissibility of coarse gridblocks
intersected by fracture can enhanced to mimic the hydraulic fracture effect. . Simple
relations are obtained by averaging of fracture plane and gridblock as beds in parallel or in
series [El-Ahmadi and Wattenbarger, 2004]. Assuming that the vertical fracture is oriented
along the x-axis, following relations are derived for gridblock permeability:

̅

(

)

(5.6)

̅

(

)

(5.7)

100

̅

(

)
(

)

(5.8)
Here

is the size of the gridblock in y-direction, b

is the width of the fracture and is the
initial gridblock pearmeability. The porosity of grid blocks representing fractures is
modified to include fracture pore volume:

̅

(

)

(5.9)
n addition eaceman’s well connection factor is modified as following:
W
(
̅

̅

)
.

(

)
(5.10)
Here,

is the wellbore radius and

is the equivalent well block radius that is given by
following relationship:
r

.
[√

̅

̅

(

)

√

̅

̅

(

)

]
.
(

̅

̅

)
.
(

̅

̅

)
.
(5.11)
Multi-phase flow in the fractured grid block is often represented by the original form the
gridblock relative permeability curves.

101

Chapter 6
Coarse-scale simulation of multi-fractured
wells

102

6 Coarse-scale simulation of multi-fractured wells

Tight and unconventional resources commonly cover large area and economic production
is attainable by increased access to the reservoir through induced fractures that counter
balances the low permeability effect on the productivity of the reservoir. Efficient
numerical flow simulation in these reservoirs requires use of large grid blocks (in the order
of hundreds of feet) where accurate representation of flow in induced hydraulic fractures
becomes challenging. Simulation of fast-moving pressure transient induced from fracture
surface, due to production, is achieved historically by locally refining the grid. This
becomes of significant importance in oil reservoirs where modeling of pressures in the
vicinity of the well and fracture affects the multi-phase flow considerably. Prediction of
dissolved gas evolution and nucleation into a continuous phase, characterized by critical
gas saturation, is also a key factor in productivity analysis as it significantly reduces the oil
flow (production).
Production from oil and gas wells is strongly affected by initial and subsequent
near-well reservoir rock and fluid properties. The pressure transient and decline in the
vicinity of the well influences reservoir rock properties such as porosity and permeability
and changes fluid properties, specifically when it falls below saturation pressure. For a
reservoir under natural depletion, a pressure reduction below saturation pressure
accompanies near-well liquid drop-out and release of dissolved gas in condensate gas and
oil reservoirs, respectively. The appearance of a new phase and further development alters
the oil or gas flow significantly and reduce or prevent the production of otherwise

103

producible far-field hydrocarbons. This becomes more complicated in fractured horizontal
wells where dynamics of fracture are involved.
In order to predict the well performance realistically, simulation models must
represent the aforementioned near-well dynamics adequately. This can be accomplished
primarily by fine-scale simulation of wells and fracture(s). Explicit simulation of hydraulic
fractures is not currently practical and coarser size gridblocks must be used to represent
fractures. Coarse scale simulation of the near-well dynamics requires special numerical
treatments that must be implemented with care. Coasre scale representation of near-well
dynamics is difficult as the gridblocks are too large to capture localized behavior
(Nakashima, 2010). This has led to the development of near-well upscaling techniques
[Nakashima and Durlofsky (2010), Rojas et al. (2013)].
In near-well upscaling, a domain around the well that is represented by one coarse
gridblock or a ring of coarse gridblocks, known as a Local Well Model (LWM), is used to
predict the observed single-phase and multi-phase flow dynamics of the underlying fine
grid. Single-phase and multi-phase flow parameters of the coarse gridblock are calculated
by minimizing the mismatch between the fine and coarse scale flow in an iterative process.
Nakashima and Durlofsky (2010) used adjoint gradient-based procedure for minimization
with boundary conditions interpolated from a coarse grid solution, similar to the ALG
method discussed in Chapter 2. These techniques are presented primarily for
heterogeneous reservoirs without the fractures. This work applies similar procedure to
fractured horizontal wells under natural depletion and describes a framework to calculate
a new set of gas-oil relative permeability curves for fractured and non-fractured gridblocks.

104

In the following sections, the governing equations of depletion-type three-phase flow
are described and multi-phase flow upscaling is presented. Explicit and coarse scale
simulation of fractured wells with local grid refinement is described and used as the
reference solution. The accuracy of coarse scale simulation with single-phase and multi-
phase upscaling is then investigated.
6.1 Governing equations
Simultaneous flow of oil, gas and water is often modeled by a three-phase three-component
flow model. Oil and water components flow merely in the designated oil and water phases
while the gas component can exist in both the oil and gas phases. Three phase flow in a
porous media is formulated by combining the mass conservation equations with the
extended form of Darcy’s equation. Following equations can be derived for oil, gas and
water flow (Aziz and Settari, 1979):

(

) .[

(

g )] q ̅

(6.1)

(

) .[

(

g )] q ̅

(6.2)

[ (

)] .(

(

g )

(

g )) q ̅

(6.3)
Here, B

g, q ̅ represent pressure, saturation, porosity, formation volume
factor, absolute permeability tensor, relative permeability, viscosity, the gravity
acceleration, volumetric source/sink flow rate per unit bulk volume and elevation,
respectively. Subscripts o w g refer to the oil, water and gas phase. Fluid properties are
functions of pressure and phase relative permeability is related to phase saturation (s).

105

Equations 6.1-3 are often solved numerically by spatial and temporal discretization of the
involved derivatives on a Cartesian grid. Well flow rates are calculated as:
q ̅

W (

) (6.4)

is the well pressure, is the gridblock volume and W is the well index often defined by
Peaceman formula (1978). For a well penetrating in the z-direction
W
√

(

)
(6.5)
Here r

is the wellbore radius and r

is the equivalent well block radius given by:
r

.
[√

( )

√

( )

]
.
(

)
(6.6)
At coarse scale, similar equations with upscaled properties are assumed. These parameters
are designated by superscript *. In this work,

and

is upscaled and other
parameters are assumed as in the original form. We solve the single phase incompressible
flow to calculate

and

. Three phase flow in the context of depletion is then solved on
a representative fine grid domain involving the explicitly represented fractures to obtain

and

relationship with coarse scale saturation. The single-phase upscaling procedure
is discussed in Chapter 2.

106

6.2 Multi-phase flow upscaling under depletion condition
Depletion of oil reservoirs accompanies gas phase appearance at the local level. Multi-
phase upscaling is required to model this process adequately in a coarse scale simulation.
Flow from a reservoir into hydraulic fractures is often linear and the use of linear boundary
conditions can be quite suitable. A cross-section of a reservoir domain involving a
horizontal well and a single hydraulic fracture, represented by 9 coarse gridblocks, is
shown in Figure ‎ 6-1.
Fracture
Well

Figure ‎ 6-1: Cross sectional view of a fractured well in a domain with 9 coarse blocks.

The phase flux at the interface of the coarse gridblock ( ) and ( ) is calculated
as the sum of fine scale flux at the interface:

∑ ∑ ∑ (

)

[

g(

)] (6.7)

107

Here (

)

and

are the phase mobility and density at the interface of the fine cell
( ) and ( ) that are calculated with upstream weighting.

is the
transmissibility between fine cell ( ) and ( ) defined as:

(

)

(6.8)
Given the average values of gridblock properties, phase fluxes can be estimated as
following:

[
̅

̅

g(

)] (6.9)
Here,

is the transmissibility at the interface of coarse gridblock ( ) and
( ) that is calculated from single-phase upscaling. Average phase pressure in a
coarse gridblock ( ) is calculated as the corresponding phase pore volume weighted
average of the underlying fine cells:

̅

∑ ∑ ∑

∑ ∑ ∑

(6.10)
Here, denote fine cell indexes in and direction and

,

and

are the number
of fine cells in a corresponding direction inside the coarse gridblock. Average viscosity is
similarly calculated as:

∑ ∑ ∑

∑ ∑ ∑

(6.11)

108

In this work, phase saturation is averaged based on pore volume:

̅

∑ ∑ ∑

∑ ∑ ∑

(6.12)
6.2.1 Workflow for multiphase flow upscaling of multi-fractured well
Selection of gridblock sizes for the simulation stage of a reservoir study depends on the
areal size of the reservoir and available computational resources. Given the size of the
coarse gridblocks, following multi-phase upscaling workflow is proposed:
1. Construct an extended local well model that involves a fracture and at least 3
gridblocks at each direction.
2. Simulate the depletion process on the extended local well model at fine scale
a. Apply no-flow boundary condition on the external boundaries
b. Reduce well pressure gradually to very low values such that all possible in-
situ pressure/saturation conditions are represented.
3. Use equation 6.9 to calculate the average saturation and gas-oil relative
permeability of the coarse grid blocks
a. Fit relative permeability “data” to a orey -type model
b. Cluster derived curves into few representative ones, e.g. fractured and non-
fractured block curves
4. Classify directional relative permeability data for fractured and non-fractured
gridblocks

109

6.3 Coarse scale simulation of horizontal well with a single fracture
A simple homogeneous model is considered to evaluate the coarse grid performance. Model
dimensions and properties are listed in Table ‎ 6-1. A hydraulic fracture extends in the Y-Z
plane at the middle of the model as shown in Figure ‎ 6-1. In the fine scale simulation, the
fracture is explicitly represented by a gridbloc si e of . ’ ’ and ’ in the X -, Y- and Z-
directions, respectively. Multi-phase flow in the matrix is represented by Corey-type
relative permeability, as shown in Figure ‎ 6-2 and Figure ‎ 6-3, and a Corey-Brooks type
capillary pressure curve (Figure ‎ 6-18). The initial saturation in the model is calculated by
the gravity-capillary equilibration and oil and gas PVT properties are reported in
Figure ‎ 6-5 and Figure ‎ 6-6. Depletion of the model is simulated with the bottom-hole
pressure trend as shown in Figure ‎ 6-4.
Table ‎ 6-1: Homogeneous test model properties

Model dimension: Lx, Ly, Lz (ft) 150,150, 150
Areal permeability (md) 1
Vertical permeability (md) 1
Porosity ( ) 0.3
Fracture width (ft) 0.1
Fracture conductivity (md.ft) 500

110

Figure ‎ 6-2: Gas-oil relative permeability in
the test model

Figure ‎ 6-3: Water-oil imbibition relative
permeability in the test model

Figure ‎ 6-4: Step-wise bottom-hole pressure reduction in the depletion simulation of
homogeneous test model

111

Figure ‎ 6-5: Oil PVT properties used in the test model

Figure ‎ 6-6: Gas PVT properties used in the test model

112

6.3.1 Explicit simulation of hydraulic fracture
Accurate modeling of comprehensively characterized multi-fractured well is achievable by
explicit representation of the fractures. This is accomplished by use of extremely small grid
sizes, on the order of the fracture width, and modification of the relative permeability and
capillary pressure curves in the fractured gridblocks. Typically, linear relative permeability
curves and zero capillary pressure are used to represent multi-phase flow in hydraulic
fractures. In the following studies, explicit simulation provides the reference solution for
the comparison purpose.
6.3.2 Coarse grid with local grid refinement
As fine-scale explicit solution is often not available, a coarse model with Local Grid
Refinement (LGR) is considered as the base solution and used to evaluate other coarse
scale simulation techniques. In this approach, a coarse grid is locally refined at the vicinity
of the hydraulic fracture and the well. Figure ‎ 6-7 and Figure ‎ 6-8 report the accuracy of the
LGR, compared to the fine-scale simulation, and the sensitivity to the smallest gridblock
size. Coarse grid simulation with LGR underestimates early production slightly and under-
predicts the gas production unless sufficiently small size gridblock in the direction
perpendicular to fracture face is used.

113

Figure ‎ 6-7: Oil production rate simulated with explicit and LGR representation of fracture

Figure ‎ 6-8: Production gas-oil ratio simulated with explicit and LGR representation of
fracture

114

6.3.3 Coarse grid with transmissibility modification
Next, we investigate the performance of uniform coarse grid performance in the simulation
of production by natural depletion. The fine grid model of the test case in section ‎ 6.3 is
coarsened to a 3x3x5 grid system. Calculation of the fractured gridblock transmissibility
was discussed in chapter 5. When multi-phase upscaling is not applied, matrix type relative
permeability is used for the fracture gridblocks. We will show in the next section ( ‎ 6.3.4) the
coarse-scale simulation with transmissibility modification overestimates the early time oil
production and gas production is under-predicted when near-well pressure falls below
saturation pressure.

6.3.4 Coarse grid with multi-phase upscaling
The multi-phase upscaling procedure, discussed in section ‎ 6.2, is applied to the depletion
process of the homogeneous model. Coarse-scale gas-oil relative permeability curves
calculated in this manner reveal interesting trends as shown in Figure ‎ 6-9 to Figure ‎ 6-13.
In the well gridblock, the relative permeability endpoint is reduced from 1 to 0.75 and has
near-linear trend. In the contrary, gas pseudo relative permeability becomes more non-
linear compared to the original form. As the fracture extends in the Y-Z plane, and flow
occurs linearly along the X-direction into the fracture, the X-direction pseudo-relative
permeability to oil and at coarse scale is systematically increased. In the Y- and Z-
directions, oil and gas pseudo-relative permeability scales differently in the fractured and
non-fractured coarse gridblocks. The Y-direction gas-oil relative permeability in non-
fractured coarse blocks scales similar to the X-direction counterparts, while it is very
scattered in the fractured blocks. In the Z-direction, due to significance of gravity forces in

115

the fracture, the spread in relative permeability data is large and it becomes very difficult
to obtain a single curve. Relative permeability and average saturation history of each
coarse gridblock is fitted to the Corey-type model and normalized to the largest of the
endpoint relative permeabilities (defined here as the scaling ratio). The scaling ratio is then
used as an additional transmissibility multiplier. Therefore, equation 6.9 reduces to the
following form:

[
̅

̅

g(

)] 6.12
Here is the directional scaling ratio. As there is significant spread in the relative
permeability data of fractured blocks, we used relative permeability history of individual
fractured coarse gridblocks to build gridblock specific curves.

Figure ‎ 6-9: Relative permeability in the well
gridblock

Figure ‎ 6-10: Relative permeability in X
direction in all gridblocks

116

Figure ‎ 6-11: Relative permeability in Y
direction in non-fractured gridblocks

Figure ‎ 6-12: Relative permeability in Y
direction in fractured gridblocks (3 layers)

Figure ‎ 6-13: Relative permeability in Z
direction in non-fractured gridblocks

Figure ‎ 6-14: Relative permeability in Z
direction in fractured gridblocks (3 layers)

117

Coarse-scale simulation of the pressure depletion scenario, discussed above, with use of the
modified transmissibility and modified (pseudo) relative permeability curves are shown in
Figure ‎ 6-15 and Figure ‎ 6-16. Oil production is simulated very accurately by the coarse grid
with multi-phase upscaling, in comparison to the coarse grid with transmissibility
modification, and prediction of gas production is significantly improved.

Figure ‎ 6-15: Oil production rate simulated with explicit and coarse representation of
fracture with transmissibility alteration (TRAN) and with multi-phase upscaling (MPU)

118

Figure ‎ 6-16: Production gas-oil ratio simulated with explicit and coarse representation of
fracture with transmissibility alteration (TRAN) and with multi-phase upscaling (MPU).

6.4 A history matching study
Reservoir models should be calibrated to available production data in order to be used for
reservoir studies (forecasting, development, etc). To this end, a reservoir model should
reasonably match oil, water and gas production and satisfy the available pressure
measurements. In the history matching phase of reservoir studies, the total or one fluid
rate is commonly enforced as a simulation constraint and simulated fluid ratios and
pressure data are compared with the observed field data.
Reservoir characterization, as an important stage of reservoir studies, was
discussed in chapter 5. Dynamics of reservoir rock and hydraulic fractures must be
modeled for accurate simulation of multi-fractured wells in stress-sensitive reservoirs. This

119

is generally the case with majority of low permeability and unconventional reservoirs.
Lolon et al. (2007) indicate that explicit simulation of a fracture is necessary to capture the
post-stimulation well clean-up. The post-fracturing well performance depends, in addition
to reservoir properties, on fracture conductivity, fracture dimensions and fracturing fluid
filtrate distribution around the fracture. Initial multi-phase fluid production predictions
can be erroneous if these are not captured correctly in by the simulation model (Lolon et
al., 2007). They noted that rapid production decline cannot be history-matched by constant
fracture permeability and that stress-dependent permeability must be used to model the
deteriorating fracture and matrix permeability.
We note that stress-dependent permeability is not sufficient to predict the rapid
production decline in the shale-oil reservoirs. This was discussed with observations from
lab measurement and transient test interpretations [Clarkson (2013), De Souza (2012),
Weaver et al. (2007), Weaver et al. (2010)]. To our knowledge, there is not any time-lapse
conductivity model or correlation for productivity modeling. In this work, we assume that
fracture conductivity varies exponentially with time:

exp [ ] , (6.13)
where is a model (adjustable) parameter and

is the initial conductivity of the fracture.
As these parameters are primarily dependent on the type, size and concentration of the
proppants used, one may attempt to find average values of these parameters by history
matching of productivity of wells that are fractured with similarly proppants and apply to
the new wells. Note that, in oil reservoirs, the well bottom-hole pressure drops quickly to a
fixed value and does not vary considerably unless artificial lift is installed. Therefore,

120

equation 6.13 assumes that loss of fracture conductivity occurs gradually during the well
life while the well pressure remains constant.
6.4.1 Multi-fractured horizontal well history-matching case study
A multi-fractured horizontal well in the Monterey formation is considered. In the reservoir,
rock and fluid properties have been adequately characterized and indicate a layered setting
with poor vertical communication. Average petrophysical properties of the reservoir are
listed in the Table ‎ 6-2.
Table ‎ 6-2: Petrophysical properties of the case study
Layer
Horizontal
permeability, K H
(md)
Porosity K V/K H
1 0.4 0.3 0.01
2 1.4 0.35 0.01
3 1.4 0.35 0.01
4 0.01 0.25 0.01
5 0.01 0.25 0.01

121

The water-oil relative permeability curve for this water-wet reservoir is shown in
Figure ‎ 6-17. Using the Leveret J-function, measured capillary pressure is scaled with
porosity and permeability and used to initialize the fluid saturation (Figure ‎ 6-18). Oil and
gas PVT properties are shown in Figure ‎ 6-19. This is an under-saturated oil reservoir with
a bubble-point pressure of 2960 psia.

Figure ‎ 6-17: Oil-water system relative
permeability

Figure ‎ 6-18: J-function from capillary
pressure data of the layers in the case study

Lab measurement indicates that the reservoir rock is considerably compressible
(Figure ‎ 6-20). This is modeled in the reservoir simulation with the use of pressure-
dependent porosity and permeability multiplier as shown in Figure ‎ 6-21.

122

Figure ‎ 6-19: Oil and gas PVT properties of the field example

Figure ‎ 6-20: Stress-dependent rock
compressibility in the case study

Figure ‎ 6-21: Stress-dependent rock
porosity and permeability input to the
simulator

10.00
12.00
14.00
16.00
18.00
20.00
22.00
24.00
26.00
0 2000 4000 6000
Pore compressibility, msip
Confining stress
0.88
0.9
0.92
0.94
0.96
0.98
1
0 2000 4000 6000
Multiplier
Pressure, psia
Porosity
Permeability

123

A 3200 ft horizontal well with 8 stages of hydraulic fracturing is considered in this
simulation study. Orientation and geometry of the hydraulic fractures are not known and is
inferred here by using the production data. We use water cut data to estimate the height of
the fractures while total liquid production is used to estimate length and conductivity of
fractures. As this is an iterative procedure, analytic solution for fractured horizontal well
can be used to estimate productivity and thus inversely calculate the fracture properties.
Due to several assumptions (limitations) of these techniques, this solution is considered as
an initial guess for the more realistic numerical simulation model. To refine the solution, a
single well simulation model with grid refinement in the vicinity of the well and hydraulic
fractures is constructed and an inverse problem is defined that minimizes the mismatch
between observed and calculated production rate to estimate the relevant model
parameters.
Production history matching of this example is performed by constraining the
simulator to the oil production rates and comparing the available measurements (water cut
and GOR) with the simulated values. We observe that dynamic fracture conductivity must
be used to match the initial high production rates and later-time low flow rates while also
reproduce the bottom-hole pressure and gas and water production.
Transmissibility of the fractured gridblocks is reduced at monthly time-steps based
on the equation 6.13. In this process, we estimated the initial conductivity value of 300 md-
ft and an value of 0.1 day
-1
. As

the fracture geometry was not known, an inverse
procedure discussed above was applied and fracture half-length of 250 ft and height of 280
ft was approximated. The height of the fractures was estimated by matching the water-cut
trend. Involved layers with varying properties are characterized in the vicinity of this well

124

and feed into the fracture based on their flow capacity (permeability × thickness). Layers 2
and 3 have larger oil saturation and contribute more to the initial production while layers 1
and 4 deplete at later-time resulting in increasing trend of water production.

Figure ‎ 6-22: Oil production rate observed
and constrained in the simulator
Figure ‎ 6-23: Producing water cut
observation and simulated values

Bottom-hole pressure measurements are not available for this well. However, it
produces with natural drive at a constant wellhead pressure of about 200 psia and an
average bottom-hole pressure can be estimated by adding the wellhead pressure to the
well column hydrostatic pressure. This calculates an average well flowing pressure of
about 3500 psi. Simulated bottom-hole pressure, as shown in Figure ‎ 6-24, is in good
agreement with the calculated average value until artificial lift is installed in 2002 that
causes the bottom-hole pressure to fall below the calculated average value.
Simulated and observed producing gas-oil ratio is reported in Figure ‎ 6-25.
Simulation of gas production is primarily controlled by the well and reservoir pressure and
the accuracy of the gas-oil relative permeability curve. Since depletion type gas-oil relative

125

permeability measurements are not available for this reservoir, we synthesized a Corey-
type relative permeability curve that approximates the gas production.

Figure ‎ 6-24: Well bottom-hole pressure
from simulation

Figure ‎ 6-25: Producing gas-oil ratio
simulated and observed on the well

6.5 Concluding remarks
In the second part of the dissertation, modeling and coarse scale simulation of a multi-
fractured horizontal well is studied. A multi-phase upscaling workflow was adapted for the
coarse scale simulation of multi-fractured horizontal wells. Coarse scale simulation with
multi-phase upscaling indicated that relative permeability curves scale differently in
different directions. In the direction perpendicular to the fracture, relative permeability in
the fractured and non-fractured gridblocks scale similarly and thereby a single curve can
be used for all gridblocks. Parallel to the fracture plane, relative permeability in the
fractured and non-fractured gridblocks are different and at least two sets of curves must be
used for each direction. We observe that vertical flow in the fractured gridblocks is not

126

consistently captured by this workflow and further studies are required to investigate the
effect of gravity on the coarse scale relative permeability scaling.
A real-life example of multi-fractured horizontal well productivity is studied. We
observe that dynamic fracture conductivity must be used to match the initial high
production rates and later-time low flow rates and simultaneously reproduce the bottom-
hole pressure and gas and water production. Based on the assumption that fracture
conductivity loss occurs gradually during well life, we used a simple two-parameter
exponential function and inferred the parameter values via history matching process.

127

Chapter 7
Future directions

128

7 Future directions

In the first part of the dissertation, transformation of fine-scale heterogeneous geomodels
into a dual-continuum model is accomplished by utilizing streamline information to
identify key flow features. However, streamline tracing on geological models with very
large number of fine cells can be computationally demanding. We adapt a dual-porosity
dual-permeability model to simulate the relevant mass transfer phenomena. Simple
calculations were used to obtain coarse scale inter-gridblock transmissibility and we
assumed that fine cells comprising the coarse gridblock primary and secondary porosity is
continuous and did not distinguish between connected and disconnected fine cells at each
level of porosity. In addition, our simulation studies indicated that dual-porosity coarse
models are significantly less sensitive to the level of upscaling, compared to conventional
single-porosity upscaling. Therefore, following studies can be performed to improve the
accuracy or highlight the features of the proposed dual-porosity upscaling technique:

1. Simple static or dynamic measures of heterogeneity may be used, as a secondary
alternative, to highlight the major flow pathways. Performance of simpler measures
of heterogeneity in the transformation of a geological model into a coarse dual-
continuum can be investigated in a separate study.
2. In this study, two levels of porosity at each coarse gridblock were considered.
Similar idea can be applied with multiple-porosity modeling. In addition, by
distinguishing between connected and disconnected fine cells at each level of

129

porosity, the sub-grid effects on local displacement efficiency can then be
represented more adequately.
3. The sensitivity of the proposed technique to the upscaling level should be
investigated by further numerical tests. Upscaling errors due to loss of
heterogeneity and numerical dispersion in single-porosity and dual-porosity
modeling is different. At this point, we do not have adequate understanding of the
interplay of these two sources of errors on the dual-porosity modeling. Investigation
of this subject can provide guidelines for the selection of an optimal upscaling level.
4. Field scale evaluation of EOR processes necessitates full three-dimensional
compositional and/or black-oil modeling and simulation. We did not include the
effect of gravity and capillarity and a separate study can investigate the
performance of proposed technique in the presence of these forces.

Modeling and coarse scale simulation of multi-fractured horizontal wells was studied in the
second part of the dissertation. We demonstrated that multi-phase upscaling can improve
the accuracy of the coarse scale simulation techniques significantly. Several aspects of the
proposed workflow can be investigated. Dynamic modeling of fracture as a key component
impacting the productivity of multi-fractured wells was studied. We noted that dynamic
fracture conductivity must be used to match the initial high production rates and later-time
low flow rates while simultaneously reproducing the bottom-hole pressure and gas and
water production. Additional research in this area could include:

130

1. Simple averaging was used to calculate the coarse grid block properties. Moreover,
average gridblock saturation was used to calculate gridblock interface properties. A
comprehensive study of calculation of effective properties can potentially improve
the accuracy of this approach.
2. Pseudo-relative permeability curves calculate for fractured gridblocks indicated a
significant spread. Improvements in the calculation of these curves can further
enhance the improvements introduced by this approach.
3. In layered reservoirs, coarse scale simulation with multi-phase upscaling can be
more appealing. We considered a homogeneous case only and further studies can
provide more insight into the applicability of this approach in heterogeneous
systems.

131

Bibliography

Audigane, P., & Blunt, M. J. (2004). Dual mesh method for upscaling in waterflood
simulation. Transport in Porous Media, 55(1), 71-89.

Aziz, K., & Settari, A. (1979). Petroleum reservoir simulation (Vol. 476). London: Applied
Science Publishers.

Babaei, M. (2013). Multiscale wavelet and upscaling-downscaling for reservoir
simulation (Doctoral dissertation, Imperial College London).

Bar er J.W. and hibeau .: “ ritical eview of the Use of seudorelative ermeabilities
for Upscaling ” (May 997) 9.

Barree, R.D., Fisher, M.K., and Woodroof, R.A.: "A Practical Guide to Hydraulic Fracture
Diagnostic Techniques," paper SPE 77442 presented at the SPE ATCE; San Antonio, TX,
USA; Oct 2002

Bonney, K., Abacioglu, Y., & Sebastian, H. (2013, January). Case History of Using Integrated
Reservoir Modeling Workflow for Tight Gas Reservoirs. In 2013 SPE Middle East
Unconventional Gas Conference & Exhibition.

Brown, M., Ozkan, E., Raghavan, R., Kazemi, H.: Practical solutions for pressure transient
responses of fractured horizontal wells in unconventional reservoirs. Paper SPE 125043
presented at the SPE Annual Technical Conference and Exhibition, New Orleans, Louisiana,
4–7 October 2009. doi:10.2118/125043-MS

Cao, H., & Aziz, K. (1999). Evaluation of pseudo functions. In SPE Western Regional
Meeting.

Castellini, A.: Flow based grids for reservoir simulation. MS Dissertation, Stanford
University (2001)

Chen, Y., Durlofsky, L.J.: Adaptive local–global upscaling for general flow scenarios in
heterogeneous formations. Transp. Porous Media 62, 157–185 (2006)

132

Chen, Y., Durlofsky, L. J., Gerristen, M. and Wen, X. H.: 2003, A coupled local–global
upscaling approach for simulating flow in highly heterogeneous formations, Adv. Water
Resour. 26, 1041–1060.

Chen, Y., Durlofsky L.J., Wen X.H.: Robust Coarse Scale Modeling of Flow and Transport in
Heterogeneous Reservoirs. 9th European Conference on the Mathematics of Oil Recovery,
Cannes, France (2004)

Chen, Yuguang, and Louis J. Durlofsky. "Efficient incorporation of global effects in upscaled
models of two-phase flow and transport in heterogeneous formations." Multiscale
Modeling & Simulation 5.2 (2006): 445-475.

Chen, Y., Li, Y., & Efendiev, Y. (2013). Time-of-flight (TOF)-based two-phase upscaling for
subsurface flow and transport. Advances in Water Resources.

Christie, M. A., and Clifford, P. J., 1997, A fast procedure for upscaling in compositional
simulation, SPE 37986 presented at the SPE Symposium on Reservoir Simulation, Dallas,
TX, June 8–11.

Christie, M. A., and Blunt, M. J., 2001, Tenth SPE comparative solution project: A comparison
of upscaling techniques: SPE Res. Eval. Eng., v. 4, no. 4, p. 308–317.

Cinco-Ley, H., & Samaniego-V, F. (1981). Transient pressure analysis for fractured
wells. Journal of petroleum technology, 33(9), 1749-1766.

Clarkson, C.R., Qanbari, F., Nobakht, M. et al. 2013. Incorporating Geomechanical and
Dynamic Hydraulic-Fracture-Property Changes Into Rate-Transient Analysis: Example
From the Haynesville Shale. SPE Res Eval & Eng 16 (3): 303-316. SPE-162526-PA.

Comisky, J., Newsham, K., Rushing, J., & Blasingame, T. (2007, November). A comparative
study of capillary-pressure-based empirical models for estimating absolute permeability in
tight gas sands. In SPE Annual Technical Conference and Exhibition.

Cvetkovic, B., Norge, B., (2009, April). Effective Wellbore Parameters and A Multifractured
Horizontal Well Productivity. In SPE Production and Operations Symposium.

Datta-Gupta, A. (2000). Streamline simulation: a technology update. Journal of petroleum
technology, 52(12), 68-73.

133

Davies, D. R., & Kuiper, T. O. H. (1988). Fracture conductivity in hydraulic fracture
stimulation. Journal of petroleum technology, 40(5), 550-552.

Diaz De Souza, O., Sharp, A., Martinez, R., Foster, R., Piekenbrock, E., Reeves Simpson, M., &
Abou-Sayed, I. (2012, June). Integrated Unconventional Shale Gas Reservoir Modeling: A
Worked Example From the Haynesville Shale, De Soto Parish, North Louisiana. In SPE
Americas Unconventional Resources Conference.
Durlofsky L. J., Numerical calculation of equivalent gridblock permeability tensors for
heterogeneous porous media, Water Resources Research 27 (1991) 699–708.

Durlofs y .J. Jones . . and Milli en W.J.: “ New Method for the Scale Up of
Displacement rocesses in Heterogeneous eservoirs ” roceedings of the
th
European
conference on the Mathematics of Oil Recovery, Roros, Norway, June 1994.

Durlofsky, L.J., Behren, R.A., Jones, R.C.: Scale up of heterogeneous three dimensional
reservoir description. SPEJ 1, 313–326 (1996)

Durlofsky L. J., Upscaling of geocellular models for reservoir flow simulation: a review of
recent progress. 7
th
International Forum on Reservoir Simulation, BÄuhl/Baden-Baden,
Germany, June 23-27, 2003.

Elgmati, M., Zhang, H., Bai, B., Flori, R., & Qu, Q. (2011, June). Submicron-pore
characterization of shale gas plays. In North American Unconventional Gas Conference and
Exhibition.

Evazi, M., & Mahani, H. (2010). Generation of Voronoi grid based on vorticity for coarse-
scale modeling of flow in heterogeneous formations. Transport in Porous Media, 83(3),
541-572.

Evazi, M., & Mahani, H. (2010). Unstructured-coarse-grid generation using background-grid
approach. SPE Journal, 15(2), 326-340.

Evazi Yadecuri, M. & Jessen, K., (2013, April). Simulation of Waterflooding with Coarse-
Scale Dual-Porosity Representation of Highly Heterogeneous Reservoirs. In SPE Western
Regional & AAPG Pacific Section Meeting, 2013 Joint Technical Conference.

134

Farmer C. L., Upscaling: a review, International Journal for Numerical Methods in Fluids 40
(2002) 63–78

Fetkovich, M. J. (1980). Decline curve analysis using type curves. Journal of Petroleum
Technology, 32(6), 1065-1077.

Fjerstad, P., Dasie, W., Sikandar, A., Cao, H., & Liu, J. (2005, December). Next generation
parallel computing for large-scale reservoir simulation. In SPE International Improved Oil
Recovery Conference in Asia Pacific.

Gautier, Y., Blunt, M. J., & Christie, M. A. (1999). Nested gridding and streamline-based
simulation for fast reservoir performance prediction.Computational Geosciences, 3(3-4),
295-320.

Garcia, M.H., Journel, A.G., Aziz, K.: Automatic grid generation for modeling reservoir
heterogeneities. SPERE 7, 278–284 (1992)

Gilman, J. R., & Kazemi, H. (1988). Improved Calculations for Viscous and Gravity
Displacement in Matrix Blocks in Dual-Porosity Simulators (includes associated papers
17851, 17921, 18017, 18018, 18939, 19038, 19361 and 20174). Journal of petroleum
technology, 40(1), 60-70.

Grader, A., Mu, Y., Toelke, J., Baldwin, C., Fang, Q., Carpio, G., ... & Kalam, M. (2010,
November). Estimation of relative permeability using the Lattice Boltzmann method for
fluid flows, Thamama Formation, Abu Dhabi. In Abu Dhabi International Petroleum
Exhibition and Conference.

Guérillot, D. R., & Verdiere, S. (1995, February). Different pressure grids for reservoir
simulation in heterogeneous reservoirs. In SPE Reservoir Simulation Symposium.

Guo, B., Yu, X., and Khoshgahdam, M. 2009. A Simple Analytical Model for Predicting
Productivity of Multifractured Horizontal Wells. SPE Reservoir Evaluation & Engineering
12(6): 879-885.

Haajizadeh M., Fayers F. J. and Cockin A. P.: 2000, Effects of phase behaviour, dispersion
and gridding on sweep patterns for nearly miscible gas displacement, SPE 62995,
presented at the SPE Annual Technical Conference, Dallas, Texas, 1–4 October, 2000

135

Hesse, M. A., Mallison, B. T., & Tchelepi, H. A. (2008). Compact multiscale finite volume
method for heterogeneous anisotropic elliptic equations. Multiscale Modeling &
Simulation, 7(2), 934-962.

Hill, A. C., & Thomas, G. W. (1985, March). A new approach for simulating complex
fractured reservoirs. In Middle East Oil Technical Conference and Exhibition.

Holden L., Nielsen B. F., Global upscaling of permeability in heterogeneous reservoirs: The
Output Least Squares (OLS) method, Transport in Porous Media 40 (2000) 115–143.

Hong, K. C. (1982, April). Lumped-component characterization of crude oils for
compositional simulation. In SPE Enhanced Oil Recovery Symposium.

Jack, H., Smith, O., & Mattax, C. C. (1973). The modeling of a three-dimensional reservoir
with a two-dimensional reservoir simulator-the use of dynamic pseudo functions. Old SPE
Journal, 13(3), 175-185.

Jessen, K., & Rastegar Moghadam, R. (2009, February). A Flow Based Lumping Approach for
Compositional Reservoir Simulation. In SPE Reservoir Simulation Symposium.

Johns . . ayers .J. and Orr .M. Jr.: “ ffect of Gas nrichment and Dispersion on Near ly
Miscible Displacements in ondensing/ apori ing Drives ” dvanced echnology
Series (April 1993) 26.

Kyte, J. R., & Berry, D. W. (1975). New pseudo functions to control numerical
dispersion. Old SPE Journal, 15(4), 269-276.

Kazemi, H., Merrill, L.S., Porterfield, K.L., and Zeman, P.R.: "Numerical Simulation of Water-
Oil Flow in Naturally Fractured Reservoirs," SPEJ (Dec. 1976) 317-326.

Lake, L. W. (1989). Enhanced oil recovery.

Li, D., Cullick, A.S., Lake, L.W.:Global scale-up of reservoirmodel permeability with local grid
refinement. J. Petrol. Sci. Eng. 14, 1–13 (1995)

Lolon, E., Shaoul, J., & Mayerhofer, M. (2007, October). Application of 3-D Reservoir
Simulator for Hydraulically Fractured Wells. In Asia Pacific Oil and Gas Conference and
Exhibition.

136

Lu, B., Alshaalan, T., & Wheeler, M. (2007, November). Iteratively coupled reservoir
simulation for multiphase flow. In SPE Annual Technical Conference and Exhibition.

Mahani, H., Muggeridge, A.H.: Improved coarse grid generation using vorticity. Paper SPE
94319 presented at the of the 14 Biennal SPE/Europec Conference, Madrid, Spain, June
2005

Medeiros, F., Ozkan, E., and Kazemi, H. 2008. Productivity and Drainage Area of Fractured
Horizontal Wells in Tight Gas Reservoirs. SPE Res Eval & Eng11 (5)

Muggeridge A.H., Generation of effective relative permeabilities from detailed simulation of
flow in heterogeneous porous media. In Reservoir Characterization II, ed. L.W. Lake, H.B.
Carroll & T.C. Wesson. Academic Press, San Diego, 1991.

Muskat, M. 1949. Physical Principles of Oil Production. McGraw-Hill, New York.

Nakashima, T., & Durlofsky, L. J. (2010). Accurate representation of near-well effects in
coarse-scale models of primary oil production. Transport in Porous Media, 83(3), 741-770.

Okouma, V., Guillot, F., Sarfare, M., Sen, V., Ilk, D., & Blasingame, T. A. (2011). Estimated
ultimate recovery (EUR) as a function of production practices in the Haynesville
shale. SPE, 147623, 2011-11.

Palish T., Duenckel R., Bazan L., Heidt H., Turk G. Determining Realistic Fracture
Conductivity and Understanding its Impact on Well Performance — Theory and Field
Examples, Society of Petroleum Engineers, SPE (2007), p. 106301

Pedrosa O.A. Jr, Pressure transient response in stress-sensitive formations, Paper SPE
15115. Presented at the 1986 California Regional Meeting, Oakland, CA, April (1986)

Pedroso, C.A and Correa, A.C.F, A New Model of a Pressure-Dependent-Conductivity
Hydraulic Fracture in a Finite Reservoir: Constant Rate Production Case, SPE 38976-MS,
Presented at the Latin American and Caribbean Petroleum Engineering Conference, Rio de
Janeiro, Brazil, September (1997).

eng D.Y. and obinson D.B.: “ New wo - onstant quation of tate ” nd. ng. hem.
Fund. (1976) 15, 59.

137

Prévost, M., Lepage, F., Durlofsky, L. J., & Mallet, J. L. (2005). Unstructured 3D gridding and
upscaling for coarse modelling of geometrically complex reservoirs. Petroleum
Geoscience, 11(4), 339-345.

Rame, M., & Killough, J. E. (1992). A new approach to flow simulation in highly
heterogeneous porous media. SPE formation evaluation, 7(3), 247-254.

Renard P. and de Marsily G., (1997), Calculating equivalent permeability: a review. Advance
in Water Resources.

Rojas, D., Li, H., Kumar, M., & Chen, Y. (2013, April). Development and Application of Near-
Well Multiphase Flow Upscaling for Forecasting of Heavy Oil Primary Production. In SPE
Western Regional & AAPG Pacific Section Meeting, 2013 Joint Technical Conference.

Saad, N., Maroongroge, V., & Kalkomey, C. T. (1996, April). Ranking geostatistical models
using tracer production data. In European 3-D Reservoir Modelling Conference.

Sablok, R., & Aziz, K. (2008). Upscaling and discretization errors in reservoir
simulation. Petroleum Science and Technology, 26(10-11), 1161-1186.

Sadrpanah, H., Charles, T., & Fulton, J. (2006, June). Explicit simulation of multiple hydraulic
fractures in horizontal wells. In SPE Europec/EAGE Annual Conference and Exhibition.

ho aei H. astegar . and Jessen : “Mixing and Mass r ansfer in Multicontact Miscible
Displacements” rans port in Porous Media, 2012, Volume 94, Number 3, Pages 837-57.

Shojaei H. and Tajer E. (2013), Analytical Solution of Transient Multiphase Flow to a
Horizontal Well with Multiple Hydraulic Fractures, SPE 165703, presented at the SPE
Eastern Regional Meeting, Pittsburgh, Pennsylvania.

Shook, G., & Mitchell, K. (2009, October). A Robust Measure of Heterogeneity for Ranking
Earth Models: The F PHI Curve and Dynamic Lorenz Coefficient. InSPE Annual Technical
Conference and Exhibition.

Stalkup, F. I. (1983), Miscible displacement. SPE Monograph series, Society of Petroleum
Engineers of AIME.

138

Stiles, W. M. (1949). Use of permeability distribution in water flood calculations.Journal of
Petroleum Technology, 1(1), 9-13.

Stone, H. L. (1991, February). Rigorous black oil pseudo functions. In SPE Symposium on
reservoir simulation.

Swanson, B. F. (1981). A simple correlation between permeabilities and mercury capillary
pressures. Journal of Petroleum Technology, 33(12), 2498-2504.

Thompson, J., Nobakht, M., & Anderson, D. (2010, October). Modeling well performance
data from overpressured shale gas reservoirs. In Canadian Unconventional Resources and
International Petroleum Conference.

Tölke, J., Baldwin, C., Mu, Y., Derzhi, N., Fang, Q., Grader, A., & Dvorkin, J. (2010). Computer
simulations of fluid flow in sediment: From images to permeability. The leading
edge, 29(1), 68-74.

Wallstrom T. C., Hou S., Christie M. A., Durlofsky L. J., Sharp D. H., Zou Q., Application of
effective flux boundary conditions to two-phase upscaling in porous media, Transport in
Porous Media 46 (2002) 155–178.

Wan, J., Penmatcha, V. R., Arbabi, S., & Aziz, K. (2000). Effects of Grid Systems on Predicting
Horizontal-Well Productivity (includes associated papers 65961 and 66620). SPE
Journal, 5(3), 309-314.

Warren, J. E., & Price, H. S. (1961). SPE-001579-G-Flow in Heterogeneous Porous
Media. Old SPE Journal, 1(3), 153-169.

Warren, J.E. and Root, P.J.: "The Behavior of Naturally Fractured Reservoirs," SPEJ (Sept.
1963) 245-255.

Weaver, J., Parker, M., van Batenburg, D., & Nguyen, P. (2007). Fracture-Related Diagenesis
May Impact Conductivity. SPE Journal, 12(3), 272-281.

139

Weaver, J., Rickman, R., and Luo, H. 2010. Fracture-Conductivity Loss Caused by
geochemical Interactions Between Man-Made Proppants and Formations. SPE J. 15 (1):
116-124. SPE-118174-PA. doi: 10.2118/118174-PA.

White C. D., Horne R. N., Computing absolute transmissibility in the presence of fine-scale
heterogeneity, Paper SPE 16011 presented at the SPE Symposium on Reservoir Simulation,
San Antonio, TX, Feb. 1-4, (1987).

Wei, Y., and Economides, M. 2005. Transverse Hydraulic Fractures from a Horizontal Well.
Paper SPE 94671 presented at the SPE Annual Technical Conference and Exhibition, Dallas,
TX, 9-12 October.

Yuan, H., and Zhou, D. 2010. A New Model for Predicting Inflow Performance of Fractured
Horizontal Wells. Paper SPE 133610 presented at the SPE Western Regional Meeting,
Anaheim, CA, 27-29 May.

Zhou, W., Banerjee, R., Poe, B., Spath, J., and Thambynayagam, M. 2012. Semi-Analytical
Production Simulation of Complex Hydraulic Fracture Network. Paper SPE 157367
presented at the SPE International Production and Operations Conference & Exhibition,
Doha, Qatar, 14-16 May.

Abstract (if available)

Abstract Oil and gas reservoirs in fluvial and deltaic depositional systems comprise a significant portion of the so-called conventional resources. In these systems, the inherent reservoir heterogeneity is often the major challenge of exploitation that complicates the application and simulation of enhanced oil recovery techniques. Despite the advancements in the computational technology, fine-scale simulation of stochastically generated geomodels is not feasible as hundreds of realizations are analyzed to quantify the impact of geological uncertainty on the reservoir performance. This has motivated the development of upscaling techniques that, for highly heterogeneous permeability fields, still remains a challenge. ❧ In the first part of this work, the effect of heterogeneity on displacement efficiency in coarse scale modeling is studied. Pore space is ranked based on flow contribution to two levels of porosity and a dual-porosity dual-permeability flow model is adapted for coarse-scale flow simulation. We use fine-scale streamline information to transform heterogeneous geomodels into a dual-continuum coarse model that preserves the global flow pathways adequately. ❧ The proposed technique is applied and tested on two heterogeneous models with different types of fluid flow modeling (compositional and black oil). In order to evaluate the performance of the proposed technique, displacement calculations are performed on the original fine grid and on a uniform coarse grid with single-porosity and dual-porosity upscaling. We demonstrate that dual-porosity coarse models predict the breakthrough time accurately and reproduce the post-breakthrough response adequately while single-porosity models overestimate arrival time and oil recovery. By preserving large scale heterogeneities, dual-porosity coarse models are demonstrated to be significantly less sensitive to the level of upscaling, compared to conventional single-porosity upscaling. This makes the proposed upscaling approach a relevant and suitable technique for upscaling of highly heterogeneous geomodels. ❧ Modeling of multi-fractured horizontal wells is the subject of the second part of this dissertation. Recovery from low and ultra-low permeability reservoirs have been unlocked by the introduction of horizontal well and hydraulic fracturing technology that has increased the world’s hydrocarbon reserve significantly. In order to maximize the productivity from these reservoirs, numerical simulation is often used to investigate the optimized well spacing and stimulation design. Fracture properties vary during the well life and must be properly modeled to enable realistic production simulation. Production history matching of a real-life multi-fractured horizontal well in a liquid-rich Monterey formation reservoir is considered. We infer from numerical simulations that, in these systems, time-dependent fracture conductivity model is required to match early-time flush production and later-time low flow rates. ❧ Accurate representation of the near-well dynamics is key to the realistic simulation of fractured wells. Due to the large size of these reservoirs, efficient numerical simulation requires use of large gridblocks that may not be adequate to capture the near-well dynamics. We adapt a multi-phase upscaling workflow for multi-fractured horizontal wells under natural depletion. Fine-scale simulation of the depletion process on a local well model is used to calculate a new set of gas-oil pseudo-relative permeability curves for fractured and non-fractured gridblocks. Simulation results demonstrate improvements in the accuracy of the coarse scale simulation over existing approaches.

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

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

PDF

Mass transfer during enhanced hydrocarbon recovery by gas injection processes

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

PDF

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

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

PDF

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

PDF

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

PDF

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

PDF

Asphaltene deposition during co₂ injection in enhanced oil recovery applications

PDF

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

PDF

Reactivation of multiple faults in oilfields with injection and production

$Integrated workflow for characterizing and modeling fracture network in unconventional reservoirs using microseismic data$

PDF

Integrated workflow for characterizing and modeling fracture network in unconventional reservoirs using microseismic data

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

PDF

An extended finite element method based modeling of hydraulic fracturing

$Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure$

PDF

Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure

PDF

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

PDF

Real-time reservoir characterization and optimization during immiscible displacement processes

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

PDF

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

PDF

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

PDF

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

PDF

Assessing induced seismicity rate increase based on deterministic and probabilistic modeling

Asset Metadata

Creator Evazi Yadecuri, Mohammad (author)

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

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Publication Date 11/26/2013

Defense Date 10/21/2013

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

Tag dual-porosity upscaling,heterogeneous reservoir,horizontal wells,hydraulic fracturing,OAI-PMH Harvest,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 evaziyad@usc.edu,mevazi@gmail.com

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

Unique identifier UC11297289

Identifier etd-EvaziYadec-2186.pdf (filename),usctheses-c3-351527 (legacy record id)

Legacy Identifier etd-EvaziYadec-2186.pdf

Dmrecord 351527

Document Type Dissertation

Format application/pdf (imt)

Rights Evazi Yadecuri, Mohammad

Type texts

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

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

Repository Name University of Southern California Digital Library

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