University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Modeling and Simulation of Multicomponent Mass Transfer in Tight Dual-Porosity Systems
(Unconventional)
by
Mohammed Samir Raslan
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirementsfor the Degree
DOCTOR OF PHILOSOPHY
(PETROLEUM ENGINEERING)
August 2024
Copyright 2024 Mohammed Samir Raslan

ii
Dedication
To my wife, kids, parents, and siblings.
To my wife, for working overtime and going above and beyond to cover for my shortcomings.
To my kids, for being the source of my joy and the spark that kept me going throughout this journey.
To my parents, for their sleepless nights worrying about me and keeping me in their prayers.
To my siblings, for having my back during my absence.
Thank you all for believing in me.

iii
Acknowledgements
I would like to express my appreciation to Professor Kristian Jessen, my dissertation advisor, for his
guidance and mentorship throughout our work together. His passion for giving back to academia and the
industry has been very inspiring to me.
I would like to extend my gratitude to my defense committee members Professor Roger Ghanem
and Associate Professor Birendra Jha for their insightful feedback on this dissertation. I would like to thank
Professor Iraj Ershaghi, who served on my PhD qualifying exam, for providing me with eye-opening
feedback that helped strengthen my research, as well for the invaluable discussions on non-research
related topics. I also want to thank Professor Donald Paul for serving on my PhD qualifying exam.
I am grateful to Saeed Alahmari for being exceptionally helpful in times of need, and to Ye Lyu,
whom I find to be a very interesting and insightful researcher.
I am indebted to Saudi Aramco for the sponsorship. Special thanks to the President and CEO of
Saudi Aramco, Mr. Amin Al-Nasser, for his support, inspiration, and words of encouragement.

iv
Table of Contents
Dedication.................................................................................................................................ii
Acknowledgements....................................................................................................................iii
List of Tables............................................................................................................................. vi
List of Figures........................................................................................................................... vii
Abstract................................................................................................................................... xii
Chapter 1. Introduction......................................................................................................... 1
1.1. Background and Motivation........................................................................................... 1
1.2. The Need for Advanced Techniques for Modeling of Tight Fractured Reservoirs..................... 3
1.3. Problem Statement ...................................................................................................... 5
1.4. Dissertation Objectives and Deliverables ......................................................................... 6
1.5. Dissertation Organization .............................................................................................. 7
Chapter 2. Literature Review.................................................................................................. 9
2.1. Modeling of Fractured Reservoirs................................................................................... 9
2.2. Multiphase Transfer Functions ......................................................................................16
2.3. Diffusive Mass Transfer................................................................................................18
Chapter 3. A Generalized Dynamic Transfer Function for Tight Dual-Porosity Systems*
...................23
3.1. Derivation of the Generalized Vermeulen Transfer Function...............................................24
3.2. Numerical Examples....................................................................................................29
3.3. Summary and Conclusion .............................................................................................37
Chapter 4. Upscaling of Discrete Fracture Characterization to Dual-Porosity Models Using an
Advanced Transfer Function*
......................................................................................................39
4.1. Importance of Capturing the Transient Behavior ..............................................................40
4.2. Upscaling of Fractured Rock..........................................................................................43
4.3. Application of the Upscaling Workflow to Irregular Fracture Geometries.............................45
4.4. Discussion..................................................................................................................56
4.5. Summary and Conclusion .............................................................................................58
Chapter 5. A Two-Phase Transfer Function for Simulation of Tight Fractured Gas-Condensate
Reservoirs* 59
5.1. A Dual-Porosity Compositional Model in MRST ................................................................60
5.2. Numerical Examples....................................................................................................62
5.3. An Improved Two-Phase Transfer Function for Condensates..............................................67
5.4. Discussion..................................................................................................................75
5.5. Application of the Adaptive Transfer Pressure – Examples 5.3.1 and 5.3.2 Revisited ..............77

v
5.6. Summary and Conclusion .............................................................................................78
Chapter 6. Diffusive Mass Transfer and Recovery During CO2 Huff-n-Puff in an Eagle Ford Shale Core*
80
6.1. Lab Experimental Materials..........................................................................................81
6.2. Lab Experimental Procedure.........................................................................................82
6.3. Experimental Observations and Evaluation......................................................................84
6.4. Numerical Modeling....................................................................................................93
6.5. Validation of The Molecular Diffusion Models .................................................................95
6.6. Modeling of CO2 HnP on the EF-3 Core ...........................................................................97
6.7. Summary, Conclusion, and Future Recommendations.....................................................112
6.8. Appendix .................................................................................................................114
Chapter 7. Summary of Dissertation and Future Work.............................................................121
7.1. Summary of Dissertation............................................................................................121
7.2. Recommendations for Future Work .............................................................................122
References.............................................................................................................................124

vi
List of Tables
Table 3-1: Model parameters for planar, cylindrical, and spherical geometries for isotropic matrix
blocks......................................................................................................................................27
Table 3-2: Generalized Vermeulen transfer function parameters for anisotropic rectangular matrix
blocks......................................................................................................................................29
Table 4-1: Generalized Vermeulen transfer function parameters after fitting the upscaled single-block
DPM to the EDFM for the different EDFM realizations ....................................................................52
Table 4-2: Transfer function parameters obtained by matching the respective DPM to the reference
EDFM solution..........................................................................................................................54
Table 4-3: Transfer function parameters obtained by matching the respective DPM (advanced and
traditional) to the reference EDFM solution ..................................................................................57
Table 5-1: EOS parameters of nC10 used in Example 5.2.1................................................................63
Table 5-2: EOS parameters of the quaternary dry-gas fluid description used in Example 5.2.2...............64
Table 5-3: EOS parameters of the quaternary condensate fluid description used in Example 5.2.3.........65
Table 5-4: Final average matrix composition for all models for Example 5.3.1. ...................................72
Table 6-1: Molecular weight and composition of the synthetic oil.....................................................82
Table 6-2: HnP injection pressure per cycle ...................................................................................82
Table 6-3: List of pressure pulse experiments ................................................................................85
Table 6-4: Oil remaining (in mole fraction) in the core after each HnP cycle........................................90
Table 6-5: EOS parameters of the synthetic oil and CO2 ...................................................................99
Table 6-6: Binary interaction parameters of the oil binaries with CO2 ................................................99
Table 6-7: Multipliers to correct for the infinite dilution diffusion coefficients obtained from Hayduk
and Minhas (1982)..................................................................................................................104

vii
List of Figures
Figure 1-1: U.S. tight oil production from selected plays (EIA 2023) and U.S. dry natural gas production
(EIA – AEO2022)........................................................................................................................ 1
Figure 1-2: Global oil production due to EOR activities with projections from 2018 to 2040 (IEA 2018)
............................................................................................................................................... 3
Figure 1-3: Average matrix pressure response of 1 mD and 0.8 D rocks............................................ 3
Figure 1-4: Average matrix pressure response during the early characteristic time of the 0.8 D rock..... 4
Figure 2-1: Warren and Root (1963) idea of the dual-continuum model with constant shape factors
that reflect the geometry of the fractured system. In later years, time-dependent transfer functions
have been introduced for tight fractured reservoirs .......................................................................10
Figure 3-1: Comparison of the exact solution and the approximate solution to obtain the generalized
Vermeulen transfer function parameters for isotropic spherical matrix blocks (three sets of fractures)
..............................................................................................................................................26
Figure 3-2: Comparison of the respective exact solutions and the approximate solution to obtain the
generalized Vermeulen transfer function parameters for isotropic planar and cylindrical matrix blocks
..............................................................................................................................................27
Figure 3-3: Comparison of the respective exact solutions and the approximate solution to obtain the
generalized Vermeulen transfer function parameters for anisotropic matrix blocks with one, two, and
three sets of fractures ...............................................................................................................28
Figure 3-4: Fluid (dry gas) properties used in Example 3.2.1: Formation volume factor (𝐵𝑔), viscosity,
and density ..............................................................................................................................30
Figure 3-5: Average matrix pressure for a 1 D matrix block: SPM calculations with 3 simulators (left
panel) and DPM calculations using the traditional and new transfer functions (right panel) .................31
Figure 3-6: Average matrix pressure for a 5 nD matrix block: SPM calculations with 3 simulators (left
panel) and DPM calculations using the traditional and new transfer functions (right panel) .................31
Figure 3-7: Average matrix pressure in the DP sector model for a 1 D rock (left panel) and 5 nD rock
(right panel).............................................................................................................................33
Figure 3-8: Schematics of an anisotropic permeability matrix block represented using two
approaches: Single-block DPM (left panel) and the equivalent fine-grid SPM (right panel) ...................33
Figure 3-9: Average matrix pressure for an anisotropic rectangular (3D) matrix block..........................34

viii
Figure 3-10: Variations in diffusivity with pressure for two different fluids (left panel) and average
matrix pressure have we repeated Example 3.1.1 using pure CH4 production (right panel)...................35
Figure 3-11: Schematics of an irregular isotropic matrix block using two different approaches: Singleblock DPM (left panel) and fine-grid SPM (right panel)....................................................................36
Figure 3-12: Average matrix pressure of the irregular matrix from SPM and two DPMs (with tuned
parameters): New transfer function (left panel) and traditional transfer function (right panel)..............37
Figure 4-1: Domain representing a 1D reservoir with a known fracture spacing (left panel); results of
the classical DPM and the coarse EDFM (middle panel); and results of the DPM (gVer) and the refined
EDFM (right panel)....................................................................................................................41
Figure 4-2: Top view of the domain which represents a 2D reservoir with a fracture network, where
all fractures are orthogonal (left panel). The middle and right panels report the average pressure of
the respective models ...............................................................................................................42
Figure 4-3: Top view of the domain that includes two sets of parallel (orthogonal) fractures (left
panel). Simulation calculations of SPM, coarse EDFM, and two DPMs (both performed using predetermined shape factors assuming cylindrical geometry) ..............................................................43
Figure 4-4: Results of the upscaling technique when applied to Example 4.1.3. Both DPMs use transfer
function parameters obtained from the reference SPM solution ......................................................45
Figure 4-5: Two intersecting impermeable fractures. Only the west boundary permits pressure
depletion, all others are no-flow .................................................................................................46
Figure 4-6: Pressure map at the end of simulation of the SPM and EDFM, respectively........................47
Figure 4-7: Average pressure for all models (left panel): Reference solution, refined EDFM, and 2
upscaled DPMs (traditional form and gVer). Error in average pressure relative to that of the reference
solution (right panel).................................................................................................................47
Figure 4-8: A synthetic fractured reservoir with 300 natural fractures in a 500 nD matrix modeled
using an EDFM (left panel: 3D view and right panel: top view).........................................................48
Figure 4-9: Average matrix pressure from the EDFM and the two upscaled DPMs (left panel). Insert,
in left panel, demonstrates early-time departure of the traditional transfer function by up to 600 psi
(see difference in the right panel)................................................................................................49
Figure 4-10: A) Schematics for setting #1 where fractures are parallel (non-intersecting); b) Settings
#2, 3, and 4 where EDFM contains different number of fractures: 200, 100, and 10, respectively; c)
EDFM with different average fracture lengths (ft): 16, 29, and 43, respectively; and d) EDFM of
different domain sizes (ft): 393 x 393 x 157 and 984 x 984 x 393 ......................................................51
Figure 4-11: Fracture conductivities that are assigned with two distinct values (left panel) and matrix
permeability describing a heterogeneous and tight reservoir (right panel).........................................53

ix
Figure 4-12: Upscaling process of a complex fractured reservoir. Left panel: EDFM with two subregions. Right panel: Equivalent upscaled dual-block DPM ..............................................................54
Figure 4-13: Boundary conditions of the sub-divided EDFM.............................................................55
Figure 4-14: Average pressure on the primary y-axis, while the secondary y-axis is the difference
between each DPM and EDFM (right panel). The right panel reports the transfer function correction
factor......................................................................................................................................55
Figure 4-15: Average matrix pressure of the two EDFM sub-regions of Example 4.3.3..........................56
Figure 4-16: Potential upscaling of the complex EDFM of Example 4.3.3............................................57
Figure 5-1: Schematic of the SPM where the red boundaries represent fractures and yellow blocks
represent the matrix (left panel). Right panel reports the comparison between the pressure response
of the analytical solution (Crank, 1975), the fine-grid SPM, and the DPM (gVer). ................................62
Figure 5-2: Average matrix pressure for 1 D permeability (left panel) and 5 nD permeability (right
panel) .....................................................................................................................................64
Figure 5-3: Phase envelope for the quaternary fluid (left panel), gas/liquid relative permeability curves
(middle panel), and CCE liquid dropout (right panel) ......................................................................66
Figure 5-4: Average matrix pressure (left panel), cumulative gas production (middle), and cumulative
liquid production (right). Note: All production streams are converted to surface condition ..................66
Figure 5-5: Results of the fine-grid SPM: Pressure profiles (left panel), liquid saturation profiles
(middle panel), and liquid saturation distribution after 2 years of production (right panel)...................67
Figure 5-6: Flowchart comparing the two ways to evaluate the transfer function. Note: “phase
equilibrium” represents all calculations that are required to evaluate the transfer function parameters
..............................................................................................................................................70
Figure 5-7: Plots in order (top row left to right): Gas production rate, cumulative gas production,
average pressure, liquid production rate, cumulative liquid production, and liquid saturation at the
matrix/fracture interface............................................................................................................71
Figure 5-8: Phase envelopes generated from the final average matrix composition of the respective
models....................................................................................................................................72
Figure 5-9: Results are in the order of: Gas production rate, cumulative gas production, average matrix
pressure, liquid production rate, and cumulative liquid production ..................................................73
Figure 5-10: Comparison of the component cumulative production for the first 12 lighter components
for all four models: reference SPM (red solid line), DPM using unmodified LA (green dashed line),
DPM using unmodified gVer (blue circles), and DPM using modified gVer (purple stars) ......................74

x
Figure 5-11: Comparison of the component cumulative production for the last 12 heavier components
for all four models: reference SPM (red solid line), DPM using unmodified LA (green dashed line),
DPM using unmodified gVer (blue circles), and DPM using modified gVer (purple stars) ......................74
Figure 5-12: Average matrix pressure, transfer pressure obtained by matching the DPM to a reference
solution, and transfer pressure obtained by an analytical solution of Examples 5.3.1 (left panel) and
5.3.2 (right panel).....................................................................................................................76
Figure 5-13: Updated results of Example 5.3.1 with the adaptive transfer pressure that is generated
using the analytical solution .......................................................................................................77
Figure 5-14: Updated results of Example 5.3.2 with the adaptive transfer pressure that is generated
using the analytical solution .......................................................................................................78
Figure 6-1: Eagle Ford shale core used in the HnP experiment .........................................................81
Figure 6-2: Pressure-pulse decay and gas flow-through experimental setup.......................................84
Figure 6-3: Upstream and downstream pressure data from the 2nd pulse (left) and its interpreted ln(pD)
vs. time plot (right) ...................................................................................................................85
Figure 6-4: Calculated permeability from the pressure-pulse decay measurements.............................86
Figure 6-5: Permeability estimated from gas flow-through measurements: black symbols represent
the upstream pressure, blue symbols depict the downstream pressure, and red symbols show the
estimated permeability..............................................................................................................87
Figure 6-6: Calculated permeability from the gas flow-through measurements...................................87
Figure 6-7: Pressure and temperature readings from the HnP experiment: Pressure is shown by blue
symbols while green symbols represent the temperature ...............................................................88
Figure 6-8: Experimental oil recovery at the end of each of the four HnP cycles..................................89
Figure 6-9: Measured produced oil composition for the four synthetic oil components .......................90
Figure 6-10: Observed pressure drop during the soaking period for the four HnP cycles ......................91
Figure 6-11: 1D domain showing the initial composition of the binary system (left panel) and CO2 mole
fraction (after 0.2 day) along the domain for the analytical and numerical solutions (right panel) .........96
Figure 6-12: Schematic of Arnold and Toor (1967) experimental setup (left panel). Average CH4 and Ar
compositions versus time in the middle and right panels, respectively. The solutions (reported in both
panels) are in the order of lab data (circles), MRST (solid lines), Hoteit (2013) in dashed lines, and
Taylor and Krishna (1993) in doted lines. ......................................................................................97

xi
Figure 6-13: Numerical grid that mimics the EF-3 HnP experiment (in MRST): 3D view of the extraction
cell and the core (left panel), 3D view of the core (middle panel), and horizontal cross-section from
the top of the core (right panel)..................................................................................................98
Figure 6-14: Comparison of oil recovery from HnP experimental data and simulation results: Lab data
(blue circles), simulation with no diffusion (solid black line), simulation using classical Fick’s law
(dotted yellow line), and simulation using generalized Fick’s law (dashed purple line) .......................100
Figure 6-15: Comparison between the interpolated infinite dilution coefficients for CO2 in oil and the
calculated values from each correlation (left panel). Interpolated values as a function of carbon
number (right panel)...............................................................................................................102
Figure 6-16: Experimental infinite dilution diffusion coefficients for oil in CO2 and the calculated values
from each correlation..............................................................................................................103
Figure 6-17: Comparison of oil recovery after applying the infinite dilution coefficients multipliers: Lab
data (blue circles), simulation with no diffusion (black solid line), simulation unmodified diffusion
model (dotted yellow line), and simulation when applying the infinite dilution multiplier (dashed
purple line) ............................................................................................................................105
Figure 6-18: Comparison of oil recovery from experiment and simulation: The classical Fick’s law (left
panel) and the generalized Fick’s law (right panel). Blue circles represent lab data, dotted-dashed
purple line represents ‘Diffusion Case – 1,’ and dashed red line represents ‘Diffusion Case – 2’ ..........106
Figure 6-19: Comparison of produced oil composition from experiment and simulation. The first
column reports simulation results where diffusion is neglected, columns 2-3 report calculations with
the classical Fick’s law, while columns 4-5 provide results from the generalized Fick’s law ..................107
Figure 6-20: Vertical cross-section of the core reporting the CO2 mole fraction during soaking (first
cycle). Note: The snapshots shown here are based on calculations with the classical Fick’s law –
Diffusion Case 1 ......................................................................................................................108
Figure 6-21: Comparison of oil recovery calculated with different relative permeability functions:
quadratic (blue columns) and straight-lines (orange columns). Note: The results reported here are
based on calculations with the classical Fick’s law – ‘Diffusion Case 1’.............................................109
Figure 6-22: Comparison of produced oil composition with different relative permeability functions:
quadratic (blue columns) and straight-lines (orange columns). Note: The results reported here are
based on calculations with the classical Fick’s law – ‘Diffusion Case 1’.............................................110
Figure 6-23: Pressure variation during soaking period of the first cycle (experimental vs simulation) ...111

xii
Abstract
The current motivation and efforts of the oil and gas industry are hinged on the fact that any future
increase in hydrocarbon production will require further development of unconventional resources. Tight
sand and shale hold vast hydrocarbon reserves that could support increasing energy demands; however,
the current low recovery factors remain a challenge. Therefore, technical studies have been ongoing to
improve existing practices and unlock the potential of these reservoirs. To this end, gas injection into tight
fractured reservoirs has been a successful enhanced oil recovery (EOR) method, making it a promising
investment and an active research area.
To support ongoing field activities and incentivize new developments, we rely on available
technology for simulation and optimization of production operations. It is important that these tools are
based on models that accurately reflect the relevant physical processes. As of today, many details still
remain to be resolved to delineate the dynamics of the tight fractured systems and their implementation
into reservoir simulation tools.
In this research, we explore fit-for-purpose models to describe multicomponent matrix/fracture
mass transferin unconventional formations. The objective is to provide a modeling platform, via an opensource code, where our findings are tested and validated. Our development includes effective algorithms
to facilitate further advancements in research and implementation of simulation at larger scale.
We start by formulating a new transfer function to represent the matrix/fracture mass transfer for
tight fractured systems via the dual-porosity modeling (DPM) approach. The formulation is inspired by the
work of Vermeulen (1953), who approximated the analytical solution for diffusion in a spherical geometry
using a closed exponential form to describe the average change of concentration inside the domain of
interest. By generalizing Vermeulen’s approximation, we introduce two parameters to represent the
geometry of the system: Ashape factor 𝜎and an exponent 𝑉. The final form, referred to as the generalized
Vermeulen (gVer) transfer function, includes an implicit time-dependent correction factor that represents

xiii
the sub-scale physics of the problem that allows for accurate prediction of the prolonged transient
behavior in tight fractured systems.
The generalized Vermeulen transfer function was implemented in the MATLAB Reservoir
Simulation Toolbox (MRST) to allow for a comparison with commercial simulators (CMG and ECLIPSE). An
initial validation of the new transfer function was performed via single-phase pressure depletion
calculations, where a fine-grid single-porosity model (SPM) is compared with an equivalent single-block
DPM for a tight fractured system with a 5 nD matrix. First, we observe that the SPM reference solution,
from the three simulation tools, agree with one another. Then, by comparing the results from DPM
calculations with MRST and commercial tools, we observe that the application of the Warren and Root
(1963) transfer function with a shape factor from Lim and Aziz (1995) provides for a reasonable agreement.
However, these results, obtained via the traditional transfer function, depart from the reference (SPM)
solution at early time but reach agreement at later time. In contrast, the gVer transfer function is
demonstrated to provide for an excellent agreement with the reference solution throughout the entire
depletion process.
Several exampleswere subsequently defined for a broader range of problems. The objectives were
to document the limitations of the traditional DPM,and to further investigate the performance of the new
gVer transfer function formulation. These calculation examples include: 1) matrix blocks with isotropic
permeability; 2) matrix blocks with anisotropic permeability; and 3) sensitivity to variable hydraulic
diffusivity (highly compressible fluids). In all the cases, the DPM using the generalized Vermeulen transfer
function is demonstrated to agree with the response of the reference solution for matrix permeability
values as low as 5 nD.
The initial testing of the gVer transfer function was performed on well-defined matrix geometries,
i.e., for problems that can be described by readily available analytical solutions. We then proceed to
demonstrate that the new transfer function is not limited to simple geometries but can also provide an

xiv
accurate description of flow behaviors for complex fracture geometries. In such settings, the relevant
matrix/fracture interaction parameters, 𝜎and 𝑉, are obtained from a reference solution, that, for example,
can be generated by calculations with an SPM or an embedded discrete fracture model (EDFM). We apply
this approach to address problems with arbitrary fracture orientations and geometries in heterogeneous
nano-darcy domains. We furthermore demonstrate the versatility of our approach by refining the DPM in
regions of higher spatial variation to capture the details of the fractured system. A workflow for converting
the EDFM to an equivalent DPM is subsequently discussed in the context of upscaling of discrete fracture
characterizations.
We then developed a DP compositional model in MRST to facilitate studies of multicomponent
transport in tight fractured systems. This model was first validated against an analytical solution for a
single-phase pressure depletion process. Following that, we extendedour study to include compositional
two-phase problems to model tight fractured gas-condensate reservoirs. However, we find that for such
problems, the implicit and coarse representation of the matrix/fracture interaction, described in a similar
manner to that of Kazemi et al. (1976), resultsin an overprediction of the liquid production: The averaging
of matrix block properties masks the compositional changes near the matrix/fracture interface, that is
critical to predict production behaviors below the dew point. To overcome the limitations of the current
models, we propose a novel approach that evaluates the gVer transfer function at conditions that
reflect/approximate the pressure and composition/saturation gradients within a matrix block. By
approximating the fluid state near the interface, we approximate the liquid dropout and buildup that
dictate the fluid mobility more accurately. We present calculation results for a single-block DPM,
representative of a tight fractured gas-condensate reservoir considering two condensate fluid examples:
1) A 4-component analog rich gas; and 2) A realistic 24-component fluid from an existing gas-condensate
reservoir. We demonstrate that calculations with the new modified gVer transfer function substantially
agree with the fine-grid SPM results in terms of the amounts and compositions of the produced phases.

xv
We subsequently directed our attention to study the role of molecular diffusion on recovery for
tight fractured reservoirs. To do that, we conducted an integrated lab experimental and numerical
modeling effortsof carbon dioxide huff-n-puff (HnP) on an Eagle Ford shale core (EF-3) with a porosity of
11.9% and permeability of 2.5 µD. The lab experimental work was performed by Saeed Alahmari and is
briefly summarized here: Four cycles of HnP were performed where the EF-3 core was initially saturated
with a synthetic oil composed of n-decane, n-dodecane, n-tetradecane, and n-hexadecane. These
experiments were carried out at pressures and temperature that mimic reservoir conditions (from 3,900
psi to 4,400 psi and 158 F/343.15 K) under first-contact miscibility conditions. Based on the pressure drop
during soaking as well as the compositional analysis of the produced oil, molecular diffusion is observed
to be the dominant recovery mechanism in the performed HnP experiments. To interpret the experimental
observations, a compositional simulation study, via the SPM approach, was performed. The governing
equations of MRST were modified to represent molecular diffusion described by either the classical Fick’s
law or the generalized Fick’s law, with composition-dependent diffusion coefficients. For both diffusion
models, our simulation calculations show very good agreement with the lab data (for all cycles) in terms
of the recovery and produced oil compositions, with the former model demonstrating higher accuracy and
computational efficiency. The results of our modeling and interpretation suggest that the application of a
simple diffusion model (classical Fick’s law) may be sufficient and certainly more computationally efficient
than more sophisticated diffusion models (i.e., the generalized Fick’s law).
In summary, this dissertation reports on our efforts to combine physical models and efficient
simulation techniques to advance the understanding of multicomponent mass transfer in tight fractured
rocks. The overarching goal was to explore and develop computationally simpler, yet accurate, approaches
to enable/support more practical and reliable forecasting for unconventional reservoirs at larger scale.

Chapter 1. Introduction
1.1. Background and Motivation
Over the past several decades, the global demand for energy has increased to support the growing
population, the increase in manufacturing, and the higher living standards. However, with the decline of
production from conventional reservoirs, it is imperative that we explore and tap into more challenging
hydrocarbon resources to circumvent the shortage of the hydrocarbon supply.
Unconventional reservoirs are an attractive target for development due to their large hydrocarbon
potential. Considerable exploration and production efforts in the U.S. began in 2007: Oil production from
the U.S. shale plays increased nearly 20-fold from 2007 to 2022 (EIA 2023) as reported in Fig. 1-1 (left
panel). The right panel of Fig. 1-1 reports the evolution of dry gas production from unconventional
reservoirs in the U.S. (EIA – AEO2022). As can be seen, the bulk of the projected future production will be
derived from fractured systems of modest permeability.
Figure 1-1: U.S. tight oil production from selected plays (EIA 2023) and U.S. dry natural gas production (EIA – AEO2022)

2
It is, however, questionable how we can deliver future production given the modest oil recovery factors
observed today. For example, recovery figures, via primary production, of less than 10% have been
reported from the Bakken field (Hawthorne et al., 2013; Hoffman, 2018). Therefore, to improve
conventional production practices and incentivize further investments, horizontal drilling and multistage
hydraulic fracturing have been successfully implemented to enhance well productivity. EOR techniques
have also gained significant attention, where gas injection into tight fractured reservoirs, for example, has
been reported to yield promising incremental recovery. At the lab scale, experiments related to carbon
dioxide (CO2) injection performed on cores, extracted from unconventional formations, have shown
significant increasesin attainable incremental recovery of more than 40% (Darvish et al., 2006; Kovscek et
al., 2008; Vega et al., 2010; Alharthy et al., 2015; Ghasemi et al., 2018; Alahmari et al., 2023).
Using CO2 as an injectant has several advantages. Environmentally, CO2 injection, and possibly
retaining a portion of it underground, can contribute to mitigating climate change by reducing carbon
emissions. In terms of EOR design and recovery, CO2’s ability to develop miscibility with hydrocarbon
mixtures at significantly lower pressures, compared to natural gas or nitrogen, helps target reservoirs
beyond primary production. Furthermore, CO2 injection results in oil swelling and hence in a reduction of
the oil viscosity, enhancing oil mobility and allowing for the mobilization of more oil from the tight matrix.
Therefore, given the potential profitability and feasibility of such methods, we expect an increasing global
oil production due to CO2 injection projects over the next decades (Fig. 1-2).

3
Figure 1-2: Global oil production due to EOR activities with projections from 2018 to 2040 (IEA 2018)
1.2. The Need for Advanced Techniques for Modeling of Tight Fractured Reservoirs
Embarking on the extraction of unconventional resources presents a fundamental shift in the time scales
of the relevant recovery mechanisms at play. To conceptualize this, we simulate dry gas production from
a fixed volume of rock and consider two permeability settings: 1 mD and 0.8 µD. The former represents a
reservoir with higher permeability, as expected in conventional developments. The latter represents the
permeability of a core extracted from an unconventional reservoir, which, as will be shown later, becomes
essential in our modeling of HnP experiments. The average pressure in the rock volume for the two
permeability scenariosis reported in Fig. 1-3.
Figure 1-3: Average matrix pressure response of 1 mD and 0.8 D rocks

4
From Fig. 1-3, it is evident that the time required to reach a stabilized pressure reading is substantially
different between the two rocks (1 day versus 2 years). This difference is expected, given that the
characteristic time (𝑡
∗
) for the pressure decay is inversely proportional to the rock permeability (𝑡
∗ =
𝐿
2𝜇𝜙𝑐⁄𝑘). The key question that we explore in this work is: Can the current DPMs, implemented in
commercial simulators, accurately forecast the production from tight fractured rocks? Most of the
commercialtools were developed to support production from rocks in the 1 mD range; however,they have
not been adaptedto unconventional developments that entailmuch tighter rocks with substantially lower
permeabilities.
Therefore, to explore the performance of commercial simulators, we continue on the above
example and compare the pressure response of the 0.8 D rock predicted from several DPMsto a finegrid SPM reference solution, with special attention to the early characteristic time. From Fig. 1-4, it can be
seen how the selection of the transfer function (and associated shape factor) dictates the accuracy of the
DPM, signifying the importance of proper representation of matrix/fracture mass transfer, especially in
tight fractured reservoirs.
Figure 1-4: Average matrix pressure response during the early characteristic time of the 0.8 D rock

5
In an unconventional dry gas reservoir, as considered in Fig. 1-4 for example, large pressure gradient can
exist near the matrix/fracture interface, with the transient time extending well beyond what is observed
in conventional reservoirs. Such distinct behavior makes it imminent to look for more advanced techniques
to provide reliable predictions of tight fractured media.
1.3. Problem Statement
Many unknowns still exist regarding the fluid flow and transport in unconventional reservoirs. The lack of
sufficient data and reliable simulation tools poses a difficulty in advancing our understanding of the physics
at play and to develop models that can support the development of these resources.
Currently, several options/approaches exist to model the matrix/fracture mass transfer in
fractured reservoirs. These options are (see Chapter 2 for further details):
1. Single-Porosity Models: Traditionally, and quite commonly, an SPM is used where the matrix and
fractures are represented explicitly and therefore no additional terms are required to model the mass
transfer between the two domains. Instead, local grid refinement and/or unstructured grids may be
used to capture the transient behavior in tight rocks. However, this can result in a large number of
grid cells making this approach very CPU intensive.
2. Dual-Porosity Models: Warren and Root (1963) revolutionized the interpretation of a fractured
reservoir where fractures are modeled implicitly, allowing for a coarserrepresentation of a fractured
domain. They suggested the use of a transfer function to represent the mass transfer between the
two continua. Despite its wide use, the accuracy of the DPM is dictated by the representation of the
transfer function, and the associated shape factor.
3. Embedded Discrete Fracture Models: In recent years, the EDFM has been used to model fractured
systems where discrete fractures are embedded in the grid as independent objects from the matrix
rock. To allow for communication between the two domains, a linear driving-force expression

6
(transfer function) is employed. Such mass transfer takes on a similar form as used in the DPM
approach requiring improved representation of the transient behavior (as demonstrated above). In
addition, the large number of control volumes used in the EDFM rendersthis option computationally
expensive, especially when simulating larger-scale models.
As indicated above, to be able to exploit options 2 or 3, the selection of the matrix/fracture
transfer function becomes of paramount importance. The transfer function serves as a fundamental
formulation that couples the matrix/fracture and allows for modeling of the exchange between the two
domains (more details provided later). The transfer function was first developed by Warren and Root
(1963) to handle mass transfer during a simple pressure depletion process. Since then, researchers have
developed variations of this transfer function to solve other problems, for example predicting the
performance of fractured reservoirs undergoing water injection, and more recently to describe diffusive
mass transfer problems. As we will demonstrate, without the proper representation of slow transientsin
unconventional reservoirs, the formulation of Warren and Root (1963), and the related variants,may result
in significant errors.
1.4. Dissertation Objectives and Deliverables
The objective of the research presented in this disseration is to explore fit-for-purpose models to describe
the multicomponent mass transfer in tight fractured reservoirs, to support interpretation, design, and
optimization of EOR operations in such settings. The research aims to develop and test effective algorithms
in an open-source reservoir simulation environment. A central component of the efforts is to investigate
the limitations of the existing commercial simulation tools, as well as to bridge the gap that exists between
simulation results and field/experimental data. Ultimately, this will support the industry and commercial
simulator development by providing practical and reliable approaches for forecasting of hydrocarbon
recovery from unconventional resources.

7
To facilitate technology transfer, the development, testing, and validation of this research was
primarily performed using MRST, an open-source reservoir simulation code developed at SINTEF Digital,
with input from several academic institutions. The reader is referred to Lie (2019) and Lie and Moyner
(2021) for more information related to the base formulations of this platform.
1.5. Dissertation Organization
The remainder of this dissertation is organized as follows:
Chapter 2: Provides an overview of the existing efforts to model mass transfer in tight DP reservoirs.
Specific topics of the discussion includes modeling of the transient state in single-phase applications as
well asrepresentation of flow and transport behavior in multiphase problems. Lastly, we provide a review
of past studies related to diffusive mass transfer in the context of EOR and related implementations in
reservoir simulation tools.
Chapter 3: Introduces the new generalized Vermeulen transfer function for DPM calculations of mass
transfer between matrix and fracturesfor tight DP systems. It was implemented in MRST and validated
against commercial simulators to serve as the reference formulation for a variety of DP problems
presented in this dissertation.
Chapter 4: Describes a workflow to allow for modeling of complex tight fractured systems. To demonstrate
its accuracy, the DPM, via the generalized Vermeulen transfer function, is shown to match the flow
behavior of a reference numerical solution (including fine-grid SPM and EDFM) describing depletion
problems from arbitrary matrix/fracture geometries.
Chapter 5: Details the development of a DP compositional model (in MRST), which is initially validated
against an analytical solution, and the modification to the generalized Vermeulen transfer function to
model compositional two-phase (gas condensate) problems in tight fractured formations. Two fluid

8
examples are presented (a 4-component analog rich gas and a 24-component realistic fluid from a wet gas
reservoir) to demonstrate the performance of the modified transfer function.
Chapter 6: Reports our experimental and modeling efforts to study diffusive mass transfer and recovery
during CO2 HnP in an Eagle Ford shale core. Experimental observations include HnP, pressure pulse decay,
and gas flow-through measurements. MRST compositional simulations are presented and analyzed to
interpret the experimental data.
Chapter 7: Provides an overview of the key findings and contributions of the research. Additionally, it offers
recommendations for future studies to further narrow the existing knowledge gaps.

9
Chapter 2. Literature Review
2.1. Modeling of Fractured Reservoirs
The mass transfer between fracture and matrix subsystems dictatesthe recovery from unconventional
reservoirs. Accordingly, an efficient representation of the two domains as well as proper modeling of the
matrix/fracture interactions are critical to provide reliable and accurate production forecasting. Below we
provide a review of three main approaches used for modeling of fractured systems.
2.1.1. Single-Porosity Models
SPMs use an explicit computational representation of fractures where each fracture segment is commonly
defined by a set of control volumes. This approach is desirable to use because no additional terms are
needed to describe the mass transfer between the fracture and matrix domains. The approach, however,
may require local grid refinement (LGR)to capture the transients in the matrix, and/or unstructured grids
to honor the orientation and geometry of the fracture segments. Despite its accuracy, the computational
cost of the SPM approach is prohibitively high due to the large number of grid blocks that are required,
with the additional complication that variations in geological parameters between adjacent cells can result
in numerical instabilities (Hoteit and Firoozabadi, 2006; Sarma and Aziz, 2006).
2.1.2. Dual-Porosity Models
Barenblatt et al. (1960) introduced the general foundation for modeling of flow behaviors in two-porosity
systems. This idea was later extended to petroleum engineering applications by Warren and Root (1963)
via a geometrical interpretation of a fractured reservoir. The DPM characterizes the fractures as the main
flow conduits, while the matrix acts as the main storage capacity (see Fig. 2-1).

10
Figure 2-1: Warren and Root (1963) idea of the dual-continuum model with constant shape factors that reflect the geometry of
the fractured system. In later years, time-dependent transfer functions have been introduced for tight fractured reservoirs
A central contribution of Warren and Root (1963) is the coupling of the two domains specified by a transfer
function. By assuming that single-phase fluid flow is described by Darcy’s law, via the pressure in each
domain, flow of a slightly compressible fluid through a homogeneous fractured medium, is governed by
two coupled diffusion equations:
Fracture: 𝜙𝑓𝑐𝑓
𝜕𝑃𝑓
(𝑥𝑓,𝑡)
𝜕𝑡 =
𝑘𝑓
𝜇
∇
2𝑃𝑓
(𝑥𝑓,𝑡) − 𝛤(𝑥𝑓,𝑡) , (2.1)
Matrix: 𝜙𝑚𝑐𝑚
𝜕𝑃𝑚(𝑥𝑚,𝑡)
𝜕𝑡 =
𝑘𝑚
𝜇
∇
2𝑃𝑚(𝑥𝑚,𝑡) . (2.2)
The subscripts 𝑓 and 𝑚 denote fracture and matrix domains, respectively. The parameter 𝜙 is the porosity,
𝑐 is the fluid compressibility, 𝑃 is the pressure, 𝑥 is location, 𝑡 is time, 𝑘 is the permeability, and 𝜇 is the
fluid viscosity. Eqs. (2.1) and (2.2) are coupled through 𝛤, that definesthe matrix/fracture interflow term,
i.e., the transfer function that acts as a source/sink term. We note that 𝛤 does not appear in Eq. (2.2) as it
serves as the boundary condition.
Warren and Root (1963) assumed that the flow inside the matrix occurs under pseudo steadystate (PSS) condition. As a result, the matrix pressure is no longer defined at each point within the matrix

11
domain, forcing the spatial coordinate, 𝑥𝑚,to be dropped from Eq. (2.2). This assumption introduces the
idea of an average matrix pressure (denoted as 𝑃̅𝑚) that is now governed by an ordinary differential
equation:
𝜙𝑚𝑐𝑚
𝑑𝑃̅𝑚(𝑡)
𝑑𝑡
= 𝛤(𝑡) . (2.3)
Finally, to represent the exchange term, Warren and Root considered a linear driving force model
according to:
𝛤 = 𝜎
𝑘
𝜇
(𝑃𝑓 − 𝑃̅𝑚) . (2.4)
Eq. (2.4) models the mass transfer rate by a shape factor, 𝜎, rock/fluid properties, and a PSS pressure drop.
The shape factor represents the geometry of the system and defines the flow between the two domains
via the simple formulation (Warren and Root, 1963):
𝜎 =
4𝑁(𝑁+2)
𝐿
2
, (2.5)
where 𝑁 is the number of sets of parallel fractures (1, 2, or 3), and 𝐿 is the fracture spacing.
A more general definition was introduced by Kazemi et al. (1992), that provides a physical meaning
of the shape factor:
𝜎 =
1
𝑉𝑚
∑
𝐴
𝐿
. (2.6)
Eq. (2.6) defines the shape factor in terms of the interface area open to flow, 𝐴, per the bulk volume of
the matrix block, 𝑉𝑚, and the characteristic length of the matrix/fracture system, 𝐿.
There have been several contributions to improve the approximation of the shape factor. Kazemi
et al. (1976) was the first to use shape factors in numerical simulation to predict the mass transfer between
the matrix and fracture domains. Their derivation also assumed that the matrix/fracture transfer takes
place under PSS conditions. For 1D problems, the shape factor is:
𝜎 =
4
𝐿
2
. (2.7)

12
Dykhuizen (1990) derived a solution to the transient inside the matrix block and suggested using two
different equations to describe the transient and PSS flow behaviors. The switch between the equations
depends on a critical time (when the diffusion front reaches the center of the block).
Chang (1993) suggested that the underlying matrix dynamics can be approximately captured by
solving the pressure diffusion equation. Using the analytical solution for a 1D problem, the resulting
expression led to the shape factor that is a function of time:
𝜎 =
𝜋
2
𝐿
2
∑ exp[−(2𝑛+1)
2𝜋
2 ∞ 𝑡𝐷] 𝑛=0
∑
1
(2𝑛+1)
2 exp[−(2𝑛+1)2𝜋2 ∞ 𝑡𝐷] 𝑛=0
, (2.8)
where the dimensionless time, 𝑡𝐷, is defined as:
𝑡𝐷 =
𝑘𝑡
𝐿
2𝜙𝜇𝑐 = 𝜎𝐷𝑡 , (2.9)
with the hydraulic diffusivity, 𝐷, given by:
𝐷 =
𝑘
𝜙𝜇𝑐 , (2.10)
and the product 𝜎𝐷 representing the inverse of the characteristic time for diffusion.
The shape factor obtained by Eq. (2.8) assumes larger values at early time and converges at later
time, allowing for a more flexible representation of the transient behavior. However, this approach is
limited to a fixed geometryand is expected to increase the computational requirement for evaluating the
series expansion for each control volume.
A similar derivation was performed independently by Lim and Aziz (1995) who proposed new
shape factors for the three simple geometries (planar, cylindrical, and spherical). However, they decided
to truncate the series solution and approximate the pressure response using only the first term. For
example, their approximation for a spherical matrix block can be written as:
𝑃𝐷 =
𝑃̅𝑚−𝑃0
𝑃𝑓−𝑃0
= 1 − 0.61exp(
−𝜋
2𝐷𝑡
𝑅2
) , (2.11)
where 𝑃𝐷 is the dimensionless pressure, 𝑃0
is the initial pressure, and 𝑅 is the radius.

13
Their approximation (seen in Eq. (2.11)) departed from the exact solution at early dimensionless
time. If we define the characteristic time, 𝑡
∗
, as 𝑡
∗ = 𝑡/𝑡𝐷, Lim and Aziz suggested that a departure from
the analytical solution is insignificant for reservoirs with a small value of 𝑡
∗
. However, as the characteristic
time depends on rock and fluid properties (𝑡
∗~𝐿
2/𝐷), substantially larger values of 𝑡
∗
should be expected
for tight reservoirs compared to conventional reservoirs, translating to a large real time departure of the
pressure response. Therefore, while Eq. (2.11) provides an approximation of the pressure transient in the
matrix, the accuracy of the truncated series expansion, at early time, depends strongly on the magnitude
of 𝐷/𝑅
2
(or the inverse of the characteristic time for pressure diffusion). For small values of 𝐷 (relevant
to shales), Zimmerman et al. (1993) demonstrated that Eq. (2.11) provides a poor approximation of the
average pressure in the matrix. Nonetheless, due to its simplicity, the use of the shape factors from Lim
and Aziz (1995) in the traditional transfer function (see Eq. (2.4)) is still popular in commercial reservoir
simulators.
Zimmerman et al. (1993) worked on resolving the transient in the matrix block in settings where
Eq. (2.4) is inaccurate at earlytime (large characteristic time for pressure diffusion). By adapting the work
of Vermeulen (1953), they suggested to use of a closed form approximation (for spherical geometry):
𝑃̅𝑚−𝑃0
𝑃𝑓−𝑃0
= [1 − exp(−
𝜋
2𝐷𝑡
𝑅2
)]
1/2
. (2.12)
To evaluate the mean pressure response in the matrix block, they differentiated Eq. (2.12) with respect to
time and eliminated time resulting in:
𝑑𝑃̅𝑚
𝑑𝑡
=
𝐷𝜋
2
2𝑅2
(𝑃𝑓−𝑃0
)
2
−(𝑃̅𝑚−𝑃0
)
2
(𝑃̅𝑚−𝑃0
)
. (2.13)
Re-arranging Eq. (2.13) results in the following form:
𝑑𝑃̅𝑚
𝑑𝑡
=
𝐷𝜋
2
2𝑅2
(𝑃𝑓 − 𝑃̅𝑚)
|𝑃𝑓−𝑃0
|+|𝑃̅𝑚−𝑃0
|
|𝑃̅𝑚−𝑃0
|
= 𝜎𝐷(𝑃𝑓 − 𝑃̅𝑚)𝜒 . (2.14)
Eq. (2.14) introduces the idea of a boost (or correction) factor, 𝜒, that accounts for the departure from
PSS in the matrix: The boost factor assumes a relatively large value at early time and approaches unity at

14
later time when the PSS is approached, and the transfer function revertsto the form of Warren and Root
(see Eq. (2.4)). Zimmerman et al. demonstrated that Eq. (2.14) performs well at both early and later times.
Ranjbar and Hassanzadeh (2011) introduced a 1D shape factor specifically for single-phase
compressible fluids. The solution of the nonlinear gas diffusivity equation is used to derive two separate
expressions to describe the early- and late-time shape factors. Similar to Dykhuizen (1990), a
dimensionless time is used to switch between the two equations.
2.1.3. Discrete Fracture Modeling
Discrete Fracture Models (DFM) were developed to explicitly account for the effect of individual fractures
on fluid flow. It relies on an unstructured mesh where the matrix discretization must conform to the
fractures, which are represented with a lower dimensionality than that of the model (𝑛-1). Although such
requirement eliminates the need for a transfer function to describe the matrix/fracture exchange (Shakiba,
2014), the resulting model consists of a large number of control volumes and complicated gridding,
limiting its application for simulation of large 3D problems (Moinfar, 2013).
Several modeling approaches have evolved from the original DFM idea, aiming to improve the
computational efficiency. The most notable contribution is the EDFM (Lee et al., 2001; Li and Lee, 2008).
The goal of this approach is to avoid conforming the mesh of the fractures and matrix domains, and to use
a structured grid to represent the matrix. The EDFM includes a hierarchical method when fractures exist
at multiple length scales: First, small-scale fractures (fully contained in a single matrix cell) are
homogenized and integrated into the matrix through an effective permeability. The remaining fractures
(large-scale features spanning more than a single matrix block) are then embedded into the main matrix
grid (Moinfar, 2013).
Similar to the DPM approach, fracture and matrix domains are discretized independently in the
EDFM approach and are coupled by a matrix/fracture transfer term, that commonly takes the form:
𝛹 = 𝐶𝐼 𝑘
𝜇
(𝑃𝑓 − 𝑃̅𝑚) , (2.15)

15
where 𝐶𝐼 is the connectivity index that reflects the geometry of the system:
𝐶𝐼 = 𝐴⁄⟨𝑑⟩ . (2.16)
In the above equation, the parameter 𝑑represents the distance of a given fracture to the matrix boundary.
We observe that Eq. (2.15) is conceptually similar to the traditional transfer function of the DPM (see Eq.
(2.4)). However, given the nature of the EDFM, 𝑑 may represent more than a single fracture embedded
within a host matrix cell. During a pre-processing step, all matrix/fracture connections are identified and
⟨𝑑⟩ is calculated.
The resemblance of Eqs. (2.4) and (2.15) is clear: Both transfer functions use a linear driving force
to express the mass transfer between the two domains and use a constant coefficient (𝐶𝐼 or 𝜎) to
represent the geometry of the system. In the context of EDFM, numerous authors have worked to improve
the representation of matrix/fracture interactions. One approach is to refine the matrix mesh, and while
that has been proven to improve the calculation accuracy, it comes at the cost of an increased
computational complexity. Another approach was examined by Vo et al. (2020), who combined mesh
refinement with the application of a constant correction factor, 𝜎𝑓, to the matrix/fracture transfer
expression, to facilitate agreement between the EDFM calculation and the reference solution. They arrive
at a value of 𝜎𝑓 = 20.3/𝐿
2
, and highlight that this provides a similar value to that of Lim and Aziz (1995).
Jiang and Younis (2016) modeled fluid flow in unconventional reservoirs by integrating an EDFM
with a multiple-interacting-continua (MINC) model. The objective of integrating the MINC model was to
resolve the transient flow in areas with large pressure gradient, and to overcome the limitation of the
EDFM approach to model matrix/fracture mass transfer at early time. Therefore, this hybrid EDFM-MINC
model eliminates the need for LGR in the traditional EDFM.
Xu et al. (2019) demonstrated that the use of MINC is computationally expensive relative to a
classical DPM with an improved transfer function. Therefore, they suggested the use of a time-dependent

16
shape factor to capture the transient flow (similar to the idea of a correction/boost factor discussed
above):
𝜎transient = 𝜎𝑚
(𝑃𝑖−𝑃𝑓
)
2
−(𝑃𝑖−𝑃𝑚)
2
2(𝑃𝑖−𝑃̅𝑚)(𝑃̅𝑚−𝑃𝑓
)
, (2.17)
and applied the shape factor of Kazemi and Gilman (1993) to represent the PSS condition:
𝜎𝑚 = 𝜋
2 (
1
𝐿𝑥
2 +
1
𝐿𝑦
2 +
1
𝐿𝑧
2
) . (2.18)
It is worth noting, however, that a central difference exists between EDFM and DPM in the modeling of
mass transfer. In addition to the matrix/fracture communication (see Eqs. (2.4) and (2.15)), the EDFM also
accounts for matrix/matrix interactions similar to that of the traditional SPM approach. In the context of
the DPM, such mass transfer is not considered, but is incorporated in the flow equations of dualporosity/dual-permeability models (DPDPM). The selection between DPM or DPDPM to study fractured
systems depends on the nature of the matrix rock and whether its permeability is large enough to
contribute to fluid flow.
2.2. Multiphase Transfer Functions
The development of multiphase transfer functions has also received some attention in the existing
literature. The first contribution by Kazemi et al. (1976) extended the formula of Warren and Root (1963)
to multiphase flow settings by accounting for relative permeability and fluid phase properties:
𝛤𝑗 = 𝜎𝜌𝑗
𝑘𝑘𝑟𝑗
𝜇𝑗
(𝑃𝑓,𝑗 − 𝑃̅𝑚,𝑗
) . (2.19)
Here,the subscript 𝑗 denotes the phase. The parameters 𝜌𝑗
, 𝜇𝑗 and 𝑘𝑟𝑗 are the fluid density, viscosity, and
relative permeability of phase 𝑗, respectively. The introduction of Eq. (2.19) enabled the evaluation of
waterflood patterns in fractured reservoirs. However, the use of a direct generalization of the single-phase
equation of Warren and Root (1963) provides a somewhat coarse approximation of the complex
multiphase flow and may lead to inaccurate forecasting. Nonetheless, Eq. (2.19) remains the default model
in most commercial reservoir simulators.

17
Dean and Lo (1988) studied two-phase (gas/oil and oil/water) flow in fractured systems using a DP
representation and introduced pseudocapillary pressure and pseudorelative permeability functions. These
pseudofunctions were then adjusted to match the recovery calculation from the DPM to a fine-grid SPM
reference solution. By introducing the pseudofunctions, the authors effectively lump rock/fluid
interactions with unresolved physics (fluid properties and/or flow physics) and use the resultant group as
an adjustable parameter (pseudocurves).
Rossen and Shen (1989) used a similar approach to that of Dean and Lo (1988) to study gas/oil
drainage and oil/water imbibition processes. However, instead of adjusting the pseudofunctions to match
the recovery of the DPM to a fine-grid model, as was done by Dean and Lo (1988), they presented a
formula to calculate the pseudocapillary pressure as a function of time-dependent fluid saturations that
were obtained from an SPM reference solution.
Coats (1989) developed a compositional model with a general cubic equation of state (EOS) to
describe the phase equilibrium calculation. The exchange term in his DPM followed the form of Eq. (2.19)
and was also based on pseudofunctions.
Penuela et al. (2002) studied flow in fractured gas-condensate reservoirs and used
pseudofunctions in their representation. In their work, however, they introduced pseudopressure and
pseudotime to model the multiphase effects. This approach simplifies the complex multiphase problem
by accounting for the reduction in relative permeability and the non-monotonic behavior of the fluid
compressibility.
Sarma and Aziz (2006) performed single-block DPM calculations for a water imbibition experiment
and compared results with a fine-grid model. They developed a two-phase transfer function that includes
pressure (and saturation) diffusion, each expressed by a separate time-dependent shape factor. The final
form of the transfer function for the water phase is:
𝑞𝑤𝑚 = 𝑉𝜌𝑤 𝜆𝑤𝜎𝑃𝐷(𝑝̅𝑤 − 𝑝𝑤𝑓)− 𝑉𝜙𝜌𝑤𝜎𝑆𝐷(𝑆̅𝑤 − 𝑆𝑤0
) , (2.20)

18
where the subscript 𝑤 denotes the water phase. They demonstrated that the shape factor attributed to a
pressure change, 𝜎𝑃𝐷 , is equivalent to the usual single-phase shape factor. On the other hand, the shape
factor representing an imbibition/saturation diffusion, 𝜎𝑆𝐷, is a function of two parameters, 𝑏 and 𝑚, that
must be obtained either experimentally or numerically:
𝜎𝑆𝐷 = 𝑏𝑡
−𝑚 . (2.21)
Lu et al. (2008) developed a DPDPM and proposed a transfer function that integrates transfer rates from
different physical mechanisms: fluid expansion, diffusion, and displacement. However, their formulation
was only validated for an imbibition process (modeled by capillary and gravitational forces) in blackoil
systems and does not provide a basis for modeling of multicomponent compositional systems (Sherafati
and Jessen, 2018).
2.3. Diffusive Mass Transfer
Molecular diffusion can play an important role in the recovery from tight fractured reservoirs during HnP
operations due to the large surface area available for diffusion. Below, we briefly review the common
methods to represent molecular diffusion and their implementation in reservoir simulation models.
2.3.1. Common Methods to Represent Molecular Diffusion
The generalized Fick’s law and the classical Fick’s law are two main approaches to model single-phase
multicomponent molecular diffusion. A central difference between these models is the representation of
the flux and the evaluation of the diffusion coefficients. The diffusion flux of a given component in the
generalized Fick’s law depends on the concentration gradients of other components and uses a full matrix
for diffusion coefficients:
𝐽𝑖 = −𝜌 ∑ 𝐷𝑖𝑘
𝑛𝑐−1
𝑘=1
∇𝑧𝑘, 𝑖 = 1,… , 𝑛𝑐 − 1 , (2.22)

19
where 𝐽𝑖 denotes the diffusive molar flux of component 𝑖 in bulk phase and 𝐷 are the generalized Fickian
diffusion coefficients (𝑛𝑐 − 1 × 𝑛𝑐 − 1 matrix). The generalized Fick’s law is based on the Maxwell-Stefan
framework to obtain the relevant diffusion coefficients.
In contrast, the classical Fick’s law uses effective diffusion coefficients(per component), and the
diffusive flux of a component is driven by the related composition gradient:
𝐽𝑖 = −𝜌𝐷𝑖
eff∇𝑧𝑖
, 𝑖 = 1,… , 𝑛𝑐
. (2.23)
Here, 𝐷⬚
eff are the effective diffusivities per component in the mixture (𝑛𝑐 vector). To facilitate the
calculation with the classical Fick’s law, it is common to use the extended Sigmund correlation to obtain
binary diffusion coefficients of the mixture (Sigmund, 1976a, 1976b; da Silva and Belery, 1989). Following
that, the formulation of Wilke (1950) is applied to obtain the effective diffusivity (see e.g.,Shi et al., 2022).
The interested reader can refer to Taylor and Krishna (1993) and Hoteit (2013) for more details on
the formulation of the above diffusion models (more information is also provided in Chapter 6). Below we
summarize how these models have been utilized to evaluate recovery performance of EOR specifically on
fractured cores.
2.3.2. Gas Injection (Enhanced Oil Recovery) in Fractured Cores
The development of fit-for-purpose models for unconventional reservoirs is often based on interpretation
of lab-scale results (lab data from cores) for benchmarking purposes, followed by some upscaling of the
findings to larger scale. Below, we review the existing efforts on modeling experimental data of gas
injection in fractured cores in the context of EOR.
Darvish et al. (2006) conducted miscible CO2 injection experiments on composite cores with a
permeability of 4 mD to study mass transfer between matrix and fracture subsystems. They concluded
that molecular diffusion contributes significantly to recovery. However, they were unable to predict the
experimental results using a commercial simulator (ECLIPSE) due to the limitations in modeling of diffusion

20
across gas/liquid interfaces. Therefore, an artificial two-phase zone was introduced wherein the
compositions were adjusted to obtain a reasonable agreement with the experimental data.
Vega et al. (2010) studied miscible CO2 injection into a siliceous shale core (1.3 mD) with 3 different
boundary conditions to mimic relevant field scenarios. They demonstrated the significance of molecular
diffusion and convection/dispersion on oil displacement and recovery. In their modeling work, they used
a commercial simulator (CMG-GEM) to predict the experimental data and to study the mass transfer
between matrix and fracture domains. The binary diffusion coefficients were first modeled as a function
of composition using the Sigmund correlation (Sigmund, 1976a, 1976b). Despite their efforts, a good
agreement with the experimental data was not obtained. They also attempted to assign constant
diffusivities as input parameters, however with limited improvement of the calculated results. Wan and
Sheng (2015), suggest that the numerical grid constructed by Vega et al. may have lacked a proper fracture
characterization.
Alharthy et al. (2015) presented laboratory and modeling work for HnP processes using CO2, a
methane/ethane (CH4/C2H6) mixture, and N2 in Bakken cores saturated with Bakken hydrocarbon fluid.
They report, for a middle Bakken core, that CO2 injection results in the highest recovery, while the lowest
recovery was observed from N2. To determine representative diffusion coefficients, they started by
comparing the approximation of several correlations to predict diffusion coefficient measurements for
C1/C2 mixtures (Sigmund, 1976). They concluded that the correlations of Hayduk and Minhas (1982) and
Leahy-Dios and Firoozabadi (2007) agree very well with the lab observations. It is unclear how the Sigmund
(1976) correlation provided a ‘poor’ agreement with the data that was used to develop the correlation.
Alharthy et al. then proceed to explore the approximation of Bakken fluid diffusion coefficients. They
highlighted that the two selected models (Hayduk and Minhas, 1982; Leahy-Dios and Firoozabadi, 2007)
gave comparable approximations, and these values formed the basis for their subsequent numerical
analysis. To performcore-scale calculations, anSPM withradial grid was used (CMG-GEM). The parameters

21
including residual oil saturation, matrix permeability, and relative permeability of the matrix were used to
tune the results and obtain a match with the experimental data.
Sofla et al. (2016) conducted lab experiments and numerical modeling on gas injection into
fractured carbonate cores. They used an SPM approach with grid refinement near the fracture. The
diffusion coefficients for the multicomponent mixture were calculated based on the correlations of Hayduk
and Minhas (1982) and Kooijman and Taylor (1991). From the simulation results, Sofla et al. highlighted
that diffusion contributed little to oil recovery and stated that such conclusion contradicts their
experimental observations. They suggest that this could be due to an inadequate diffusion model being
used in the simulator.
Li et al. (2017) performed HnP experiments on Wolfcamp cores using: CO2, CH4, and N2 at 2,000
psi and 104 F (313.15 K). CO2 injection provided for the highest recovery followed by N2 and CH4. The
experimental conditions were simulated with a commercial simulator (CMG-GEM) however did not
provide for a consistent agreement with the experimental oil recovery: A mismatch was observed for N2
and CH4 during the last cycles, while a reasonable agreement was obtained for CO2. The reported
calculation results were limited to oil recovery while produced oil compositions were not provided.
Ghasemi et al. (2018) studied CO2 injection in a composite chalk core where ECLIPSE was used as
the basis for their numerical simulation efforts. To model diffusive mass transfer, diffusion coefficients
were used as input parameters and were assumed constant throughout the simulation calculations.
Moreover, and as pointed out by Ghasemi et al., the long diffusion time (450 hours) presented a challenge
in matching the oil production using diffusivities that do not reflect the mixture composition. Therefore,
to match the experimental data, the authors divided the experiment into 3 time-segments, each with a
separate set of adjusted diffusion coefficients.
Fu et al. (2021) investigated the role of molecular diffusion during hydrocarbon gas HnP in Eagle
Ford cores (permeability in the range of 1.3-2 D). The experiment was conducted below minimum

22
miscibility pressure; hence, both liquid and gas hydrocarbon phases were expected. They used a DPDPM
approach (ECLIPSE) to represent the core and the diffusion model included inter-phase and intra-phase
flux calculations. To match the pressure-drop during the soaking period andthe produced oil compositions
(of the first cycle), several model inputs were used as adjustable parametersincluding relative permeability
curves, fracture porosity, fracture permeability, and tortuosity.
Alahmari et al. (2023) conducted an integrated experimental and modeling study of CO2 HnP in
Middle East shale cores with a permeability of about 0.85 D. They used an SPM (CMG-GEM) to model
the experimental conditions and to interpret the results. They demonstrated that diffusion (using
composition-dependent coefficients) contributed significantly to the calculated oil recovery; however, the
simulation calculations overestimated the recovery factors observed from the experimental data. To
achieve an improved agreement with the experimental observations, constant effective diffusion
coefficients, selected manually as input parameters into CMG, were used to calculate the diffusionalfluxes.

23
Chapter 3. A Generalized Dynamic Transfer Function for Tight DualPorosity Systems*
Warren and Root (1963) introduced the transfer function and laid the foundation for describing the mass
transfer between the matrix and fracture blocks in DP reservoir simulation. However, the PSS assumption
embedded in the approach of Warren and Root is no longer applicable when the duration of the transient
state is prominent (tight oil or shale gas reservoirs). In later years, Lim and Aziz (1995) derived new shape
factors that avoid the PSS assumption. However, truncating the series solution in the approximation of Lim
and Aziz may also lead to limited accuracy in capturing the pressure gradients inside tight matrix blocks.
In this chapter, we introduce a new generalized dynamic transfer function that can accurately
predict the mass transfer in tight DP formations. The formulation is first derived by generalizing an
approximation of the analytical solution for spherical domains, and then extended to a broader range of
geometries and rock/fluid settings. We demonstrate that the new formulation is similar in form to the
Warren and Root (1963) model, with the addition of a correction term that reflects on the physical
behavior of the fractured system and allows for representation of the transient in the matrix block.
We highlight that given the new transfer function was derived from the analytical solution, the
formulation is equally useful to describe other diffusion-like processes (e.g., heat and concentration). The
ability to model mass transfer insimilar representationsfurther endorses the simplicity of the DP modeling
approach.
*The results in this chapter were presented in: Zhang, J., Raslan, M., Wu, C., et al. 2022. A Generalized Dynamic
Transfer Function for Ultra-Tight Dual-Porosity Systems. Paper presented at the SPE Western Regional Meeting,
Bakersfield, California, USA, April 2022. https://doi.org/10.2118/209324-MS

24
3.1. Derivation of the Generalized Vermeulen Transfer Function
Inspired by the work of Zimmerman et al. (1993), we seek to generalize the closed-form exponential
solution of Vermeulen (1953) to approximate the exact analytical solution in relevant geometries. We then
determine the parameters of the new transfer function for a range of geometries of the matrix segments
(planar, cylindrical, and spherical) and report the performance of the formulation for isotropic and
anisotropic matrix segments. We then extend the application to irregularly shaped matrix blocks. Finally,
we demonstrate the accuracy of the formulation for a highly compressible fluid where hydraulic diffusivity
is expected to vary considerably with pressure.
3.1.1. Isotropic Matrix Blocks
We start by considering a matrix block that is surrounded by three sets of perpendicular fractures with a
fracture spacing of 𝐿, and assume that the matrix block can be approximated by a sphere with an
equivalent radius, 𝑅, providing the same volume as a cube with a length 𝐿 on all sides. The exact pressure
distribution within the matrix block can be found by solving the differential equation:
𝜕𝑃
𝜕𝑡 = 𝐷(
𝜕
2𝑃
𝜕𝑟
2 +
2
𝑟
𝜕𝑃
𝜕𝑟 ) , (3.1)
subjected to the initial and boundary conditions:
𝑝 = 𝑝0
, 0 ≤ 𝑟 ≤ 𝑅 𝑡 = 0
𝑝 = 𝑝𝑓
, 𝑟 = 𝑅 𝑡 > 0
. (3.2)
Following the work of Crank (1975), the exact solution to the variation in the average pressure of this
problem can be written as:
𝑃𝐷 =
𝑃̅𝑚−𝑃0
𝑃𝑓−𝑃0
= 1 −
6
𝜋2
∑
1
𝑛2
exp(−𝐷
𝑛
2𝜋
2
𝑡
𝑅2
)
∞
𝑛=1
. (3.3)
Vermeulen (1953) provided a simple exponential form as an approximation to Eq. (3.3), which is expressed
as (here in terms of pressure):
𝑃𝐷 =
𝑃̅𝑚−𝑃0
𝑃𝑓−𝑃0
= [1 − exp(−
𝜋
2𝐷𝑡
𝑅2
)]
1/2
. (3.4)

25
To extend the application to any matrix geometry, we generalize Vermeulen’s approximation, Eq. (3.4), by
introducing two parameters, 𝜎 and 𝑉:
𝑃𝐷 =
𝑃̅𝑚−𝑃0
𝑃𝑓−𝑃0
= [1 − exp(−𝜎𝐷𝑡)]
1/𝑉
. (3.5)
To arrive at a functional form that is similar to the traditional transfer function of Warren and Root (1963),
Eq. (2.4): 𝛤 = 𝜎
𝑘
𝜇
(𝑃𝑓 − 𝑃̅𝑚), we first differentiate Eq. (3.5) with respect to time and then eliminate time
on both sides of the equation. This results in the following final expression for the generalized Vermeulen
(gVer) transfer function:
𝑑𝑃̅m
𝑑𝑡
= 𝜎𝐷(𝑃𝑓–𝑃̅𝑚) [
1
𝑉
𝑃𝑓−𝑃0
𝑃𝑓−𝑃̅𝑚
((
𝑃̅𝑚−𝑃0
𝑃𝑓−𝑃0
)
1−𝑉
−
𝑃̅𝑚−𝑃0
𝑃𝑓−𝑃0
)] = 𝜎𝐷(𝑃𝑓–𝑃̅𝑚)𝛽 . (3.6)
The expression in the square brackets represents the correction factor, 𝛽, which is analogous to the boost
factor of Zimmerman et al. (1993) (see Eq. (2.14)): 𝛽 assumes values > 1 at early times, and approaches
unity at later times. Further, Eq. (3.6) includes the parameter 𝜎 referring to the concept of the traditional
shape factor and to arrive at the new transfer function that is consistent, in the functional form, to the
classical transfer function. Further, the product 𝜎𝐷 denotes the inverse of the characteristic time of the
relevant diffusion process inside the matrix block. We note that the term 1/𝐿
2
is embedded into the shape
factor and that in this representation of the matrix geometry, an explicit calculation of the characteristic
lengthsis no longer needed. This becomes particularly useful when working with irregularly shaped matrix
blocks and complex fracture settings.
The model parameters, 𝜎 and 𝑉, are constants to be determined by fitting the transfer function to
a reference solution which can be in the form of anexact analytical solution or a numerical solution (more
on this in Chapter 4). To demonstrate the application of the workflow for the case of isotropic spherical
matrix blocks, the approximate solution seen in Eq. (3.5) is fitted to the exact solution given in Eq. (3.3)
using nonlinear regression. From Fig. 3-1, we observe an excellent agreement between the reference
solution and the approximate solution.

26
Figure 3-1: Comparison of the exact solution and the approximate solution to obtain the generalized Vermeulen transfer function
parametersfor isotropic spherical matrix blocks(three sets of fractures)
A similar nonlinear regression analysis is carried out for two other matrix/fracture geometries: 1) One set
of fractures (planar); and 2) Two sets of fractures (cylindrical). In the case of the planar geometry, we
consider a matrix block surrounded by a set of parallel fractures of length 𝐿. The exact solution can be
found by solving the pressure diffusion equation:
𝜕𝑃
𝜕𝑡 = 𝐷
𝜕
2𝑃
𝜕𝑥
2
, (3.7)
subjected to the initial and boundary conditions:
𝑝 = 𝑝0
, − 𝐿/2 ≤ 𝑥 ≤ 𝐿/2 𝑡 = 0
𝑝 = 𝑝𝑓, 𝑥 = −𝐿/2 𝑡 > 0
𝑝 = 𝑝𝑓, 𝑥 = +𝐿/2 𝑡 > 0
. (3.8)
The analytical solution can be expressed as follows (Crank, 1975):
𝑃𝐷 =
𝑃̅𝑚−𝑃𝑖
𝑃𝑓−𝑃𝑖
= 1 − ∑
8
(2𝑛+1)2𝜋2
exp(𝐷
−(2𝑛+1)
2𝜋
2
𝑡
𝐿
2
)
∞
𝑛=0
. (3.9)
In the case of the cylindrical geometry, we consider a matrix block that is surrounded by two sets of
perpendicular fractures, each set having a fracture spacing of 𝐿. The pressure diffusion can be
approximated by the diffusion of a cylindrical matrix block with an equivalent radius, 𝑅, providing the same
volume as a bar with a 𝐿 × 𝐿 cross-section. The partial differential equation is accordingly:

27
𝜕𝑃
𝜕𝑡 =
1
𝑟
𝜕
𝜕𝑟 (𝑟𝐷𝜕𝑃
𝜕𝑟 ) , (3.10)
subjected to the initial and boundary conditions:
𝑝 = 𝑝0
, 0 ≤ 𝑟 ≤ 𝑅 𝑡 = 0
𝑝 = 𝑝𝑓
, 𝑟 = 𝑅 𝑡 > 0
. (3.11)
It follows that the analytical solution can be expressed as (Crank, 1975):
𝑃𝐷 =
𝑃̅𝑚−𝑃0
𝑃𝑓−𝑃0
= 1 − ∑
4
𝑅2𝛼𝑛
2
exp(−𝛼𝑛
2 ∞ 𝑡𝐷) 𝑛=1
, (3.12)
where 𝛼𝑛’s are the roots of the Bessel function of the first kind of order zero: 𝐽0
(𝑅𝛼𝑛
) = 0.
Once again, the approximate solution (see Eq. (3.5)) is fitted to the respective exact solutionsgiven
in Eqs. (3.9) and (3.12) using nonlinear regression (see Fig. 3-2).
Figure 3-2: Comparison of the respective exact solutions and the approximate solution to obtain the generalized Vermeulen
transfer function parameters for isotropic planar and cylindrical matrix blocks
Finally, the resulting values for the model parameters, 𝜎 and 𝑉, from the regression analysis pertaining to
one, two, and three sets of fractures are tabulated in Table 3-1 below.
Table 3-1: Model parameters for planar, cylindrical, and spherical geometries for isotropic matrix blocks.
Set(s) of fractures Geometry 𝜎𝐿
2 𝑉
1 Planar 8.76 1.54
2 Cylindrical 15.97 1.87
3 Spherical 23.16 2.06

28
3.1.2. Anisotropic Rectangular Matrix Blocks
Next, we consider a rectangular(prism)matrix block with the dimensions 𝐿𝑥
, 𝐿𝑦 , and 𝐿𝑧 andpermeabilities
of 𝑘𝑥
,𝑘𝑦, and 𝑘𝑧
in the corresponding directions. For a depletion process with an initial pressure of 𝑃0
subjected to a constant fracture pressure of 𝑃𝑓
(𝑡 > 0), the analytical solution for one set of fractures
surrounding an anisotropic matrix block is given by (Holman, 1990):
𝑃𝐷 =
𝑃̅𝑚−𝑃0
𝑃𝑓−𝑃0
= 1 − ∑ (
8
𝜋2
)
exp[
−𝜋
2𝑡
𝜙𝜇𝑐𝑡
{
𝑘𝑥
𝐿𝑥
2
(2𝛼+1)
2}]
(2𝛼+1)2
∞
𝛼=0
. (3.13)
The solution for two sets of fractures surrounding an anisotropic matrix block is given below:
𝑃𝐷 =
𝑃̅𝑚−𝑃0
𝑃𝑓−𝑃0
= 1 − ∑ ∑ (
8
𝜋2
)
2
exp[
−𝜋
2𝑡
𝜙𝜇𝑐𝑡
{
𝑘𝑥
𝐿𝑥
2
(2𝛼+1)
2+
𝑘𝑦
𝐿𝑦
2
(2𝛽+1)
2
}]
(2𝛼+1)2(2𝛽+1)2
∞
𝛽=0
∞
𝛼=0
. (3.14)
Lastly, the solution for three sets of fractures surrounding an anisotropic matrix block is:
𝑃𝐷 =
𝑃̅𝑚−𝑃0
𝑃𝑓−𝑃0
= 1 − ∑ ∑ ∑ (
8
𝜋
2
)
3
exp[
−𝜋
2𝑡
𝜙𝜇𝑐𝑡
{
𝑘𝑥
𝐿𝑥
2
(2𝛼+1)2+
𝑘𝑦
𝐿𝑦
2
(2𝛽+1)2+
𝑘𝑧
𝐿𝑧
2
(2𝛾+1)2}]
(2𝛼+1)2(2𝛽+1)2(2𝛾+1)2
∞
𝛾=0
∞
𝛽=0
∞
𝛼=0
. (3.15)
Similarly, we use nonlinear regression to fit our approximate solution, Eq. (3.6), to the respective analytical
solutions, Eqs. (3.13)-(3.15), to find the relevant model parameters. A comparison of the approximate and
exact solutions is provided in Fig. 3-3.
Figure 3-3: Comparison of the respective exact solutions and the approximate solution to obtain the generalized Vermeulen
transfer function parameters for anisotropic matrix blocks with one, two, and three sets of fractures

29
The values of themodel parametersof the gVer transfer function for anisotropic matrix blocksare reported
in Table 3-2 below.
Table 3-2: Generalized Vermeulen transfer function parameters for anisotropic rectangular matrix blocks.
Set(s) of fractures 𝜎𝐿
2 𝑉
1 8.76 1.54
2 17.44 1.97
3 27.49 2.27
When comparing the values from Tables3-1 and 3-2, we note that 𝜎 and 𝑉 vary according to the geometry
and permeability anisotropy.
3.2. Numerical Examples
The generalized Vermeulen transfer function was implemented in the DP module in MRST for validation
and testing. In this section, we report the results of such efforts using single-block and sector model
calculations for fractured systems. The objective is to demonstrate, for tight systems, the accuracy of the
gVer transfer function,as well as the limitations of the traditional transfer function with shape factors from
Lim and Aziz (1995).
3.2.1. Single-block Calculations for an Isotropic Matrix Block
We proceed to validate the new generalized Vermeulen transfer function in a two-step procedure. First,
the results of MRST are validated against commercial simulators (CMG and ECLIPSE) using an SPM
representation of a tight fractured matrix block. Then, single-block DPM calculations are performed and
compared with the fine-grid SPM solution. The DPM is modeled via 1) the traditional transfer function
(from Eq. (2.4): 𝛤 = 𝜎
𝑘
𝜇
(𝑃𝑓 − 𝑃̅𝑚)) using the shape factor fromLim and Aziz (1995); and 2) the generalized
Vermeulen transfer function (see Eq. (3.6)).

30
We simulate the production of a dry gas (single-phase system) from a tight fractured matrix block.
We consider a fine-grid SPM with 50 x 50 x 1 cells of 20 ft x 20 ft x 20 ft, and a single-block DPM of the
same size. A grid refinement study was initially conducted to ensure convergence of the SPM reference
solution. We consider a matrix block with a porosity of 20% that is in contact with2 sets of perpendicular
fractures: Fracture permeability = 10 mD and fracture porosity = 2%. We evaluate the pressure response
of two scenarios for the matrix permeability: 1 D and 5 nD. The fluid used in this calculation example
represents a dry natural gas mixture with properties extracted from real reservoir data (fluid properties
are detailedin Fig. 3-4). The initial model pressure is set at 8,000 psi and the fluid is produced at a constant
fracture pressure of 1,000 psi.
Figure 3-4: Fluid (dry gas) properties used in Example 3.2.1: Formation volume factor (𝐵𝑔
), viscosity, and density
The calculation results, in terms of the average matrix pressure over time, are reportedin Figs. 3-5 (for a
1 D matrix block) and 3-6 (for a 5 nD matrix block). InFig. 3-5 (left panel), we compare calculation results
for the SPMs and we observe that MRST is in excellent agreement with the commercial simulators.
Furthermore, for the DPMs (right panel of Fig. 3-5), when using the traditional transfer function and the
shape factor from Lim and Aziz, we observe a good agreement between the three simulators. However,
we observe a large deviation from the reference solution at the early time with better agreement after 10
days. In contrast, the generalized Vermeulen transfer function agrees with the fine-grid solution at both
the early and late production times. Figure 3-6 reports a similar comparison for a 5 nD matrix block.

31
Figure 3-5: Average matrix pressure for a 1 D matrix block: SPM calculations with 3 simulators(left panel) and DPM
calculations using the traditional and new transfer functions(right panel)
Figure 3-6: Average matrix pressure for a 5 nD matrix block: SPM calculations with 3 simulators (left panel) and DPM
calculations using the traditional and new transfer functions (right panel)
We observe that the traditional transfer function with the shape factor from Lim and Aziz approaches the
reference solution only after 1,500 days of production, which is significantly later than for the 1 D rock
example above. Again, the results from the new transfer function outperforms the traditional approach
and agrees substantially with the reference SPM solution for all times.

32
This initial example validates the newly developed generalized Vermeulen transfer function and
shows its accuracy relative to a fine-grid SPM reference solution. The calculation results also emphasize
the problem with the truncated series solution used to derive the shape factors by Lim and Aziz (1995),
where errors at early dimensionless time may be acceptable for higher matrix permeability (Fig. 3-5) but
result in a substantial overestimation of the average reservoir pressure for tight fractured rocks (Fig. 3-6).
3.2.2. DP Sector Model Calculations
Next, we consider the production of a dry gas from a fractured reservoir represented by a DP sector model
using: 1) the traditional transfer function with shape factor from Lim and Aziz (using MRST, CMG, and
ECLIPSE); and 2) the generalized Vermeulen transfer function.
We assume that the matrix blocks are uniform and surrounded by two sets of perpendicular
fractures. The computational domain represents a 1,000 ft x 1,000 ft x 50 ft section of the subsurface and
is represented by 50 x 50 x 1 DPM cells. Production calculations were performed with the fluid and rock
properties of Example 3.2.1. The initial reservoir pressure is set at 8,000 psi with one production well
completed in the center of the domain and operating at a constant bottomhole pressure of 1,000 psi.
Figure 3-7 reports a comparison of the average pressure, as predicted by the three simulators, for
two scenarios: 1) A 1 D matrix permeability (left panel); and 2) A 5 nD matrix permeability (right panel).
For the larger matrix permeability, the average pressure calculated via the traditional transfer function
agrees well with the new transfer function. However, for the matrix permeability of 5 nD, the difference in
the average pressure predicted by the two methods can reach values exceeding 1,000 psi during early
time. The overestimation of the average pressure translates into a large difference in the production
behavior for approximately 1,500 days (~4 years).

33
Figure 3-7: Average matrix pressure in the DP sector model for a 1 D rock (left panel) and 5 nD rock (right panel)
3.2.3. Single-block Calculations for an Anisotropic Matrix Block
In the previous two examples, we have studied cases with an isotropic permeability. Here, we extend the
test to investigate the performance of the new transfer function for an anisotropic permeability problem.
We consider depletion of a single DP matrix block surrounded by three sets of fractures. The dimensions
of the matrix block are 20 ft x 20 ft x 20 ft while the permeabilities are 100 nD, 50 nD, and 10 nD, in the
respective directions (left panel of Fig. 3.8). An equivalent fine-grid SPM was constructed with 50 x 50 x 50
cells (right panel of Fig. 3-8) and we consider the same fluid description as in the previous examples.
Figure 3-8: Schematics of an anisotropic permeability matrix block represented using two approaches: Single-block DPM (left
panel) and the equivalent fine-grid SPM (right panel)

34
Figure 3-9 compares the average matrix block pressure as calculated by the SPM and the two DPMs with
traditional and new transfer functions. We find that the results confirm the observations made in the
previous examples where the new transfer function demonstrates an excellent agreement with the finegrid model, while the traditional DPM approach overestimates the average pressure at early times.
Figure 3-9: Average matrix pressure for an anisotropic rectangular (3D) matrix block
3.2.4. Impact of Diffusivity
In this section, we investigate the impact of the diffusivity on the performance of the generalized
Vermeulen transfer function.
In the derivation of the new transfer function (section 3.1), 𝐷 refersto the hydraulic diffusivity
during pressure-driven processes. Theoretically, 𝐷 should be constant in the partial differential equation
(e.g., in Eq. (3.1)); however, in reality, 𝐷 changes with time as the fluid properties change with pressure.
From Example 3.2.1, the diffusivity was only changing marginally with time during the depletion process,
as shown by the red solid line in the left panel of Fig. 3-10 and hence, the example satisfies the assumption
embedded in Eq. (3.1).

35
Figure 3-10: Variations in diffusivity with pressure for two different fluids(left panel) and average matrix pressure have we
repeated Example 3.1.1 using pure CH4 production (right panel)
To test this assumption and its potential limitation on the performance of the new transfer function, we
performed simulation calculations with pure CH4 (properties extracted from NIST webbook), where we
observe a substantially larger variation in 𝐷 than for the base fluid, as shown by the blue dashed line in
the left panel of Fig. 3-10. The right panel of Fig. 3-10 compares the calculation results with the new
transfer functionwith the reference solution for pure methane. We observe that the new transfer function
still agrees well with the SPM results even as 𝐷 changes distinctly during the depletion process.
Accordingly, we find that the accuracy of the generalized Vermeulen transfer function is largely
independent of changes in diffusivity.
3.2.5. Irregular Matrix Blocks
Due to the complex geologic features of the subsurface, it is common to use corner-point geometry or an
unstructured mesh in reservoir simulation. Therefore, a test of our new transfer function for irregular
matrix blocks is warranted to demonstrate the broader application of this work. However, as the shape
factor in the transfer function depends on the geometry of the matrix block, the model parameters

36
presented earlier in Tables 3-1 and 3-2 are not applicable anymore. Therefore, we need to fit our
approximate solution to a reference solution that describes this DP system.
In this example, we use an SPM as our reference solution to determine a new set of 𝜎 and 𝑉 of
the new transfer function, Eq. (3.6), for an irregular matrix block defined by corner point geometry. The
left panel of Fig. 3-11 shows the model setup where the DPM represents the domain as a single block. We
then divide the block into 50 × 50 × 50 cells to define the SPM (right panel of Fig. 3-11) and allow the
production from the outer boundaries of the grid. The fluid properties are the same as used in Example
3.2.1 while the matrix permeability is assumed to be isotropic and is equal to 5 nD.
For both approaches, the initial reservoir pressure is set at 8,000 psi and a constant production
pressure of 1,000 psi is set at the boundaries.
Figure 3-11: Schematics of an irregular isotropic matrix block using two different approaches: Single-block DPM (left panel) and
fine-grid SPM (right panel)
The average pressure from the SPM was then used to determine the parameters 𝜎and 𝑉 of the generalized
Vermeulen transfer function via nonlinear regression. The resulting transfer function parameters are 𝜎 =
0.94 ft -2 and 𝑉 = 2.30. From Fig. 3-12 (left panel), we observe a good agreement between the two models.

37
Figure 3-12: Average matrix pressure ofthe irregular matrix from SPM and two DPMs (with tuned parameters): New transfer
function (left panel) and traditional transfer function (right panel)
To further justify the application of the proposed transfer function, we repeated the above workflow for
the traditional transfer function by estimating the optimal shape factor from the SPM calculation results.
From Fig. 3-12 (right panel), we observe that since the traditional transfer function is restricted to a fixed
(time invariant) shape factor, it is impossible to match both early and late time behaviors in tight fractured
formations: During the early time, the traditional transfer function departs substantially from the
reference solution.
3.3. Summary and Conclusion
In this chapter, a new dynamic DPM transfer function was presented for modeling of fluid flow between
matrix and fracture segments with emphasis on tight DP systems. The new transfer function is derived by
generalizing Vermeulen’s closed-form approximation and introducing 𝜎 and 𝑉. Nonlinear regression is
used to match the approximate solution to a reference solution and estimate the two model parameters.
The new formulation was validated by comparing DPM results with fine-grid SPM calculations
using MRST, CMG, and ECLIPSE. The validation was initially performed for an isotropic permeability and

38
subsequently extended to an anisotropic permeability setting. These observations were also
demonstrated to be consistent for a depletion process with large variations in the hydraulic diffusivity.
A simple workflow was presented to allow for the application of the new transfer function when
the formation is represented by an irregular mesh: Each matrix block is meshed and simulated as an SPM,
and the results are used to estimate the model parameters of the new transfer function.
The generalized Vermeulen transfer function has a simple functional form that is similar to the
traditional Warren and Root (1963) type transfer function. This provides for an easy implementation in
existing reservoir simulators like MRST.
The proposed improvement of the transfer function in DP modeling has significant practical values.
For instance, if the initial reservoir pressure is above the bubble or dew point, once the production begins,
the depletion strategy may be hinged on the prediction of the average pressure in the formation. If a
traditional transfer function is used for tight formations, the depletion rate may be misinterpreted and
resultsin unexpected production behavior (more on this in Chapter 5).

39
Chapter 4. Upscaling of Discrete Fracture Characterization to DualPorosity Models Using an Advanced Transfer Function*
In Chapter 3, we developed a new transfer function, the generalized Vermeulen transfer function, which
is capable of accurately capturing the early- and late-time behavior of tight DP formations. We validated
the new formulation for simple geometries (geometries associated with analytical solutions), as well as
irregular matrix blocks. In this chapter, we examine the performance of the new transfer function in more
complex fracture applications.
EDFMs are commonly used to study fluid flow in unconventional reservoirs. The EDFM, however,
requires extensive pre-processing of non-neighboring connections. More importantly, it becomes further
demanding when a refined grid is required to accurately model mass transfer between the fractures and
the tight matrix rock (to resolve transients in the matrix). This renders the EDFM an expensive option to
model large-scale fractured reservoirs where the number of fracture segments can be substantial.
We propose an upscaling procedure to construct a DPM from a discrete fracture characterization.
The implicit representation of the fractures in the DPM approach provides for an improved computational
efficiency. While the DPMs are often perceived as simple sugar-cube representations of complex fracture
networks, the upscaling technique presented here, demonstrates the capability of DPMs to provide
accurate and efficient solutions for a broad range of complex fractured systems.
*The results in this chapter were presented in: Raslan, M. S., Gupta, A., Ates, H., et al. 2023. Upscaling of Discrete
Fracture Characterization to Dual-Porosity Models Using an Advanced Transfer Function. Paper presented at the
Middle East Oil, Gas and Geosciences Show, Manama, Bahrain, February 2023. https://doi.org/10.2118/213208-MS

40
The potential of the DPM is unlocked via the application of the gVer transfer function, which as
mentioned in Chapter 3, incorporates a time-dependent correction factor to resolve the transient within
the matrix without the need for sub-gridding. The use of an advanced DPM becomes key to our proposed
upscaling workflow which also includes the evaluation of effective rock properties and a unique set of
matrix/fracture interaction parameters.
We present a set of example calculations for tight fractured systems to demonstrate the need for
advanced modeling techniques to accurately capture the matrix transient, as observed from refined SPMs.
For these examples, the results from explicit models are compared to EDFM and DPM (with traditional
transfer function and the gVer transfer function) and serve as an introduction to the proposed upscaling
technique. We then proceed by reporting on the performance of the upscaling approach for tight systems
with large and complex fracture networks.
4.1. Importance of Capturing the Transient Behavior
We start by examining several test cases, representing tight fractured systems, and establish the need for
advanced modeling techniques to accurately capture the matrix transient, as observed from fine-grid
SPMs. In the context of EDFM, this can be accomplished by refining the mesh of the matrix grid at an
increased computational cost, while the use of gVer transfer function to describe the matrix/fracture mass
transfer in the DPM is demonstrated to accomplish the same results at a marginal additional
computational cost.
4.1.1. One Set of Parallel Fractures
We start by considering a 1D matrix segment with a 20% porosity and a 5 nD permeability, bounded by a
set of highly permeable fractures that actsas the boundary condition on the matrix block. We construct a
30 ft x 30 ft x 5 ft fine-grid SPM (47 x 1 x 1 cells) with fractures at each end and use LGR to represent the
details of the pressure response (see left panel of Fig. 4-1). The EDFM is constructed with two embedded

41
orthogonal fractures, and we consider two different grid resolutions: A coarse model with 4 x 1 x 1 cells
and a finer model with 10 x 1 x 1 cells, both with uniform cell dimensions. The DPM uses a single block
and we consider two matrix/fracture interaction models: 1) the traditional transfer function (from Eq.
(2.4): 𝛤 = 𝜎
𝑘
𝜇
(𝑃𝑓 − 𝑃̅𝑚)) with a shape factorfrom Lim and Aziz (1995); and 2) the gVer transfer function
with model parameters reported in Table 3-2.
The system is initialized at 7,000 psi and then subjected to a drawdown of 5,500 psi by imposing
a constant pressure boundary condition in the fractures.
Figure 4-1: Domain representing a 1D reservoir with a known fracture spacing (left panel); results of the classical DPM and the
coarse EDFM (middle panel); and results of the DPM (gVer) and the refined EDFM (right panel)
The middle panel of Fig. 4-1 reports the average matrix pressure response of the fine-grid SPM along with
the traditional DPM and the coarse EDFM. We observe that both solutions overestimate the average
pressure for a significant period of time relative to the reference solution. The right panel of Fig. 4-1 reports
the solutions using the gVer transfer function (DPM) and the refined EDFM: We observe a good agreement
between all models. Accordingly, this example emphasizes that additional efforts are required (advanced
transfer function in the DPM and refinement in the EDFM) to match the reference solution for tight
fractured rocks.

42
4.1.2. 2D Flow with Intersecting Fractures
Next, we extend the previous example and introduce a third (orthogonal) fracture plane (fracture 3) with
an aperture of 0.066 ft that intersects fractures 1 and 2 on the boundaries (see left panel of Fig. 4-2). As a
result, this generates a fracture network that dictates a 2D flow setting. Here, we use anSPM with 65 x 66
x 1 cells and apply LGR near the fractures. Similar toExample 4.1.1, we construct two EDFMswith different
resolutions: 15 x 15 x 1 cells and 65 x 66 x 1 cells, both with uniform cell sizes. Rock properties and
boundary conditions are the same as used in Example 4.1.1.
Figure 4-2: Top view of the domain which represents a 2D reservoir with a fracture network, where all fractures are orthogonal
(left panel). The middle and right panels report the average pressure of the respective models
For the DPM calculations, we consider the 2D problem as equivalent to a flow in a cylinder with radius of
21 ft and evaluate the shape factors accordingly. The results from the different models are compared and
reported in Fig. 4-2. We observe that the DPM results, using the shape factor from Lim and Aziz (1995) of
about 0.023 ft-2
, as well as the results from the coarse EDFM do not agree with the reference solution. In
contrast, the refined EDFM and DPM + gVer (𝜎 = 0.023 ft-2 and 𝑉 = 1.97) are found in good agreement
with the fine-grid SPM.
4.1.3. Two Sets of Parallel Intersecting Fractures
In this example, we examine the behavior of a larger fracture network as predicted by EDFM and DPM.
Unlike in Example 4.1.2where the fractures, due to symmetry, resemble a cylindrical shape, the fracture

43
spacing between fractures 3 and 4 are now distanced sufficiently close to generating larger pressure
gradientsin this region (see left panel of Fig. 4-3). In a first set of calculations, we reuse the shape factors
from Example 4.1.2 (pre-calculated from a cylindrical approximation) and compare the results from DPMs
(traditional and gVer) and the coarse EDFM (20 x 20 x 1 cells with uniform cell sizes) with the fine-grid
solution (63 x 63 x 1 cells with LGR).
Figure 4-3: Top view of the domain that includestwo sets of parallel (orthogonal) fractures (left panel). Simulation calculations
of SPM, coarse EDFM, and two DPMs (both performed using pre-determined shape factors assuming cylindrical geometry)
The calculated average pressures, reported in the right panel of Fig. 4-3, demonstrate that all modeling
approaches (both DPMs with pre-calculated transfer function parameters and the coarse EDFM) fail in
reproducing the reference solution. In other words, a better representation of the system geometry is
required to model the matrix/fracture interactions in this setting. We address this challenge in the
following section.
4.2. Upscaling of Fractured Rock
Dershowitz et al. (2000) explored the idea of integrating a discrete fracture network (explicit
representation of fracture) with a DPM representation. They used a geometrical method combined with a

44
statistical method to evaluate the fracture spacing and to calculate the relevant shape factors. They,
furthermore, performed fracture connectivity analysis to obtain an optimal grid size of the upscaled
model.
Karimi-Fard et al. (2006) developed a flow-based upscaling technique, the Multiple Subregion
Model, to construct a DPM from a given fracture characterization. However, the determination of
boundary conditions, and hence the number of subregions, remains unclear and may become
computationally expensive for large models (Li et al., 2015; Ding, 2019).
The goal in this chapteris to introduce a workflow that allows for DPMsto provide the same flow
behavior as a detailed reference solution for a tight fractured system. This can be accomplished by using
an advanced mass transfer function, which as discussed above, incorporates a time-dependent correction
factor to address the transient flow. The generalized Vermeulen transfer function offers such advantage
by using physical parameters (𝑃𝑓,𝑃0
, and 𝑃̅𝑚), combined with 𝜎 and 𝑉 that describe the matrix/fracture
interaction. In the proposed workflow, the relevant matrix/fracture interaction parameters are obtained
from matching a reference solution, which can be in the form of an SPM or a refined EDFM. We note that
this does not impose any restriction on the type of the transfer function: For the traditional transfer
function, the shape factor is the only adjustable parameter that is determined/estimated from a reference
solution.
To lump an EDFM into a single-block DPM entails the calculation of an effective fracture porosity,
the volume of fractures embedded in a cell divided by the cell bulk volume. Additionally, an effective
permeability is calculated by performing single-phase flow simulations on the high-resolution model. This
is done by imposing constant pressure boundary conditions and perform a steady-state flow calculation
on the domain to evaluate an equivalent permeability.
The approach provides for flexibility in the selection of the number of subregions needed to
represent the fracture network. We demonstrate, in the following examples, that a fracture

45
characterization can be approximated by one coarse DP block, or by as many as needed to capture the
spatial variability of the reference model.
4.3. Application of the Upscaling Workflow to Irregular Fracture Geometries
Here we apply the upscaling approach to solve complex problems where we simulate the production from
nano-darcy fractured formations with fracture networks exhibiting arbitrary fracture orientations and
geometries. We furthermore demonstrate the versatility of our approach by refining the DPM in regions
of higher spatial variation to capture the details of the fractured rock as dictated by the original EDFM
representation.
Following the discussion above, we revisit Example 4.1.3 and apply the upscaling approach using
a single-block DPM. This is done by estimating the transfer function parameters by matching the pressure
calculations from the respective DPM to that of the SPM and get the following values: 1) 𝜎 =0.22 ft-2
for
the traditional transfer function; and 2) 𝜎 = 0.11 ft-2 and 𝑉 = 2.16 for the gVer transfer function. A
comparison of the upscaled DPMs and a refined EDFM (63 x 63 x 1 with uniform cell sizes) with the
reference SPM solution is provided in Fig. 4-4.
Figure 4-4: Results of the upscaling technique when applied to Example 4.1.3. Both DPMs use transfer function parameters
obtained from the reference SPM solution

46
We observe that the traditional transfer function, with an estimated 𝜎, closely follows the pressure
behavior, while splitting the departure (error) between early and late times. The gVer transfer function,
with 𝜎and 𝑉 estimated from the SPM, agrees very well with the refined EDFM and the reference solution.
We attribute the upscaled DPM agreement, over the entire time range, to the use of the correction factor
that facilitates the representation of the matrix transient.
4.3.1. Intersecting and Impermeable Fractures
Dictated by depositional environments, one may have to model fluid flow in a fractured reservoir with a
low fracture conductivity. Figure 4-5 depicts such setting in terms of a “+” shaped fracture network located
at the center of a matrix domain (20 ft x 14 ft x 10 ft). No-flow boundary conditions are applied at all sides
except for at the west boundary allowing for gas flow. The intersecting (orthogonal) fractures have a
permeability of 2 nD, while the matrix rock is assigned a permeability of 1 D. This permeability contrast
causes the pressure diffusion to be controlled by the matrix rock, while the fractures act as flow barriers.
Figure 4-5: Two intersecting impermeable fractures. Only the west boundary permits pressure depletion, all others are no-flow
This is a classic example inspired by the work of Tene et al. (2017), who demonstrated that the EDFM
approach can fail when the fracture permeability is lower than that of the matrix. Figure 4-6 reports the
pressure distribution as calculated by an SPM (58 x 58 x 1 with LGR) and a refined EDFM (58 x 58 x 1 cells

47
with uniform cell sizes). The results demonstrate the limitation of the EDFM in capturing the response of
the reference SPM solution: The EDFM allows for leakage through the sealing fractures, implying a
homogeneous domain. Tene et al. (2017) proposed a novel approach to resolve this limitation, however,
at the cost of creating additional matrix/fracture connections.
Figure 4-6: Pressure map at the end of simulation of the SPM and EDFM, respectively
Figure 4-7: Average pressure for all models(left panel): Reference solution, refined EDFM, and 2 upscaled DPMs (traditional form
and gVer). Error in average pressure relative to that of the reference solution (right panel)
The relevant matrix/fracture interaction parameters are: 1) 𝜎 = 6.58 ft-2
for the traditional transfer
function; and 2) 𝜎 = 2.03 ft-2 and 𝑉 = 6.34 for the gVer transfer function. Figure 4-7 reports the average

48
pressure in the domain and the error percentage as calculated by the DPMs (traditional and gVer with
parameters estimated using the upscaling technique), the EDFM and the reference fine-grid solution. We
find that the upscaled DPM using the gVer transfer function agrees well with the reference solution,
whereas the EDFM (independent of the grid resolution) and the upscaled DPM with the traditional transfer
function depart from the fine-grid SPM result.
4.3.2. Complex Fractured Rock
Next, we test the upscaling approach for a problem containing multiple fractures with arbitrary
orientations in a 3D model. The fracture system used in this example was generated as one realization of
a fracture network (Alghalandis, 2017). The model represents a synthetic fractured reservoir with 20 x 20
x 5 matrix cells of uniform grid cell sizes (16 ft x 16 ft x 26 ft) including 300 natural fractures, resulting in a
total of 5,164 active cells. The fracture network has a mean length of 26 ft, and a maximum length of 33
ft. The fractures are represented with 45o or -45o dip angles. The conductivity of the fractures was set to
27 mD-ft, and the matrix rock permeability to 500 nD. The resulting EDFM grid is shown in Fig. 4-8.
The model was initialized at 7,000 psi and a constant pressure boundary condition of 500 psi is
imposed on all four sides in the vertical plane.
Figure 4-8: A synthetic fractured reservoir with 300 natural fractures in a 500 nD matrix modeled using an EDFM (left panel: 3D
view and right panel: top view)

49
Gas production was simulated with EDFM and DPMs for this complex model: We applied the upscaling
approach for a single-block DPM and estimated the transfer function parameters from the EDFM as the
reference solution. The relevant model parameters are: 1) 𝜎 = 3.76e-4 ft-2
for the traditional transfer
function; and 2) 𝜎 = 2.73e-4 ft-2 and 𝑉 = 1.62 for the gVer transfer function. A comparison of calculation
results from the two upscaled DPMs and the EDFM is provided in Fig. 4-9.
Figure 4-9: Average matrix pressure from the EDFM and the two upscaled DPMs (left panel). Insert, in left panel, demonstrates
early-time departure of the traditional transfer function by up to 600 psi (see difference in the right panel)
The results demonstrate that the upscaling technique, using the gVer transfer function, can successfully
reproduce the results of the EDFM with very high accuracy. The traditional transfer function, however,
approaches the reference solution after 4 years of production with an observed overestimation of the
pressure by 600 psi at early times. At later times the traditional transfer function is observed to
underpredict the pressure: A single and constant value of the shape factor cannot represent the pressure
transient accurately at both early and late times in tight fractured formations.
We extend Example 4.3.2 by carrying out additional simulation calculations based on several
fracture realizations. The objective is to demonstrate the impact of the generalized Vermeulen transfer
function parameters on different matrix/fracture settings. Recalling that the original problem (Example

50
4.3.2) considers 300 natural fractures with an average length of 26 ft that are intersecting at a 90o angle
and contained in a 328 ft x 328 ft x 131 ft domain (seen earlier inFig. 4-8 above). We conducteda sensitivity
analysis of EDFM input parameters including the number of natural fractures,dip angle, average length of
natural fractures, and domain size. The EDFM schematicsof the additional settings are shown in Fig. 4-10.
Furthermore, Table 4-1 summarizes the resulting 𝜎 and 𝑉 values after fitting the equivalent upscaled
single-block DPM, with the generalized Vermeulen transfer function, to the reference EDFM solution for
the respective setting.

51
Figure 4-10: a) Schematicsfor setting #1 where fractures are parallel (non-intersecting); b) Settings #2, 3, and 4 where EDFM
contains different number of fractures: 200, 100, and 10, respectively; c) EDFM with different average fracture lengths(ft): 16,
29, and 43, respectively; and d) EDFM of different domain sizes (ft): 393 x 393 x 157 and 984 x 984 x 393

52
Table 4-1: Generalized Vermeulen transfer function parameters after fitting the upscaled single-block DPM to the EDFM for the
different EDFM realizations
Sensitivity study Generalized Vermeulen
Transfer Function
# Input Parameter Input Value 𝜎, 10-4 ft-2 𝑉
Original problem
(Example 4.3.2)
• Orientations of 45o and -45o
• 300 natural fractures
• Average length 26 ft
• Domain 328 ft x 328 ft x 131 ft
2.73 1.62
1 Orientation of fractures All 45o
(parallel fractures) 2.52 1.66
2
Number of fractures
200 2.29 1.70
3 100 1.97 1.76
4 10 1.71 1.80
5
Average fracture length, ft
16 2.13 1.72
6 29 3.84 1.47
7 43 7.71 1.20
8
Domain size, ft
393 x 393 x 157 1.54 1.70
9 984 x 984 x 393 0.22 1.69
From Table 4-1, a number of observations can be made regarding the behavior of the matrix/fracture
interaction for the selected fracture realization. For example, the impact of specifying parallel natural
fractures (non-intersecting) (see panel “a” in Fig. 4-10) introduces a very moderate change in 𝜎 and 𝑉
relative to the original problem (see Fig. 4-8). In addition, we observe that removing 100 natural fractures
results in a 16% reduction of the shape factor and a have negligible impact on 𝑉 (settings #2 and 3 in Table
4-1). Further, if we only consider 10 natural fractures (setting #4), we observe a 37% decrease in the
resulting shape factor, while the value of 𝑉 is similar in magnitude to settings #2 ad 3 (see Fig. 4-10 – b).
Based on the above sensitivity analysis, we can conclude that the EDFM input parameters, to
describe the fracture segments (or fracture network), have varying degrees of impact on the resulting
transfer function parameters. Certainly, the observations made from this study is consistent with our
understanding of the shape factor: From Eq. (2.6), 𝜎 = ∑ 𝐴⁄(𝐿𝑉𝑚), the larger the matrix/fracture surface

53
area, per unit bulk volume and characteristic length, the larger the value of the shape factor, and vice
versa. Additionally, from the above cases, we observe that the value of 𝑉 remains greater than 1.
4.3.3. Application to Unconventional Reservoir
All calculation examples presented above considered a simple set of petrophysical properties describing
the fractured rock. In this example, we assume a more realistic fractured formation with a heterogeneous
permeability distribution. This allows us to further examine the performance of the proposed upscaling
technique in a reservoir exhibiting spatial variability of the rock and fracture properties.
We utilize the fracture network from Example 4.3.2 but assign two different fracture conductivities
of 27 mD-ft and 53 mD-ft (left panel of Fig. 4-11). Furthermore, the histogram (right panel of Fig. 4-11)
illustrates the bimodal nature of the matrix permeability distribution. The model has 5 layers with one
being much tighter than the remaining with a permeability of around 60 nD.
Figure 4-11: Fracture conductivities that are assigned with two distinct values (left panel) and matrix permeability describing a
heterogeneous and tight reservoir (right panel)
To upscale the EDFM, we assume two sub-regions and select a dual-block DPM grid that serves as ‘one
realization’ of the fractured reservoir (see Fig. 4-12).

54
Figure 4-12: Upscaling process of a complex fractured reservoir. Left panel: EDFM with two sub-regions. Right panel: Equivalent
upscaled dual-block DPM
To construct the upscaled DPM, we start by conducting a single-phase flow simulation of the EDFM to
obtain the effective permeability per sub-region. Then, simulation of each grid block of the DPM was
performed, separately, to obtain the transfer function parameters by matching the results to the
corresponding EDFM sub-region (reported in Table 4-2).
Table 4-2: Transfer function parameters obtained by matching the respective DPM to the reference EDFM solution
EDFM Sub-region Advanced Transfer Function Traditional Transfer Function
𝜎, 10-4 ft-2 𝑉 𝜎, 10-4 ft-2
Lower 3.77 1.63 5.19
Upper 2.19 1.59 3.1
The constant pressure boundary conditions applied in this example are shown in Fig. 4-13. With all the
relevant input parameters available, the two blocks were combined in calculations for the upscaled 2-block
DPM representation of the EDFM.

55
Figure 4-13: Boundary conditions of the sub-divided EDFM
Figure 4-14 compares the upscaled DPMs with the reference EDFM in terms of the average matrix
pressure. The results are consistent with the observations made from previous examples: The upscaled
DPM using the gVer transfer function provides for an excellent agreement with the complex EDFM. The
right panel of Fig. 4-14 reports the correction term used in the gVer transfer function (per block) to capture
the early-time response. At early times, the correction factor assumes a value of ~20, and declines until it
reaches PSS (~1) and hence converges to the traditional transfer function at later times.
Figure 4-14: Average pressure on the primary y-axis, while the secondary y-axis is the difference between each DPM and EDFM
(right panel). The right panel reportsthe transfer function correction factor

56
4.4. Discussion
The objective of this study was to eliminate the need for simulation of complex and massive fracture
networks when, in reality, there exists no technology today that determines the exact orientation,
geometry, and location of all fractures in the subsurface. In fact, one can argue that a high-resolution
fracture model could contain more details than can be justified by the available data. Gong and Rossen
(2017) have demonstrated that, for certain applications, even for a well-connected fracture network, the
flow path is preferential to a sub-network of the fractures. They suggested that by removing unimportant
fractures, the system can still behave like the original fracture network. Therefore, it may be tempting to
believe that with the increase in computing power one will seamlessly be able to run high-resolution
models with complex fractures, which would in turn justify the use of discrete representation of fracture
networks.
In Example 4.3.3, the selection of the two sub-regions were made so that each region is exposed
to a different pressure regime (Fig. 4-15). This reflects the upscaling technique’s ability to honor the spatial
variability by assigning distinct model parameters for each region. Based on an independent study of rock
quality regions, one could extend this exercise by selecting different grid block sizes in the DPM to further
reflect the heterogeneity of the reservoir.
Figure 4-15: Average matrix pressure of the two EDFM sub-regions of Example 4.3.3

57
To this extent, Example 4.3.3 was re-run by upscaling the detailed fracture characterization of the EDFM
to a 4-block DPM. The workflow follows the same steps as presented earlier (relevant model parameters
are reported in Table 4-3).
Table 4-3: Transfer function parameters obtained by matching the respective DPM (advanced and traditional) to the reference
EDFM solution
EDFM Sub-region Advanced Transfer Function Traditional Transfer Function
𝜎, 10-4 ft-2 𝑉 𝜎, 10-4 ft-2
Lower Left 2.56 1.71 3.37
Lower Right 2.65 1.63 3.72
Upper Left 2.39 1.60 3.31
Upper Right 2.86 1.56 3.89
The results reported in Fig. 4-16 indicate that the proposed method still works and that the advanced
transfer function accurately predicts the pressure behavior, per sub-region. An important point from this
exercise is the inherent reduction in CPU time between simulating a 4-block DPM relative to the full EDFM
with 5,164 active cells: This, in turn, provides for efficient upscaling of larger EDFMs, where the full model
can be divided into any number of sub-regions for fast evaluation of transfer functions parameters.
Figure 4-16: Potential upscaling of the complex EDFM of Example 4.3.3

58
As pointed out previously, the application of the gVer transfer function consistently results in non-unit
exponents, 𝑉. This supports the earlier discussion regarding the role of this parameter to describe the
matrix transient that is critical for unconventional reservoirs. In other words, this implies that a traditional
transfer function will never be able to capture early-time and late-time behavior in tight fractured rocks
when only using a fixed value of 𝜎.
4.5. Summary and Conclusion
We have proposed a technique to upscale complex discrete fracture characterization to a coarse DPM
representation. The use of an advanced transfer function, which honors the longer transient time of
unconventional reservoirs and hence has the advantage of capturing early and late times at a minimal
additional computational cost, is key to the accuracy of this method.
The advanced transfer function was first validated against fine-grid SPM for idealized fracture
geometries. The upscaling technique was then tested for complex fracture representations by comparing
pressure responses to calculations using the EDFM approach. It was demonstrated that the proposed
approach adapts very well to settings with variability in rock and fracture properties by considering
different resolutions without losing accuracy or efficiency.

59
Chapter 5. A Two-Phase Transfer Function for Simulation of Tight
Fractured Gas-Condensate Reservoirs*
In modeling of tight fractured gas-condensate reservoirs, it is particularly important to represent liquid
dropout, and its impact on fluid flow behavior, accurately. This, in turn, allows for design and optimization
of reservoir production strategies, to mitigate potential production loss below the dewpoint pressure.
Unfortunately, the traditional transfer functions, used to describe the mass transfer between
matrix/fracture segments in DPMs, do not explicitly represent the complex physics of the inherent
multiphase problem, which may complicate the application of the existing simulators. Therefore, the
objective of this work is to develop a two-phase transfer function that improves the modeling of gascondensate systems via a DP representation.
Commercial simulators typically offer a simple matrix/fracture mass transfer model to perform
calculations in DP systems: The default shape factors are based on a PSS approximation, and for multiphase
problems, the transfer function of Kazemi et al. (1976) is commonly used. In the coarse DPM approach,
the constant shape factors cannot represent the extended transient behavior observed in tight fractured
reservoirs. Furthermore, the traditional transfer functions average block properties to evaluate the mass
transfer: In compositional modeling of gas-condensate reservoirs, such averaging masks the compositional
changes near the matrix/fracture interface, which is critical to prediction of production behavior below
the dewpoint.
*The results in this chapter were presented in: Raslan, M. S., Gupta, A., Ates, H., et al. 2023. A Two-Phase Transfer
Function for Simulation of Fractured Ultratight Gas-Condensate Reservoirs. Paper presented at the Middle East Oil,
Gas and Geosciences Show, Manama, Bahrain, February 2023. https://doi.org/10.2118/213205-MS

60
To overcome the limitations of the current models, we introduce a new matrix/fracture transfer
function that includes:
1) A time-dependent correction factor (from the generalized Vermeulen transfer function) that
effectively allows the shape factors to vary with time.
2) A modified representation of two-phase flow between matrix and fracture segments. Instead
of evaluating the mass transfer at an average pressure and composition in the matrix, we
propose a novel approach that evaluates the transfer function at conditions that
reflect/approximate the pressure and composition/saturation gradients within the matrix
block. By approximating the fluid state near the matrix/fracture interface, we represent the
liquid buildup that dictates the fluid mobility more accurately.
The proposed method is motivated by the physics of the problem and requires no-subgridding or
pseudofunctions. The enhanced simulation of gas-condensate systems enables improved decision making
in reservoir management and supports optimization of field development plans. In addition, the improved
calculation of the surface streams, obtained from the new method, provides for a more accurate design
of field surface facilities.
5.1. A Dual-Porosity Compositional Model in MRST
In this work, we continue to use MRST as a modeling and simulation platform for testing and
implementation of new developments. The base code supports fully implicit compositional simulation for
SPMs and was previously validated against other simulators (ECLIPSE and AD-GPRS) for a test case that
describes a CO2 injection process (Moyner, 2021). Building on this model, we developed a two-phase
(hydrocarbon vapor/liquid) DPM to facilitate studies relevant of mass transfer in tight fractured reservoirs.

61
The isothermal multicomponent and multiphase flow in a fractured porous medium can be
described by a set of component material balance equations. We formulate the mass conservation of all
hydrocarbon components (in residual form) as:
Fracture: 𝜙𝑓
𝜕𝑀𝑖,𝑓
𝜕𝑡 + 𝛻 ∙ 𝐹𝑖 + 𝛤𝑖 = 0 , 𝑖 = 1,… , 𝑛𝑐
, (5.1)
Matrix: 𝜙𝑚
𝜕𝑀𝑖,𝑚
𝜕𝑡 − 𝛤𝑖 = 0 , 𝑖 = 1, …, 𝑛𝑐
, (5.2)
where 𝑀 and 𝐹 denote the total molar density and convective flux (described by Darcy’s law) of
component 𝑖 in all the available phases, respectively. 𝑛𝑐
is the number of components in the mixture.
To model mass transfer, either the traditional transfer function or gVer transfer function can be
extended to compositional problems by including the mole fraction, 𝑧𝑖𝑗 ,of a given component 𝑖 in phase
𝑗. To account for two phases, the generalized Vermeulen transfer function was initially extended following
the approach of Kazemi et al. (1976):
𝛤𝑖𝑗 = 𝛽𝜎𝜌𝑗
𝑘𝑘𝑟𝑗
𝜇𝑗
𝑧𝑖𝑗(𝑃𝑓,𝑗 − 𝑃̅𝑚,𝑗
) , 𝑖 = 1,… , 𝑛𝑐
, (5.3)
while the traditional transfer function is expressed as follows:
𝛤𝑖𝑗 = 𝜎𝜌𝑗
𝑘𝑘𝑟𝑗
𝜇𝑗
𝑧𝑖𝑗 (𝑃𝑓,𝑗 − 𝑃̅𝑚,𝑗
) , 𝑖 = 1,… , 𝑛𝑐
. (5.4)
The overall transfer term, 𝛤𝑖
, per component, is evaluated by summing over the relevant phases. In this
chapter, we focus on primary recovery from tight fractured condensate reservoirs and, will accordingly,
ignore the effects of gravity drainage on the mass transfer between matrix and fracture domains.

62
5.2. Numerical Examples
5.2.1. Validation of the Dual-Porosity Compositional Model in MRST
We start by validating the new DP compositional model against an analytical solution to the pressure
diffusion equation and use the gVer transfer function that has been demonstrated to accurately represent
the pressure response.
We consider a 10-ft long single-block DPM saturated with a single-phase fluid. For further
reference, we also construct a 1D SPM with 100 x 1 x 1 cells, including grid refinement near the fractures,
to capture the transient in the matrix (left panel of Fig. 5-1).
Figure 5-1: Schematic of the SPM where the red boundaries represent fractures and yellow blocks represent the matrix (left
panel). Right panel reports the comparison between the pressure response of the analytical solution (Crank, 1975), the fine-grid
SPM, and the DPM (gVer).
We model primary production of pure decane (nC10) (PR EOS (1978) properties reported in Table 5-1) at a
temperature of 125 F (325 K) from a 500 nD rock with a 20% porosity. The initial model pressure is set to
7,000 psi, and the fluid is produced at a constant fracture pressure of 1,500 psi. The relevant shape factor
and 𝑉 parameters are reported in Table 3-1.

63
Table 5-1: EOS parameters of nC10 used in Example 5.2.1
Component 𝑃𝑐
(psi) 𝑇𝑐
(F/K) 𝜔 𝑀𝑤 (lb/mol) 𝑧feed (mol%)
nC10 305 652.19/617.70 0.4884 0.3137 100
𝑃𝑐
= Critical Pressure, 𝑇𝑐
= Critical Temperature, 𝜔 = Acentric Factor, 𝑀𝑤 = Molecular weight
The right panel of Fig. 5-1 reports the dimensionless pressure versus dimensionless time for the three
solutions (analytical, SPM and DPM). We observe a good agreement between all calculations supporting
the accuracy of the implemented DP compositional model in MRST.
5.2.2. Multicomponent Dry-Gas System
Prior to the discussion of gas-condensate problems, we assess the performance of the traditional transfer
function, that is typically available in commercial simulators, for modeling of multicomponent single-phase
gas production. We consider the same SPM grid as discussed in Example 5.2.1, and an equivalent singleblock DPM using one of two transfer functions: 1) Warren and Root (1963) transfer function with a shape
factor from Lim and Aziz (1995); and 2) Generalized Vermeulen transfer function with parameters 𝜎 and 𝑉
from Table 3-1.
For the 10-ft domain, the calculated shape factors from Lim and Aziz and gVer are 9.87e-2 ft-2 and
8.76e-2 ft-2
, respectively. The parameter 𝑉 of gVer is set to 1.54. The matrix porosity is 20%, while we
consider two values of the matrix permeability: 1 D and 5 nD. We model the production of a
multicomponent dry gas with a molar composition of 94% CH4, 4.7% C2H6, 0.8% n-Propane (C3H8), and
0.5% n-Butane (nC4) at a temperature of 170 F (350 K). Table 5-2 summarizes the input parameters for the
EOS calculations with the quaternary fluid system used in this example. The initial pressure is set to 7,000
psi and the gas is produced at a constant pressure boundary condition of 1,500 psi.

64
Table 5-2: EOS parameters of the quaternary dry-gas fluid description used in Example 5.2.2
Component 𝑃𝑐
(psi) 𝑇𝑐
(F/K) 𝜔 𝑀𝑤 (lb/mol) 𝑧feed (mol%)
CH4 667 -116.6/190.6 0.0114 0.0354 94.0
C2H6 707 89.90/305.3 0.0990 0.0663 4.70
C3H8 616 206.2/369.9 0.1521 0.0972 0.80
nC4 550 305.5/425.1 0.2010 0.1281 0.50
Figure 5-2 reports the comparison of the average pressure from the three models (reference SPM and two
DPMs). In both scenarios (1 D in the left panel and 5 nD in the right panel), the traditional transfer
function with a shape factor from Lim and Aziz (1995) departs from the reference solution at early-time,
with early-time ranging from 2 days to 2 years depending on the matrix permeability: The characteristic
time for mass transfer is inversely proportional to the rock permeability.
Figure 5-2: Average matrix pressure for 1 D permeability (left panel) and 5 nD permeability (right panel)
In contrast, we observe that gVer agrees well with the SPM reference solution for both matrix permeability
settings for all times. This is attributed to the time-dependent correction factor, despite the use of a lower
value of the shape factor relative to that of Lim and Aziz (1995).

65
5.2.3. 4-Component Gas-Condensate System
Next, we investigate the accuracy of DPM calculations in settings where two phases may form during the
depletion process. More specifically, we examine the Kazemi-type transfer function to describe the twophase mass transfer problem in a tight fractured condensate system. This allows for testing of the gVer
transfer function in the context of two-phase systems, recalling that the multiphase formulation was
initially extended in a manner similar to that of Kazemi et al. (1976), as seen previously in Eq. (5.3).
We compare calculation results from a single-block DPM to a fine-grid SPM (reference solution),
using the domain of Example 5.2.1, with a matrix permeability of 5 nD and a porosity of 20%. We use a 4-
component gas-condensate fluid description with a molar composition of 0.01% CO2, 87.49% CH4, 10%
nC4, and 2.5% nC10 (the relevant fluid parameters are reported in Table 5-3).
Table 5-3: EOS parameters of the quaternary condensate fluid description used in Example 5.2.3
Component 𝑃𝑐
(psi) 𝑇𝑐
(F/K) 𝜔 𝑀𝑤 (lb/mol) 𝑧feed (mol%)
CO2 1070 87.70/304.1 0.2240 0.0970 0.01
CH4 667 -116.6/190.6 0.0114 0.0354 87.49
nC4 550 305.5/425.1 0.2010 0.1281 10.0
nC10 305 652.2/617.7 0.4880 0.3137 2.50
The model temperature is set at 125 F (325 K). We utilize quadratic gas/liquid relative permeability
functions in the matrix with a residual liquid saturation of 20%, and straight-line relative permeability
functions in the fractures.
The condensate phase envelope, the matrix gas/liquid relative permeability curves, and the liquid
dropout curve from a constant composition expansion (CCE) are shown in Fig. 5-3. At the given
temperature, the fluid has a saturation pressure of 3,252 psi and exhibits a maximum liquid dropout of
10.5% (relative to the saturated volume). With an initial model pressure set at 7,000 psi, we start at singlephase gas conditions. The fracture pressure is then reduced and maintained at 1,500 psi to initiate
production.

66
Figure 5-3: Phase envelope for the quaternary fluid (left panel), gas/liquid relative permeability curves (middle panel), and CCE
liquid dropout (right panel)
Calculation results from the three models in terms of the average matrix pressure and the cumulative
gas/liquid production are reported in Fig. 5-4.
Figure 5-4: Average matrix pressure (left panel), cumulative gas production (middle), and cumulative liquid production (right).
Note: All production streams are converted to surface condition
We observe a substantial difference in the surface liquid production calculated by both DPMs, relative to
the reference solution: After 4 years, the cumulative liquid production of both DPMs is ~4 times larger
than that of the SPM. This departure mainly stems from the use of the average matrix pressure to calculate
the fluid state (see further discussion below) that can render predictions of a field performance in the
development of tight fractured gas-condensate resources inaccurate.

67
5.3. An Improved Two-Phase Transfer Function for Condensates
As demonstrated in Example 5.2.3, the multiphase transfer function introduced by Kazemi et al. (1976)
provides limited accuracy in the representation of a condensate problem. To develop an improved transfer
function, we take a closer look at the physics involved and detail the flaw of the Kazemi-type transfer
functions (when applied without the use of pseudofunctions). Following that, we introduce a new
approach to address the limitations without the need for pseudofunctions.
In the reference solution, the initial drawdown causes the pressure at the matrix/fracture interface
to rapidly drop below saturation pressure, forcing the initial reservoir fluid to split into a lighter gas and a
heavier liquid phase. As a result, the gas flow is impaired by the reduced fluid mobility (relative
permeability) and concurrently, a leaner gas is produced to surface. As the depletion progresses, pressure
and saturation gradients develop inside the matrix block, resulting in liquid buildup near the
matrix/fracture interface as demonstrated in Fig. 5-5.
Figure 5-5: Results of the fine-grid SPM: Pressure profiles(left panel), liquid saturation profiles (middle panel), and liquid
saturation distribution after 2 years of production (right panel)
In contrast to the SPM, the DPM averages the entire fluid and transport behavior into a single block
representation, with an average matrix pressure that remains above the saturation pressure and,
accordingly, predicts flow of a single-phase gas. In other words, the transfer function does not consider
the conditions that prevail at the matrix/fracture interface and sustains production of a richer gas. As a
result, the DPM with the traditional multiphase transfer function (Kazemi-type seen in either Eq. (5.3) or

68
Eq. (5.4)) overestimates the liquid production at the surface conditions. Therefore, the objective of the
proposed method is to improve the evaluation of the transfer function by implicitly avoiding the intrinsic
averaging of the DPM approach. This can be realized by performing equilibrium calculation at conditions
that reflect/approximate the pressure variation that exists inside the matrix. Given the composition
gradient, this will allow for modeling of the liquid buildup in the transfer function, which is essential to
overcome the excessive production of (average) rich gas. To achieve that, we need to define a
representative fluid composition and pressure that will serve as input to the calculation of the transfer
function. These new input parameters must reflect the dynamic behavior at and near the matrix/fracture
interface, where the largest gradients in pressure and composition/saturation reside (as shown in Fig. 5-
5).
The traditional mass transfer calculation entails determining the fluid properties and the rock-fluid
interactions (𝜌𝑗
, 𝑧𝑖𝑗, 𝜇𝑗
, and 𝑘𝑟𝑗) that are evaluated via an equilibrium calculation for the average fluid
composition at the average matrix pressure. The proposed method leverages the computational simplicity
of the transfer function of Kazemi et al. (1976) but evaluates the fluid properties and fluid mobility at a
different/representative fluid composition and conditions. In the following, we refer to these as the
transfer pressure, 𝑃transfer, and the transfer composition, 𝑧transfer.
The main contribution of the modified transfer function is the determination of a new set of
transfer pressures and transfer compositions. We start by using a pre-processed transfer pressure (see
more details below) and integrate itin the DPM to perform flash calculations of the transfer composition
(see more details below) during simulation. The results of this equilibrium calculation are then used to
evaluate the transfer function. To illustrate this, we rewrite Eq. (5.3) in the modified form:
𝛤𝑖𝑗 = 𝛽𝜎𝜌𝑗
∗
𝑘𝑘𝑟𝑗
∗
𝜇𝑗
∗ 𝑧𝑖𝑗
∗ (𝑃𝑓,𝑗 − 𝑃̅𝑚,𝑗
) , 𝑖 = 1,… , 𝑛𝑐
, (5.3a)
where the parameters labeled with an ‘∗’ superscript are obtained from the new equilibrium calculation.

69
We aim to determining a transfer composition that describes the mixture that is formed near the
matrix/fracture interface during a depletion process. In this work, we use the following approach to
calculate the interface composition:
𝑁𝑖 = [(1 − 𝑆𝑙
𝑛−1)(𝜌𝑔𝑧𝑖
)
𝑛
+ (𝜌𝑙𝑥𝑖𝑆𝑙
)
𝑛−1] 𝑀𝑤𝑖
⁄ , 𝑖 = 1, … , 𝑛𝑐
, (5.5)
where the superscript 𝑛 denotes the time level, and subscripts 𝑔 and 𝑙 denote the gas and liquid phases.
𝑁𝑖
is the number of moles of component 𝑖 (per unit bulk volume), 𝑆𝑗
is the saturation of phase 𝑗, 𝑥𝑖
is the
composition of the liquid (from the previous time level), and 𝑧𝑖
is the average composition. In Eq. (5.5),
the properties that are evaluated at the current time level describe the average composition, while the
properties that are evaluated at the previous time level (second termin the square bracket) represent the
liquid due to condensation. Any liquid dropout that occurs at a given time level is accordingly carried
forward to the next time level. Finally, the results of Eq. (5.5) are normalized by the total moles of all fluid
(at a given timestep) to obtain the new transfer composition (𝑧transfer).
To investigate the best possible approximation of 𝑃transfer, at which 𝑧transfer is equilibrated, we
start by matching the liquid production of the DPM to that of the reference solution: This approach
controls/models the amount of liquid dropout in order to mimic the behavior of the reference solution. At
this stage, we have an estimate of the transfer pressure versus time that provides for the best possible
agreement between the DPM and the SPM, subject to the assumptions imbedded in Eq. (5.5). We start
with this definition of 𝑃transfer to evaluate the performance of the modified approach to modeling
matrix/fracture transfer, before providing a discussion and example of possible simplifications that
eliminates the need for matching the DPM to a reference solution.
A comparison of the workflow used to evaluate the traditional and the modified transfer function
is provided in Fig. 5-6. As discussed earlier, the traditional approach uses the average fluid composition
and flashes it at the average matrix pressure. The resulting parameters from the equilibrium calculation
are then used to evaluate the mass transfer rate.

70
The proposed method however (right panel of Fig. 5-6), requires additional input in terms of a
transfer pressure and a transfer composition that are either pre-processed or calculated during DPM
simulation. In the first timestep, we utilize the transfer pressure to flash the average block composition
and use the equilibrium state to obtain the parameters required for evaluation of the gas and liquid mass
transfer rates. Subsequent timesteps entail the calculation of a transfer composition from Eq. (5.5) that is
then equilibrated at 𝑃transfer of the respective timestep to determine the mass transfer rates.
Figure 5-6: Flowchart comparing the two ways to evaluate the transfer function. Note: “phase equilibrium” represents all
calculations that are required to evaluate the transfer function parameters
5.3.1. 4-Component Gas-Condensate System – Revisited
To test/demonstrate the performance of the proposed method, we revisit Example 5.2.3. A comparison of
the reference solution (SPM) with the calculations from the DPM using the unmodified as well as the
modified multiphase version of the gVer transfer function (see Eqs. (5.3) and (5.3a)) is provided in Fig. 5-
7. We observe that the results of the new approach substantially agree with the fine-grid SPM solution.

71
The bottom right panel reports the liquid saturation at the interface for the three calculation
approaches. We observe that the SPM predicts a liquid buildup that reaches ~35% at early time, in contrast
to the unmodified transfer function (blue circles) that predicts formation of two phases at a later time
(almost after a year of production) and never reaches the critical saturation of 20%. However, the liquid
saturation of the new method (purple crosses), resulting from flashing the transfer composition, and
therefore representative of the conditions near the matrix/fracture interface, improves the representation
of the liquid buildup within the matrix block.
Figure 5-7: Plots in order (top row left to right): Gas production rate, cumulative gas production, average pressure, liquid
production rate, cumulative liquid production, and liquid saturation at the matrix/fracture interface
We further report on the final average matrix composition, for all models, in Table 5-4 and by comparing
their phase envelopes in Fig. 5-8. It can be observed that the DPM, with the unmodified transfer function,
departs from the new method which, in turn, agrees well with the reference solution.

72
Table 5-4: Final average matrix composition for all models for Example 5.3.1.
Approach Composition (mol%)
CO2 CH4 nC4 nC10
SPM 0.01 82.06 12.73 5.20
DPM: gVer Original (Eq. 5.3) 0.01 85.39 11.00 3.60
DPM: gVer New Method (Eq. 5.3a) 0.01 82.83 12.24 4.92
Figure 5-8: Phase envelopes generated from the final average matrix composition of the respective models
5.3.2. 24-Component Gas-Condensate
To investigate how the modified gVer transfer function performs for a realistic fluid system, we report on
a test for a 24-component fluid description from a wet gas reservoir. The reader is referred to Whitson and
Kuntadi (2005) for additional details on the compositional description.
We consider again the SPM and DPM domains from Examples 5.2.2 and 5.2.3. The gas is described
by a condensate-gas-ratio of 33.5 STB/MMSCF and a molar composition of 83% CH4, 1.98% C7+, while the
non-hydrocarbon content includes N2, CO2, and hydrogen sulfide. The initial reservoir pressure is set at
5,315 psi and a temperature of 220 F (377 K). At this temperature, the fluid has a saturation pressure of
5,135 psi and the model is initialized at a single-phase gas condition. Production is then initiated with a
constant pressure boundary condition of 1,500 psi.

73
The calculation results are presented for four models: 1) reference SPM; and 2) DPM with the
traditional transfer function using a shape factor from Lim and Aziz; 3) DPM with the unmodified gVer
transfer function; and 4) DPM using the modified gVer transfer function (new approach). Figure 5-9 reports
the average pressure and the gas/liquid volumetrics. In addition, the cumulative production by component
is provided in Figs. 5-10 (first 12 lighter components) and 5-11 (the remaining 12 heavier components).
Figure 5-9: Results are in the order of: Gas production rate, cumulative gas production, average matrix pressure, liquid
production rate, and cumulative liquid production

74
Figure 5-10: Comparison of the component cumulative production for the first 12 lighter components for all four models:
reference SPM (red solid line), DPM using unmodified LA (green dashed line), DPM using unmodified gVer (blue circles), and DPM
using modified gVer (purple stars)
Figure 5-11: Comparison of the component cumulative production for the last 12 heavier components for all four models:
reference SPM (red solid line), DPM using unmodified LA (green dashed line), DPM using unmodified gVer (blue circles), and DPM
using modified gVer (purple stars)

75
As observed in Fig. 5.9, the DPMs with the unmodified transfer functions depart from the reference
solution, particularly in terms of the cumulative liquid production. In contrast, the modified DPM
demonstrates a significant improvement in the forecasting of liquid production. Furthermore, from Figs.
5-10 and 5-11, it can be noticed that the lighter components have a good agreement between all the
approaches. However, the modified approach provides for a significant improvement in the prediction of
the heavier components where it matches the reference solution (for all components).
5.4. Discussion
As highlighted in the condensate examples above, the classical DPM approach results in a delayed twophase effect and extended production of rich gas. To overcome this limitation and enable capturing the
phase transition near the matrix/fracture interface, a transfer pressure and a transfer composition were
introduced and used to evaluate the modified transfer function. This method retains the simple DPM
representation of the Warren and Root (1963) model, as well as the Kazemi et al. (1976) multiphase
formulation, providing for an easy implementation in reservoir simulators.
So far, the improved representation of the mass transfer entails the pre-processing of the transfer
pressure by matching the DPM to a reference solution. As indicated earlier, the general idea of matching
the recovery of a DPM to a fine-grid model has frequently been used in the development of DPMs.
However, a reference solution may not always be available. Therefore, we seek an alternative approach
that eliminates the need for the evaluation/matching of a reference solution and aim to develop a
simpler/adaptive approach for the estimation of the transfer pressure.
To guide the development of an adaptive approach for evaluating 𝑃transfer, we first review the
average matrix pressure and the transfer pressure of Examples 5.3.1 and 5.3.2 (see Fig. 5-12).

76
Figure 5-12: Average matrix pressure, transfer pressure obtained by matching the DPM to a reference solution, and transfer
pressure obtained by an analytical solution of Examples 5.3.1 (left panel) and 5.3.2 (right panel)
To generate a similar transfer pressure, we propose approximating the local pressure of the system by
solving the analytical solution of the pressure diffusion equation. From Crank (1975), we have the
following solution (for the current geometry):
𝑃𝑚−𝑃0
𝑃𝑓−𝑃0
= 1 −
4
𝜋
∑
(−1)
𝑘
2𝑘+1
exp{−𝐷
(2𝑘+1)
2𝜋
2
𝑡
4𝑙
2
} cos (
(2𝑘+1)𝜋𝑥
2𝑙
)
∞
𝑘=0
, (5.6)
where 𝑃𝑚 is the matrix pressure at a given distance 𝑥 from the center of the block, and 𝑙 is the half-width
of the fracture spacing. To approximate the local pressure, we consider the flow of pure CH4 and evaluate
the pressure at a location equal to 95% of the matrix length (𝑥 = 0.95 × 𝑙 in Eq. (5.6)).
The estimated transfer pressure from Eq. (5.6) is compared to the one obtained from matching
the DPM to the reference solution in Fig. 5-12. While not identical, both curves mimic the expected
pressure behavior near the matrix/fracture interface. We test the simple/adaptive approach for evaluation
of the transfer pressure next.

77
5.5. Application of the Adaptive Transfer Pressure – Examples 5.3.1 and 5.3.2 Revisited
In Examples 5.3.1 and 5.3.2 discussed above, the transfer pressure was obtained by matching the liquid
recovery of the respective DPM to its corresponding fine-grid SPM reference solution. Here, we evaluate
the performance of the modified transfer function when 𝑃transfer is approximated from an analytical
solution. To this end, Eq. (5.6) is applied and incorporated into the DPM. The updated results (average
pressure and gas/liquid production) of Examples 5.3.1 and 5.3.2 are reported in Figs. 5-13 and 5-14. We
observe that the new method, using the analytical solution to estimate the transfer pressure, provides for
a good agreement with the fine-grid SPM.
Figure 5-13: Updated results of Example 5.3.1 with the adaptive transfer pressure that is generated using the analytical solution

78
Figure 5-14: Updated results of Example 5.3.2 with the adaptive transfer pressure that is generated using the analytical solution
5.6. Summary and Conclusion
In this chapter, we propose a novel method that improves the modeling and simulation of mass transfer
in tight fractured gas-condensate reservoirs via a DPM representation. The new method tracks an
equilibrium calculation at the matrix/fracture interface to capture the rapid changes of the fluid
behavior/composition and integrate representative fluid properties and fluid mobilities in the transfer
function. As a result, the approach overcomes the lag in the liquid dropout, observed in the traditional
DPM formulation, and provides for improved production forecasts.
Key to the modified transfer function is the calculation of the transfer composition and transfer
pressure. The transfer composition is estimated by accounting for the average composition and the liquid
dropout observed at/near the interface. Furthermore, depending on the resources available, one option
to obtain the transfer pressure is to match the recovery (liquid production) of the DPM to a reference
solution. A more general (automated) workflow approximates the transfer pressure by means of an
analytical solution to the pressure diffusion equation.

79
To demonstrate the performance of the new method, examples of gas-condensate (a 4-
component analog fluid and a realistic 24-component fluid) production were presented and the results
show an excellent agreement with the fine-grid SPM. While the use of an analytical solution to
approximate the transfer pressure is less accurate than when the one obtained from matching the DPM
to the reference solution, it renders the overall approach more suitable for larger models where a
reference solution may be very costly to evaluate.

80
Chapter 6. Diffusive Mass Transfer and Recovery During CO2 Huff-n-Puff
in an Eagle Ford Shale Core*
EOR is a promising approach to unlock the untapped potential of unconventional reservoirs. One of the
emerging interests is CO2 injection into tight fracturedformationsto stimulate additional production, and
to raise the recovery factor above what is achievable via primary production. Lab experiments and
observations from pilot tests have demonstrated the potential of CO2 EOR in shale in terms of hydrocarbon
recovery (e.g., Darvish et al., 2006; Hoffman, 2018; Alahmari et al., 2023). However, a systematic
framework for modeling and simulation of such settings, to support design of optimal field development
strategies,is still lacking. Accordingly, a further improvement of our understanding of the physics at lab-,
pilot-, and field-scales is critical to guide industrial efforts. Reservoir modeling and simulation commonly
serve as a link between lab and pilot scales by providing a path to upscale relevant processes to a field
level, as well as provide means to justify the cost of operations.
A key enabler to define a robust/physics-based model is the understanding of the recovery
mechanisms during CO2 injection into tight fractured reservoirs, which is different than our previous
observations from conventional reservoirs. Generally, advective flow, molecular diffusion, and gravity
drainage are the three main transport mechanisms. In higher permeability conventional reservoirs,
advective flow and gravity drainage often dictate the recovery of fluids. Molecular diffusion becomes
relevant in unconventional formations where a mobile fluid (CO2 in the fracture) contactsan oil phase that
is largely stagnant in the shale matrix. CO2 then permeates from the fracture into the matrix, via diffusive
mass transfer, causing the oil to swell and start mobilizing. Operation of EOR in unconventional formations
differ therefore from CO2 flooding, used in conventional reservoirs, and is often performed via a cyclic
injection/soak/production process referred to as huff-and-puff (HnP).
*Joint work with Saeed Alahmari and Ye Lyu which is currently in preparation for publication

81
In this chapter, we present an integrated experimental and modeling study of CO2 HnP in an Eagle
Ford shale core. The objective is to demonstrate/delineate the importance of diffusive mass transfer on
recovery in fractured media during EOR operations and its implementation in simulation tools. The
experimental component of this work was performed by Saeed Alahmari and Ye Lyu and is summarized
here to contextualize the presentation/discussion of our modeling/simulation efforts. The reader is
referred to Alahmari (2024) for additional information regarding the experimental work.
6.1. Lab Experimental Materials
We start by presenting the experimental work including a description of the core, the fluid system, and
the preparation of the core.
We used a shale core from the Eagle Ford outcrop (see Fig. 6-1) with a diameter of 1.00 in., a
length of 3.278 in, and a dry core weight of 91.90 g. Based on initial visuals of the EF-3 core, laminations
are observed across the core entire length (more on this later).
Figure 6-1: Eagle Ford shale core used in the HnP experiment
A synthetic oil consisting of decane (nC10), dodecane (nC12), tetradecane (nC14), and hexadecane (nC16)was
used in the HnP experiment. The molecular weights and initial compositions of the fluid system are
reported in Table 6-1.

82
Table 6-1: Molecular weight and composition of the synthetic oil
Component 𝑀𝑤, g/mole Composition, mole %
nC10 142.3 7.9
nC12 170.3 19.9
nC14 198.4 30.8
nC16 226.5 41.4
The core was evacuated at 230 F (383.15 K) for approximately three weeks, during which the change in
weight was monitored. When the observed change was less than 0.01 g, the dry core weight was recorded.
The core wasthen saturated in a high-pressure vessel for 3 months at 3,000 psi and 122 F (323.15 K). The
saturated weight of the core was then recorded to provide an estimate of the porosity.
6.2. Lab Experimental Procedure
In this section, we present the experimental setup and procedures for the HnP, pressure-pulse decay (PPD)
and gas flow-through measurements.
6.2.1. Huff-n-Puff
For the HnP experiment, we used the same setup and procedure as was detailed in Alahmari et al. (2023).
We conducted four HnP cycles on the EF-3 core where each cycle starts with CO2 injection followed by a
soaking period and then pressure depletion while temperature was kept constant at 158 F (343.15 K).
Table 6-2 reports the injection pressure for each cycle. We maintained a soaking time of approximately
100 hours for all cycles. The pressure was gradually decreased during the depletion in three stages until
reaching a final pressure of 14.7 psi.
Table 6-2: HnP injection pressure per cycle
Cycle # Injection pressure, psi
1 4,160
2 4,400
3 3,990
4 3,900

83
We planned to maintain an injection pressure of 4,000 psi for each cycle and were able to achieve that
within +/- 200 psi, except for the 2nd cycle where we had some technical issues to stay within the desired
range. It is worth noting that the selection of the injection pressure was made to operate under firstcontact-miscible (FCM) conditions during all cycles. Although the calculated FCM pressure (FCMP) for CO2-
synthetic oil, at 158 F (343.15 K), is around 2,600 psi, it was reported by Alahmari et al. (2023) that a
change in miscibility conditions from one cycle to another is expected. Therefore, to avoid this, we
assumed an extreme scenario where the EF-3 core is fully saturated with pure nC16 and calculated the
FCMP accordingly. For the given scenario, the FCMP was calculated to 3,700 psi; hence, we expect all the
HnP cycles to be within FCM conditions.
6.2.2. Pressure-Pulse Decay and Gas Flow-Through Experiments
To provide for a more detailed characterization of the core permeability, we performed PPD and gas flowthrough measurements using helium (He) to delineate the effects of the observed laminations (seen in
Fig. 6-1 above).
As for the PPD measurements, the setup used is shown in Fig. 6-2 and consists of a gas source, a
hydraulic pump to control the confining pressure, a constant temperature oven, two cells (upstream and
downstream), a core holder, a vacuum pump, two pressure gauges and a data acquisition system.

84
Figure 6-2: Pressure-pulse decay and gas flow-through experimental setup
We used the same setup for the gas flow-through measurements with the addition of mass flow meters
placed after the downstream cell to measure the flow rate. It is worth noting that both PPD and gas flowthrough measurements were conducted on the EF-3 core after the HnP experiments. Refer to Appendix-A
for more details on the PPD measurements procedure and analysis to calculate the core permeability.
6.3. Experimental Observations and Evaluation
In this section, we present the experimental observations and related analysis.
6.3.1. Core Porosity and Permeability Estimates
The core porosity was estimated from the weight of the core pre- and post- saturation (gravimetrical
approach). The dry core weight is 91.90 g, the saturated weight is 95.64 g, and the saturating fluid density
is 0.7616 g/cc resulting in a pore volume of 4.9 cc and a porosity of 11.9%. The permeability (2.5 D) was
measured using the Gas Research Institute (GRI) technique (Guidry et al., 1996) on a ground sample
obtained from a different plug that was cut from the same core.
Furthermore, a total of six PPD measurements were performed and analyzed to confirm the
Downstream Cell

85
permeability measurements. The corresponding upstream and downstream pressures that were used in
the PPD measurements are listed in Table 6-3.
Table 6-3: List of pressure pulse experiments
Experiment # 𝑃1
, psia 𝑃2
, psia
1 500 0
2 485 435
3 498 450
4 550 493
5 615 555
6 670 608
We present an example of the experimental data for the 2
nd pulse in Fig. 6-3 (left panel) where the blue
circles represent the upstream pressure and the green circles correspond to the downstream pressure.
The data from each pulse was interpreted by plotting ln(𝑝𝐷
) versus time (see example with 2nd pulse in in
the right panel of Fig. 6-3) to evaluate the matrix permeability by the method of Brace et al. (1968), Dicker
and Smits (1988), and Joens (1997). The resulting values of permeability, from all pulses, are plotted versus
the inverse of average core pressure, and a corrected permeability of 2.79 μD is estimated. Compared to
the GRI method, the two estimates of the matrix permeability are consistent.
Figure 6-3: Upstream and downstream pressure data from the 2
nd pulse (left) and its interpreted ln(pD) vs. time plot (right)

86
Figure 6-4: Calculated permeability from the pressure-pulse decay measurements
After obtaining an estimate of the matrix permeability from PPD, gas flow-through measurements along
the core’s bedding plane were performed to analyze the effects of lamination. A total of four experiments
were performed where the results from one experiment are reported in Fig. 6-5: The black symbols
represent the upstream pressure; the blue symbols depict the downstream pressure; and the red symbols
show the estimated permeability.
The permeability estimates of the four experiments are shown in Fig. 6-6 where the core
permeability is plotted versus the inverse of average core pressure. We observe that the permeability in
the axial direction is approximately 3.88 D. In comparison to the PPD and GRI estimates (2.79/2.50 D),
that primarily probe the matrix permeability, the difference indicates an increased connectivity in the core
along the bedding plane (laminations) by a factor of 1.39.
y = 685.19x + 2.7926
R² = 0.9978
1
2
3
4
5
0 0.0005 0.001 0.0015 0.002 0.0025 0.003
Permeability (uD)
1/Pressure (psia-1)

87
Figure 6-5: Permeability estimated from gas flow-through measurements: black symbols represent the upstream pressure, blue
symbols depict the downstream pressure, and red symbols show the estimated permeability
Figure 6-6: Calculated permeability from the gas flow-through measurements
6.3.2. Recovery, Produced Composition, and Pressure Behavior
The pressure and temperature profiles, as observed during the four HnP cycles, are reported in Fig. 6-7:
The pressure is shown by the blue symbols while the green symbols represent the temperature.

88
Figure 6-7: Pressure and temperature readings from the HnP experiment: Pressure is shown by blue symbols while green symbols
represent the temperature
We start analyzing the HnP experimental results by evaluating the oil recovery, followed by the produced
oil composition and finally the pressure behavior during the soaking phase.
Oil recovery (by mass) was evaluated from a material balance calculation using measurements of
the core weight after each cycle (see Fig. 6-8). From Fig. 6-8, we observe that the largest incremental oil
recovery was during the 1st cycle (58.2%) with an addition of 19.5% during the 2nd and nearly complete
recovery after the 3rd cycle (96%). The additional recovery from the 4th cycle was only ~2%.

89
Figure 6-8: Experimental oil recovery at the end of each of the four HnP cycles
The produced oil composition, reported alongside the initial synthetic oil composition in Fig. 6-9 indicates
that diffusion plays a central role in the recovery process: We observe from the 1st cycle that the produced
oil composition is very close to the synthetic oil composition. After the 1st cycle, we notice that the mole
fractions of nC10 and nC12 continuously decrease while nC14 and nC16 continue to increase. Furthermore,
we do not observe any nC10 production from the 4th cycle. Given that the estimated additional recovery
from the 4th cycle is comparatively small (~2%), we used the available GC compositional data of the
produced liquids to estimate the remaining oil composition residing in the core after the 3rd cycle: The
combined mole fraction of nC14 and nC16 remaining in the core is approximately 0.89 representing a
significant change relative to the original synthetic oil (0.72).

90
Figure 6-9: Measured produced oil composition for the four synthetic oil components
Further, the estimated oil composition left in the core after each cycle is reported in Table 6-4: The oil
remaining in the core becomes heavier after each cycle which highlights the role of diffusion that
contributes to producing lighter component first and heavier components in subsequent cycles.
Table 6-4: Oil remaining (in mole fraction) in the core after each HnP cycle
Post 1st Cycle Post 2nd Cycle Post 3rd Cycle Post 4th Cycle
nC10 0.0708 0.0598 0.0428 0.0119
nC12 0.1914 0.1798 0.1584 0.1016
nC14 0.3102 0.3152 0.3320 0.3560
nC16 0.4276 0.4452 0.4669 0.5306
Next, we explore the pressure variation during the soaking phase to identify the relevant modes of mass
transfer. The pressure-drop vs time (and square root of time) during all soaking stages are reported in Fig.
6-10. During the injection phase, the oil initially in the core is compressed from atmospheric pressure to
the injection pressure. From the calculations with bulk fluid properties, this compression amounts to
approximately 3% of the fluid volume. Subsequently, during the soaking phase, CO2 invades the core
further and causes the extraction cell pressure to decrease over time. The pressure decline should
eventually stop upon reaching an equilibrium state (had the soaking period been long enough).

91
Figure 6-10: Observed pressure drop during the soaking period for the four HnP cycles

92
For the first HnP cycle (upper row of right panel of Fig. 6-10), we observe a linear relationship between
the pressure-drop and √𝑡, suggesting a diffusion-dominated process. We note a deviation from linearity
at early times, a behavior that has been discussed by many researchers (Reamer et al., 1956; Caskey and
Michelsen, 1973; Chukwuma, 1983; Renner, 1988; Tan and Thorpe, 1992): e.g., Chukwuma attributed this
to an “incubation” period for equilibrium to be established at the gas-oil contact area. We also observe
from Fig. 6-10 (upper row of the left panel) that the pressure during soaking does not equilibrate during
the experiment.
Subsequent cycles (2-4) depart from an exponential pressure decay with respect to time,
indicating the presence of a minor leak and/or the impact of adsorption of CO2 in the core. The potential
for adsorption is discussed in more detail in Alahmari (2024). While any leak from the extraction cell would
invalidate the interpretation of the pressure responses, the overall pressure drop during soaking is
relatively small (< 60 psi) and is not expected to affect the extraction process. Accordingly, we restrict
further analysis of the pressure decline to consider only the first cycle.
To further investigate the pressure variation during the soaking of the 1st cycle, we compare the
total pressure-drop to the equilibrium pressure obtained from a Volume-Temperature (VT) flash
calculation (equilibrium calculation for a fixed 𝑉, 𝑇 and molar feed composition). We used the PR EOS
(1978), with tuned binary interaction coefficients and volume shift parameters (further is discussed in
Section 6.6.1), to calculate the equilibrium pressure of a mixture of CO2 and synthetic oil that corresponds
to the moles of initial oil in place, and the moles of CO2 injected into the cell. The CO2 moles were estimated
using CO2 density (at pressure and temperature conditions relevant to the HnP) and the total system
volume (volume of the extraction cell minus the core’s bulk volume). Recalling that the final pressure
recorded at the end of soaking during the 1st cycle is 4,107 psi, the equilibrium pressure from the VT flash
calculation is 4,150 psi: The difference might be attributed to other mechanisms such as CO2 adsorption

93
that causes the pressure during soaking to drop below the calculated equilibrium pressure (when sorption
is not included).
6.4. Numerical Modeling
We utilize MRST for the numerical modeling studies of this chapter. The base code of MRST’s compositional
simulation, which supports convective flow via an SPM representation, was previously validated against
commercial and research simulators (Moyner, 2021). Building on the existing MRST modules, in this
chapter, we develop a molecular diffusion model to facilitate studies relevant of mass transfer during HnP
operations in tight fractured reservoirs.
As highlighted earlier, it is common to use commercial tools to perform calculations towards the
interpretation of HnP experiments. The open-source environment of MRST is the novelty of this research
that allows us to construct the necessary simulation tool that provides the flexibility to develop fit-forpurpose algorithms and investigate the complex characteristics of unconventional resources. The central
advantage of using MRST is realized when we allow for modifications in two main aspects of the modeling
of diffusive mass transfer including: 1) evaluation of multicomponent diffusion coefficients; and 2)
selection of different representations of the diffusive fluxes.
In all simulation examples, we used the PR EOS (1978) to compute equilibrium conditions, and the
LBC correlation (Lohrenz et al., 1964) for calculation of phase viscosities. Next, we present the
mathematical model, validation examples, and finally our modeling efforts for interpreting the EF-3 core
HnP experiment.
6.4.1. Mathematical Model
The transport equations for 𝑛𝑐
components, existing in both oil and gas phases, can be written in a set of
multicomponent mass-balance equations:
𝜕
𝜕𝑡
(∑𝑗𝜙𝜌𝑗𝑆𝑗𝑧𝑖𝑗) ⏟
accumulation
+ ∇ ∙ (∑ 𝜌𝑗𝑧𝑖𝑗𝑢⃗ ⏟𝑗 𝑗
Convection
+ ∑ 𝐽 ⏟𝑗 𝑖,𝑗
diffusion
) = 0, 𝑖 = 1,… , 𝑛𝑐
, (6.1)

94
where 𝑡 denotes time and 𝜙 representsthe porosity. The parameters 𝜌𝑗
, 𝑆𝑗
, and 𝑢𝑗
represent the fluid
molar density, saturation, and the velocity of phase 𝑗, respectively. The parameter 𝑧𝑖𝑗 representsthe mole
fraction of component 𝑖 in phase 𝑗.
The second term in Eq. (6.1) represents the total flux due to convection and diffusion: The former
is described by Darcy’s law with an upstream weighting method; and 𝐽𝑖,𝑗
is the molecular diffusion flux of
component 𝑖 in phase 𝑗 (more on this next).
6.4.2. Modeling of Molecular Diffusion
The diffusive molar flux of component 𝑖 in a bulk phase, 𝐽𝑖
, is commonly described using one of two
common molecular diffusion models. The first is the generalized Fick’s law:
𝐽𝑖 = −𝜌 ∑ 𝐷𝑖𝑘
𝑛𝑐−1
𝑘=1
∇𝑧𝑘
, 𝑖 = 1,… , 𝑛𝑐 − 1 , (6.2)
and the second is the classical Fick’s law:
𝐽𝑖 = −𝜌𝐷𝑖
eff∇𝑧𝑖
, 𝑖 = 1,… , 𝑛𝑐
. (6.3)
We note that 𝐷 in Eq. (6.2) are the generalized Fickian diffusion coefficients represented by a
(𝑛𝑐 − 1 × 𝑛𝑐 − 1) matrix. To obtain 𝐷, we follow the Maxwell-Stefan (MS) framework, where the MS
binary diffusivities are related to the infinite dilution coefficients (𝐷
∞) using, for example, the correlation
of Kooijman and Taylor (1991). Common correlationsto estimate the infinite dilution coefficients include
that of Hayduk and Minhas (1982) or Leahy-Dios and Firoozabadi (2007). Finally, to convert the MS
diffusivities to the generalized Fickian diffusion coefficients, the following relationship is employed:
[𝐷] = [𝐵]
−1[Γ] , (6.4)
where 𝐵 is a function of the MS diffusivities and Γ is the thermodynamic matrix that accounts for nonideal behavior. Refer to Appendix-B for more details on the MS framework used in this chapter.
In the classical Fick’s law (see Eq. (6.3)), 𝐷𝑖
eff is a column vector with elements representing the 𝑛𝑐
effective diffusion coefficient in the mixture. To evaluate the binary diffusion coefficients, we use the
extended Sigmund correlation (Sigmund, 1976a, 1976b; da Silva and Belery, 1989). Following that, the

95
formulation of Wilke (1950) is applied to obtain the effective diffusivity per component in the mixture. The
details of this workflow are provided in Appendix-C.
Finally, to express the diffusive flux in porous media, we multiply Eqs. (6.2) and (6.3) by the ratio
𝜙𝑆𝑗/𝜏, where 𝜏 is the tortuosity factor.
6.5. Validation of The Molecular Diffusion Models
In this section, we present two examples to validate the implementation of the diffusion models (the
classical Fick’s law and the generalized Fick’s law) in MRST.
6.5.1. 1D Analytical Solution
A synthetic example is constructed to validate the diffusive flux model in MRST by comparing the
calculated response to an analytical solution (Li, 1972). The problem considers a 1D linear diffusion process
of a binary fluid containing CO2 and CH4 in a domain with length, 𝐿, of 3.28 ft (1 m). The domain is
discretized by 20 equally spaced grid blocks, with the left 10 cells initially containing CH4 and the right 10
cells initially containing CO2 (left panel of Fig. 6-11). This createsa sharp composition change at the center
of the domain where diffusion is expected to initiate in response to a constant and mutual diffusion
coefficient of 1.5e-5 ft2/s (1.5e-6 m2/s). The system is initialized at standard conditions of 14.7 psi (101
kPa) and 60 F (288 K).
The diffusion process can be simulated using the classical Fick’s law assuming that the diffusion
coefficient remains constant throughout. As indicated in the right panel of Fig. 6-11, which reportsthe CO2
mole fraction along the domain after 0.2 days, the simulation model agrees well with the analytical
solution.

96
Figure 6-11: 1D domain showing the initial composition of the binary system (left panel) and CO2 mole fraction (after 0.2 day)
along the domain for the analytical and numerical solutions(right panel)
6.5.2. Lab Experimental Data of Arnold and Toor (1967)
The objective here is to validate our implementation of the generalized Fick’s law, via the SPM, using data
from the lab experiment of Arnold and Toor (1967). The study involves diffusion of a ternary mixture
composed of CH4, argon (Ar), and hydrogen (H2) in a Loschmidt tube with initial compositions shown in
Fig. 6-12 (left panel). The pressure and temperature of the setup were 14.7 psi (101 kPa) and 93 F (307 K).
An important observation made by the experimentalists is that the diffusion process did not follow a
Fickian behavior. Accordingly, the dragging effect is significant and the application of the generalized Fick’s
law is warranted.
We construct a 1D grid with a total length of 2.66 ft (0.81 m) of 20 equally spaced grid blocks. The
grid is shown in the left panel of Fig. 6-12 where the upper 10 grid blocks representing the upper bulb of
the Loschmidt tube, and the bottom 10 grid blocks are used to represent the lower bulb. According to the
given conditions, single-phase gas diffusion is expected to initiate at the middle of the domain.

97
Figure 6-12: Schematic of Arnold and Toor (1967) experimental setup (left panel). Average CH4 and Ar compositions versus time
in the middle and right panels, respectively. The solutions (reported in both panels) are in the order of lab data (circles), MRST
(solid lines), Hoteit (2013) in dashed lines, and Taylor and Krishna (1993) in doted lines.
Figure 6-12 reports the average composition, in each bulb, of CH4 (middle panel) and Ar (right panel) for
three sets of calculations: Our work with MRST, Hoteit (2013), and Taylor and Krishna (1993). Starting with
Taylor and Krishna (1993), represented by dotted lines, they approached the multicomponent problem
analytically by solving the coupled diffusion equations following the method of Toor (1964), which
assumes a constant diffusivity (refer to Sherafati and Jessen (2018) for more details on this solution).
Hoteit and our work in MRST rely on numerical simulation to find the concentration variation of
this multicomponent problem, where the generalized Fick’s law is used in both solutions. In order to obtain
a match with the experimental data, Hoteit applied 2 multipliers to tune the calculated diffusion
coefficients. For the same objective, and following a similar approach of Hoteit’s, we multiply the diagonal
elements by 1.5 and the off-diagonal elements by 2.52. As can be seen from the calculation results, all
three solution approaches provide a reasonable agreement with the experimental data.
6.6. Modeling of CO2 HnP on the EF-3 Core
In this section, we present our efforts to interpret the experimental observations for the EF-3 core via
compositional simulation. We compare MRST calculations with the lab data in terms of oil recovery and

98
produced oil composition (for all cycles), while modeling of the pressure drop during soaking period is
reported in the discussion section.
6.6.1. Model Setup
A compositional SPM was built (in MRST) using a radial grid to represent the dimensions of the HnP
experimental setup (see Fig. 6-13).
Figure 6-13: Numerical grid that mimics the EF-3 HnP experiment (in MRST): 3D view of the extraction cell and the core (left
panel), 3D view of the core (middle panel), and horizontal cross-section from the top of the core (right panel)
The extraction cell, mimicking the behavior of a larger fracture, is represented by the yellow outer grid
blocks, and the core is shown with red grid blocks. Grid refinement, in radial and vertical directions, was
performed to reduce numerical diffusion. The porosity of the extraction cell was set to 99% to mimic the
void space and a ‘large’ value of permeability equal to 0.500 Darcy was assigned. The matrix porosity
(11.9%) and matrix permeability (𝑘𝑥 = 𝑘𝑦 = 2.79 μD and 𝑘𝑧 = 3.88 μD) were assigned according to the
data from the core characterization (presented in previous sections). As discussed in section 6.4.2, a
tortuosity parameteris required for the description of diffusion in the porous material. A tortuosity value

99
of 6 was used here based on literature review for the Eagle Ford formation (Zhang et al., 2019; Davudov
et al., 2020).
The system was initialized at 14.7 psi and a temperature of 343.15 K. The dimensionless volume
shift parameters, 𝑆𝑖
, and the binary interaction coefficients were obtained by matching the relevant
density and solubility data from the existing literature. Such fluid analysis was previously presented by
Alahmari et al. (2023) and the EOS input parameters are summarized in Table 6-5 (for 𝑃𝑐
, 𝑇𝑐
, 𝜔, and 𝑆𝑖
)
and Table 6-6 (for binary interaction).
Table 6-5: EOS parameters of the synthetic oil and CO2
Component 𝑃𝑐
(psi) 𝑇𝑐
(K) 𝜔 𝑆𝑖
(dimensionless)
CO2 1069.87 304.2 0.23 0.12279
nC10 305.68 617.6 0.49 0.06478
nC12 264.53 658.2 0.56 0.10308
nC14 235.14 694.0 0.68 0.15413
nC16 205.74 717.0 0.74 0.17669
Table 6-6: Binary interaction parameters of the oil binaries with CO2
Binary Binary Interaction
CO2-nC10 0.09694
CO2-nC12 0.08947
CO2-nC14 0.08777
CO2-nC16 0.07676
Initially, the core is saturated with 100% synthetic oil and the annulus with 100% CO2. Further, injection
and production are represented by two wells completed at the top and bottom of the (extraction cell) grid,
respectively. Quadratic gas/liquid relative permeability functions were assumed for the matrix and
straight-line relative permeability functions for the fracture (more on this later).
6.6.2. Representation of the Diffusive Flux
As the role of diffusion is at the center of our study, we focus on proper representation of the diffusive
mass transfer during the HnP experiment. As such, we examine the performance of the two main diffusion

100
models: 1) The generalizedFick’s law (see Eq. 6.2),and 2) the classical Fick’s law (see Eq. 6.3). In both cases,
the diffusion coefficients are calculated as a function of PVT conditions and compositions (refer to
Appendices-B and C for additional details).
6.6.3. Oil Recovery- Prediction
Starting with the oil recovery, Fig. 6-14 reports the experimental data versus the simulation results: Blue
circles represent the lab data, the solid black line indicate simulation without diffusion, the dotted yellow
line denote calculation with the classical Fick’s law (Sigmund, 1976; Wilke, 1950), and the purple dashed
line illustrates the results using the generalized Fick’s law (Kooijman and Taylor, 1991; Leahy-Dios and
Firoozabadi, 2007). We observe that the difference between calculated (10%) and experimental recovery
(60%) amounts to a factor of 6 when diffusion is not included in the calculations and only convective flow
is activated, indicating the importance of diffusion during the recovery process. Figure 6-14 also reveals
that for both diffusion models, diffusion contributes substantially to the recovery of oil, however, both
models are also observed to overpredict the experimental data.
Figure 6-14: Comparison of oil recovery from HnP experimental data and simulation results: Lab data (blue circles), simulation
with no diffusion (solid black line), simulation using classical Fick’s law (dotted yellow line), and simulation using generalized
Fick’s law (dashed purple line)

101
From such initial calculation results, we find that the current (default) methods/models do not provide an
accurate representation of the impact of molecular diffusion during hydrocarbon recovery. To this end, we
present and examine three approaches to address the disagreement between the calculated and observed
experimental data.
6.6.4. Improved Representation of Diffusive Mass Transfer
To ascertain the sensitivity of the diffusion models parameters(coefficients), we consider three different
approaches and report their performance relative to the experimental observations.
It is common to use multipliers to minimize errors that may arise in the default diffusion models,
for example,Hoteit (2013) applied two multipliers (for diagonal- and off-diagonal diffusion coefficients) to
match the results from the generalized Fick’s law to the observations of Arnold and Toor (1967), see
Example 6.5.2. Similarly, Shi et al. (2022) demonstrated through experimental and modeling work that a
single correction factor is required for each binary when applying Sigmund’s correlation for (Sigmund,
1976) to model CH4 diffusing in hydrocarbon mixtures.
In our first approach to improve the agreement between the calculations and the lab data (see
section 6.6.3), we focus on the application of the generalized Fick’s law, where the calculated infinite
dilution coefficients are adjusted based on available measurementsfrom the existing literature. As such,
we rely on data sets from Cadogan et al. (2016) for diffusion coefficients pertaining to CO2 infinitely diluted
in hydrocarbons, and data from Umezawa & Nagashima (1992) for hydrocarbons diluted in CO2, as they
represent measurements conducted relatively close to the conditions of our HnP experiment.
The procedure first entails interpolating/extrapolating the infinite dilution coefficients from the
relevant studies to 158 F (343.15 K) and 4,000 psi, which then serve as the reference values for the
respective CO2-hydrocarbon binaries. These reference values are then compared to the calculated infinite
dilution coefficients (at the same conditions) from the available correlations to obtain multipliers for each
binary.

102
The left panel of Fig. 6-15 reports the comparison of the infinite dilution coefficients for CO2 in
nC10, nC12, and nC16 obtained from interpolating the experimental data of Cadogan et al. (2016), which
serve asreference values, with the calculated values obtained from five popular correlations, including:
A. Wilke and Chang (1955)
B. Hayduk and Minhas (1982)
C. Siddiqi and Lucas (1986)
D. Wong and Hayduk (1990)
E. Leahy-Dios and Firoozabadi (2007)
We note that nC14 was not included in the experimental work of Cadogan et al. (2016). Therefore, we
developed a correlation between the infinite dilution coefficient and the carbon number using data from
Cadogan et al. (2016) and then interpolated the infinite dilution coefficientsfor CO2 in nC14 as depicted in
Fig. 6-15 (right panel). Among the applied correlations, we observe that Wong and Hayduk (1990) and
Hayduk and Minhas (1982) provide the best prediction of the experimental data with average errors of
3.87% and 10.69%, respectively.
Figure 6-15: Comparison between the interpolated infinite dilution coefficients for CO2 in oil and the calculated values from each
correlation (left panel). Interpolated values as a function of carbon number (right panel)

103
We attribute the accuracy of the prediction from the correlation of Wong and Hayduk (1990) to the fact
that their correlation was developed using experimental data for CO2 dissolved in organic solvents
including normal alkanes. Similarly, Hayduk and Minhas (1982) utilized experimental data for normal
alkanes as solvents in developing the correlation. The other correlations either did not use experimental
data for CO2 in normal alkanes or used various types of mixture and consequently deviate from specific
mixtures such as CO2 in normal alkanes.
A similar analysis was performed for the other concentration endpoint, i.e., hydrocarbon infinitely
diluted in CO2, with data obtained from Uzemawa and Nagashima (1992), see Fig. 6-16. However, for this
data set, extrapolation to our HnP experiment was required given that 𝐷
∞ were obtained at lower 𝑇, 𝑝.
We note from Fig. 6-16 that all correlations deviate more than 40% from the reference values,withHayduk
and Minhas (1982) providing the best prediction: We attribute this to the fact that their correlation was
developed using experimental data with normal alkanes as solutes.
Figure 6-16: Experimental infinite dilution diffusion coefficients for oil in CO2 and the calculated values from each correlation

104
Based on the above analysis for the two infinite dilution limits (per binary), we note that the correlation
of Hayduk and Minhas (1982) provides the best overall prediction of the experimental data. Therefore, we
proceeded to evaluate multipliers that correct for the infinite dilution coefficients, computed using the
correlation of Hayduk and Minhas (1982), to the experimental reference values (see Table 6-7).
Table 6-7: Multipliers to correct for the infinite dilution diffusion coefficients obtained from Hayduk and Minhas (1982)
nC10 nC12 nC14 nC16
CO2 in Oil 0.881 0.896 0.916 0.935
Oil in CO2 2.611 3.097 3.223 3.588
At this point, we have determined correction factors that can be applied during numerical calculations
when the correlation of Hayduk and Minhas (1982) is used to provide infinite dilution diffusion coefficients
in the generalized Fickian diffusion model. We refer to this approach as the “Infinite Dilution Multiplier”
approach. Figure 6-17 comparesthe experimental oil recovery from EF-3 with the calculations including:
1) no diffusion, 2) the unmodified diffusion model; and 3) the infinite dilution multiplier approach. We
observe that application of the correction factors does not provide for a better agreement between
calculation and the experimental data, with calculated recovery departing further from the experimental
observations than in the previous calculations approach.

105
Figure 6-17: Comparison of oil recovery after applying the infinite dilution coefficients multipliers: Lab data (blue circles),
simulation with no diffusion (black solid line), simulation unmodified diffusion model (dotted yellow line), and simulation when
applying the infinite dilution multiplier (dashed purple line)
To further investigate the gap between the experimental data and the modeling results (presented in
section 6.6.3), we examine two modifications to the formulations of the diffusive mass transfer (see
section 6.4.2) that we are applying in our HnP modeling effort.
The first approach considers the uncertainty in the calculation of the multicomponent diffusion
coefficients. It is well known that the measurement of multicomponent diffusion coefficients is
challenging, particularly at relevant reservoir conditions: Most published experiments are limited to binary
and ternary mixtures at lower temperatures and pressures. The second approach explores the role of the
anisotropy in the core fabric, as evident from the observed laminations in Fig. 6-1, on the modeling of the
HnP process. The heterogeneity of the core likely resultsin a characteristic diffusion time that depends on
the direction, which is not representedinthe model, thus far. Based on the depositional environment (and
related lamination), it is reasonable to assume that diffusive mass transfer across the laminations will be
slower than transport along the laminations. Accordingly, we expect that the overall diffusive flux in the
axial direction will exceed the diffusion flux in the radial direction.

106
To address such limitations in the diffusion models, a single multiplier was applied to reduce the
overall flux into the core. This can be done by one of two approaches: 1) A common correction factor that
is applied to the calculated diffusion coefficients (labeled ‘Diffusion Case 1’ in Fig. 6-18); and 2) A common
multiplier is applied to the radial diffusion flux (‘Diffusion Case 2’). To obtain a good agreement with the
experimental recovery data, multipliers of 0.35 (Case 1) and 0.25 (Case 2) were applied in calculations with
the classical Fick’s law, while multipliers of 0.10 (Case 1) and 0.03 (Case 2) were used in calculations with
the generalized Fick’s law (𝐷
∞ from Hayduk and Minhas (1982)). The calculated recovery obtained by this
approach is compared to the experimental observations in Fig. 6-18. We observe that both approaches
provide for a good agreement with the experimental observations.
Figure 6-18: Comparison of oil recovery from experiment and simulation: The classical Fick’s law (left panel) and the generalized
Fick’s law (right panel). Blue circles represent lab data, dotted-dashed purple line represents ‘Diffusion Case – 1,’ and dashed red
line represents ‘Diffusion Case – 2’
Next, we compare the produced oil composition (per cycle) from the lab data and modeling results in Fig.
6-19. From the figure, we observe that for the simulation results where diffusion is neglected, the
compositions remain generally unchanged throughout the cycles due to the coherence of species in
transport via advection.

107
Figure 6-19: Comparison of produced oil composition from experiment and simulation. The first column reports simulation
results where diffusion is neglected, columns 2-3 report calculations with the classical Fick’s law, while columns 4-5 provide
results from the generalized Fick’s law
On the other hand, for the remaining four columns shown in Fig. 6-19, we find that the general trend in
the variation of produced compositions between HnP cycles is captured well by the calculations. More
importantly, we observe that both modeling approaches, generalized or classical Ficks law, provide a very
good agreement with the measurements. Furthermore, the application of multipliers (overall flux or in the
radial direction only), largely provides for the same result. Finally, we find that the classical Fick’s law
provides an improved accuracy in the prediction of produced nC12 and nC16 compared to the generalized
Fick’s law.

108
6.6.5. Discussion
It is evident from the results and analysis provided above, that diffusion plays a central role in the oil
recovery during cyclic CO2 HnP process. To quantify the contribution of the diffusion mechanism, it was
demonstrated how the recovery differs by a factor of 6 between a base case calculations (with no diffusion)
and the experimental observations (also from the modified diffusion modeling). In fact, the production
period (puff) during the HnP process contributes less than 10% to the recovery during the first cycle, while
the remaining is attributed to the soaking period, and accordingly to mass transfer driven by diffusion.
Figure 6-20 illustrate how CO2 invadesthe core by means of molecular diffusion during the soaking stage.
Figure 6-20: Vertical cross-section of the core reporting the CO2 mole fraction during soaking (first cycle). Note: The snapshots
shown here are based on calculations with the classical Fick’s law –Diffusion Case 1
For the calculation results presented above, quadratic gas/liquid relative permeability functions (𝑛 = 2)
were assumed in the core due to a lack of relative permeability data for shales in the existing literature. To
investigate the role of relative permeability functions, calculations were performed also with straight-line

109
relative permeability functions (𝑛 = 1) assigned to the core and a modest difference of only 1.5% in
recovery after the completion of the first cycle is observed (see Figs. 6-21). This modest variation in the
calculated recovery further illustrates the dominant role of diffusion in the recovery process (50% recovery
during soaking at FCM conditions versus an additional 10% during the production phase where relative
permeability functions can become important below the FCMP).
Figure 6-21: Comparison of oil recovery calculated with different relative permeability functions: quadratic (blue columns) and
straight-lines (orange columns). Note: The results reported here are based on calculations with the classical Fick’s law –
‘Diffusion Case 1’
The impact of relative permeability on the produced oil composition is illustrated in Figure 6-22. Similar to
the observations made for the calculated oil recovery, we also observe a very modest impact on the
calculated oil composition as the relative permeability functions are changed from quadratic to linear.

110
Figure 6-22: Comparison of produced oil composition with different relative permeability functions: quadratic (blue columns) and
straight-lines (orange columns). Note: The results reported here are based on calculations with the classical Fick’s law –
‘Diffusion Case 1’
During the construction of the numerical model, all available core and fluid characterization information
was included. While directional permeability was observed from PPD experiments, no information is
readily available to infer the directional aspect of diffusional mass transfer. To evaluate the impact of a
directional diffusion flux, a multiplier was applied that allowed us to match the experimental recovery
data. We considered two cases, one where the overall flux was reduced by a common multiplier and one
where the radial flux was adjusted by a single multiplier. While the latter approach may be more
reasonable given the visible lamination of the core, both cases provide for an excellent agreement with
the recovery (matched) and the produced compositions (predicted) across all cycles. We note that analysis
of the produced compositions was not included in the selection of multipliers. In both cases, a multiplier
(less than 1) reducesthe impact of the diffusion process and therefore resultsin a reduced recovery.

111
Contrary to the findings of some researchers (e.g., Hoteit, 2013; Sofla et al., 2016; Sistan et al.,
2022), we observe that the classical Fick’s law provides for an accurate calculation of the produced
compositions and amounts. In fact, the classical Fick’s model performs slightly better than the generalized
model. This is promising given the favorable (reduced) computational requirement associated with this
approach compared to the generalized Fick’s law, which entails additional layers of calculations.
As part of this work, the correlation of Leahy-Dios and Firoozabadi (2007) was also investigated in
place of the correlation of Hayduk and Minhas (1982). Both correlations provide for similar trends in the
calculated results. However, the multipliers required to match the model using the correlation of LeahyDios and Firoozabadi to the lab data were less aggressive (multipliers = 18% and 10% for Diffusion Cases 1
and 2, respectively).
The focus of the modeling efforts, presented above, was on the evaluation of recovery factors and
produced compositions. Based on the calibrated model, we subsequently explored the calculated pressure
drop as observed during the soaking period. A comparison of the calculated response and the lab data for
the first cycle of HnP is reported in Fig. 6-23.
Figure 6-23: Pressure variation during soaking period of the first cycle (experimental vs simulation)

112
The simulation results are consistent with conclusions made from the results presented in Figs. 6-17 and
6-18, i.e., the largerthe impact of the diffusion, the faster the pressure drop during soaking, and the higher
the recoveryper HnP cycle. Two additional observations canbe made from Fig. 6-23. The calculation with
the unmodified version of the generalized Fick’s law stabilizes at a pressure of ~4,150 psi. As highlighted
earlier, this represents the equilibrium pressure of the system. It is also noticed how the simulation
calculations underestimate the pressure drop relative to the experimental data. This departure could be
attributed to the adsorption of the CO2 molecules onto the shale core, which was not considered in this
simulation study. Additional work is underway to further delineate the impacts of adsorption during HnP
operations.
6.7. Summary, Conclusion, and Future Recommendations
In the previous sections, we presented an integrated experimental and modeling study to investigate the
role of diffusive mass transfer during CO2 HnP in an Eagle Ford shale core. The experiment wasperformed
at 3,900-4,400 psi and 158 F (343.15 K) on a core saturated with a well-characterized normal alkane
mixture. The system pressure, recovery, and produced oil compositions were collected and interpreted.
Furthermore, PPD and gas flow-through experiments were conducted, following the 4th cycle, that allowed
for an estimation of critical information including the effective core permeability and porosity.
Then, an open-source simulation tool (MRST) was modified (to include the diffusion models) and
applied to interpret the experimental observations. Validation examples were initially presented followed
by modeling of the EF-3 core HnP experiment. Two diffusion models were used to calculate the diffusion
fluxes: The classical Fick’s law and the generalized Fick’s law, with composition- and pressure-dependent
diffusion coefficients (isothermal experiments).

113
Key findings and recommendations of this chapter can be summarized as follows:
1) The interpretation of the data collected from the EF-3 core HnP experiment indicate that diffusion
plays a central role in the oil recovery: A larger fraction of the lighter components is produced
initially, while larger fractions of the heavier components are produced during the later cycles.
Furthermore, the pressure behavior during the soaking phase (of the first cycle) exhibits a linear
relationship with respect to √𝑡 in support of a diffusion-dominated process.
2) The mainstream diffusion models cannot predict the experimental observations without
adjustments to represent the physics of the diffusion process. As such, the open-source code was
applied and modified to address the uncertainty of the calculated diffusion coefficients and to
account for the directionality of the diffusion flux. Themodifications result in a good agreement with
the recovery and produced compositions as observed from the experiment.
3) The pressure drop during the soaking period, calculated from the calibrated model, does not agree
well with the experimental observations and further model adjustments were not attempted. Work
is currently underway to incorporate an adsorption model into the open-source MRST simulator to
further explore the role of adsorption during HnP processes.

114
6.8. Appendix
Appendix-A: PPD Measurements Procedure and Analysis
In this section, we report the PPD measurements procedure as well as the associated analysis to calculate
the core permeability.
1. The core was initially wrapped in shrink-tubing and then placed in the biaxial core holder (as
shown in the bottom left corner of Fig. 6-2).
2. The oven temperature wasset at 122 F (323.15 K) and the system is monitored until thermal
equilibrium was established.
3. A confining pressure was applied on the system and was evacuated for 2 days.
4. After the vacuuming process was completed, the reference cell wasisolated from the core holder
and charged with gas at a selected pressure.
5. Gas was allowed to flow from the reference cell (upstream) through the core to the downstream
cell while the cell pressures (𝑃1 and 𝑃2
) were recorded.
A total of six PPD measurements were performed and the observed pressure decay was used to
evaluate the core permeability via a modified version of the method developed by Brace et al. (1968) and
later extended by Dicker and Smits (1988) and Jones (1997): The gas, helium in our measurements, in the
pore space of the core is initially at the same pressure as in the downstream vessel (vacuum for first pulse).
The gas pressure in the upstream vessel is set slightly higher than the core pressure and the magnitude of
the pressure in the upstream vessel is set slightly higher than the core pressure and the magnitude of the
pressure in the upstream cell should not exceed 10% of the downstream pressure (Brace et al., 1968).
Upon opening the valve connecting the upstream cell to the core holder, the upstream pressure decays as
the gas flows from the upstream cell through the core and into the downstream cell. The observed
pressure response versus time is then used to calculate the core permeability as detailed below.

115
A general solution for the dimensionless pressure difference, as a function of the dimensionless
time, during a PPD experiment was presented by Dicker and Smits (1988):
∆𝑃𝐷 = 2 ∑ exp(−𝜃𝑚
2
𝑡𝐷
)
𝑎(𝑏
2+𝜃𝑚
2 )—(−1)𝑚𝑏√(𝑎2+𝜃𝑚
2 )(𝑏
2+𝜃𝑚
2 )
[𝜃𝑚
4 +(𝜃𝑚
2 (𝑎+𝑎2+𝑏+𝑏
2) +𝑎𝑏(𝑎+𝑏+𝑎𝑏)]
∞
𝑚=1
, (A.1)
where 𝜃𝑚 represents the roots of Eq. (A.2):
tan𝜃 =
(𝑎+𝑏)𝜃
𝜃2−𝑎𝑏
, (A.2)
where 𝑎 and 𝑏 are the ratios of the core pore volume to upstream and downstream vessels, respectively:
𝑎 =
𝑃𝑉
𝑉1
, (A.3)
𝑏 =
𝑃𝑉
𝑉2
, (A.4)
where ∆𝑃𝐷 is the differential dimensionless pressure for gases:
∆𝑃𝐷 =
𝑃1
(𝑡𝐷)
2−𝑃2
(𝑡𝐷)
2
𝑃1
(0)2−𝑃2
(0)2
. (A.5)
𝑡𝐷 is the dimensionless time:
𝑡𝐷 =
6.805𝑒−5 𝑘 𝑡
𝜇 𝜙 𝐿
2
. (A.6)
The subscripts 1 and 2 denote upstream and downstream, respectively. The parameter 𝑃𝑉 is the
core pore volume, 𝑉 is the vessel volume, 𝑃 is the pressure, 𝑡 is the time, 𝑘 is the core permeability, 𝜙 is
the core porosity, and 𝐿 is the core length.
If early time dynamics are ignored, Eq. (A.1) reduces to the below equation:
ln[∆𝑃𝐷
] = ln[𝑓0
] − [
6.805𝑒−5 𝑘 𝜃1
2
𝜇𝑔𝑐𝑔 𝜙 𝐿
2
]𝑡 , (A.7)
where 𝜇𝑔 is the gas viscosity, 𝑐𝑔 is the gas compressibility, 𝜃1
2
can be approximated as shown in Eq. (A.8)
(Dicker and Smits, 1988), while 𝑓0
is obtained from Eq. (A.9):
𝜃1
2 = (𝑎 + 𝑏 + 𝑎𝑏) −
1
3
(𝑎 + 𝑏 + 0.4132𝑎𝑏)
2 + 0.0744 (𝑎 + 𝑏 + 0.0578𝑎𝑏)
3
, (A-8)
𝑓0 =
2[𝑎(𝑏
2+𝜃1
2 )+ 𝑏√(𝑎2+𝜃1
2 )(𝑏
2+𝜃1
2)]
𝜃1
2 (𝜃1
2+𝑎+𝑎2+𝑏+𝑏
2)+𝑎𝑏(𝑎+𝑏+𝑎𝑏)
. (A.9)

116
The core compressibility is negligible when compared to the gas compressibility; hence, it was ignored in
the interpretation. We can also see from Eq. (A.7) why the pressure pulse must be small: Gas viscosity and
compressibility vary with pressure and if the upstream and downstream pressures were substantially
different, at time zero, then Eq. (A.7) would provide a coarse approximation.
Appendix-B: Calculation of the Generalized Fickian Diffusion Coefficients (Generalized Fick’s law)
Provided in this appendix the framework used to calculate the generalized Fickian diffusion coefficients
used to describe the flux in the generalized Fick’s law.
The generalizedFickian diffusion coefficients, denoted as𝐷, are obtained from MS diffusivities via
the following relationship:
[𝐷] = [𝐵]
−1[Γ] , (B.1)
where
𝐵𝑖𝑗 = {
𝑥𝑖
Ð𝑖𝑛𝑐
+ ∑
𝑥𝑘
Ð𝑖𝑘
𝑛𝑐
𝑘=1
𝑘≠𝑖
𝑖 = 𝑗
−𝑥𝑖
(
1
Ð𝑖𝑗
−
1
Ð𝑖𝑛𝑐
) 𝑖 ≠ 𝑗
, (B.2)
and
𝛤𝑖𝑗 = 𝑥𝑖
𝜕 ln(𝑓𝑖
)
𝜕𝑥𝑗
, (B.3)
where 𝐵 is a function of the MS diffusivities (denoted as Ð), Γ is the thermodynamic matrix that accounts
for non-ideal behavior, and 𝑓 is the fugacity. The fugacity is obtained using the EOS calculations.
In this work, The MS diffusivities for multicomponent mixtures are estimatedusing the generalized
Vignes equation proposed by Kooijman and Taylor (1991) that are based on the infinite dilution diffusion
coefficients as such:

117
Ð𝑖𝑗 = (𝐷𝑖𝑗
∞)
𝑥𝑗 (𝐷𝑗𝑖
∞)
𝑥𝑖 ∏ (𝐷𝑖𝑘
∞𝐷𝑗𝑘
∞)
𝑛 𝑥𝑘/2
𝑐
𝑘=1
𝑘≠𝑖,𝑗
, 𝑖,𝑗 = 1, …, 𝑛𝑐
(𝑖 ≠ 𝑗) , (B.4)
where 𝐷𝑖𝑘
∞ is the diffusion coefficient of component 𝑖 infinitely diluted in component 𝑘. The MS coefficients
matrix, Ð𝑖𝑗, is symmetric and its diagonal elements do not exist. Therefore, Eq. (B.4) calculates 𝑛𝑐
(𝑛𝑐 −
1)/2 MS coefficients.
Two correlations, Hayduk and Minhas (1982) and Leahy-Dios and Firoozabadi (2007), were
developed in MRST to calculate the infinite dilution diffusion coefficients. The empirical equation by
Hayduk and Minhas (1982) is expressed as the following:
𝐷𝑖𝑗
∞ = 13.3× 10−8𝑇
1.47𝜇𝑗
(10.2/𝑉𝑏𝑖−0.791)
𝑉𝑏𝑖
−0.71
, 𝑖,𝑗 = 1, …, 𝑛𝑐
(𝑖 ≠ 𝑗) , (B.5)
where 𝑇 is the temperature and 𝑉𝑏𝑖 is the partial molar volume of component 𝑖 at the boiling point and
is obtained from:
𝑉𝑏𝑖 = 0.285𝑉𝑐𝑖
1.048
. (B.6)
In Eq. (B.6), 𝑉𝑐𝑖 is the critical volume of component 𝑖.
The correlation of Leahy-Dios and Firoozabadi (2007)relates the infinite dilutions, at a given 𝑝 and
𝑇, to the diffusivity at low pressure, 𝐷𝑖𝑗
𝑜
:
𝜌𝐷21
∞
(𝜌𝐷)0 = 𝐴0
(
𝑇𝑟,1𝑃𝑟,2
𝑇𝑟,2𝑃𝑟,1
)
𝐴1
(
𝜇
𝜇0
)
[𝐴2
(𝜔1,𝜔2
)+𝐴3
(𝑃𝑟,𝑇𝑟
)]
. (B.7)
The constants 𝐴0
to 𝐴3 are given by:
𝐴0 = 𝑒
𝑎1
𝐴1 = 10𝑎2
𝐴2 = 𝑎3
(1 + 10𝜔1 − 𝜔2 + 10𝜔1𝜔2
)
𝐴3 = 𝑎4
(𝑃𝑟,1
3𝑎5 + 6𝑃𝑟,2
𝑎5 + 6𝑃𝑟,1
10𝑎6 )+ 𝑎7𝑃𝑟,2
−𝑎5 + 𝑎2
(
𝑇𝑟,1𝑃𝑟,2
𝑇𝑟,2𝑃𝑟,1
)
, (B.8)

118
where 𝑇𝑟,𝑖 𝑎𝑛𝑑 𝑃𝑟,𝑖 are the reduced temperature and pressure of component 𝑖 and defined as 𝑇𝑟,𝑖 =
𝑇 𝑇𝑐,𝑖
⁄ and 𝑃𝑟,𝑖 = 𝑃 𝑃𝑐,𝑖
⁄ , respectively. The parameter 𝜔𝑖
is the acentric factor of component 𝑖. Finally, the
constant are given as: 𝑎1 = −0.0472, 𝑎2 = 0.0103, 𝑎3 = −0.0147, 𝑎4 = −0.0053, 𝑎5 = −0.3370, 𝑎6 =
−0.1852, and 𝑎7 = −0.1914.
The binary diffusion coefficient at low pressure, (𝜌𝐷)
0
, is evaluated using Fuller et al. (1969):
(𝜌𝐷)
0 = 1.01 × 10−2𝑇
0.75
(
1
𝑀𝑤1
+
1
𝑀𝑤2
)
0.5
𝑅[(∑𝑣1
)1⁄3+(∑ 𝑣2
)1⁄3]
2 . (B.9)
In Eq. (B.9), 𝑀𝑤 is the molecular weight, 𝑅 is the universal gas constant, ∑ 𝑣𝑖
is the diffusion volume
increments of component 𝑖 which is calculated by summing the atomic diffusion volumes (see Poling et
al., 2001). The low-pressure viscosity for a given component 𝑖 in the mixture (see Eq. (B.7)) is obtained
using the correlation by Stiel and Thodos (1961):
𝜇𝑖
0𝜉𝑖 = 34 × 10−8(𝑇𝑟,𝑖
)
0.94
for 𝑇𝑟,𝑖 < 1.5 , (B.10)
𝜇𝑖
0𝜉𝑖 = 17.78 × 10−8(4.58𝑇𝑟,𝑖 − 1.67)
5/8
for 𝑇𝑟,𝑖 > 1.5 , (B.11)
where
𝜉𝑖 =
𝑇𝑐,𝑖
1/6
𝑀𝑤𝑖
1 /2
(0.987×10−5𝑃𝑐,𝑖
)
2/3 , (B.12)
Finally, the dilute gas viscosity of the mixture is computed using a weighted average of the dilute
gas viscosity of the components:
𝜇
0 =
𝜇1
0𝑀𝑤1
1/2
+𝜇2
0𝑀𝑤2
1/2
𝑀𝑤1
1/2
+𝑀𝑤2
1/2 , (B.13)

119
Appendix-C: Calculation of the Effective Diffusion Coefficients (classical Fick’s law)
Calculation of effective diffusion coefficients, denoted as 𝐷𝑖
eff
, are required to define the diffusive flux that
is described by the classical Fick’s law.
The correlation of Sigmund (1976a, 1976b) for a binary system is given by:
𝐷𝑖𝑗 =
(𝜌𝑚𝐷𝑖𝑗 )
𝑜
𝜌𝑚
× (0.99589 + 0.096016𝜌𝑚𝑟 − 0.22035𝜌𝑚𝑟
2 + 0.032874𝜌𝑚𝑟
3 ) , (C.1)
where 𝜌𝑚𝑟 is the reduced mixture density defined by:
𝜌𝑚𝑟 = 𝜌𝑚
∑ 𝑧𝑖𝑉𝑐𝑖
𝑛 2/3
𝑖=1
∑ 𝑧𝑖𝑉𝑐𝑖
𝑛 5/3
𝑖=1
. (C.2)
In Eq. (C.1), the low-pressure limit of the density-diffusivity product is given by the Hirschfelder et
al. equation (Reid et al., 1987):
(𝜌𝑚𝐷𝑖𝑗 )
𝑜
=
2.2648𝑥10−5
𝜎𝑖𝑗
2Ω𝑖𝑗
(
1
𝑀𝑤𝑖
+
1
𝑀𝑤𝑗
)
1/2
𝑇
1/2
, (C.3)
In the above equation, 𝜎𝑖𝑗 is the collision diameter and Ω𝑖𝑗 is the collision integral and are related
to the component critical properties through the following equations:
𝜎𝑖 = 0.1866𝑉𝑐𝑖
1/3
𝑍𝑐𝑖
−6/5
, (C.4)
𝜀𝑖 = 65.3𝑇𝑐𝑖𝑍𝑐𝑖
18/5
, (C.5)
𝜎𝑖𝑗 = 0.5(𝜎𝑖 + 𝜎𝑗
) , (C.6)
𝜀𝑖𝑗 = (𝜀𝑖𝜀𝑗
)
1/2
, (C.7)
𝑇𝑖𝑗 = 𝑇/𝜀𝑖𝑗 , (C.8)
Ω𝑖𝑗 =
1.06036
𝑇𝑖𝑗
0.15610 + 0.193exp(−0.4635𝑇𝑖𝑗)+ 1.03587exp(−1.52996𝑇𝑖𝑗) + 1.76474exp(−3.89411𝑇𝑖𝑗) , (C.9)

120
where 𝜀𝑖𝑗 is the minimum energy of attraction.
The effective diffusion coefficient of component 𝑖 in the mixture is:
𝐷𝑖
eff =
1−𝑧𝑖
∑ 𝑧𝑗𝐷𝑖𝑗
𝑛𝑐 −1
𝑗=1
𝑗≠𝑖
, 𝑖 = 1,… , 𝑛𝑐
. (C.9)
da Silva and Belery (1989) indicated that the correlation of Sigmund (1976) (see Eq. (C.1)) is not applicable
to very dense gases and liquid systems. Therefore, Eq. (C.1) is assigned for fluid mixtures of 𝜌𝑚𝑟 ≤ 3. For
mixtures of 𝜌𝑚𝑟 > 3, da Silva and Belery (1989) proposed the following:
𝐷𝑖𝑗 =
(𝜌𝑚𝐷𝑖𝑗 )
𝑜
𝜌𝑚
× 0.18839exp (3− 𝜌𝑚𝑟) . (C.10)

121
Chapter 7. Summary of Dissertation and Future Work
7.1. Summary of Dissertation
We began by revisiting the classical formulation of Warren and Root (1963) and developed/introduced an
improved transfer function termed the generalized Vermeulen transfer function. The advanced modeling
of the matrix/fracture mass transfer includes a time-dependent correction factor that provides for
accurate prediction of early-time behavior in tight fractured systems. The new transfer function was
implemented in MRST and validated initially against commercial simulators using single-block calculations.
Further validationsincluded isotropic and anisotropic permeability settingsfor a single-block DPM and an
equivalent fine-grid SPM. We also demonstrated the performance of the generalized Vermeulen for a
depletion process with large variations in the hydraulic diffusivity.
We expanded the application of the new transfer function to include the modeling of complex
fracture settings wherein the coarse DPM was utilized to upscale discrete fracture characterization
represented by an EDFM. For random fracture geometries and petrophysical properties, we demonstrate
the DPM’s ability to reproduce the flow behavior of the reference solution at a lower computational cost.
Advancing our studies to two-phase flow scenarios, we modified the new transfer function to
model tight fractured gas-condensate reservoirs. To overcome the limitation of the traditional DPM, which
incorporates the Kazemi-type multiphase transfer function, we proposed a formulation that aims to
predict the correct fluid state (i.e., liquid dropout) at the matrix/fracture interface,which dominates the
physics of the fluid flow and transport. We documented the accuracy of this method in terms of recovery
and produced compositions (for all phases and per hydrocarbon component) for two condensate fluids: A
4-component analog dry gas and a 24 component realistic wet gas.
We then demonstrated the importance of molecular diffusion by conducting experimental and
modeling studies of HnP on an Eagle Ford shale core (EF-3). The experimental work included HnP, PPD,

122
and gas flow-through measurements. Diffusion was observed to be the dominant recovery mechanism
based on the analysis of the oil recovery, pressure behavior during soaking, and produced oil composition.
In our modeling efforts, we implemented a diffusion model (represented by either classical Fick’s law or
generalized Fick’s law) in MRST, via the SPM approach,to interpret the lab data. We proposed two methods
to adjust the standard diffusion models based on the physics of the problem to be able to match the
calculation results to the experimental data.
7.2. Recommendations for Future Work
To further advance the findings of this dissertation, we propose the following research directions:
1. In Chapter 4, we proposed a workflow to use the generalized Vermeulen transfer function in
modeling of the matrix/fracture mass transfer for a random fractured system. This workflow can
be tested to model more realistic scenarios. In the context of unconventional reservoirs, an
example can be defined for a pressure depletion problem from a multistage hydraulic fracture well
completed in a tight fractured reservoir. In such settings, multiple porosity domains can be
identified and upscaled via the proposed approach. An important step will be to identify
relevant/optimal subsystems to accurately represent how they contribute to the mass transfer
process. This recommendation can, for example, support studies related to dry gas production
from tight fractured formations through complex fracture networks in order to optimize well
stimulation designs.
2. We presented our modeling efforts to interpret the experimental data for the HnP experiment on
the EF-3 core using currently available diffusion models (see Chapter 6). In this experiment, CO2 is
at supercritical conditions and first contact miscible conditions prevail such that single-phase
diffusion is initiated. It is relevant also to study and understand how to represent diffusive mass
transfer below the miscibility pressure where gas/liquid interfaces will develop.

123
3. Finally, in Chapter 6, we highlighted the importance of investigating the impact of adsorption on
the pressure response during HnP operations, and how to integrate this in compositional
simulation. The affinity of the species to adsorb on the rock surfaces may have a significant impact
on the pressure behavior. Accordingly, the integration of sorption behavior may be required to
represent the soaking period of the HnP experiment in related calculations.

124
References
Alahmari, S. 2024. A Study of Diffusive Mass Transfer in Tight Dual-Porosity Systems (Unconventional).
Ph.D. dissertation, University of Southern California, Los Angeles, CA (August 2024).
Alahmari, S., Raslan, M. S., Khodaparast, P., et al. 2023. CO2 Huff-n-Puff: An Experimental and Modeling
Approach to Delineate Mass Transfer and Recovery from Shale Cores. Paper presented at Middle East
Oil, Gas and Geosciences Show, Manama, Bahrain, February 2023. https://doi.org/10.2118/213400-
MS.
Alghalandis, Y. F. 2017. ADFNE: Open Source Software for Discrete Fracture Network Engineering, Two
and Three Dimensional Applications. J. of Computers and Geosciences 102: 1-11.
http://dx.doi.org/10.1016/j.cageo.2017.02.002.
Alharthy, N., Teklu, T., Kazemi, H., et al. 2015. Enhanced Oil Recovery in Liquid-Rich Shale Reservoirs:
Laboratory to Field. SPE Res. Eval. & Eng. 21 (01): 137-159. https://doi.org/10.2118/175034-PA.
Barenblatt G.I., Zheltov, I.P., and Kochina, I.N. 1960. Basic Concepts in the Theory of Seepage of
Homogeneous Liquids in Fissured Rocks. Prikladnaya Matematikai Mekhanika, Akad. Nauk, SSSR, 24
(5): 852-864. https://doi.org/10.1016/0021-8928(60)90107-6.
Brace, W., Walsh, J. B., & Frangos, W. T. 1968. Permeability of granite under high pressure. J. of
Geophysical Research 73 (6): 2225-2236. https://doi.org/10.1029/JB073i006p02225.
Cadogan, S. P., Mistry, B., Wong, Y., Maitland, G. C., Trusler, J. P., 2016. Diffusion Coefficients of Carbon
Dioxide in Eight Hydrocarbon Liquids at Temperatures between (298.15 and 423.15) K at Pressures up
to 69 MPa. J. of Chem and Eng. Data 61 (11): 3922-3932. 10.1021/acs.jced.6b00691.
Caskey, J.A., Michelsen, D.L., and To, Y.P. 1973. The Effect of Surfactant Hydrophilic Group on Gas
Adsorption Rates. J. Colloid and Interface Sci. 42 (1): 62-69. https://doi.org/10.1016/0021-
9797(73)90007-6.
Chang, M. 1993. Deriving the shape factor of a Fractured Rock Matrix. Technical Report NIPER-696
(DE93000170). Bartlesville, Oklahoma: NIPER.
Chukwuma, F. O. 1983. Mass Transfer of Gaseous Carbon Dioxide into a Quiescent Liquid Hydrocarbon.
Ph.D. dissertation, The University of Tulsa, Tulsa, Oklahoma.
Coats, K.H. 1989. Implicit Compositional Simulation of Single-Porosity and Dual-Porosity Reservoirs. Paper
presented at the SPE Symposium on Reservoir Simulation, Houston, TX, February 1989.
https://doi.org/10.2118/18427-MS.
Computer Modelling Group LTD. 2018. Computer Modelling Group-IMEX. https://www.cmgl.ca/imex.
Crank, J. The Mathematics of Diffusion. 2
nd ed. Claredon Press, Oxford: 1975.

125
da Silva, F. V., and Belery, P. 1989. Molecular Diffusion in Naturally Fractured Reservoirs: A Decisive
Recovery Mechanism. Paper presented at SPE Annual Technical Conference and Exhibition, San
Antonio, Texas, October 1989. https://doi.org/10.2118/SPE-19672-MS.
Darvish, G. R., Lindeberg, E., Utne, S. A., et al. 2006. Reservoirs-Conditions Laboratory Experiments of CO2
Injection into Fractured Cores. Paper was prepared for presentation at SPE Europec/EAGE Annual
Conference and Exhibition held in Vienna, Austria, June 2006. https://doi.org/10.2118/99650-MS.
Davudov, D., Moghanloo, R. G., & Zhang, Y. 2020. Interplay between pore connectivity and permeability
in shale sample. International J. of Coal Geology 220 (1): 103427.
https://doi.org/10.1016/j.coal.2020.103427.
Dean, R.H., Lo, L.L. 1988. Simulations of Naturally Fractured Reservoirs. SPE Res. Eng. 3 (2): 638-648.
https://doi.org/10.2118/14110-PA.
Dershowitz, B., LaPointe, P., Eiben, T., et al. 2000. Integration of Discrete Feature Network Methods with
Conventional Simulator Approaches. SPE Res. Eval. & Eng. 3 (2): 165-170.
https://doi.org/10.2118/62498-PA.
Dicker, A. I., and Smits, R. M. 1988. A Practical Approach for Determining Permeability from Laboratory
Pressure-Pulse Decay Measurements. Paper presented at International Meeting on Petroleum
Engineering, Tianjin, China, November 1988. https://doi.org/10.2118/17578-MS.
Ding, Y. D. 2019. Modeling of Matrix/Fracture Transfer with Nonuniform-Block Distributions in LowPermeability Fractured Reservoirs. SPE J. 24 (6): 2653-2670. https://doi.org/10.2118/191811-PA.
Dykhuizen, R.C. 1990. A New Coupling Term for Dual-Porosity Models. Water Resour. Res. 26 (2): 351-356.
https://doi.org/10.1029/WR026i002p00351.
Fu Q., Cudjoe, S., Ghahfarokhi, R. B. 2021. Investigating The Role of Diffusion in Hydrocarbon Gas Huff-nPuff Injection - an Eagle Ford Study. J. Pet. Sci. & Eng., 198, 108146.
https://doi.org/10.1016/j.petrol.2020.108146.
Fuller, E.N., Ensley, K., and Griddings, J.C. 1969. Diffusion of Halogenated Hydrocarbons in Helium: The
Effect of Structure on Collision Cross Sections. J. Phys. Chem. 73 (11): 3679-3685.
http://10.1021/j100845a020.
Ghasemi, M., Astutik, W., Alavian, S. A., et al. 2018. Laboratory Tests and Modeling of Carbon Dioxide
Injection in Chalk with Fracture/Matrix-Transport Mechanisms. SPE Res. Eval. & Eng. 21 (01): 122–136.
https://doi.org/10.2118/180102-PA.
Gong, J. and Rossen, W. R. 2017. Modeling Flow in Naturally Fractured Reservoirs: Effect of Fracture
Aperture Distribution Sub-Network for Flow. Pet Sci. 14 138-154. https://doi.org/10.1007/s12182-016-
0132-3.
Guidry, K., Luffel, D., and Curtis, J. 1996. Development of Laboratory and Petrophysical Techniques for
Evaluating Shale Reservoirs. Final technical report, October 1986-September 1993. ResTech Houston,
Inc., Texas, United States.

126
Hawthorn, S.B., Gorecki, C.D., Sorensen, J.A. 2013. Hydrocarbon Mobilization Mechanisms from Upper,
Middle, and Lower Bakken Reservoir Rocks Exposed to CO2. Paper presented at SPE Canada
Unconventional Resources Conference, Calgary, Alberta, Canada, November 2013.
https://doi.org/10.2118/167200-MS.
Hayduk, W. and Minhas, B. S. 1982. Correlations for Prediction of Molecular Diffusivities in Liquids. Can.
J. Chem. Eng. 60 (2): 295-299. https://doi.org/10.1002/cjce.5450600213.
Hoffman, B. T. 2018. Huff-n-Puff Gas Injection Pilot Projects in the Eagle Ford. Paper was prepared for
presentation at SPE Canada Unconventional Resources Conference in Calgary, Alberta, Canada, March
2018. https://doi.org/10.2118/189816-MS.
Holman, J.P. 1990. Heat Transfer, 7th ed. McGraw-Hill, New York, N.Y., pp. 160-164.
Hoteit, H. 2013. Modeling Diffusion and Gas-Oil Mass Transfer in Fractured Reservoirs. J. Pet. Sci. & Eng.
105: 1-17. https://doi.org/10.1016/j.petrol.2013.03.007.
Hoteit, H. and Firoozabadi, A. 2006. Compositional Modeling of Discrete-Fractured Media Without
Transfer Functions by the Discontinuous Galerkin and Mixed Methods. SPE J. 11 (3): 341–352.
https://doi.org/10.2118/90277-PA.
Jia, B., Tsau, J.-S., Barati,R. 2018. Role ofMolecular Diffusion in Heterogeneous, Naturally FracturedShale
Reservoirs During CO2 Huff-n-Puff. J. Petrol. Sci. & Eng. 164: 31–42. 10.1016/j.petrol.2018.01.032.
Jiang, J., and Younis, R. 2016. Hybrid Coupled Discrete-Fracture/Matrix and Multicontinuum Models for
Unconventional-Reservoir Simulation. SPE J. 21 (3): 1009–1027. https://doi.org/10.2118/178430-PA.
Jones, S.C. 1997. A Technique for Faster Pulse-Decay Permeability Measurements in Tight Rocks. SPE Form
Eval 12 (01): 19-25. https://doi.org/10.2118/28450-PA.
Karimi-Fard, M., Gong, B., Durlofsky, L. J. 2006. Generation of Coarse-Scale Continuum Flow Models from
Detailed Fractured Characterizations. Water Resour. Res. 42: W10423.
https://doi.org/10.2118/178430-PA.
Kazemi, H., Merrill L.S. Jr., Porterfield K.L. et al. 1976. Numerical Simulation of Water-Oil Flow in Naturally
Fractured Reservoirs. SPE J. 16 (6): 317-326. https://doi.org/10.2118/5719-PA.
Kazemi, H., Gilman, J.R., Elsharkawy, A.M. 1992. Analytical and Numerical Solution of Oil Recovery From
Fractured Reservoirs With Empirical Transfer Functions. SPE Res. Eng. 7 (2): 219-227.
https://doi.org/10.2118/19849-PA.
Kazemi, H. and Gilman, J.R. 1993. Multiphase Flow in Fractured Petroleum Reservoirs. In Flow and
Containment Transport in Fractured Rock, ed. J. Bear, C.-F. Tsang, G. de Marsily, 267-323. San Diego,
California: Academic Press.
Kooijman, H. A., and Taylor, R. 1991. Estimation of Diffusion Coefficients in Multicomponent Liquid
Systems. Ind. Eng. Chem. Res. 30 (6): 1217-1222. https://doi.org/10.1021/ie00054a023.

127
Leahy-Dios, A. and Firoozabadi, A. 2007. Unified Model for Nonideal Multicomponent Molecular Diffusion
Coefficients. AIChe J. 53 (11): 2932-2939. https://doi.org/10.1002/aic.11279
Lee, S.H., Jensen, C.L., and Lough, M.F. 2001. Hierarchical Modeling of Flow in Naturally Fractured
Formations with Multiple Length Scales. Water Resour. Res. 37 (3): 443-455.
https://doi.org/10.1029/2000WR900340.
Lemmon, E.W., McLinden, M.O., Friend, D.G., 2005. In: Thermophysical Properties of Fluid Systems. NIST
Chemistry Webbook, vol. 69. NIST standard reference database. https://doi.org/10.18434/T4D303.
Li, W., 1972. Differential Equations of Hydraulic Transients, Dispersion, and Groundwater Flow. PrenticeHall, Englewood Cliffs.
Li, J., Lei, Z., Qin, G., et al. 2015. Effective Local-Global Upscaling of Fractured Reservoirs Under Discrete
Fractured Discretization. Energies 8 (9): 10178-10197. http:/dx.doi.org/10.3390/en80910178.
Li, L., and Lee, S.H. 2008. Efficient Field-Scale Simulation of Black Oil in a Naturally Fractured Reservoir
through Discrete Fracture Networks and Homogenized Media. SPE Res. Eval. & Eng. 11 (4): 750-758.
http://doi.org/10.2118/103901-PA.
Li, L., Sheng, J. J., and Xu, J. 2017. Gas Selection for Huff-n-Puff EOR in Shale Oil Reservoirs Based Upon
Experimental and Numerical Study. Paper was prepared for presentation at SPE Unconventional
Resources Conference held in Calgary, Alberta, Canada, 15-16 February 2017.
https://doi.org/10.2118/185066-MS.
Lie, K.-A. 2019. An Introduction to Reservoir Simulation Using MATLAB/GNU Octave: User Guide for the
MATLAB Reservoir Simulation Toolbox (MRST). Cambridge, UK: Cambridge University Press.
Lie, K.-A., and Moyner, O. 2021. Advanced Modeling with the MATLAB Reservoir Simulation Toolbox.
Cambridge, UK: Cambridge University Press. http://doi.org/10.1017/9781009019781.
Lim, K.T., and Aziz, K. 1995. Matrix-Fracture Transfer Shape Factors for Dual-Porosity Simulators. J Pet. Sci.
& Eng. 13 (3): 169-178. https://doi.org/10.1016/0920-4105(95)00010-F.
Lohrenz, J., Bray, B. G., and Clark, C. R. 1964. Calculating Viscosities of Reservoir Fluids from Their
Compositions. J of Pet Technol 16 (10): 1-171. 10.2118/915-PA.
Lu, H., Di Donato, G., Blunt, M.J. 2008. General Transfer Functions for Multiphase Flow in Fractured
Reservoirs. SPE J. 13 (3): 289-297. https://doi.org/10.2118/102542-PA.
MATLAB Reservoir Simulation Toolbox (MRST), version 2020a, SINTEF Applied Mathematics, December
2019. https://www.sintef.no/projectweb/mrst/.
Moinfar, A. 2013. Development of an Efficient Embedded Discrete Fracture Model for 3D Compositional
Reservoir Simulation in Fractured Reservoirs. Ph.D. dissertation, University of Texas, Austin, Texas
(August 2013).

128
Moyner, O. 2021. Compositional Simulation with the AD-OO Framework. In Advanced Modeling with the
MATLAB Reservoir Simulation, K.A. Lie and O. Moyner, Chap. 8, 324-371. Cambridge: Cambridge
University Press.
Penuela, G., Civan, F., Hughes, R.G. et al. 2002. Time-Dependent Shape Factors for Interporosity Flow in
Naturally Fractured Gas-Condensate Reservoirs. Paper was prepared for presentation at SPE Gas
Technology Symposium held in Calgary, Alberta, Canada, April 2002. https://doi.org/10.2118/75524-
MS.
Poling, B.E., Prausnitz, J.M., O’Connell, J.P. The Properties of Gases and Liquids, 5
th ed. McGraw-Hill, New
York (2001).
Ranjbar, E., Hassanzadeh, H. 2011. Matrix-Fracture Transfer Shape Factor for Modeling Flow of a
Compressible Fluid in Dual-Porosity Media. Adv. Water Resour. 34: 627-639.
https://doi.org/10.1016/j.advwatres.2011.02.012.
Reamer, H.H, Sage, B.H. 1956. Diffusion Coefficients in Hydrocarbon Systems. Methane-n-ButaneMethane in Liquid Phase. Chem. & Eng. Data Series 1 (01): 71-77. 10.1021/i460001a014.
Reid, R. C., Prausnitz, J. M., and Poling, B.E. 1987. The Properties of Gases and Liquids, 4th ed. McGrawHill, New York.
Renner, T.A. 1988. Measurement and Correlation of Diffusion Coefficients for CO2 and Rich-Gas
Applications. SPE Res. Eng. 3 (02): 517–523. https://doi.org/10.2118/15391-PA
Robinson, D. B. and Peng, D. Y. (1978). The Characterization of the Heptanes and Heavier Fractions for the
GPA Peng-Robinson Programs. Gas Processors Association.
Rossen, R.H., and Shen, E.I.C. 1989. Simulation of Gas/Oil Drainage and Water/Oil Imbibition in Naturally
Fractured Reservoirs. SPE Res. Eng. 4 (4): 464-470. https://doi.org/10.2118/16982-PA.
Sarma, P., Aziz, K. 2006. New Transfer Functions for Simulation of Naturally Fractured Reservoirs with
Dual-Porosity Models. SPE J. 11 (3): 328-340. https://doi.org/10.2118/90231-PA.
Schlumberger. 2011. ECLIPSE (300) Reservoir Engineering Software,
https://www.software.slb.com/products/eclipse.
Shakiba, Mahmood. 2014. Modeling and Simulation of Fluid Flow in Naturally and Hydraulically Fractured
Reservoirs Using Embedded Discrete Fracture Model (EDFM).
Sherafati, M. Jessen, K. 2018. Coarse-Scale Modeling of Multicomponent Diffusive Mass Transfer in DualPorosity Models. Paper presented at the SPE Annual Technical Conference and Exhibition, Dallas,
Texas, USA, September 2018. https://doi.org/10.2118/191568-MS.
Shi, Z., Khodaparast, P., Hu, S., et al. 2022. Measurement and modeling of methane diffusion in
hydrocarbon mixtures. Fuel 324, 124740. https://doi.org/10.1016/j.fuel.2022.124740.

129
Shojaei, H., Jessen, K. 2014. Diffusion and Matrix-Fracture Interactions during Gas Injection in Fractured
Reservoirs. Paper presented at the SPE Improved Oil Recovery Symposium, Tulsa, Oklahoma, USA, April
2014. 10.2118/169152-MS.
Siddiqi, M. A., andLucas, K. (1986). Correlations for Prediction of Diffusion in Liquids. Can. J. of Chem. Eng.
64 (5): 839-843. 10.1002/cjce.5450640519.
Sigmund, P. M. 1976a. Prediction of Molecular Diffusion at Reservoir Conditions. Part 1- Measurement
And Prediction of Binary Dense Gas Diffusion Coefficients. J. Can. Pet. Technol. 15 (2). 10.2118/76-02-
05.
Sigmund, P. M. 1976b. Prediction of Molecular Diffusion at Reservoir Conditions. Part II - Estimating the
Effects Of Molecular Diffusion And Convective Mixing In Multicomponent Systems. J Can Pet Technol
15 (3). 10.2118/76-03-07.
Sistan, M., Ghazanfari, M., Jamshidi, S. 2022. The Impact of the Double Cross-Phase Diffusion Mechanism
on Oil Recovery during CO2 Injection into Fractured Rocks: A Simulation Study. ACS Omega 7 (25):
21630-21642. 10.1021/acsomega.2c01503.
Sofla, S., Pouladi, B., Sharifi, M., et al. 2016. Experimental and Simulation Study of Gas Diffusion Effect
During Gas Injection into Naturally Fractured Reservoirs. J. of Nat. Gas Sci. & Eng. 33: 438-447.
10.1016/j.jngse.2016.05.035.
Stiel, L. I. and Thodos, G. 1961. The Viscosity of Nonpolar Gases at Normal Pressures. AiChE J. 7 (4): 611-
615. 10.1002/aic.690070416.
Tan, K. K., Thrope, R. B. 1992. Gas Diffusion into Viscous and Non-Newtonian Liquids. Chem. Eng. Sci. 47
(13-14): 3565-3572. 10.1016/0009-2509(92)85071-I.
Taylor, R. and Krishna, R. 1993. Multicomponent Mass Transfer, first ed. John Wily & Sons, New York.
Tene, M., Bosma, S. B. M., Al Kobaisi, M. S., et al. 2017. Projection-Based Embedded Discrete Fracture
Model (pEDFM). Adv. Water Resour. 105: 205-216. 10.1016/j.advwatres.2017.05.009.
Tian, Y., Zhang, C., Lei, Z., et al. 2021. An Improved Multicomponent Diffusion Model for Compositional
Simulation of Fractured Unconventional Reservoirs. SPE J. 26 (05): 3316–3341. 10.2118/204010-PA.
Toor, H. L. 1964. Solution of the Linearized Equations of Multicomponent Mass Transfer: 1. AIChe J. 10 (4)
448-455. 10.1002/aic.690100408.
Umezawa, S., Nagashima, A., 1992. Measurement of The Diffusion Coefficients of Acetone, Benzene, and
Alkane in Superctitical CO2 by the Taylor Dispersion Method. J. of Supercritical Fluids 5 (4): 242-250.
10.1016/0896-8446(92)90014-B.
Vermeulen, T. 1953. Theory for Irreversible and Constant-Pattern Solid Diffusion. Industrial & Engineering
Chemistry 45 (8), 1664-1670. 10.1021/ie50524a025.

130
Vo, H., Mallison, B., Kamath, J., et al. 2020. Study of Oil Recovery Mechanisms in Complex Natural Fracture
Systems using Embedded Discrete Fracture Models. Paper presented at the SPE Annual Technical
Conference and Exhibition, Virtual, October 2020. 10.2118/201463-MS.
Wan, T. and Sheng, J. 2015. Compositional Modelling of the Diffusion Effect on EOR Process in Fractured
Shale-Oil Reservoirs by Gasflooding. J. Can. Pet. Technol. 54 (02): 107–115. 10.2118/2014-1891403-
PA.
Warren, J.E., and Root, P.J. 1963. The Behavior of Naturally Fractured Reservoirs. SPE J. 3 (3): 245-255.
10.2118/426-PA.
Wilke, C. R. 1950. Diffusion Properties of Multicomponent Gases. Chem. Eng. Prog. 46: 95-104.
Wilke, C. R., and Chang, P. 1955. Correlation of Diffusion Coefficients in Dilute Solutions. AIChE J. 1 (2):
264-270. 10.1002/aic.690010222.
Whitson, C.H. and Kuntadi, A. 2005. Khuff Gas Condensate Development. Paper presented at the
International Petroleum Technology Conference, Doha, Qatar, 21-23 November 2005. 10.2523/IPTC10692-MS.
Wong, C. F., and Hayduk, W. 1990. Correlations for Prediction of Molecular Diffusivities in Liquids at
Infinite Dilution. Can. J. of Chem. Eng. 68 (5): 849-859. 10.1002/cjce.5450680516.
Wu, C., Jessen, K. 2020. Consistent and Efficient Representation of Diffusive Mass Transfer in Fractured
Reservoirs. Paper presented at the SPE Improved Oil Recovery Conference, Virtual, August 2020.
10.2118/200453-MS
Xu, J., Sun, B., Chen, B. 2019. A Hybrid Embedded Discrete Fracture Model for Simulating Tight Porous
Media with Complex Fracture Systems. J. Pet. Sci. & Eng. 174: 131-143.
https://doi.org/10.1016/j.petrol.2018.10.094.
Yu, W., Zhang, Y., Varavei, A., et al. 2019. Compositional Simulation of CO2 Huff ’n’ Puff in Eagle Ford Tight
Oil Reservoirs with CO2 Molecular Diffusion, Nanopore Confinement, and Complex Natural Fractures.
SPE Res. Eval. & Eng. 22 (02): 492–508. 10.2118/190325-PA.
Zimmerman, R.W., Chen, G., Hadgu, T. et al. 1993. A Numerical Dual-Porosity Model with Semianalytical
Treatment of Fracture/Matrix Flow. Water Resour. Res. 29 (7): 2127-2137. 10.1029/93WR00749.
Zhang, Y., Ghanbarnezhad Moghanloo, R., & Davudov, D. 2019. Pore Structure Characterization of a Shale
Sample Using SEM Images. Paper presented at the SPE Western Regional Meeting, San Jose, California,
USA, April 2019. 10.2118/195352-MS.
Zuloaga-Molero, P., Yu, W., Xu, Y., et al. 2016. Simulation Study of CO2-EOR in Tight Oil Reservoirs with
Complex Fracture Geometries. Sci. Rep. 6 (1): 33445-33445. 10.1038/srep33445.

Abstract (if available)

Abstract The current motivation and efforts of the oil and gas industry are hinged on the fact that any future increase in hydrocarbon production will require further development of unconventional resources. Tight sand and shale hold vast hydrocarbon reserves that could support increasing energy demands; however, the current low recovery factors remain a challenge. Therefore, technical studies have been ongoing to improve existing practices and unlock the potential of these reservoirs. To this end, gas injection into tight fractured reservoirs has been a successful enhanced oil recovery (EOR) method, making it a promising investment and an active research area.

To support ongoing field activities and incentivize new developments, we rely on available technology for simulation and optimization of production operations. It is important that these tools are based on models that accurately reflect the relevant physical processes. As of today, many details still remain to be resolved to delineate the dynamics of the tight fractured systems and their implementation into reservoir simulation tools.

In this research, we explore fit-for-purpose models to describe multicomponent matrix/fracture mass transfer in unconventional formations. The objective is to provide a modeling platform, via an open-source code, where our findings are tested and validated. Our development includes effective algorithms to facilitate further advancements in research and implementation of simulation at larger scale.

We start by formulating a new transfer function to represent the matrix/fracture mass transfer for tight fractured systems via the dual-porosity modeling (DPM) approach. The formulation is inspired by the work of Vermeulen (1953), who approximated the analytical solution for diffusion in a spherical geometry using a closed exponential form to describe the average change of concentration inside the domain of interest. By generalizing Vermeulen’s approximation, we introduce two parameters to represent the geometry of the system: A shape factor σ and an exponent V. The final form, referred to as the generalized Vermeulen (gVer) transfer function, includes an implicit time-dependent correction factor that represents the sub-scale physics of the problem that allows for accurate prediction of the prolonged transient behavior in tight fractured systems.

The generalized Vermeulen transfer function was implemented in the MATLAB Reservoir Simulation Toolbox (MRST) to allow for a comparison with commercial simulators (CMG and ECLIPSE). An initial validation of the new transfer function was performed via single-phase pressure depletion calculations, where a fine-grid single-porosity model (SPM) is compared with an equivalent single-block DPM for a tight fractured system with a 5 nD matrix. First, we observe that the SPM reference solution, from the three simulation tools, agree with one another. Then, by comparing the results from DPM calculations with MRST and commercial tools, we observe that the application of the Warren and Root (1963) transfer function with a shape factor from Lim and Aziz (1995) provides for a reasonable agreement. However, these results, obtained via the traditional transfer function, depart from the reference (SPM) solution at early time but reach agreement at later time. In contrast, the gVer transfer function is demonstrated to provide for an excellent agreement with the reference solution throughout the entire depletion process.

Several examples were subsequently defined for a broader range of problems. The objectives were to document the limitations of the traditional DPM, and to further investigate the performance of the new gVer transfer function formulation. These calculation examples include: 1) matrix blocks with isotropic permeability; 2) matrix blocks with anisotropic permeability; and 3) sensitivity to variable hydraulic diffusivity (highly compressible fluids). In all the cases, the DPM using the generalized Vermeulen transfer function is demonstrated to agree with the response of the reference solution for matrix permeability values as low as 5 nD.

The initial testing of the gVer transfer function was performed on well-defined matrix geometries, i.e., for problems that can be described by readily available analytical solutions. We then proceed to demonstrate that the new transfer function is not limited to simple geometries but can also provide an accurate description of flow behaviors for complex fracture geometries. In such settings, the relevant matrix/fracture interaction parameters, σ and V, are obtained from a reference solution, that, for example, can be generated by calculations with an SPM or an embedded discrete fracture model (EDFM). We apply this approach to address problems with arbitrary fracture orientations and geometries in heterogeneous nano-darcy domains. We furthermore demonstrate the versatility of our approach by refining the DPM in regions of higher spatial variation to capture the details of the fractured system. A workflow for converting the EDFM to an equivalent DPM is subsequently discussed in the context of upscaling of discrete fracture characterizations.

We then developed a DP compositional model in MRST to facilitate studies of multicomponent transport in tight fractured systems. This model was first validated against an analytical solution for a single-phase pressure depletion process. Following that, we extended our study to include compositional two-phase problems to model tight fractured gas-condensate reservoirs. However, we find that for such problems, the implicit and coarse representation of the matrix/fracture interaction, described in a similar manner to that of Kazemi et al. (1976), results in an overprediction of the liquid production: The averaging of matrix block properties masks the compositional changes near the matrix/fracture interface, that is critical to predict production behaviors below the dew point. To overcome the limitations of the current models, we propose a novel approach that evaluates the gVer transfer function at conditions that reflect/approximate the pressure and composition/saturation gradients within a matrix block. We present calculation results for a single-block DPM, representative of a tight fractured gas-condensate reservoir considering two condensate fluid examples: 1) A 4-component analog rich gas; and 2) A realistic 24-component fluid from an existing gas-condensate reservoir. We demonstrate that calculations with the new modified gVer transfer function substantially agree with the fine-grid SPM results in terms of the amounts and compositions of the produced phases.

We subsequently directed our attention to study the role of molecular diffusion on recovery for tight fractured reservoirs. To do that, we conducted an integrated lab experimental and numerical modeling efforts of carbon dioxide huff-n-puff (HnP) on an Eagle Ford shale core (EF-3) with a porosity of 11.9% and permeability of 2.5 µD. Four cycles of HnP were performed where the EF-3 core was initially saturated with a synthetic oil composed of n-decane, n-dodecane, n-tetradecane, and n-hexadecane. These experiments were carried out at pressures and temperature that mimic reservoir conditions (from 3,900 psi to 4,400 psi and 158 F/343.15 K) under first-contact miscibility conditions. Based on the pressure drop during soaking as well as the compositional analysis of the produced oil, molecular diffusion is observed to be the dominant recovery mechanism in the performed HnP experiments. To interpret the experimental observations, a compositional simulation study, via the SPM approach, was performed. The governing equations of MRST were modified to represent molecular diffusion described by either the classical Fick’s law or the generalized Fick’s law, with composition-dependent diffusion coefficients. For both diffusion models, our simulation calculations show very good agreement with the lab data (for all cycles) in terms of the recovery and produced oil compositions, with the former model demonstrating higher accuracy and computational efficiency.

In summary, this dissertation reports on our efforts to combine physical models and efficient simulation techniques to advance the understanding of multicomponent mass transfer in tight fractured rocks. The overarching goal was to explore and develop computationally simpler, yet accurate, approaches to enable/support more practical and reliable forecasting for unconventional reservoirs at larger scale.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

PDF

Modeling and simulation of complex recovery processes

PDF

Mass transfer during enhanced hydrocarbon recovery by gas injection processes

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

PDF

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

PDF

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

$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

Investigation of gas transport and sorption in shales

PDF

Reactivation of multiple faults in oilfields with injection and production

PDF

Investigation of gas sorption and mass transfer in gas shales

PDF

Modeling and simulation of circulating tumor cells in flow. Part I: Low-dimensional deformation models for circulating tumor cells in flow; Part II: Procoagulant circulating tumor cells in flow

PDF

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

$Isotopic insights to human-impacted hydrologic systems: from urban watersheds to hydraulically fractured basins$

PDF

Isotopic insights to human-impacted hydrologic systems: from urban watersheds to hydraulically fractured basins

Asset Metadata

Creator Raslan, Mohammed Samir (author)

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

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Degree Conferral Date 2024-08

Publication Date 09/11/2024

Defense Date 05/20/2024

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

Tag compositional modeling,diffusive mass transfer,dual-porosity model,embedded discrete fracture model,enhanced oil recovery,HnP,Huff-n-puff,matrix fracture transfer function,molecular diffusion,multiphase transfer function,OAI-PMH Harvest,shale,single-porosity model,tight fractured systems,transient behavior,unconventional reservoirs,Warren and Root (1963)

Format theses (aat)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Jessen, Kristian (committee chair), Ghanem, Roger (committee member), Jha, Birendra (committee member)

Creator Email mraslan@usc.edu

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

Unique identifier UC11399ANKZ

Identifier etd-RaslanMoha-13497.pdf (filename)

Legacy Identifier etd-RaslanMoha-13497

Document Type Dissertation

Format theses (aat)

Rights Raslan, Mohammed Samir

Internet Media Type application/pdf

Type texts

Source 20240912-usctheses-batch-1211 (batch), 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 author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright.

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

Repository Email cisadmin@lib.usc.edu