University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content UNIVERSITY OF SOUTHERN CALIFORNIA

Pore-Scale Characterization and Fluid-Flow Simulation:
Application of the Lattice Boltzmann Method to Conventional
and Unconventional Rocks

Shahram Farhadi Nia

PHD DISSERTATION
FACULTY OF THE USC GRADUATE SCHOOL
MORK FAMILY DEPARTMENT OF CHEMICAL ENGINEERING AND
MATERIALS SCIENCE
LOS ANGELES, CALIFORNIA

December 2014
© SHAHRAM FARHADI NIA 2014
ii

Abstract
Flow and transport in porous rocks is controlled by the pore structures and surface chemistry
that, in part, defines a porous material. Pores are commonly in the micron range for conventional
oil and gas formations while often in the nano-meter range for unconventional resources.
Traditionally, core experiments are performed for characterization of reservoir rocks. Recently,
with advances in imaging and computational power, we are able to image the pores in most
relevant rocks and to construct a digital rendering of the rock that can subsequently be used for
modeling of flow and transport processes.
In this manuscript, we introduce a workflow that takes the rock images as input. In the first
stage of the workflow, a 3D image is analyzed to find effective and total porosities. Single-phase
flow simulation, using the Lattice Boltzmann Method (LBM), is then performed to evaluate
directional (e.g. x, y and z) permeability of the rock sample. Finally, two-phase flow simulations
are performed on the pore geometries dictated by the sample image to determine multiphase
characteristics, including relative permeability of the rock sample.
To place imaging technology in context, we compare and discuss the insight and
understanding about the pore structure characteristics of rocks, specifically for Monterey shale
samples that can be obtained from various characterization techniques. The characterization
techniques considered here include: (i) Mercury injection capillary pressure experiments
(MICP), (ii) Brunauer, Emmett and Teller (BET) – Nitrogen adsorption experiments, (iii) High-
Resolution X-ray CT (HRXCT), and (iv) Focused Ion Beam – Scanning Electron Microscopy
(FIB-SEM). The HRXCT and FIB-SEM characterization data are utilized to reconstruct
(visualize) the 3D structures of the rocks; which are then used to calculate the permeabilities of
these shale samples via LBM flow simulations.
iii

The LBM approach for flow simulation is a specific discretized form of the Boltzmann
equation. The Boltzmann equation, the cornerstone of Kinetic theory, is discussed and the
derivation of different 2D and 3D LBM schemes is presented. The 3D single-phase LBM
approach is then used for calculating the flow field at the pore scale and allows us to calculate
directional permeabilities. We demonstrate the application of LBM in this context for sandstone
as well as shale images.
Shan and Chen multiphase LBM is used for two-phase flow modeling in this work.
Benchmark studies for two-phase systems, including surface area minimization and Laplace law,
are presented. Wetting phenomena in 2D and 3D is subsequently discussed. With all components
of LBM validated and in place, we then proceed to calculate the relative permeability curves for
2D examples of porous rocks. The effect of viscous coupling on relative permeability curve is
then discussed.
In the last section, we study the dynamics of capillary-driven fluid invasion in three
different settings including: a) a single capillary tube, b) a porous column and c) a fractured
porous media. A Lambert W functional form is proposed to describe the invasion dynamics in a
single capillary tube, that predicts both early time Washburn-type behavior ( ) and late time
exponential decay behavior. We extend the formulation to describe homogenous porous media
and to include viscosity, pressure and gravity effects in both advancing and defending fluids.
Solutions for closed systems, where the advancing fluid compresses the defending fluid, are then
developed. Finally, we extend the theory to describe fractured systems, and propose a
convolution integral formulation and a new explicit solution.

iv

Dedicated to My Parents,
Kobra and Hassan

v

Acknowledgements
I would like to express my gratitude to everyone who inspired me to start, encouraged me to
continue and helped me to finish my doctoral studies.
In particular, I would like to thank my PhD advisor, Professor Kristian Jessen, for his
guidance through my education and research at the University of Southern California. His
exceptional enthusiasm for discovery led us to many interesting discussions on different research
subjects, some of which became parts of this thesis.
My appreciation is extended to the committee members of my dissertation, Professor Iraj
Ershaghi and Professor Aiichiro Nakano. I am especially thankful to Professor Ershaghi, director
of the Petroleum Engineering program, for being a great mentor through my academic and
industrial journeys.
I would like to thank Occidental Petroleum Corporation for supporting this research at
USC. I owe OXY a debt of gratitude for providing internship opportunities that helped me gain
invaluable industrial experience. I am grateful to all my mentors, managers and colleagues that
helped me with projects. I am extremely thankful to Bryan Bell for his insightful comments that
improved my work tremendously.
My gratitude goes to all the great friends that I have been fortunate to have in my life. I
cannot name all of them but I wish to thank Mohammad Evazi, Mehran Hosseini, Arman
Khodabakhshnezhad, Mohsen Taghavifar, Hamed Darabi, Cyrus Ashayeri and Azarang
Golmohammadi.
vi

In the end, my deepest gratitude goes to my parents, my brother, Hadi, and my lovely
sisters, Mojgan and Masoumeh, for their lifetime support and encouragement.

vii

Table of Contents
Abstract ........................................................................................................................................... ii
Acknowledgements ......................................................................................................................... v
List of Tables ................................................................................................................................. ix
List of Figures ................................................................................................................................. x
CHAPTER ONE: INTRODUCTION ............................................................................................. 1
1.1 Heterogeneity and Scale in Reservoir Rocks ................................................................... 1
1.2 Pore-Scale Modeling and Simulation ............................................................................... 3
1.2.1 Pore Network Modeling ............................................................................................ 3
1.2.2 Direct Calculations.................................................................................................... 5
1.3 Pore Scale Characterization Workflow ............................................................................ 7
1.3.1 Rock Imaging ............................................................................................................ 8
1.4 Areas of Application ...................................................................................................... 13
1.4.1 Unconventional Resources (Physics at Play) .......................................................... 13
1.4.2 Conventional Resources (Co- and Counter-Current Flows) ................................... 15
1.5 Organization of the Manuscript...................................................................................... 16
CHAPTER TWO: CHARACTERIZATION OF SHALE SAMPLES (FROM ROCKS TO
NUMBERS) .......................................................................................................................... 18
2.1 Introduction .................................................................................................................... 18
2.2 Porosity Characterization Techniques ............................................................................ 21
2.2.1 Mercury Porosimetry (MICP) ................................................................................. 21
2.2.2 BET Analysis .......................................................................................................... 22
2.2.3 Imaging and Image Analysis .................................................................................. 24
2.3 Permeability from Flow Calculations ............................................................................ 31
2.4 Spontaneous Imbibition Experiments and Simulations ................................................. 32
2.5 Summary and Conclusions ............................................................................................. 36
CHAPTER THREE: THE LATTICE BOLTZMANN METHOD ............................................... 39
3.1 The Boltzmann Equation ................................................................................................ 39
3.2 Relaxation Time Approximation (RTA) for Collision Term ......................................... 41
3.3 From the Boltzmann Equation to the Lattice Boltzmann Equation: .............................. 42
3.4 Numerical Implementation and Boundary Conditions................................................... 46
3.5 Directional Permeability from LBM Simulations .......................................................... 48
CHAPTER FOUR: LATTICE BOLTZMANN SIMULATIONS OF TWO PHASE FLOW IN
POROUS MEDIA ................................................................................................................. 54
4.1 Introduction .................................................................................................................... 54
4.2 Multiphase LBM schemes .............................................................................................. 54
4.3 Surface Minimization Test (2D and 3D) ........................................................................ 56
4.4 Laplace Law Validation ................................................................................................. 59
4.5 Implementation of Adhesion Forces .............................................................................. 60
4.6 Wetting in 2D and 3D .................................................................................................... 61
4.7 Relative Permeability and its Functional Dependencies ................................................ 65
4.8 Viscous coupling effects (analytic and LBM)................................................................ 70
viii

CHAPTER FIVE: THEORETICAL ANALYSIS OF CAPILLARY RISE IN THE POROUS
MEDIA .................................................................................................................................. 79
5.1 Introduction .................................................................................................................... 79
5.2 Single Capillary Tube..................................................................................................... 81
5.3 Homogenous Porous Media ........................................................................................... 85
5.4 Heterogeneous Porous Media ........................................................................................ 93
5.5 Discussions and Conclusions ......................................................................................... 98
CHAPTER SIX: POTENTIAL FUTURE WORK AND RESEARCH DIRECTIONS ............. 101
REFERENCES ........................................................................................................................... 104
APPENDIX-A: Maxwell-Boltzmann Equilibrium Distribution Function ................................. 113
APPENDIX-B: Derivation of D
2
Q
9
Scheme .............................................................................. 119

ix

List of Tables
Table 2-1. Parameters of the imbibition model ………………………………………….35
Table 3-1. Directional Permeabilities ……………………………………………………53
Table 4-1. Different simulation cases in 2D channel ……………………………………72

x

List of Figures
Figure 1.1. Pore network for a Brea sandstone (Piri, 2003) ..............................................................4
Figure 1.2. Velocity profile is disturbed by solid particles, Re=10 (Rahimain and
Pourshaghaghy, 2002)................................................................................................................5
Figure 1.3. Pore Scale Characterization Workflow ..........................................................................8
Figure 1.4. Gray scale image of a carbonate rock. The L.H.S. matrix is the number
representation of the highlighted portion of the picture. 0 is black (pore) and 255 is white
(matrix). (Courtesy of Dr. Masa Prodanovic) ..........................................................................10
Figure 1.5. Gray scale image with binary images at different thresholds with the corresponding
calculated porosities. The threshold value is adjusted from porosity experiments. .................11
Figure 2.1. Mercury intrusion data for a quartz-phase Monterey sample. .......................................22
Figure 2.2. BJH cumulative pore volume data for the shale sample. ..............................................23
Figure 2.3. HK cumulative pore volume data for the shale sample. ...............................................24
Figure 2.4. Vugs/pores (dark features) and high-density regions (white features). .........................25
Figure 2.5. Spatial distribution of pyrite from imaging at 2-µm resolution. ...................................26
Figure 2.6. Pore size distribution from 2-µm HRXCT imaging. .....................................................27
Figure 2.7. Pore size distribution from the 0.5-µm HRXCT imaging. ............................................28
Figure 2.8. High-resolution surface SEM image with visible nano pores. ......................................29
Figure 2.9. 2D slices of the FIB-SEM image at a resolution of 30 nm. ..........................................30
Figure 2.10. Pore size distribution from FIB-SEM imaging. ..........................................................31
Figure 2.11. (a) Solid phase, (b) pore space and (c) velocity in the x direction. The domain size
is 1003 voxel (for illustration). ...............................................................................................32
Figure 2.12. Water uptake (cc) vs. square root of time (hours) over the 23 hours of the
spontaneous imbibition experiment. ........................................................................................33
Figure 2.13. Summary of porosity characterization data. ................................................................37
Figure 3.1. N molecules away from equilibrium state .....................................................................40
Figure 3.2. The D2Q9 scheme, two dimensional with nine velocities. ...........................................44
Figure 3.3. Inlet velocity boundary condition ..................................................................................47
xi

Figure 3.4. Pore structure in a Fontainebleau sandstone (125^3) ....................................................48
Figure 3.5. Simulation setup ............................................................................................................49
Figure 3.6. Contour plot of velocity in x direction. Flow is in x direction. .....................................50
Figure 3.7. Contour plot of velocity in y direction. Flow is in x direction. .....................................50
Figure 3.9. Inlet and outlet pressure, flow is in the x direction. ......................................................51
Figure 3.8. Contour plot of velocity in z direction. Flow is in x direction. .....................................51
Figure 3.11. Convergence behavior of the permeability .................................................................52
Figure 4.1. 2D Surface minimization for two immiscible phases, domain size is 128x128 and
periodic boundary conditions are applied in all directions. Blue is water and red is oil. ........57
Figure 4.2. 3D Surface minimization for two immiscible phases, domain size is 100x100x100
and periodic boundary conditions are applied in all directions ...............................................58
Figure 4.3. Simulation setup for droplets of different radii. Periodic boundaries are applied in
all directions. ............................................................................................................................59
Figure 4.4. Validation of the Multi Component Multi Phase (MCMP) LBM code via Laplace
law: A linear relationship between the pressure difference and inverse radius. For the two
fluids in this system, the IFT is 0.0051 mu/ts^2. mu (mass unit) and ts (time step) are
LBM units for mass and time, respectively. ............................................................................60
Figure 4.5. Initial conditions for wettability modeling in 2D: A 30x30 droplet of oil in a
120x80 domain. The lower boundary condition defines the solid wall and all other
boundary conditions are periodic. ............................................................................................61
Figure 4.6. Different wetting properties: Go and Gw are water and oil adsorption coefficients.
A complete range from non-wetting to completely water wet cases is shown. Gow is set
to 0.2. .......................................................................................................................................63
Figure 4.7. Wetting in 3D, starting from water droplet. This droplet remains spherical for the
oil-wet system and spreads on the surface for water-wet system (bottom row). Go and
Gw are water and oil adsorption coefficients. Gow (oil-water interface constant) is set to
0.2.............................................................................................................................................64
Figure 4.8. Distribution of two fluids at Sw=0.4. Water (blue) is the wetting phase. Black is
the solid and red is the non-wetting phase. The saturations start from an initial random
distribution. ..............................................................................................................................66
Figure 4.9. Vertical (flow direction) component of velocity at t=75000. ........................................67
xii

Figure 4.10. Distribution of wetting (blue) and non-wetting (orange) phases at different
wetting phase saturations. Densities and viscosities are the same and flow is in the
horizontal direction. .................................................................................................................69
Figure 4.11. Relative permeability curves from the simulated flow field at different saturations ..69
Figure 4.12. Schematic for layered flow of two immiscible phases in a channel ...........................70
Figure 4.13. Comparison of analytical solution and LBM (a=25, b=50, Gw=Gn=g) .....................73
Figure 4.14. Comparison of analytical solution and LBM (a=25, b=50, Gw=-g, Gn=g) ................74
Figure 4.15. Comparison of analytical solution and LBM (a=25, b=50, Gw=g, Gn=-g) ................75
Figure 4.16. Relative permeability curves for co and counter-current flow (Gw=g, Gn=-g).
Viscosities and densities are the same for both phases. ...........................................................77
Figure 4.17. Relative permeability curves for co-current flow at different viscosity ratios ............77
Figure 4.18. Relative permeability curves for counter-current flow (wetting phase flows in the
negative direction) at different viscosity ratios ........................................................................78
Figure 5.1. Schematic of capillary rise in a single tube ...................................................................82
Figure 5.2. Capillary rise in a tube with radius 10-4 m including inertial effects (Equation (5.5)
- dots) and without inertial effects (Equation (5.8) - solid line). ............................................84
Figure 5.3. Dynamics of capillary rise in nano- and micron-tubes. ................................................85
Figure 5.4. Schematic of a porous cylinder. L is the length of the system while Z denotes the
frontal position of the advancing fluid. The pore spaces are initially saturated 100% with
a defending fluid. .....................................................................................................................87
Figure 5.5. The effect of viscosity ratio on rise dynamics for a vertical homogeneous porous
medium (L = 15m) with a permeability of 1 Darcy (9.869*10-13 m2) and a porosity of
30% .........................................................................................................................................90
Figure 5.6. The dampening effect of pressure in a water-air closed system. The pressure
increases with the rise in the closed system and as a result, the fluid imbibition stops
earlier. ......................................................................................................................................93
Figure 5.7. Sugar-Cube representation of a fractured porous medium. The fracture system has
a porosity φf, a permeability kf and an average weighted pore size of E[rf], and the
matrix has a porosity φm, a permeabil ity km and an average weighted pore size of E[rm]. ..94
Figure 5.8. Fracture, matrix and total volumetric uptake for a fractured system with φt=0.3,
φm=0.25, φf=0.5, Kf =1 D, Km=1e -8 D, E[rf]=10 μm and E[rm]=10 nm. Total pore
volume is 3 cc and the fracture opening (w) is 0.5 mm. ..........................................................98
1

Chapter One: Introduction

1.1 Heterogeneity and Scale in Reservoir Rocks
Flow and transport in porous media is a common phenomenon in nature and many areas of
science and engineering. Water transport form tree roots to its leaves and nutrient or contaminate
transport in soil are examples of transport in natural porous systems. For oil to migrate from a
source rock to the, so-called, reservoir rock it has to transport through porous rocks. Primary and
secondary oil recovery processes are also examples of flow and transport in porous media. The
main challenge in reservoir characterization is high spatial heterogeneity which makes this class
of problems very scale dependent. In other words, depending on the size of a subsample under
study a new picture of the system is commonly extracted. Heterogeneities can be classified into
four different categories (Haldorsen and Lake, 1984):
1. Microscopic heterogeneities that exist at the level of the pore surfaces and grains
and can be characterized using Scanning Electron Microscopy (SEM), X-ray
computed tomography (CT scan), thin sections, etc. Characterization and
simulations at this scale, also called pore scale modeling, is the focus of this thesis.
2. Macroscopic heterogeneities that are observed at the core scale where, for example,
two cores that are taken from the same well in close proximity might show
completely different flow properties. There might be heterogeneities in a single
core itself. Different parts of the core may belong to different lithologies with
different minerals and transport characteristics. To understand this variability,
however, we should use microscopic scale characterization tools. Lab experiments
are traditionally performed on core samples.
2

3. Megascopic heterogeneities that are at the reservoir scale and include faults as well
as large fractures. A reservoir is an infinitely large number of cores and its behavior
depends on flow properties of cores and their spatial distribution. In addition, there
are heterogeneities that only appear at this scale and can change the overall
behavior of the flow system drastically. For example, a fault between injection and
production wells affects the water-flood performance by reducing the sweep
efficiency and recovery factor. This scale is the reservoir engineering scale where
models of the whole reservoir are built and oil in place, well placement, recovery
factors, etc., are studied.
4. Gigascopic heterogeneities that appear at the landscape scale, where reservoirs are
studied in the geologic settings relative to rivers, mountains, etc. The
characterization at this scale is considered during the exploration phase where
seismic techniques are utilized to image the underground rocks.
In macroscopic (and larger) scales, fluid flow is described by Darcy’s law (Scheidegger,
1957). In a one-dimensional horizontal system, Darcy’s law stat es that the average fluid velocity
is linearly proportional to the pressure gradient (

). The constant of proportionality is
called permeability ( ) and is, in fact, a measure of the effective (open to flow) area. Darcy’s
law, originally an empirical equation, proposed based on experimental results of laminar flow in
sand packs (Darcy, 1856). However, for simple geometries with very low Reynolds numbers
(creeping flow), this relationship can be derived from Navier-Stokes equations (Sahimi, 2012)
and (Bear, 2013).
3

Like many other empirical equations (e.g. Newton’s law of cooling), Darcy’s law is a
simple equation with an inherently complicated constant, k. Permeability, especially for two
phase systems, is a complex function of the rock pore structure and specific flow settings.
In addition to experimental measurements, macroscopic properties can be calculated using
simulations at the microscopic level (Blunt and King, 1990). For pore-scale simulations, the
detailed pore structure is required. The necessary imaging resolution and hence imaging
technique varies for different rock types. Differential equations for momentum, mass and energy
conservation, with appropriate initial and boundary conditions can then be solved on the rock
image. Macroscopic properties are calculated by averaging the microscopic simulation results.
The domain size in pore-scale simulation should be statistically significant for the averaging to
be meaningful. The domain size where variability in macroscopic properties disappears is called
the Representative Elementary Volume (REV) (Bear, 213). However, the main challenge to this
end is the pore size dependency of the governing physics in the nanometer pores that are
observed in unconventional systems (Darabi et al., 2012).

1.2 Pore-Scale Modeling and Simulation
During the last 50 years, different techniques have been developed for simulation of flow and
transport at the pore scale. A review of these techniques and our proposed workflow is presented
in this section.

1.2.1 Pore Network Modeling
Researchers from different disciplines have used simple geometrical configurations as
approximate and simple models for flow in porous media. The first model of this kind was the
4

‘bundle of capillary tubes’ model that was used to calculate macroscale flow coefficients. These
simple models do not capture the full interconnectivity and tortuous nature of porous materials
and thus more complex 2D and 3D networks were introduced to better represent the porous
structures (Fatt, 1956). In these models, porous media is approximated by a collection of pores
and throats. Pores are bodies with large volumes and are connected through throats (with very
small volume). Pores control the porosity and residual saturations while throats control the flow
conductivity of the system. Figure 1.1 shows the pore network approximation of a Brea
sandstone (Piri, 2003).

Figure 0.1. Pore network for a Brea sandstone (Piri, 2003)

Rules are developed to determine the two-phase configurations in these idealized
geometries. These rules are then applied to the pore network model to study its two-phase flow
and capillary pressure characteristics (Blunt et al., 2002). Pore network models are
computationally less intensive than direct calculations on rock images. However, finding a pore
network that is equivalent to the actual rock images is a challenging task.

5

1.2.2 Direct Calculations
In direct calculations, the rock image is the input geometry to flow simulations and each voxel is
a computational grid. This problem can be addressed in computational fluid dynamics (CFD)
framework. Conventional (Navier-Stokes equations) and modern (Lattice Boltzmann Method)
CFD techniques are discussed in the following sections.

1.2.2.1 Conventional CFD Techniques
Rahimain and Pourshaghaghy (2002) solved Navier-Stokes equations (using marker and cell
algorithm) on a random porous media with a specified porosity (Figure 1.2). They divided the
channel into three parts and placed stationary solid particles in the second part, no-slippage (zero
velocity) boundary condition were applied in the fluid-solid boundaries.

Figure 0.2. Velocity profile is disturbed by solid particles, Re=10 (Rahimain and
Pourshaghaghy, 2002)

Flow in natural porous media is surface dominated and Reynolds numbers are usually very
low (

). In this case the non-linear terms in the original Navier-Stokes equation can
be neglected and the equation reduces to Stokes form. Stokes equation is linear and its solution
does not require iterations. Adler (1990) used this method to calculate permebility in
Fontainebleau sandstone.
6

1.2.2.2 The Lattice Boltzmann Method
The Lattice Boltzmann Method (Qian, Humiéres and Lallemand, 1992; Benzi, Succi and
Vergassola, 1992; Chen and Doolen, 1998; Succi, 2001 and Sukop and Throne, 2005) is a
mesoscopic method in computational fluid dynamics and represents a discrete form of the
Boltzmann equation near equilibrium. In this method, the fluid is modeled with a distribution
function of particles that can move on specified lattice sites. During each lattice time step,
particles move to their neighboring sites and exchange momentum through collisions. At each
time step external forces and boundary conditions are applied. LBM can be derived directly from
the Boltzmann equation (He and Lou, 1997), however historically it is the evolved version of the
Lattice Gas Automata (LGA).
The LGA became popular in fluid dynamics when researchers realized that fluids could be
modeled with discrete particles (Pomeau and Frisch, 1986). In this method, the fluid is modeled
as particles that can rest and move on uniform grids. Different LGA schemes for single phase
(Rothman and Zaleski, 2004) and multiphase (Rothman and Zaleski, 1994) systems were
developed. There were some problems associated with the initial LGA methods. Most of these
problems were solved by introducing non-linear (McNamara and Zanetti, 1988) and linear
(Higuera and Jimenez, 1989) Boltzmann collision operators. In addition, the Galilean invariance
1

of the macroscopic behavior was satisfied by defining zero velocities (rest particles) and the
Maxwellian distribution for equilibrium (Higuera, Succi and Benzi, 1989). The simplest and
most widely used collision operator is BGK (after Bhatnagar, Gross and Krook), which is a
linear operator that is based on a single time relaxation parameter towards equilibrium (Qian,

1
“Galilean invariance or Galilean relativity states that the laws of motion are the same in all inertial frames.”
7

d’Humiéres and Lallemand , 1992). Multiple relaxation time LBM schemes were introduced
(d’Humiéres, 2002) to improve the numerical stability.
LGA (and later LBM) methods found to be very appealing for studying flow in porous
media. Balasubramanian, Hayot and Saam (1987) studied the validity of Darcy’s law for a 2D
random porous media. Succi et al., (1989) published one of the first studies on the applicability
of LBM for flow modeling in 3D complex structures. Later, researchers used single phase LBM
to study flow in spherical packs (Cancelliere et al., 1990), 3D fractal geometry (Chen et al.,
1991) and Fontainebleau sandstone (Ferréol and Rothman, 1990). After introducing two-phase
models in LBM (Shan and Chen, 1993), the approach was used to model two-phase phenomena
in porous media (Pan et al., 2004).

1.3 Pore Scale Characterization Workflow
In this section, we propose a workflow that takes 3D rock images as input and calculates
properties at different levels (Figure 1.3).
8

Figure 0.3. Pore Scale Characterization Workflow

As we will explain below, zero level properties refer to the properties that can be extracted
from the binary images (e.g. porosity) while 1
st
level properties are the ones that can be
calculated from single phase flow simulations (e.g. permeability). 2
nd
level properties are the
results of two-phase flow simulations (e.g. relative permeability curves).

1.3.1 Rock Imaging
Three dimensional imaging techniques can be categorized into two groups. Techniques that scan
the sample by propagating energy wave through its volume, for example X-ray CT. These
techniques are non-destructive. However, the sample needs to be perfectly symmetrical and a
3D-Imaging
Image Analysis

Multiphase
(Minerals-Pore)
Two-Phase
(Solid-Pore)
Zero Level Properties
(Image Analysis)
1
st
Level Properties
(Single Phase Flow Simulation)
2
nd
Level Properties
(Two-Phase Flow Simulation)
9

high resolution imaging requires a high energy instrument. The images, in this case, have the
same resolution in all three directions. In the second group of imaging techniques, 3D images are
generated from a stack of 2D images. Focused Ion Beam - Scanning Electron Microscopy (FIB-
SEM) belong to this category where 3D images are constructed from surface SEM images (2D)
that are exposed by cutting/milling the sample using the ion beam. The distance between
different 2D layers can be as small as a few nanometers and the resulting 3D image can
accordingly have the same resolution in all three directions.
Rock imaging is specially employed in characterization of unconventional resources where
most of the conventional flow-through experiments are not applicable. Lemmens and Butcher
(2011) presented a workflow to characterize shale from micron to nano scale using SEM/EDX
and FIB-SEM. Curtis et al. (2010) used FIB-SEM imaging to conceive the pore size distribution
and their results were in good agreement with the data obtained from conventional
characterization techniques like Mercury Injection Capillary Pressure (MICP) and Nuclear
Magnetic Resonance (NMR). They also adopted Scanning Transmission Electron Microscopy
(STEM) imaging to calculate the porosity from a 2D image with pore sizes that are smaller than
those captured by the FIB-SEM resolution. Riepe et al. (2011) used X-ray μCT imaging and pore
network modeling (PNM) to calculate porosity and permeability and the results were in close to
the routine core analysis techniques. Elgmati et al. (2011) found good agreement between
calculated total organic content (from 2D SEM images) and measurements in shale gas samples.
A 3D image can be imagined as a stack of 2D images. An image is a matrix with an
assigned value for each pixel. For grayscale digital images, this value is between 0 and 255.
Figure 1.4 shows the image of a carbonate rock.
10

Figure 0.4. Gray scale image of a carbonate rock. The L.H.S. matrix is the number
representation of the highlighted portion of the picture. 0 is black (pore) and 255 is white
(matrix). (Courtesy of Dr. Masa Prodanovic)

Image segmentation refers to the task of defining different classes of values and assigning
a phase to each class. The phase can be a mineral and its corresponding range of values is
obtained from images of the pure mineral. For our purpose, we are mainly interested in two
phase, matrix or pore, segmentation. For this purpose, only one threshold value is required. The
pixels with values above that threshold are defined as solid and the ones below the threshold are
assigned to be pores. The porosity observed from experimental measurements, is commonly used
to select the threshold value. Figure 1.5 illustrates this concept.

11

I- Gray Scale

II- Threshold = 159 (Porosity = 0.27)

III- Threshold = 149
Porosity = 0.18

IV- Threshold = 169
Porosity = 0.39
Figure 0.5. Gray scale image with binary images at different thresholds with the corresponding
calculated porosities. The threshold value is adjusted from porosity experiments.

Skeletenization is a technique to extract the skeleton of the pore space. We use this
technique to determine if there is a flow path connecting the inflow and outflow faces. If there is
no connected porosity in the porous sample, the permeability will be zero. At this stage we have
a digital rock sample (a 3D matrix with 0 for pores and 1 for solids). From the binary image,
properties at different levels can be calculated (Figure 1.3).
12

Zero level properties are calculated from the binary image and can be categorized into
porosity classification, specific surface area calculation, pore size distribution and connectivity
determination (zero or non-zero permeability). The total porosity is, simply, calculated by
counting the pore voxels. The pores that are not in the cluster of the connected skeleton are no-
flow porosity. If there are micro-fractures in the image, they can be quantified based on their
pore sizes and micro-fracture porosity is calculated.
Specific Surface Area (SSA) is calculated by counting the areas of the porous material
where the pores are in contact with solid voxels. SSA has effective and ineffective components.
Effective SSA is the surface area of the connected pores that contribute to flow while ineffective
SSA is due to ineffective (isolated) porosity. Effective SSA is a measure of the resistance of the
porous material to flow. Generally speaking, the higher the effective SSA is, the lower the
permeability will be. SSA is an important factor in surface controlled phenomena like acidizing
for carbonate rocks and adsorption/desorption in shale gas systems. The pore-throat data can also
be extracted from 3D binary images using pore-network algorithms (Dong and Blunt, 2009). If
the skeleton of pores has at least one cluster connecting the inflow and outflow, permeability is
non-zero and one can proceed with the flow simulations.
1
st
level properties are calculated from the single phase flow field. Equation (1.1) describes
Darcy’s law in the vector rom.

is the velocity vector (average velocities of the pore scale) and

is the flow potential. The permeability is a tensor and is defined by Equation (1.2) in three
dimensional Cartesian coordinates. As it is explained in chapter 3, to mimic the flow
experiments, we close four faces of the digital sample and apply a pressure gradient to calculate
permeability in a specific direction. Due to compaction in sedimentary layers, usually the vertical
permeability is less than horizontal permeability.
13

(0.1)

(0.2)

The pores with very small velocity (compared to the average velocity in the domain) can
also be counted as ineffective porosity. At this stage, sensitivity analysis of permeability to
different pore sizes can be performed. For example, one can close all the pores except the micro-
fractures and calculate fracture permeability.
2
nd
level properties are two-phase properties. Steady-state relative permeability curves are
the result of two-phase flow simulations at different saturations. This workflow provides an
opportunity for the microscopic engineering of reservoir rocks. Recovery factor of a rock sample
under a specific process can be calculated and with the information of saturation distribution one
can decide on the most feasible EOR technique.

1.4 Areas of Application
In this section, we discuss different relevant geological settings and processes that can be studied
in more detail using the pore scale analysis and modeling.

1.4.1 Unconventional Resources (Physics at Play)
In recent years, both shale gas and shale oil have been predicted to potentially fill the increasing
void between the growing demand and declining supply of conventional oil and gas resources,
and are thus attracting a significant economic interest from the oil and gas industry. In
unconventional oil-shales, the oil is thought to be trapped in the pores of these tight rocks at the
14

sub-micron level. In order to extract this oil in a more economic and efficient manner, it is
important to understand the pore-size distribution and the pore connectivity that dictate fluid
flow and mass transfer at the sub-micron scale. Standard core analysis techniques are not always
applicable to unconventional systems. Measurement of liquid permeability in tight pore systems
is a very time consuming task, and when the aim is to extract e.g. oil-water relative permeability
data, direct measurements become even more tedious. For characterizing pore throats as small as
a few nanometers, mercury injection capillary pressure (MICP) experiments require pressures of
up to 60,000 psi (or higher); material integrity under these high injection pressures is a key
concern.
In shale oil reservoirs, due to low permeability of the rock (usually ),
measurement of liquid permeability is extremely time-consuming and gas permeability
experiments are conducted instead. Transforming gas permeability in sub-micron pore systems
into a representative liquid permeability is not well established yet. In the context of the
proposed workflow, we simulate the flow of an incompressible liquid on the images of Monterey
shale oil sample and directly calculate liquid permeability.
There are many settings where oil and water coexist in the shale systems. During hydraulic
fracturing, large volumes of water, with some additives, are injected to create fractures in the
rock. This water can then be imbibed into the shale rock and block the flow path for the
hydrocarbon in the subsequent production stage. Furthermore, since we are producing from the
source rock, the chance of having partially charged reservoirs is higher. In this case, oil and
water will be both flowing and we might have high water cuts even in the early time of
production. During secondary recovery, like in conventional reservoirs, water injection will be
15

an interesting EOR option for unconventional systems. In this case, water will be displacing oil
in the sub-micron pores.
The simulation of two-phase flow in these rocks at the pore scale will allow us to better
understand imbibition of these systems at the basic pore level. In chapter 5, we discuss modeling
of imbibition at pore and core scale. This model is used for matching spontaneous imbibition
experiments in a shale oil sample.

1.4.2 Conventional Resources (Co- and Counter-Current Flows)
In the extension of Darcy’s law to two-phase flow settings (Equation (1.3)); it is assumed that the
only effect from the presence of one phase on the flow of the other phase is decreasing
permeability by introducing relative permeability (

). In this notation, the two phases are
hydro-dynamically decoupled and

does not move phase .

(0.3)

However, in some cases viscous coupling between the fluids becomes important. Javaheri
and Jessen (2013) investigated the significance of considering counter-current effects (fluids
flow in two different directions) during CO
2
storage in aquifers. They used two different relative
permeability models for co- and counter-current flow regimes and chose the appropriate rel-perm
model based on the flow directions in each block at each time step. Yiotis et al. (2006)
considered the effect of viscosity ratio on the relative permeability curves in 2D and 3D simple
porous geometry using LBM while Hunag and Haibu (2009) used Shan and Chen two-phase
LBM to study the effect of flow direction on relative permeability of 2D simple structures.
16

For 2D porous materials, we calculate the relative permeability curves for both co- and
counter-current settings. Flow in the opposite directions reduces the relative permeability.

1.5 Organization of the Manuscript
A summary on different chapters of this manuscript is discussed in this section.
Chapter 2 discusses the rock digitization process. Different characterization techniques are
applied on samples from the oil-rich Monterey shale to understand the volumes associated with
each pore size. MICP, BET and HRXRCTs (0.5 µm and 2 µm) are performed. FIB-SEM is
carried out at a resolution of 30 nm and is used for calculating permeability in the shale oil
sample. Spontaneous imbibition experiment in Monterey shale sample and its interpretation is
presented.
Chapter 3 provides the theory of LBM method that we use for single- and two-phase flow
simulations. The Boltzmann equation with the appropriate collision operator is introduced. The
Maxwell-Boltzmann equilibrium distribution function is discussed and a two dimensional LBM
scheme is derived from the Boltzmann equation with the weights, etc. LBM is then used for
calculating directional permeability on a 3D rock image.
Chapter 4 discusses the extension of LBM for multiphase systems. Different multiphase
LBM schemes are introduced and the Shan and Chen formulation is discussed in more detail.
Surface area minimization and Laplace law are studied as benchmarking problems. The solid-
fluid force implementation is then used for modeling wetting phenomena in 2D and 3D. Viscous
coupling in layered two phase flow in 2D channels is solved analytically and the results are
compared with LBM simulation results for different viscosity ratios. Viscous coupling in the
context of relative permeability in a 2D porous media is discussed.
17

In chapter 5, we present a theoretical analysis of capillary rise dynamics in a single
capillary tube. We extend the formulation to describe homogenous porous media and to include
viscosity, pressure and gravity effects in both advancing and defending fluids. Finally, we extend
the theory to describe fractured systems, and propose a convolution integral formulation and a
new explicit solution.
Chapter 6 suggests the potential future work and research directions in this area.

18

2. Chapter Two: Characterization of Shale Samples (from Rocks to Numbers)

2.1 Introduction
In this chapter, we compare and discuss the insight and understanding of the pore structure
characteristics of Monterey shale samples using various characterization techniques. The
techniques considered here include: (i) Mercury injection capillary pressure (MICP)
experiments, (ii) Brunauer, Emmett and Teller (BET) – Nitrogen adsorption experiments, (iii)
High-Resolution X-ray CT (HRXCT), and (iv) Focused Ion Beam – Scanning Electron
Microscopy (FIB-SEM). The HRXCT and FIB-SEM images are utilized to calculate the
permeabilities of these shale samples via Lattice Boltzmann flow simulations. HRXCT data
indicates that pores larger than 0.5 µm represent about one third of the sample’s porosity.
Furthermore, a large portion of the pores, mainly in the clay structure, appear to be isolated and
not to contribute to the flow. Calculated permeability values, based on image analysis, are in
good agreement with experimental observations.
Production of hydrocarbons from unconventional source rocks (oil/gas shales) has
increased significantly over the past decade due to a shift in the focus of the oil and gas industry
from conventional systems towards unconventional oil/gas resources. The total domestic crude
oil production has increased from 5.66 million barrels per day (bpd) in 2011 to 8.46 million bpd
in 2014, and this is largely due to the expansion of shale oil operations (EIA 2014). As noted in a
2012 report from the Energy Information Administration (EIA), there are presently 33.2 billion
barrels of technically recoverable crude oil in tight formations in the USA, the recovery of which
will require the drilling of ~220,000 new wells. The Monterey formation holds about 40% of
these resources (13.7 billion barrels) in highly heterogeneous and laminated rocks. An improved
19

understanding of the heterogeneities at the micro-scale in these rocks will enable us to
understand and optimize the production behavior from this vast oil resource.
The main distinguishing feature of shale oil rock samples is their small (sub-micron) pore
sizes. The distribution and connectivity of pores in a rock determines its overall fluid flow
behavior, for example, its single-phase and multi-phase permeabilities. The small pore sizes
make most of the standard core analyses and tests, which were developed for conventional
systems with much larger pores, difficult to apply for shale systems. Measurements of liquid
permeability and especially of two-phase relative permeability become tedious (at best) in these
low permeability systems. Mercury injection capillary pressure (MICP) tests are routinely
applied to measure the pore-size distribution (PSD) of conventional rock samples. However, for
mercury to invade the small pores in an unconventional sample requires a very high injection
pressure (up to 60,000 psi), which may potentially cause the rock pore structure to deform and
change.
More recently, researchers are using pore-scale imaging of rocks as a new core analysis
technique (Knackstedt et al., 2009). This technique, which is often referred to as digital rock
physics, is specifically useful for shale oil and shale gas samples. Curtis et al. (2010) utilized
Focused Ion Beam – Scanning Electron Microscopy (FIB-SEM) imaging to extract the pore sizes
from shale samples, and found them to be dimensionally consistent with the pore sizes obtained
from MICP and NMR measurements. In addition, they used high resolution scanning
transmission electron microscopy (STEM) imaging to capture the details of porosities in clay and
kerogen. Riepe et al. (2011) used 3D micro-tomograms at different scales combined with pore-
network modeling to compute petrophysical properties (e.g., porosity, permeability, formation
resistivity factor, and drainage capillary pressure curves), and demonstrated a good agreement
20

between calculations and experimental observations. Elgmati et al. (2011) applied MICP and
FIB-SEM on shale gas samples to investigate porosity and PSD. They found the MICP porosities
to be 10.34% and 14.56% with median pore diameters of 6.5 nm and 30 nm, respectively, from
two separate tests. They also calculated a porosity of 28.22% with a mean pore diameter of 40
nm from the FIB-SEM 3D pore modeling. Their SEM image revealed a kerogen porosity of 40-
50%. In addition, they used X-ray diffraction (XRD) analysis to quantify the clay content of their
samples, which were found to contain a high fraction of illite. However, they were not able to
perform comparison studies due to sample limitations.
In this chapter, we have used different techniques in order to take a closer look at a
collection of shale oil rock samples from a specific depth in the Monterey formation and to study
their pore structure characteristics. We use the Brunauer, Emmett and Teller (BET) nitrogen
adsorption technique to characterize pores as small as 1 nm. High-Resolution X-ray CT
(HRXCT) images, at two different resolutions (0.5 µm and 2 µm), are then used to investigate
pore connectivity for the relevant range of pore sizes. FIB-SEM is used for 3D imaging of rock
samples at a 30 nm resolution. Pore-connectivity data are then extracted from these 3D images
and are compared to mercury injection capillary pressure (MICP) data. The 3D binary images
from the 0.5 micron HRXRCT and the FIB-SEM analyses are used to calculate the sample
permeability using the Lattice Boltzmann Method (LBM) (Guodong et al., 2004). Finally, we use
the experimental data from spontaneous imbibition experiments to estimate porosity and
permeability values.

21

2.2 Porosity Characterization Techniques

2.2.1 Mercury Porosimetry (MICP)
The mercury injection capillary pressure method is a well-established technique for measurement
of volumes associated with different pore sizes. It is a relatively fast and inexpensive approach to
study petrophysical attributes of reservoir rocks (Shafer and Neasham, 2000). The total volume
of sample invaded provides a measure of its porosity, while a cross-plot of saturation vs.
injection pressure provides valuable information related to the PSD. Several empirical equations
have been proposed to date to estimate permeability from MICP data in both conventional
(Swanson, 1981) and tight-pore systems (Comisky et al., 2007). Furthermore, Li and Horne
(2001) have used MICP data to generate (two-phase) relative permeability curves for steam-
water systems.
During the MICP test, a clean sample is placed in a mercury container whose injection
pressure is gradually increased. At each pressure step, the mercury volume that has invaded the
rock sample is recorded. As the pressure increases, smaller and smaller size pores are
progressively invaded. We note here, that during the MICP test mercury is the non-wetting phase
and, hence, the MICP data reflect a drainage process. The mercury pressure and the
corresponding invaded pore size are related through the Young-Laplace equation:
(2.1)
where r=pore radius, σ=surface tension, θ=contact angle between mercury and the rock surface,
and p=absolute injection pressure. To invade pores as small as a few nanometers, the above
equation indicates that injection pressure must be increased to about 60,000 psi. Our MICP test
results for a Monterey siliceous shale sample are shown in Figure 2.1.
22

Figure 2.1. Mercury intrusion data for a quartz-phase Monterey sample.

From Figure 2.1, one observes that the technique detects some larger pores that are usually
due to surface roughness and should be excluded from the data. After correcting for this, the
maximum pore size is detected as 0.86 µm. The smallest pores that are probed at 60,000 psi of
mercury pressure are 3.6 nm in diameter, and we find that half of the porosity is due to pores
with diameters smaller than 0.22 µm. From the data reported in Figure 2.1, one can evaluate the
Swanson permeability (Swanson, 1981) to be 0.21 mD.

2.2.2 BET Analysis
Another method that is commonly used to characterize porous materials, including shale, is the
BET technique. In contrast to the MICP method, BET is a low pressure gas adsorption technique
in which the surface area, pore volume and the average pore size are evaluated based on the
physical adsorption and desorption of probe gas molecules, commonly N
2
and CO
2
, on the pore
surface (Adesida et al., 2011 and Clarkson et al., 2013). In this work, liquid N
2
(temperature=77
23

K) sorption was carried out on an intact sample of shale, and the data are used to estimate the
surface area and to characterize the PSD via BET analysis (Brunauer et al., 1938).
Several approaches, including the Barrett-Joyner-Halenda (BJH), Horvath-Kawazoe (HK),
the t-plot methods, and density functional theory (DFT), are all commonly used to analyze the
experimental isotherm data in order to determine the sample’s pore structure characteristics. In
this work, the BJH method was used to calculate the PSD in the mesoporous region (2-50 nm)
assuming a cylindrical pore geometry (Barrett et al., 1951), while the HK method was used to
characterize the pores in the microporous region (<2nm), assuming a slit-like pore geometry
(Kawazoe, 1983).
The total BET specific surface area for the sample was calculated to be 6.624 m
2
/g. Figure
2.2 shows the pore diameter as a function of the fractional cumulative pore volume evaluated
from the BJH analysis of the experimental isotherm. The total mesoporous volume is calculated
to be 0.027cm
3
/g, with an average pore diameter of 17 nm.

Figure 2.2. BJH cumulative pore volume data for the shale sample.
24

Figure 2.3 shows the average pore diameter as a function of the cumulative pore volume as
obtained from HK analysis of the experimental isotherm. The total microporous volume and the
average pore diameter are calculated to be 2.3 e-3 cm
3
/g and 1.15 nm, respectively.

Figure 2.3. HK cumulative pore volume data for the shale sample.

Based on the sample bulk volume of 0.52 cc, the sample porosity was calculated to be
5.1% in the micro and the mesoporous regions (<45 nm).

2.2.3 Imaging and Image Analysis
Digital rock evaluation or characterization of shale samples using various imaging techniques is
today becoming a routine approach to capture and analyze the submicron features of rocks.
Properties like pore connectivity, liquid permeability, and porosity can all be calculated from a
3D image of a given volume of the rock. One such nondestructive imaging technique is high-
resolution X-ray computed tomography. It uses a high-energy beam to detect the sub-micron
25

features of dense materials such as shales, which is a substantially greater resolution than that
attained by conventional medical CAT instruments. For this study, a shale sample was imaged at
two different resolutions of 2 μm and 0.5 μm, respectively, using a Scanco Medical μCT50
instrument.

2.2.3.1 2-μm HRXRCT Imaging
For this study, 945 2D slices were imaged at a 2-μm resolution using a sample from the same
depth as the sample used in the MICP and the BET experiments. For illustration purposes, Figure
2.4 shows one of the 2D slices. The two main features observed in this image are: (1) The darker
or black features that represent the pores, and (2) the lighter or white features that represent high-
density minerals. Based on XRD data for samples from the same depth (not discussed in this
paper, but available to our team), the lighter features are believed to be pyrite. The 3D structure
constructed from 2D slices was segmented in order to quantify the features discussed above.

Figure 2.4. Vugs/pores (dark features) and high-density regions (white features).

As seen in Figure 2.5, the pyrite structures visible at a 2-µm resolution, do not appear to
form a connected network in the reconstructed 3D image. Pyrites may, however, still be
26

connected through smaller features not detected by HRXRCT at this resolution (it should be
noted, that due to its high conductivity an interconnected network of pyrite can be
misinterpreted, from conventional well-logging methods, as water/brine saturation). After the
segmentation, the volume of pyrite inclusions of size 2 µm or larger were calculated, using
computations analysis, to be 0.5% of the total sample volume.

Figure 2.5. Spatial distribution of pyrite from imaging at 2-µm resolution.

Image segmentation algorithms were also applied to identify pores from the 2-µm
HRXRCT images. We use the maximal-ball algorithm to extract the pore network from the
images (Dong, 2009). Figure 2.6 shows the results from the analysis of the 2-µm resolution
images. We note, that pores with sizes <2 µm are artifacts of the analysis due to unsmoothed
edges in the images.

27

Figure 2.6. Pore size distribution from 2-µm HRXCT imaging.

The porosity from this image is approximately 4%. The pore network, at this resolution,
does not percolate in any direction and as a result the sample permeability is zero. Higher
resolution images, however, presented in the next section, reveal that these large pores (>2 µm)
are connected through smaller pores.

2.2.3.2 0.5-µm HRXRCT Imaging
The HRXCT imaging at a 2-µm resolution, presented in the previous section, captures all the
pores that are larger than 2μm in size. The same sample was subsequently modified (to fit the
sample holder) and imaged at a 0.5-µm resolution. The outcome of such imaging analysis was
157 slices (images), each corresponding to an area of 0.81 mm
2
. An “object” was then defined as
a set of pores that form a connected flow path in 3D space; based on this definition, we detected
256,696 such “objects” in the 0.5-µm 3D image (composed of the 157 2D images, as noted
28

previously). However, most of these “objects” are relatively small in size and isolated, i.e., not
fully spanning any of the sample’s directions. From all 256,696 such “objects”, only eight were
connected (i.e., spanning) in the z-direction, while only one “object” was found to be connected
(spanning) in all three (xyz) directions. The fact that more “objects” span the entire thickness in
the z-direction is most likely due to the fact that the 3D image is thinner in the z-direction as
compared to the x- and y-directions.
The total porosity of the 0.5-µm resolution image was calculated to be 13.7%, while the
connected porosity was 10.7%. Approximately 50% of the pores detected by the 0.5-μm
resolution imaging lie in the size-range below 2μm, as shown in Figure 2.7.

Figure 2.7. Pore size distribution from the 0.5-µm HRXCT imaging.

2.2.3.3 SEM and FIB-SEM
The HRXCT imaging discussed in the above section can achieve a resolution of 0.5 µm. In order
to capture the pores and other rock features in the size range of a few nm (micropores and
mesopores), the shale sample was imaged using a SEM instrument. A sample with an
29

approximate cross-sectional area of (1 cm x 1 cm) and a thickness of 2 mm was prepared for
SEM imaging. Rock samples, such as those studied here, are known to be insulating and, as a
result, charging effects may introduce imaging artifacts. The sample’s surface was, therefore,
coated with gold to reduce the potential charging effects. Imaging was performed using a JEOL
JSM-7001F-LV instrument at various magnifications. One such image of the sample’s surface at
65,000x magnification can be seen below in Figure 2.8. The sample displays randomly scattered
porous areas throughout its surface, with pore diameters <50 nm clearly visible in Figure 2.8.

Figure 2.8. High-resolution surface SEM image with visible nano pores.

To capture the connectivity of such pores, a 3D imaging of the porous material is required.
The technique we have used to accomplish this is FIB-SEM that produces a stack of SEM
images, ideally spaced with the same distance as the SEM resolution to represent the rock in 3D
at high resolution (Bera et al., 2011). Figure 2.9 shows a subset (i.e., several 2D slices) of a FIB-
SEM 3D image from a neighboring sample. The spatial resolution in this image is 30 nm in all
30

directions. The size of this sample is 49 μm x 27 μm x 17 μm. Sample porosity is 40% and the
image percolates in all three directions. We note that the FIB-SEM porosity is high compared to
the porosity (~30%) that we have measured with the same sample using He pycnometry. The
difference can be explained by the fact that usually more porous regions of the rock are milled
for FIB-SEM imaging.

Figure 2.9. 2D slices of the FIB-SEM image at a resolution of 30 nm.

Figure 2.10 reports the pore diameters and their contribution to the total pore volume, as
extracted from the 3D image analysis. The plot shows a good agreement with the PSD extracted
from the 0.5-µm HRXCT imaging (Figure 2.7). The same is also true when comparing calculated
permeabilities using these two different images, see discussion in the next section.
31

Figure 2.10. Pore size distribution from FIB-SEM imaging.

2.3 Permeability from Flow Calculations
The binary images from the 0.5-µm HRXRCT and the FIB-SEM analyses are used for flow
simulations using the Lattice Boltzmann Method (LBM). The Lattice Boltzmann Method (as will
be explained in more detail in chapter 3) is a mesoscopic computational fluid dynamics method
and represents a discrete form of the Boltzmann equation near equilibrium.
In this application, the lattice is often selected with its points being identical to voxels in
the rock image (Guodong, 2004). To calculate the permeability in the transport (x) direction, the
inlet and outlet faces are open and all the other faces are closed with no-flow boundaries (this
closely approximates the experimental configuration). A periodic boundary condition with a
body force is then applied in the flow direction. The LBM simulation results in a flow field
(pressure and velocity) that we use (average) to calculate the permeability from Darcy’s
equation. Figure 2.11 shows the rock binary image and the x direction component of velocity
calculated using the LBM. For the 0.5-µm HRXRCT image, the simulation was performed with
32

a sub-sample of (500 x 500 x 157) voxels and the permeability was calculated to be 0.961 mD.
For the FIB-SEM image, the whole sample was used in the flow calculations, and the
permeability was found to be equal to 5.2 mD. The higher permeability in the FIB-SEM sample
(30 nm resolution) can be attributed to the presence of the pores smaller than 0.5 µm which were
not captured in the HRXRCT image (and thus not accounted for in the permeability simulations
using that image).

Figure 2.11. (a) Solid phase, (b) pore space and (c) velocity in the x direction. The domain size
is 100
3
voxel (for illustration).

2.4 Spontaneous Imbibition Experiments and Simulations
Spontaneous imbibition experiments have long been used to characterize various rock properties
such as the rock’s wettability, capillary pressure and the porosities (fracture and matrix)
accessible to fluids (Ma et al., 1999 and Roychaudhuri et al., 2013 and Takahashi and Kovscek,
33

2010). A water-air spontaneous imbibition experiment was performed on a ~1.4 cc shale sample,
at the same depth as the samples used in previous sections, with all sides open during the
imbibition process. The total porosity (

) of this sample was calculated to be 31%, based on the
grain density and the sample weight. During the experiment, the sample was completely
submerged in water. Water (the predominant wetting-phase) invasion from all sides leaves no
room for the air to escape, thus at the end of the imbibition process some of the air will be most
likely trapped inside the matrix. We have included the air pressure (an average between initial
and anticipated final) in the modeling effort reported in Figure 2.12.

Figure 2.12. Water uptake (cc) vs. square root of time (hours) over the 23 hours of the
spontaneous imbibition experiment.

The experimental data is bilinear in a square root time plot, which is characteristic of
fractured systems. The first line, faster imbibition, is attributed to the larger pores in the siliceous
structure while the second line is due to clay porosity. This observation confirms our findings in
the previous sections that we have pores as small as few nanometers (BET analysis) as well as
34

pores in the few hundred nanometers range (from FIBSEM analysis for example). From the BET
measurements, we concluded that about 17% of the pore volume was due to clay (<50 nm)
porosity which matches the clay contribution in the imbibition model (see Figure 2.12).
.To match the experimental observations, we use a model of capillary rise in fractured
porous media (Chapter 5). In this model, the siliceous structure is assumed to be fractures while
the clay content represents the matrix portion of the porous system. The model parameters
include total porosity of the system (ϕ
t
), ratio of matrix volume to total volume (R
m
) , the fracture
opening (w), fracture and matrix porosity (ϕ
f
and

ϕ
m
), fracture and matrix permeability (K
f
and
K
m
) and an average weighted pore size of E[r] for both matrix and fracture systems. Note that
the model considers both the fracture and the matrix as porous structures. The value of
parameters that is used for the match in Figure 2.12 is summarized in Table 2.1.

35

Table 2.1 - Parameters of the imbibition model.
Parameter Value Source
Fixed model parameters:
K
f
0.4 mD Lab measurement (gas flow)
E[r
f
] 0.39 μm FIB-SEM
E[r
m
] 12.5 nm BET
P
air
1.45 atm Average trapped air pressure

0.31 Experimentally measured
Adjusted model parameters:

0.38 Matched to exp. observations
K
m
2.0e-9 mD Matched to exp. observations
W 0.027 cm Matched to exp. observations
Θ (55
0
) 66
0
Lab measurement, adjusted from static to dynamic, and matched

The value and source of parameters that are used in the model to represent the experimental
observations (Figure 2.12) are summarized in Table 2.1. We define the values of five of the
parameters: P
air
, is taken from the experiments, as an average between the initial and the final
pressure, while the total porosity ϕ
t
is measured experimentally. Fracture permeability (K
f
) is
assumed to be equal to the experimental gas permeability, measured from an adjacent sample.
The average pore size for the fracture and matrix systems are calculated, respectively, from the
36

FIB-SEM and BET pore size distributions. R
m
is the ratio of matrix volume to the total sample
volume and can be calculated from the matrix porosity (from BET) and the absolute matrix
porosity

(which is here an adjustable parameter determined from fitting the experimental
data) from the simple correlation R
m

. The absolute fracture porosity (

) can
be calculated from R
m
,

and

(experimentally measured) as

. Model
adjustable parameters, which are calculated using the experimental data (see Table 2.1), include,
in addition to

the fracture opening (w) and the matrix permeability K
m
.
We have performed static contact angle measurements on these rocks. The average value for
static angle is 55
0
that is corresponding to advancing angle between 65
0
-70
0
(Figure 3.8 Lake,
1989). We found that a contact angle of 66
0
matches the data best. Note that the driving force for
the fracture portion (

) is very sensitive to contact angle values.
The calculated matrix permeability for the clay system is very low, consistent with values of clay
permeability in the literature (Neuzil, 1994) and as a result controls the long-term imbibition
behavior

2.5 Summary and Conclusions
Different techniques have been applied to characterize the pore size distribution and connectivity
of quartz-phase, Monterey shale samples. Figure 2.13 provides a summary of the findings from
the various characterization techniques in terms of the range of pore sizes and the value of
porosity measured.
37

Figure 2.13. Summary of porosity characterization data.

HRXRCT was used to image the samples at 2-μm and 0.5-μm resolution. The data indicate
that pores larger than 2 μm make up ~4% of the total pore volume, but do not form a connected
flow path. This then suggests that large pores are connected through smaller pores with
diameters less than 2 μm. This, in turn, suggests that the connecting smaller pores (throats) are
potentially those captured by MICP, and that the volume of larger pores is, incorrectly, counted
by this technique towards the volume of the smaller pores. The total porosity calculated based on
the 0.5-μm resolution images is ~14%, and the associated pores constitute a connected flow path;
however the porosity is still low compared to the experimental values using He pycnometry,
~30%, and thus there are still pores that cannot be seen at this resolution. FIB-SEM imaging at a
resolution of 60 nm indicates a porosity of ~40%. The average pore diameters from the 0.5
HRXRCT and FIB-SEM images are 1.27 μm and 0.79 μm, respectively. This difference is likely
due to the fact that smaller pores are detected in the FIB-SEM, which “weigh” the average pore
38

size towards the smaller values. BET was performed to characterize the pores in the range of
diameters from 1-45 nm, and the porosity of these pores was found to be ~5%. These pores
belong to the clay portion of the matrix, and in the presence of the siliceous flow backbone, most
likely, do not contribute significantly to the liquid flow. However, they can be quite significant
for gas flow in shale gas samples; BET, therefore, can improve the characterization of such shale
gas samples.
Our analysis of HRXCT data indicates that about one third of the sample’s porosity is
represented by pores that are larger than 0.5 µm. LBM flow simulations on the inter-connected
pore structures indicate a permeability of ~1 mD, which is in good agreement with other
experimental results (gas permeability measurements). The pore size distribution calculated from
the FIB-SEM images and the MICP experiments do not completely agree, however, which is
likely due to the assumptions made in the analysis of MICP data and the limitations imposed by
the image resolution. The BET analysis shows that the porosity from pores that are smaller than
45 nm (lower range of FIB-SEM resolution) corresponds to ~13% of the total porosity of ~40%.
Since the rock matrix is shaly-silica, this porosity is believed to be clay porosity, which does not
significantly contribute to flow in the presence of a connected siliceous pore structure.
Furthermore, the calculated permeability from FIB-SEM images indicates that, for this type of
rock, a resolution of 30 nm is good enough to capture the flow backbone of the shale sample.

39

3. Chapter Three: The Lattice Boltzmann Method

3.1 The Boltzmann Equation
The Boltzmann Equation is the cornerstone of Kinetic theory. Kinetic theory is a branch of
statistical mechanics that studies the dynamics of non-equilibrium systems and their relaxation
towards thermodynamic equilibrium (Huang, 1963).
Consider molecules in a container of volume and temperature (Figure 3.1). If the
size of molecule is much smaller than the mean interparticle separation

then the
molecules can be treated as point-like structureless particles and classical Newtonian equations
govern the dynamic of the system
1

(3.1)

(3.2)

Here, m is the particle mass,

is the position vector of particle ,

is the momentum and

is the force exerted on particle by other particles (inter-particle interaction) and external
fields (gravity for example). These equations can be solved for given initial and boundary
conditions. The only problem being that is extremely large e.g. 10
15
for a sample as small as a
salt grain. With the computational advances, these equations can be solved for

(billion
atom molecular dynamic simulations) (Kadau et al., 2006).

1
If this assumption is violated then quantum effects should be taken into account and we will be dealing with a Bose
gas or a Fermi gas

40

Figure 3.1. N molecules away from equilibrium state

However, the behavior of macroscopic systems is the statistical average of microscopic
dynamics and hence a probabilistic approach with less computation will suffice in most cases. In
statistical mechanics, the distribution function is introduced where is the probability
of finding a particle with momentum in position at time . The six
dimensional space is called phase space. The main goal of kinetic theory is
to find a distribution function for a specific molecular interaction.
Now, let’s write the equation of motion for the distribution function . Assume an element

in the phase space , the number of particles in this element at
time is

. At a later time , the number of molecules in the element

will be

, where

is the velocity of the
particle and

is the force exerted on a particle. Neglecting the effect of collisions:

(3.3)

41

This equality is called Liouville’s theorem and states that distribution function behaves as
an incompressible fluid (Kardar, 2007):
Assuming that is small:

and

(3.4)

Inserting the effect of collisions:

(3.5)

First order expansion of the left hand side yields:

(3.6)

3.2 Relaxation Time Approximation (RTA) for Collision Term
Finding a general form for the collision operator,

, in Equation (3.6) is challenging. In
the Lattice Boltzmann formulation, a Relaxation Time Approximation (RTA) is used which is
valid for the systems close to equilibrium (Bhatnagar et al., 1954).

(3.7)

Replacing in Equation (3.6)

(3.8)

Assuming that the system is at equilibrium, no net forces are applied and the distribution is
the same throughout the entire domain. Now if we perturb the system by adding some particles
then Equation (3.8) reduces to:
42

(3.9)

Solving this equation results in:

(3.10)

The solution indicates that using RTA as the collision operator, any small perturbation in
the system goes away exponentially with a characteristic (dimensionless) time .
In the final form of the Boltzmann equation the equilibrium distribution term (

)
appears. The distribution function at equilibrium is independent of time and location and is a
function of velocities only (

). Since

are independent variables
then

. Because of the importance of the equilibrium
distribution in the derivation of different LBM schemes, a summary of Maxwell-Boltzmann
derivation is presented in Appendix A.

(3.11)

3.3 From the Boltzmann Equation to the Lattice Boltzmann Equation:
Macroscopic properties of a system are the moments of the probability distribution function:

(3.12 a)

(3.12 b)

(3.12 c)

43

Equations (3.8), (3.11) and (3.12) form a closed set and can be solved to find macroscopic
properties of the system. However, even at this stage, the problem is computationally expensive
and the probability distribution of all possible velocities is not required to describe the macro-
state of the system. Instead, the system is discretized in a way that the particles can only occupy
lattice sites and the velocities are defined such that particles can only travel to their neighboring
lattice sites at each time step. This specific discretization of Equation (3.8) is called the Lattice
Boltzmann Method (LBM). In LBM, the equilibrium distribution function is obtained by Taylor
series expansion of Equation (3.11) with a small velocity (low Mach number) approximation:

(3.13)
Writing Equation (3.12) in a general from and introducing a numerical approximation of
the integral we find

. (3.14)

is a polynomial of and

is the weight associated with each discrete velocity for
estimating the continuous integral. Equation (3.12) can be rewritten as:
44

(3.15 a)

(3.15 b)

(3.15 c)

(3.16)

Depending on the number of discrete velocities, different LBM schemes are defined. LBM
schemes are generally denoted by D
m
Q
n
, where m is the dimensions of the system and n is the
number of discrete velocities. D
2
Q
9
scheme is discussed in this part. The complete derivation of
this scheme is presented in Appendix B. The nine velocity vectors are:

Figure 3.2. The D2Q9 scheme, two dimensional with nine velocities.
45

The equilibrium distribution function for this scheme is as the following:

(3.17)

At this step, Equation (3.8) (Boltzmann Equation) can be written as:

(3.18)

Equation (3.18) and Equation (3.15) are numerically solved for distribution functions and
consequently macroscopic properties.
Conservation equations (continuity and Navier-Stokes) are derived from LBM using
Chapman-Enskog expansions (He and Lou, 1997). During this derivation, one can find:

. is the dimensionless relaxation time and is the kinematic viscosity of the fluid.

is the lattice spacing and

is the lattice time step. Normalized pressure is related to density
via

.
(3.19)

In the porous media, 2D is a different world form 3D. To model flow and transport in
rocks 3D LBM schemes should be utilized. Here we introduce D
3
Q
15
scheme, it can be derived
in the same way as D
2
Q
9
(Appendix B). There are 15 discrete velocities in this model (

which are the columns of the matrix E.

46

The density and macroscopic flow velocity are defined in the same way as Equations 3.15a
and 3.15b. The equilibrium distribution function is defined as:

(3.20)

The pressure is related to the density (

).

3.4 Numerical Implementation and Boundary Conditions
There are two steps involved in LBM. The first step is the streaming where the particles move to
the neighboring sites and the second step is the collision

.
At each time step, macroscopic properties and equilibrium distributions are calculated (Sukop,
2007).
The distribution functions start from an initial value and iterate to reach equilibrium. To
initialize the system, the initial macroscopic properties (density and velocity) are assigned and
equilibrium distribution function is calculated (from Equation (3.17)) and used as the initial
distribution function. Usually, initializations are uniform. However, if we can initialize the
system closer to equilibrium it converges faster.
One of the main advantages of LBM over conventional computational fluid dynamics is its
strength in modeling fluid flow in complex geometries. In the case of solid boundaries the
particles are bounced back with an equal in length and opposite in direction velocity, thus the net
47

velocity on the wall is zero and no slip boundary condition is satisfied. For pressure and velocity
boundary conditions (for example at the inlet), one can find expressions for the unknown
incoming in terms of known outgoing and the specified macroscopic property (flux or
pressure) at that site. Periodic boundary conditions are also very appealing, because they do not
need any programming effort and the streaming step, in fact, satisfies this type of boundary
condition.
Here we provide an example of the von Neumann (flux) boundary condition in 2D.
Implementation of different boundary conditions can be found in the literature (Chen and
Doolen, 1998; Zou and He, 1996).

Figure 3.3. Inlet velocity boundary condition

The input velocity is known (

). Some components of the distribution
function

are known because they are the result of streaming from nodes
inside the domain. The incoming distribution functions (

and , however, are unknowns.
To solve for these unknowns we need four equations. Three of them come from Equations (3.15
48

a, b, c). The fourth equation can be written by assuming that the bounce back condition holds in
the direction normal to the boundary

(Zou and He 1996). The four
equations can then be summarized as:

(3.21 I)

(3.21 II)

(3.21 III)

(3.21 IV)
3.5 Directional Permeability from LBM Simulations
Next, we consider the results of single-phase flow simulation on Fontainebleau sandstone. Figure
3-4 shows a 125
3
sample image. This figure shows the structure of the pore space which is the
images at isosurfaces of pore values (for example, if 1 stands for solid and zero for pores this
shows the zero value isosurfaces).

Figure 3.4. Pore structure in a Fontainebleau sandstone (125^3)
49

The computational domain is 175x125x125, for and the cells are fluid.
Constant pressure boundary conditions have been applied at the inlet and outlet flow, for the
plane , pressure is set to and for the plane , the pressure is set to .
Figure 3.5 shows the simulation set up. All other four sides are no-flow boundaries, defined as
solid points.

Figure 3.5. Simulation setup

50

Figure 3.6. Contour plot of velocity in x direction. Flow is in x direction.

Figure 3.7. Contour plot of velocity in y direction. Flow is in x direction.

51

Figure 3.9. Inlet and outlet pressure, flow is in the x direction.

Note that the variables in the LBM have the Lattice units and the conversion to physical
units is straightforward. Since we are interested in permeability values, only the conversion of
permeability units is discussed. Permeability has the dimension of length squared, and in the
LBM the length dimension is called lu (length unit) which is, in fact, the resolution of the input
Figure 3.8. Contour plot of velocity in z direction. Flow is in x direction.

52

image. . From Darcy’s law , the permeability can
be calculated as the following:

(3.22)

Here, is the sum of the x-direction components of velocity along a constant x-plane,
which due to mass conservation for incompressible systems is the same for all planes
perpendicular to the x axis. is the viscosity of the fluid, is the length, and is the pressure
difference. The typical CPU time for this calculation is roughly 10 hours on a single. Figure 3.10
shows the convergence behavior of the calculated permeability as the simulation progresses.
Note that the negative permeability values are related to the initial state of the system, away from
equilibrium (steady state).

Figure 3.10. Convergence behavior of the permeability

53

To calculate the Y and Z directional permeability one needs to rotate the image and
calculate the flow field for the new input image. These simulations were run and the results are
summarized in Table 3.1.
Table 3.1. Directional Permeabilities

0.186 1.456 1.730 1.469

54

4. Chapter Four: Lattice Boltzmann Simulations of Two Phase Flow in Porous
Media

4.1 Introduction
The true beauty of hydrodynamics appears in multiphase systems. Magical interfaces form as
two immiscible phases are separated through the minimization of interfacial force. The
interfacial energy arises from the molecular force difference between two phases acting on the
interfacial molecules. Modeling of two-phase phenomena has long been studied in the realm of
computational fluid dynamics (Kolev, 2007). Since its introduction, LBM has been widely used
for studying two-phase systems and a range of different multiphase LBM techniques have been
developed (Kazmin, 2010).
In this chapter, we review different multiphase LBM schemes and discuss the Shan and
Chen formulation (Shan and Chen, 1994). We use this model for two phase segregation in two
and three dimensional settings. The application of Laplace law to estimate interfacial tension for
specified two phase settings is discussed. Next, we introduce fluid-solid forces to include
wettability effects in the presence of solid phase. With all components of LBM validated and in
place, we then proceed to calculate the relative permeability curves for 2D examples of porous
rocks. The effect of viscous coupling on relative permeability curve is then discussed.

4.2 Multiphase LBM schemes
The multiphase system can consist of one or more components. For each component a
distribution function is defined. In both cases the implementation of interfacial forces is critical
for two-phase modeling. Gunstensen et al. (1991) introduced the color model where the two
phases are flagged as red and blue. The interfacial force is applied at each time step and the
55

interface locations are calculated. This method is computationally demanding and becomes
unstable for large density ratios. The free energy model was proposed by Swift et al. (1995). The
original idea was to incorporate the surface energy in the collision operator (Succi, 2001). This
model, however, is not stable for systems with high density and viscosity ratios. Shan and Chen
(1994) introduced a pseudo-potential formulation for fluid-fluid forces where the potential of
component is given by:

(4.1)

It means that at each lattice site the force on component is the sum of forces from all
other components on the neighboring sites. The binary interaction potential is defined as:

(4.2)

In this formulation

is a form of the free energy and

is zero
for non-neighbors.

is a matrix determining the strength of interaction between different
components. For example, for a three component system:

(4.3)

Diagonal values are zero, as there is no force applied on a phase by the phase itself. Due to
symmetrical nature of forces (Newton’s 3 rd law of motion), the matrix is symmetric and

. For a two component system the notation can be furthermore simplified

, and the fluid-fluid force for a two-component system are written as:

(4.4)

56

controls the strength of the fluid-fluid forces and below a certain value the two phases become
miscible. For multiphase simulations, we have used Taxila LBM (a Los Alamos National
Laboratory software). Shan and Chen implementation is used in Taxila.

4.3 Surface Minimization Test (2D and 3D)
Here, we provide examples of phase separation and interface minimization. The system starts
with a random distribution of phases (50% saturation). The two phases have equal densities and
is selected to ensure immiscible conditions. Figure 4.1 shows the initial states of the
simulation along with the phase distributions at later times approaching equilibrium (minimum
interfacial length). The simulation starts with a random initialization of two immiscible phases in
a 128x128 domain. Fluid-fluid forces decrease the interfacial surface area (length in 2D) until it
reaches a minimum in equilibrium.

t=0

t=100

57

t=200 t=500

t=1000

t=3000

t=10000

t=20000
Figure 4.1. 2D Surface minimization for two immiscible phases, domain size is 128x128 and
periodic boundary conditions are applied in all directions. Blue is water and red is oil.

The same problem is solved for a 3D system (Figure 4.2). Initial and boundary conditions
are similar to those used for the 2D example. The physical equivalent of this numerical
experiment is when you shake two immiscible phases (oil and water for example) and let them
rest. Isosurfaces are shown which can be considered as interface between the two phases. The
system starts form a complete emulsion with large interfacial areas and the fluid-fluid forces act
to minimize these surfaces toward equilibrium.
58

t=0

t=200

t=400

t=800

t=1000

t=2000

t=4000

t=20000
Figure 4.2. 3D Surface minimization for two immiscible phases, domain size is 100x100x100
and periodic boundary conditions are applied in all directions
59

4.4 Laplace Law Validation
Laplace law states that the pressure difference between the interfaces of a circular droplet (with
radius r) of one phase surrounded by another phase is inversely proportional to the radius of and
the constant of proportionality is called the interfacial tension (σ).

(4.5)

As discussed in the Shan and Chen formulation, controls the fluid-fluid force strength.
However, there is no explicit equation to relate and . For a constant , simulations for drops
with different radii are carried out and the pressure difference at the equilibrium is calculated.
The slope of the linear plot ( vs

) is . Figure 4.3 shows the experimental setup and phase
distributions for the test. The results of simulation on 8 different radii are shown in Figure 4.4
and the interfacial tension is equal to the slope of the curve.

Figure 4.3. Simulation setup for droplets of different radii. Periodic boundaries are applied in all
directions.

60

Figure 4.4. Validation of the Multi Component Multi Phase (MCMP) LBM code via Laplace
law: A linear relationship between the pressure difference and inverse radius. For the two fluids
in this system, the IFT is 0.0051 mu/ts^2. mu (mass unit) and ts (time step) are LBM units for
mass and time, respectively.

4.5 Implementation of Adhesion Forces
Reservoir rocks have been exposed to oil and/or water for millions of years and in most cases
have different tendency towards different fluids. This effect is caused by surface chemistry of the
rock and is known as wettability. In this section, we include solid-fluid interaction forces that
account for the wettability. Wettability is an input to the simulation of two phase flow in rocks
and should be determined from other techniques (e.g. imbibition tests or contact angle
measurements). These tests, however, provide the average wettability of the core and no
measurement of the actual pore surface wettability is reported in the literature. Hence, it remains
as the main uncertainty in the context of pore scale modeling.
To model this effect, we introduce a force that acts at the solid-fluid interfaces. The surface
force has the following form:

(4.6)

y = 0.0051x + 0.0003
0
0.0001
0.0002
0.0003
0.0004
0.0005
0.0006
0.0007
0.0008
0.0009
0.001
0 0.02 0.04 0.06 0.08 0.1 0.12
ΔP
1/r
61

where the subscript denotes the phases,

is the adsorption strength and is a multiplier
that is set to 1 if a point in the calculation domain represents a solid and is set to zero otherwise.

4.6 Wetting in 2D and 3D
The final form of a static droplet, surrounded by another fluid, on a surface is the result of
equilibrium between fluid-fluid and fluid-solid forces. In this part we show the results of
wettability modeling for a wide range of wetting properties in 2D and 3D.

Figure 4.5. Initial conditions for wettability modeling in 2D: A 30x30 droplet of oil in a 120x80
domain. The lower boundary condition defines the solid wall and all other boundary conditions
are periodic.

Wettability is defined by contact angles and can be calculated by the force balance or by
geometrical considerations for each case in Figure 4.6. The initial condition configuration for all
cases is given in Figure 4.5.

(4.7)

62

63

Figure 4.6. Different wetting properties: G
o
and G
w
are water and oil adsorption coefficients. A
complete range from non-wetting to completely water wet cases is shown. G
ow
is set to 0.2.

The wetting phenomena in 3D are shown in Figure 4.7. The initial condition is a spherical
droplet on the solid plate and two cases of oil wet and water wet systems are simulated. When
the rock is water-wet, droplet spreads on the surface and for oil-wet surface the droplet remains
the same. The interfacial forces smooth the initial droplet.

64

Initial Condition

, Oil Wet Surface

, Oil Wet Surface
Figure 4.7. Wetting in 3D, starting from water droplet. This droplet remains spherical for the oil-
wet system and spreads on the surface for water-wet system (bottom row). G
o
and G
w
are water
and oil adsorption coefficients. G
ow
(oil-water interface constant) is set to 0.2
65

4.7 Relative Permeability and its Functional Dependencies
Due to the viscosity effects, fluids lose energy as they flow: the higher the viscosity is the higher
the pressure drop will be for a constant flow rate. In porous media, due to the presence of
tortuous paths, this pressure drop can be significantly higher than pipe flow. This force is called
the viscous force and is the only force in single phase flow settings. In two-phase situations,
interfacial forces (fluid-fluid) become important as well. Rocks have surface areas as large as
few meter squares per each centimeter-cube volume. The surface, depending on its chemistry,
may prefer one phase to another. This force is called wetting, as described in the previous
section, and appears at the interface of rock surface and two fluids. The dynamic behavior of two
phase flow in rock results from the combination of these three forces.
Figure 4.8 shows an example of two phase flow in a 2D porous media. The porous media
has a high porosity (0.77) and is used for validation of concepts. The domain size is 400x400 and
solid particles (black) are 20x20. Le+ft and right boundaries are walls and upper and lower
boundaries are periodic. Wetting phase saturation is 0.4 and is distributed randomly at the initial
state. The surface tension is 0.0051 and a gravity force of 0.00001 is applied (all in LBM units).
The surface is completely water wet. Phase distribution at different time steps is shown. Blue
phase (wetting) surrounds the solids. There configuration of two phases is the same at t=50000
and 75000 which indicates the equilibrium state of the system.

66

t=0

t=500

t=1500

t=3000

t=10000

t=30000

t=50000

t=75000
Figure 4.8. Distribution of two fluids at Sw=0.4. Water (blue) is the wetting phase. Black is the
solid and red is the non-wetting phase. The saturations start from an initial random distribution.
67

Figure 4.9. Vertical (flow direction) component of velocity at t=75000.

The vertical component of velocity is shown in Figure 4.9. Note that most of flow happens
through the connected oil saturation, at equilibrium, that runs from inlet to the outlet.
A series of two-phase simulations at different saturations is used for calculating relative
permeability curves. Figure 4.10 shows the final configuration of the two phases at different
saturations. Note that as we add the wetting phase to the system it sticks to the surface. The flow
is in horizontal direction and upper and lower boundaries are walls. Saturations are shown after
150,000 time steps. The simulation time is ~3 hours for each saturation using 64 processors.

68

Sw=0.0

Sw=0.1

Sw=0.2

Sw=0.3

Sw=0.4

Sw=0.5

Sw=0.6

Sw=0.7
69

Sw=0.8

Sw=0.9
Figure 4.10. Distribution of wetting (blue) and non-wetting (orange) phases at different wetting
phase saturations. Densities and viscosities are the same and flow is in the horizontal direction.

Figure 4.11 shows the relative permeability curves resulted from two-phase simulations
shown in Figure 4.10. Note that there are saturations that none of the phases are flowing. The
reason is that the slug of one phase formed in the middle of the domain closes the flow path to
the other phase.

Figure 4.11. Relative permeability curves from the simulated flow field at different saturations

70

4.8 Viscous coupling effects (analytic and LBM)
In the extension of Darcy’s law to two phase flow (

), it is
assumed that the two phases are hydrodynamically decoupled. To investigate the importance of
viscous coupling in two phase flow and to study the accuracy of LBM in capturing this effect,
analytical solution for layered flow in a 2D channel (Figure 4.12) is presented in this part.

Figure 4.12. Schematic for layered flow of two immiscible phases in a channel

Two dimensional incompressible Steady State Navier-Stokes Equation:

(4.8)

In this case

and Equation (4.8) reduces to the following form.

(4.9)

The boundary conditions are summarized below.
71

(4.10)

Solving for velocities in both phases:

(4.11)

Where:

,

(4.12)

is the gravity constant and

is the viscosity ratio. As it can be seen from the
analytical solution, the flow of two phases is coupled. For example if, if there is no gravity force
on the wetting phase (

), the wetting phase still moves due to the flow of non-wetting
phase.
To investigate the accuracy of LBM in modeling viscous coupling phenomena, velocity
profiles for different cases are simulated and the numerical results are compared to the analytic
solution (Figures 4.13, 4.14 and 4.15). Table 4.1 summarizes nine different cases with variable
viscosity ratios and flow directions. In all cases a=25, b=50 and the two phases have equal
densities.

72

Table 4.1. Different simulation cases in 2D channel
Case 1 2 3 4 5 6 7 8 9
Figure 4-13 a 4-13 b 4-13 c 4-14 a 4-14 b 4-14 c 4-15 a 4-15 b 4-15 c
1 3 0.33 1 3 0.33 1 3 0.33

g G g g g g -g -g -g

g G g -g -g -g G g g

Considering that Shan and Chen multi-phase scheme is not completely immiscible, the
interface thickness is not zero in LBM simulation while in the analytical solution no interface
layer is considered. This leads to a slight mismatch between the numerical and analytical
solutions at the interfaces.

73

a) M=1

b) M=3

c) M=0.33
Figure 4.13. Comparison of analytical solution and LBM (a=25, b=50, Gw=Gn=g)

74

a) M=1

b) M=3

c) M=0.33
Figure 4.14. Comparison of analytical solution and LBM (a=25, b=50, Gw=-g, Gn=g)

75

a) M=1

b) M=3

c) M=0.33
Figure 4.15. Comparison of analytical solution and LBM (a=25, b=50, Gw=g, Gn=-g)

76

In the extension of single phase Darcy’s law to two -phase systems, it is assumed that the
two phases act independently and viscous coupling is neglected in many settings. For the 2D
channel problem, however, the two phases are coupled and there are non-zero off diagonal terms
in the two phase formulation of Darcy’s equation.

(4.13)

In this formulation

is the relative permeability of wetting phase due to the flow of
non-wetting phase. Here we derive the explicit equations for the four rel-perms in the 2D
channel. Obviously,

. In two dimensions the rate is and the
relative-permeabilities are calculated (Equation (4.14)).

(4.14)

As it is seen the off-diagonal terms (

) are non-zero and flow of two phases
is coupled. To investigate this effect in porous media we have calculated the relative
permeability curves for the 2D porous media in Figure 4.10. Cases 1-6 in Table 4.1are simulated.
The effect of counter-currency is shown in Figure 4.16 where negative driving force on the
wetting phase shifts the rel-perm curve to smaller values.

77

Figure 4.16. Relative permeability curves for co and counter-current flow (Gw=g, Gn=-g).
Viscosities and densities are the same for both phases.

Figures 4.17 and 4.18 demonstrate the effect of viscosity ratios on rel-perm curves in co-
and counter-current settings.

Figure 4.17. Relative permeability curves for co-current flow at different viscosity ratios
78

Figure 4.18. Relative permeability curves for counter-current flow (wetting phase flows in the
negative direction) at different viscosity ratios

79

5. Chapter Five: Theoretical Analysis of Capillary Rise in Porous Media

5.1 Introduction
In this chapter, we study the dynamics of capillary-driven fluid invasion in three different
settings including: a) a single capillary tube, b) a porous column and c) a fractured porous media.
A Lambert W functional form is proposed to describe the invasion dynamics in a single capillary
tube, that predicts both early time Washburn-type behavior ( ) and late time exponential decay
behavior. We extend the formulation to describe homogenous porous media and to include
viscosity, pressure and gravity effects in both advancing and defending fluids. Solutions for
closed systems, where the advancing fluid compresses the defending fluid, are then developed.
Finally, we extend the theory to describe fractured systems, and propose a convolution integral
formulation and a new explicit solution.
Capillary assisted liquid flow is the main transport mechanism in many natural systems. It
dictates the maximum height of trees on earth (Koch et al., 2004) and contributes to the
regulation of blood flow in our brain (MacVicar and Salter, 2006). Oil extraction from reservoir
rocks and immiscible contaminant dissipation in subsurface water resources are examples where
capillary dominated flows play a central role (Bear, 2013). Furthermore, a clear understanding of
capillarity in micro and nano systems is pivotal to the advance in design of micro/nano electrical
mechanical systems (MEMS and NEMS) (Pelesko and Bernstein, 2002; Gad-el-Hak, 2010; Ho et
al. 1998).
Capillarity was first noted by the immaculate eye of Leonardo da Vinci (McCurdy, 1906).
In the early 18th century, Jurin (De Gennes, 2004) proposed an equation for the final rise, he, of
a fluid with density, ρ, in a capillary tube with radius, R, in the form
80

(5.1)

where σ is the interfacial tension, θ is the contact angle, and g is the gravitational acceleration.
This equation predicts an equilibrium rise of water in the range of 15 km for a capillary tube with
a radius of 1nm (10^-9m). Caupin et al. (2008) investigated the maximum possible capillary rise
and report equilibrium heights, for water, exceeding 100 km in cylindrical graphitic pores.
Washburn (1921) first studied the dynamics of capillary driven flow in the 1920’s and derived a
time dependence of the imbibing fluid height in the form . Dimitrov et al. (2007) later
verified Washburn’s equation in nanopores using molecular dynamics simulation. Ponomarenko
et al. (2011) demonstrated that capillary rise in corners of different geometries obeys the
relationship h

, while Weislogel (2012) reported equations to describe capillary rise in
conduits with multiple corners in low-gravity environments. Weislogel’s closed form solutions
follow Washburn’s h dependence. Capillary rise in non-uniform tubes was investigated via
numerical calculations by Erickson et al. (2002) and analytically by Liou et al. (2009).
Imbibition in porous media is often denoted as the cumulative fluid uptake of a porous
sample and has been studied extensively both experimentally and numerically (1997). Delker et
al. (1996) report results from experiments on capillary rise in a column of glass beads with
different diameters. Handy (1960) derived an equation, similar to Washburn’s equation, for
imbibition in porous rocks to estimate the effective capillary pressure in sandstones. At larger
scale, during water injection into fractured reservoirs, imbibition of water into the water-wet
matrix block is the main oil production mechanism (1978). During hydraulic fracturing in shale
systems, large volumes of water are injected to break the rock with the purpose to increase the
surface area of the contact area between reservoir and the wellbore. The injected water will then
imbibe into the shale rock and potentially block the flow paths of the oil/gas during the
81

subsequent production stage. Roychaudhuri et al. (2013) addressed this challenge experimentally
and demonstrated that imbibition processes can account for the observed fluid loss (~70%)
during stimulation of shale gas wells. At the pore scale, Hyväluoma et al. (2006) applied the
lattice Boltzmann method to model liquid penetration in paperboard and concluded that
multiphase lattice Boltzmann methods accurately captures the essential physics of capillary
penetration.
We present a generalized form of Handy’s equation (1960) for capillary rise in porous
media that includes viscosity, pressure and gravity for both advancing and defending fluids. For
a capped porous cylinder, the coupling between defending fluid’s pressure and rise height is
introduced and the rise dynamics is compared for open and closed systems. Finally, we discuss
imbibition in fractured systems and propose a new explicit solution.

5.2 Single Capillary Tube
The driving force for liquid rise in a capillary tube, as shown in Figure 5.1, is a combination of
cohesive (fluid-fluid interactions) and adhesive (fluid-solid interactions) forces while the
opposing forces include gravity and viscous forces.
82

Figure 5.1. Schematic of capillary rise in a single tube

For the single capillary tube, we assume initially that the defending fluid is invicid with
negligible mass (e.g. air displaced by water) and write the momentum balance in the tube as
1

(5.2)
where

is the shear stress (=μ

for a Newtonian fluid). The entry length (Le) is the length that
the flow velocity fully develops in the tube and is approximated as (1978)

(5.3)
where

is the Reynolds number. For the majority of the flow L
e
is small compared to
the capillary height h and we assume that fully developed Poiseuille flow prevails along the
height (h). By this assumption, the radial variation in the fluid velocity is given by

1
Note that the contact angle θ is a function of velocity

and static contact angle (Barcke et al., 1989). Here, we
assume that θ is equal to the static contact angle. However, one can replace θ with an appropriate function.
83

(5.4)
where

is the average velocity of the liquid column. By substituting

and

into Equation (5.2), we find

(5.5)
Equation (5.5) is an ordinary nonlinear second-order differential equation that can be
solved numerically. Initially, h is small and the capillary force controls the dynamics. During this
early stage, the fluid height is proportional to t
2
and then t, depending on the fluid viscosity,
while at later times when

, is the kinematic viscosity, the dynamics become viscous
dominated and Equation (5.5) reduces to

(5.6)
The differential Equation (5.6) has the analytical solution

(5.7)

where

and

with the equilibrium height h
e
given by Equation (5.1). For

, one finds that

, which is the well-known Washburn equation. As the capillary
height approaches equilibrium, (

), Equation (5.7) can be written as

which is an exponential decay. To express

as a function of

, Equation (5.7) is inverted to

(5.8)

where is the Lambert W function (2000).
84

As suggested by Zhmud et al. (2000) the L.H.S. term in Equation (5.2) can result in an
oscillatory behavior at the late time behavior. Masoodi et al. (2013) conclude that for the cases
with

the late time oscillations disappear and solutions to Equation (5.5) and
Equation (5.8) become identical. As an example, for a capillary tube with R = 10
-4
m Equation
(5.5) and Equation (5.8) are compared for a water-air fluid system with a contact angle of 45
o
(Figure 5.2). We observe that the early inertial effects disappear quickly and that Equation (5.8)
is accurate for viscous systems. There is no late time oscillation in this case ( = 1.47e-4). This
finding is used later to develop new expressions for capillary rise in porous media.

Figure 5.2. Capillary rise in a tube with radius 10
-4
m including inertial effects (Equation (5.5) -
dots) and without inertial effects (Equation (5.8) - solid line).

For smaller radii, the equilibrium height increases while the rise rate reduces. The
dynamics of capillary rise is demonstrated for two tubes with R=1nm and R=1µm in Figure 5.3.
We observe that it will take 125 years for water to climb 10m in a nanotube, while for the
85

microtube a similar advance takes only 0.6 years. As observed, the imbibition process can be
extremely slow and it can take millions of years to travel distances of relevance to geological
settings.
In most natural systems, capillarity is driven by porous structures that can be modeled as
3D networks of capillary tubes with different radii (Dong and Blunt, 2009). For homogenous
systems (e.g. a packed bed of spheres), the pore size variation is relatively small. However, for
heterogenous systems(such as subsurface rocks), pore sizes can vary by orders of magnitude.
This heterogenity, in turn, will have a significant impact on the imbibition characteristics of the
system. We consider both settings in the following sections.

Figure 5.3. Dynamics of capillary rise in nano- and micron-tubes.

5.3 Homogenous Porous Media
To develop a description of capillary rise in porous media, we start by assuming that the porous
media is homogenous. In other words, the pores are similar in diameter and dispersion is small
86

(narrow pore size distribution). The purpose is to find an equation for the frontal position of the
imbibing fluid (advancing fluid). To find a general solution, we do not make any assumptions
regarding the viscosity of the advancing and defending fluids. In addition, the advancing and
defending fluids can be pressurized to represent settings of both spontaneous and assisted
capillary rise dynamics. Figure 5.4 shows a schematic of the system under study. The bottom and
top faces of the cylindrical porous medium are initially open while all other sides are closed. We
note that the formulation is not limited to circular cross sections.
The driving forces for fluid displacement are capillary forces and the pressure of the
advancing fluid. Capillarity appears where a solid is in contact with two fluids of different
affinity to the solid surface(s) and is the only driving force during spontaneous imbibition.

87

Figure 5.4. Schematic of a porous cylinder. L is the length of the system while Z denotes the
frontal position of the advancing fluid. The pore spaces are initially saturated 100% with a
defending fluid.

For forced imbibition processes, the advancing fluid pressure P
A
must also be considered.
The total capillary force at the displacement front is

(5.9)

where σ is the interfacial tension between advancing and defending fluids and

is the pore
radius. n is the average number of pores at each cross section. Next, we approximate the
summation in Equation (5.9) by

where
88

(5.10)

is the expected value of r and f is the distribution function for pore radii in the porous system.
Based on this approximation, we write the capillary force as

(5.11)

To find an equation for n, one can argue that the total porous area in each cross section is
equal to porosity φ multiplied by cross-sectional area A to find

(5.12)

The driving force due to the pressure of the rising fluid is P
A
φA. The forces of resistance
acting against capillary rise are the weights of the advancing ρ
A
gZφA and the defending ρ
D
g(L-
Z)φA fluids, the pressure of the defending fluid P
D
φA and the viscous forces of both advancing
and defending fluids. Here, we calculate the viscous pressure drop from Darcy’s equation

(5.13)

where

is the pressure drop due to viscous effects to arrive at

(5.14)

In Equation (5.13) and Equation (5.14), V
A
is the Darcy velocity of the advancing fluid, k is
the permeability, µ
A
is the viscosity and Z is the frontal location of the advancing fluid. S
A
is the
saturation of the advancing fluid behind the front (1-S
R
) and takes into account any residual
saturation (S
R
) of the defending fluid. Similarly, we write the viscous resistance of the defending
fluid as

(5.15)

89

where L is the total height of the porous column and µ
D
is the viscosity of the defending fluid.
Equivalent to the single tube problem discussed above, the inertial terms can be neglected and
the force balance for the rising fluid takes form

(5.16)

At the equilibrium t =

and the final height is found directly from Equation (5.16)

(5.17)

By separation of variables in Equation (5.16) we find

(5.18)

When µ
A
S
A
- µ
D
≠ 0, Equation (5.18) simplifies to

(5.19)

where

is the dimensionless advancing front and the dimensionless time t* and the
constant c are given by

(5.20)

The solution to this ODE with the initial condition, Z
*
(0) = 0, is

(5.21)

The inverse form of this function gives an explicit solution to the front location as a
function of time

(5.22)

Equation (5.22) takes the same functional form as the single capillary problem given in Equation
(5.8), and the volume of the imbibed fluid is found from
90

(5.23)

For the case where µ
A
S
A
- µ
D
= 0, Equation (5.18) reduces to

(5.24)

and the advancing front obeys

(5.25)

Figure 5.5 illustrates the effect of the defending fluid viscosity on the rise dynamics. For an
invicid defending fluid, c in Equation (5.22) is zero and the solution becomes independent of the
system length (L) For the case where µ
A
S
A
= µ
D
, Equation (5.25) is applied and the behavior is
exponential from the beginning (t=0).

Figure 5.5. The effect of viscosity ratio on rise dynamics for a vertical homogeneous porous
medium (L = 15m) with a permeability of 1 Darcy (9.869*10
-13
m
2
) and a porosity of 30%

91

Equation (5.19) provides the general solution for early and late times where both fluids can
be pressurized and viscid. A special case is that of free imbibition in the absence of a defending
fluid (ρ
D
= μ
D
= P
D
= P
A
= 0). For this case, Equation (5.21) is simplified to

(5.26)

From a Taylor series expansion of Equation (5.26) at early time

we find

(5.27)

By substituting for t* and Z* we find

(5.28)

Here

is the capillary pressure. The imbibed volume is derived to be:

(5.29)

which is identical to the expression derived by Handy (1960).
If the homogeneous system (see Figure 5.4) is closed at the top, the defending fluid is
trapped and compressed by the advancing fluid. Accordingly, for an incompressible defending
fluid, the advancing fluid cannot invade the closed system. However, in such cases, slight
heterogeneities will promote a counter-current flow where the advancing fluid can expel the
defending fluid. For a compressible defending fluid, P
D
will increase with time due to the
imbibition of the advancing fluid. The compression of the defending fluid has a dampening
effect on imbibition dynamics. Here we assume that the defending fluid is compressible and
cannot escape through the imbibed zone. This assumption can be justified for a fractal structure
associated with induced fracture networks. The change in P
D
is

(5.30)

92

where is the compressibility factor of the defending fluid, that is a function of pressure and
temperature through an equation of state.

is the compressibility factor at the initial P
D
. When
P
D
is variable, Equation (5.22) and Equation (5.30) are coupled and must be solved iteratively.
To evaluate the equilibrium height in this closed system we insert

into Equation
(5.16).

+

(5.31)

The second-degree polynomial given in Equation (5.31) has two roots and the smallest root
is the equilibrium height for a closed system. The imbibition dynamics of a porous column with
a length of 15m, a permeability of 1 Darcy and a porosity of 30% is compared for open and
closed settings in Figure 5.6. The average pore size is 1µm and the pore space is initially filled
with air at atmospheric pressure. The equilibrium rise for the open system is Z
eq
= 4.35m, while
for the closed system Z
eq
= 2.40m. The pressure in the closed system increases from the initial
1atm to 1.2atm at equilibrium.
93

Figure 5.6. The dampening effect of pressure in a water-air closed system. The pressure
increases with the rise in the closed system and as a result, the fluid imbibition stops earlier.

5.4 Heterogeneous Porous Media
Many natural systems are not homogenous as sketched in Figure 5.4. Fractured rocks are an
example of heterogeneous porous media. Imagine a core sample, centimeters in length and
diameter, which has a network of micro-fractures embedded in a matrix with pores at the
nanometer range. This setting is representative of unconventional oil and gas resources in the
context of shale reservoirs. The overall system (matrix and micro-fracture network), with a total
porosity φ
t
, is commonly referred to as a dual porosity system (Warren and Root, 1963). The
fractures represents one porous system with a porosity φ
f
, a permeability k
f
and an average
weighted pore size of E[r
f
], while the matrix represents a second porous system with a porosity
φ
m
, a permeability k
m
and an average weighted pore size of E[r
m
]. We note that the fracture
porosity φ
f
is defined relative to the fracture volume. The two porous systems are in
94

communication across the matrix-fracture interface areas. In fractured systems, the advancing
fluid rises, initially, through the fractures and imbibition into the matrix starts when the wetting
fluid spreads on the matrix surfaces. Accordingly, the matrix-fracture interface area A
mf
controls
the volume imbibed into the matrix. We assume here a simplified representation of a fractured
rock “sugar -cube model” as shown in Fig ure 5.7 to estimate the matrix-fracture area:

(5.32)

where N is the number of cubes with length l and Z
f
(t) is the frontal location in fractures.

Figure 5.7. Sugar-Cube representation of a fractured porous medium. The fracture system has a
porosity φ
f
, a permeability k
f
and an average weighted pore size of E[r
f
], and the matrix has a
porosity φ
m
, a permeability k
m
and an average weighted pore size of E[r
m
].

The fracture opening is w. Here, we estimate the imbibition in a 3D matrix cube with an
equivalent 1D system that has a base area of 6l
2
and a height of l/6. This assumption increases
the rate of imbibition, as in a 1D system there are no resistances imposed by fluids moving from
95

other directions. For a sphere with volume l
3
, we find that an equivalent 1D system will have a
base area of 4.86l
2
and a height of l/4.86. In the following, we focus on the simpler formulation.
Next, we introduce a new variable, R
m

(5.33)

which defines the volume of matrix blocks relative to the total system volume. A and L are
sample area and length of the overall system. R
m
is also related to fracture porosity, matrix
porosity and total system porosity

(5.34)

and the average fracture opening (w) and R
m
are related to the matrix cube length (l)

(5.35)

To simplify the equations, it is assumed that L >> l to arrive at an expression for the matrix-
fracture interface area

(5.36)

We note that this formulation applies to both porous fractures and open (empty) fractures
via φ
f
=1 and w=2E(r
f
). During an imbibition process, the total uptake of the wetting phase (Q) is
the result of imbibition into both matrix and fractures

. (5.37)

We find the fluid invasion into the fracture system from

(5.38)

where S
A,f
is the saturation of the advancing fluid in the fractures (accounts for residual
nonwetting phase in the fractures) and Z
f
(t) is calculated from Equation (5.22). Evaluation of Q
m
96

is more complex due to the time lag in the arrival of the advancing fluid to the matrix surfaces,
and Q
m
is hence expressed in the form of the convolution integral

(5.39)

Here, T=min(t,t
F
), where t
F
is the time that it takes for the advancing fluid to rise to the end of
the fracture network (Z
f
= L). By inserting Equation (5.36) into Equation (5.39), we establish

(5.40)

where Z
m
is the front location of the advancing fluid in matrixcube, that is calculated from
Equation (5.22) based on φ
m
, k
m
and E[r
m
]. The imbibition into a matrix cube stops when it is
saturated with the advancing fluid. This dictates a maximum value of Z
m
equal to l/6. Equation
(5.40) can be evaluated numerically using the full solution for Z
m
from Equation (5.22).
However, for small systems (

), one can use Equation (5.28) for Z
m
and

Z
f
to find

(5.41)

In this formulation, C
f
, C
m
and t
F
are given by

(5.42)

(5.43)

.
(5.44)

Finally, we define a characteristic matrix time as
97

(5.45)

and a closed form solution for the matrix imbibition is obtained as

(5.46)

(5.47)

We note that this solution, similar to the integrand in Equation (5.41), has an asymptote at
zero. We can now calculate the total uptake of the wetting phase from the sum of Equation (5.38)
and Equation (5.47)

.
(5.48)

The imbibition characteristics of a fractured porous system, as predicted by Equation
(5.48), are shown in Figure 5.8 for w=0.5 mm, φ
t
=0.3, φ
m
=0.25, φ
f
=0.5, K
m
=10
-8

D, K
f
=1

D,
E[r
f
]=10 μm and E[r
m
]=10 nm. R
m,
evaluated from Equation (5.34), is 0.2 and l
,
as evaluated
from Equation (5.35), is 6 mm. For an overall system length (L) of 10 cm, the total pore volume
is 3 cc. The advancing fluid is water and the system is initially at vacuum conditions. The
observed bi-linear imbibition dynamics (Figure 5.8) are consistent with the spontaneous
imbibition experiments on gas shales reported by Roychaudhuri et al. (2013).

98

Figure 5.8.

Fracture, matrix and total volumetric uptake for a fractured system with φ
t
=0.3,
φ
m
=0.25, φ
f
=0.5, K
f
=1

D, K
m
=1e-8

D, E[r
f
]=10 μm and E[r
m
]=10 nm. Total pore volume is 3 cc
and the fracture opening (w) is 0.5 mm.

5.5 Discussions and Conclusions
Previous sections of this chapter present a theoretical analysis of capillary flow in a single
capillary tube, a homogeneous porous medium and a heterogeneous (two-scale) porous material.
For a single tube, we show that the initial inertial regime disappears quickly for a water-air
system. As a result, the momentum conservation equation is simplified and a Lambert W
analytical equation applies, that describes early and late time behavior accurately. This solution
predicts that a capillary rise of 10m in a nanometer tube takes approximately 125 years and
clearly demonstrates the potential of very slow moving fluid fronts in unconventional reservoirs.
We present a general equation for homogenous porous systems that includes pressure,
viscosity and density of both advancing and defending fluids. Input parameters for this model
99

include porosity, permeability and the weighted average of the pore radii. We show that by
increasing the viscosity ratio (µ
D
/µ
A
), the early time behavior deviates from the well-known
dynamics, while for a unit mobility ratio (µ
D
/µ
A
=1), the uptake dynamics are represented by an
exponential decay. This observation helps, for example, to better understand the imbibition
dynamics of fracking fluid during stimulation of shale oil wells, where both advancing and
defending fluids are pressurized and viscous.
An interesting component of this work is the imbibition dynamics in a closed system. As
an example for a porous material, we show that the equilibrium rise in closed settings is
approximately the half of that in an equivalent open system. We observe the same effect for a
single capillary tube. Considering the fact that in any porous material, dead-end pores represent a
fraction of the pore space; this observation can further help to estimate the residual saturation of
the defending fluid during imbibition processes.
In a closed pore, as the advancing fluid rises, the pressure of the trapped defending fluid
increases. If the pressure increases sufficiently to balance the capillary force, the imbibition
process stops and equilibrium is established. However, if pores are small enough, the capillary
forces can drive the advancing fluid to a level where the defending fluid pressure exceeds the
mechanical strength of the material. In such case, the pores will act as small natural pumps to
create and potentially grow micro/nano fracture networks.
For a fractured system, liquid rises faster in the fractures and imbibition in the matrix
segments start as the wetting fluid spreads on the matrix surfaces. We use a simple sugar-cube
model to represent a fracture system and find an explicit formulation to predict the liquid uptake
as a function of time. The input parameters for this model are total porosity of the system,
average fracture opening in addition to porosity, permeability and weighted average of the pore
100

radii for both matrix and fracture systems. While this work considers porous systems with two
discrete levels of porosity, an extension to n-scale porous systems is straightforward.
For a closed system, coupling between the increasing defending pressure and rise dynamics
is discussed for a homogenous system. The same discussion is relevant for the heterogeneous
porous media. As the fluid rises in the fracture system, gas is trapped in the matrix cubes and the
capillary pressure should overcome the (increasing) gas pressure for the matrix to be saturated.
For an atmospheric gas pressure, this effect becomes critical to capillary dynamics in micron
range porous systems. However, it can be neglected for nano scale (shale gas) porous rocks.
Further theoretical and experimental work is required to better quantify the closed system effect.
Imbibition experiments are widely used as a characterization technique in petrophysical
evaluations. The work presented here provides a theoretical background to promote
interpretation of experimental observations. One approach will be to match the model to
imbibition data to find relevant characteristic rock properties. To promote the characterization
further, multiple imbibition tests with a different viscosity of the defending fluid (e.g. water-air
and water-oil) will allow for a better characterization of porous materials. Finally, the theory and
models discussed in this work are directly applicable to the characterization of fractured systems.
Carefully designed imbibition experiments will allow us to differentiate between matrix and
fracture systems and to delineate their characteristic properties, individually, from the imbibition
data.

101

6. Chapter Six: Potential Future Work and Research Directions

In this work, we have studied shale oil samples using image based characterization as well as
pore and core scale modeling efforts. We have applied Lattice Boltzmann Method for pore scale
modeling of single and two-phase flow phenomena. On the core scale, we have developed new
equations for imbibition dynamics in homogeneous and heterogeneous porous rocks.
We have addressed different problems in each of these areas (characterization, pore and
core scale physics). However, there are components that require further studies to improve
fundamental knowledge of transport phenomena in tight porous systems and enhance outputs of
image-based characterizations.
On the characterization part, we concluded that FIB-SEM images capture the main pore
structure for shale oil samples and can be used as inputs for the digital rock physics workflow.
Currently, the imaging resolution is the same in all directions which requires a large number of
SEM images. Any attempts, e.g. geostatistical interpolation, to reduce the number of slices
required for rendering the 3D image can highly decrease the imaging cost and make the
workflow more appealing for industrial applications.
In pore scale modeling of two-phase flow, wettability information is required at each pore
surface voxel. However, currently wettability is measured at core scale and pore surface wetting
property is the main uncertainty in the context of pore scale modeling. Developing methods to
understand wettability variations at the pore scale is a big step in characterization of materials
with complex surface chemistry. A possible approach is to design nano particles with a
preference to deposit on oil-wet (or water-wet) surfaces. A gas containing such nano particles is
102

injected into the porous media. The porous media can later be imaged to map the sites with
deposited particles and hence characterize wettability.
We have calculated relative permeability curves in two-phase settings using LBM
simulations. However, extension to three, immiscible, phases is straightforward in the current
LBM framework. Three-phase relative permeability curves, critical gas saturation etc can be
studied in this context.
In shale gas systems, permeability is often in the nano-Darcy range. Flow through
experiments, even with gas, become increasingly tedious. Furthermore, given that a large
percentage of pores are in the nano-meter range and the fact that at each pore size range,
different physics governs the gas flow complicates the interpretation of experiments. In large
pores (order of tens of microns) Knudsen number is very small, , and the Navier-Stokes
equation with no slip boundary condition is valid. In medium pores (from ~100nm to few
microns) Knudsen number is ~0.1 and the Navier-Stokes equation with slip boundary condition
can be used to model gas flow. In small pores (order of a few nanometers) the Knudsen number
is large, , and the gas flow cannot by represented accurately by the continuum
hydrodynamics (Navier-Stokes). In this case, the mean free path of the gas molecules is
comparable to the pore diameter and the number of collisions between molecules and pore walls
dominates the number of intermolecular collisions. In the extreme case of very small pores at
low pressure, the intermolecular collisions can be neglected altogether and Knudsen diffusion
dominates mass transfer.
The original Boltzmann equation is derived with the assumption of no wall effect (large
container). However, one can formulate a probabilistic approach for the systems with high wall
effect. Adsorption and desorption effects can be included in such a model.
103

In the last part of chapter 5, we derive an equation to describe imbibition in fractured
systems. We use a simple sugar-cube model to represent the fracture system and find an explicit
formulation to predict the liquid uptake as a function of time. While this work considers porous
systems with two discrete levels of porosity, an extension to n-scale porous systems is
straightforward.
In most parts of this thesis, we have discussed methods to estimate transport properties of
submicron rocks from their images. The inverse of this problem, i.e. designing porous
geometries with desired transport properties, is interesting and can have a wide range of
applications including micro-fluidics, biomaterials etc.

104

7. References
Adesida, A., Akkutlu, I. Y., Resasco, D., &Rai, C. (2011, October). Characterization of Barnett
Shale Kerogen Pore Size Distribution using DFT Analysis and Grand Canonical Monte
Carlo Simulations. In SPE Annual Technical Conference and Exhibition.
Adler, P. M., Jacquin, C. G., & Quiblier, J. A. (1990). Flow in simulated porous
media. International Journal of Multiphase Flow, 16(4), 691-712.
Balasubramanian, K., Hayot, F., & Saam, W. F. (1987). Darcy’s law from lattice -gas
hydrodynamics. Physical Review A, 36(5), 2248.
Barrett, E. P., Joyner, L. G., &Halenda, P. P. (1951). The determination of pore volume and area
distributions in porous substances. I. Computations from nitrogen isotherms. Journal of the
American Chemical Society, 73(1), 373-380.
Barrett, E. P., Joyner, L. G., &Halenda, P. P. (1951). The determination of pore volume and area
distributions in porous substances. I. Computations from nitrogen isotherms. Journal of the
American Chemical society, 73(1), 373-380.
Bear, J. (2013). Dynamics of fluids in porous media. DoverPublications. com.
Benzi, R., Succi, S., & Vergassola, M. (1992). The lattice Boltzmann equation: theory and
applications. Physics Reports, 222(3), 145-197.
Bera, B., Mitra, S. K., & Vick, D. (2011). Understanding the micro structure of Berea Sandstone
by the simultaneous use of micro-computed tomography (micro-CT) and focused ion beam-
scanning electron microscopy (FIB-SEM). Micron, 42(5), 412-418.
Bhatnagar, P. L., Gross, E. P., & Krook, M. (1954). A model for collision processes in gases. I.
Small amplitude processes in charged and neutral one-component systems. Physical
review, 94(3), 511.
Blunt, M. J., Jackson, M. D., Piri, M., & Valvatne, P. H. (2002). Detailed physics, predictive
capabilities and macroscopic consequences for pore-network models of multiphase
105

flow. Advances in Water Resources, 25(8), 1069-1089.
Blunt, M., & King, P. (1990). Macroscopic parameters from simulations of pore scale
flow. Physical Review A, 42(8), 4780.
Bracke, M., De Voeght, F., & Joos, P. (1989). The kinetics of wetting: the dynamic contact
angle. In Trends in Colloid and Interface Science III (pp. 142-149). Steinkopff.
Brunauer, S., Emmett, P. H., & Teller, E. (1938).Adsorption of gases in multimolecular
layers. Journal of the American Chemical Society, 60(2), 309-319.
Cancelliere, A., Chang, C., Foti, E., Rothman, D. H., & Succi, S. (1990). The permeability of a
random medium: Comparison of simulation with theory.Physics of Fluids A Fluid
Dynamics, 2(12), 2085.
Caupin, Frédéric, et al. "Absolute limit for the capillary rise of a fluid." EPL (Europhysics
Letters) 82.5 (2008): 56004.
Chen, S., & Doolen, G. D. (1998). Lattice Boltzmann method for fluid flows.Annual review of
fluid mechanics, 30(1), 329-364.
Chen, S., Diemer, K., Doolen, G. D., Eggert, K., Fu, C., Gutman, S., & Travis, B. J. (1991).
Lattice gas automata for flow through porous media. Physica D: Nonlinear
Phenomena, 47(1), 72-84.
Clarkson, C. R., Solano, N., Bustin, R. M., Bustin, A. M. M., Chalmers, G. R. L., He, L., ... &
Blach, T. P. (2013). Pore structure characterization of North American shale gas reservoirs
using USANS/SANS, gas adsorption, and mercury intrusion. Fuel, 103, 606-616.
Comisky, J., Newsham, K., Rushing, J., & Blasingame, T. (2007, November). A comparative
study of capillary-pressure-based empirical models for estimating absolute permeability in
tight gas sands. SPE Annual Technical Conference and Exhibition.
Comisky, J., Newsham, K., Rushing, J., & Blasingame, T. (2007, November). A comparative
106

study of capillary-pressure-based empirical models for estimating absolute permeability in
tight gas sands. SPE Annual Technical Conference and Exhibition.
Curtis, M., Ambrose, R., Sondergeld, C., & Sondergeld, C. (2010, October). Structural
characterization of gas shales on the micro-and nano-scales. In Canadian Unconventional
Resources and International Petroleum Conference.
Darabi, H., Ettehad, A., Javadpour, F., & Sepehrnoori, K. (2012). Gas flow in ultra-tight shale
strata. Journal of Fluid Mechanics, 710, 641.
Darcy, H. (1856). Les fontaines publiques de la ville de Dijon: exposition et application... Victor
Dalmont.
De Gennes, P. G., Brochard-Wyart, F., & Quéré, D. (2004). Capillarity and wetting phenomena:
drops, bubbles, pearls, waves. Springer.
de Swaan, A. (1978). Theory of waterflooding in fractured reservoirs. Society of Petroleum
Engineers Journal, 18(02), 117-122.
Delker, T., Pengra, D. B., & Wong, P. Z. (1996). Interface pinning and the dynamics of capillary
rise in porous media. Physical review letters, 76(16), 2902.
d'Humières, D. (2002). Multiple–relaxation–time lattice Boltzmann models in three
dimensions. Philosophical Transactions of the Royal Society of London. Series A:
Mathematical, Physical and Engineering Sciences, 360(1792), 437-451.
Dimitrov, D. I., Milchev, A., & Binder, K. (2007). Capillary rise in nanopores: molecular
dynamics evidence for the Lucas-Washburn equation. Physical review letters, 99(5),
054501.
Dong, H., & Blunt, M. J. (2009). Pore-network extraction from micro-computerized-tomography
images. Physical review E, 80(3), 036307.
Elgmati, M., Zhang, H., Bai, B., Flori, R., &Qu, Q. (2011, June). Submicron-pore
107

characterization of shale gas plays. In North American Unconventional Gas Conference and
Exhibition.
Erickson, D., Li, D., & Park, C. B. (2002). Numerical simulations of capillary-driven flows in
nonuniform cross-sectional capillaries. Journal of colloid and interface science, 250(2),
422-430.
Fatt, I. (1956). The network model of porous media I. Capillary pressure characteristics. Trans.
AIME, 207(7), 144-159.
Ferreol, B., & Rothman, D. H. (1995). Lattice-Boltzmann simulations of flow through
Fontainebleau sandstone. Transport in porous media, 20(1-2), 3-20.
Gad-el-Hak, M. (Ed.). (2010). MEMS: introduction and fundamentals. CRC press.
Grad, H. (1949). Note on N-dimensional hermite polynomials. Communications on Pure and
Applied Mathematics, 2(4), 325-330.
Gunstensen, A. K., Rothman, D. H., Zaleski, S., & Zanetti, G. (1991). Lattice Boltzmann model
of immiscible fluids. Physical Review A, 43(8), 4320.
Guodong, J., Patzek, T., &Silin, D. (2004, September). Direct prediction of the absolute
permeability of unconsolidated and consolidated reservoir rock. In SPE Annual Technical
Conference and Exhibition.
Haldorsen, H., & Lake, L. (1984). A new approach to shale management in field-scale
models. Old SPE Journal, 24(4), 447-457.
Handy, L. L. (1960). Determination of effective capillary pressures for porous media from
imbibition data.
He, X., & Luo, L. S. (1997). Lattice Boltzmann model for the incompressible Navier–Stokes
equation. Journal of Statistical Physics, 88(3-4), 927-944.
Higuera, F. J., & Jimenez, J. (1989). Boltzmann approach to lattice gas simulations. EPL
108

(Europhysics Letters), 9(7), 663.
Higuera, F. J., Succi, S., & Benzi, R. (1989). Lattice gas dynamics with enhanced collisions. EPL
(Europhysics Letters), 9(4), 345.
Ho, C. M., & Tai, Y. C. (1998). Micro-electro-mechanical-systems (MEMS) and fluid
flows. Annual Review of Fluid Mechanics, 30(1), 579-612.
Huang, H., & Lu, X. Y. (2009). Relative permeabilities and coupling effects in steady-state gas-
liquid flow in porous media: A lattice Boltzmann study.Physics of Fluids, 21, 092104.
Huang, K. (1963), Statistical Mechanics, John Wiley & Sons, Inc.
Hyväluoma, J., Raiskinmäki, P., Jäsberg, A., Koponen, A., Kataja, M., & Timonen, J. (2006).
Simulation of liquid penetration in paper. Physical Review E, 73(3), 036705.
Javaheri, M., & Jessen, K. (2013, February). CO Mobility and Transitions Between Co-current
and Counter-Current Flows. In 2013 SPE Reservoir Simulation Symposium.
Jin G., Patzek T.W. and Silin D.B. (2004). Direct Prediction of Absolute Permeability of
Unconsolidated and Consolidated Reservoir Rock, SPE 90084
Kadau, K., Germann, T. C., & Lomdahl, P. S. (2006). Molecular dynamics comes of age: 320
billion atom simulation on BlueGene/L. International Journal of Modern Physics C, 17(12),
1755-1761.
Kardar, M. (2007). Statistical physics of particles. Cambridge University Press.
Kawazoe, H. G. (1983). Method for the calculation of effective pore size distribution in
molecular sieve carbon. Journal of Chemical Engineering of Japan, 16(6), 470-475.
Knackstedt, M. A., Latham, S., Madadi, M., Sheppard, A., Varslot, T., & Arns, C. (2009).
Digital rock physics: 3D imaging of core material and correlations to acoustic and flow
properties. The Leading Edge, 28(1), 28-33.
109

Koch, G. W., Sillett, S. C., Jennings, G. M., & Davis, S. D. (2004). The limits to tree
height. Nature, 428(6985), 851-854.
Kolev, N. I. (2007). Multiphase Flow Dynamics 2: Mechanical Interactions (Vol. 2). Springer.
Kuzmin, A. (2010). Multiphase simulations with lattice boltzmann scheme.
Lemmens, H., & Butcher, A. (2011, October). A Workflow Concept to Characterize Core
Samples From The Microscale To The Nanoscale. In SPE Annual Technical Conference
and Exhibition.
Lake, L. W. (1989). Enhanced oil recovery.
Li, K., & Horne, R. N. (2001, September). Steam-water relative permeability by the capillary
pressure method. In proceedings of the International Symposium of the Society of Core
Analysts, Edinburgh, UK.
Liou, W. W., Peng, Y., & Parker, P. E. (2009). Analytical modeling of capillary flow in tubes of
nonuniform cross section. Journal of colloid and interface science, 333(1), 389-399.
Ma, S., Zhang, X., Morrow, N. R., & Zhou, X. (1999). Characterization of wettability from
spontaneous imbibition measurements. Journal of Canadian Petroleum Technology, 38(13),
1-8.
MacVicar, B. A., & Salter, M. W. (2006). Neuroscience: controlled
capillaries.Nature, 443(7112), 642-643.
Masoodi, R., Languri, E., & Ostadhossein, A. (2013). Dynamics of liquid rise in a vertical
capillary tube. Journal of colloid and interface science, 389(1), 268-272.
McCurdy, E. (1906). Leonardo da Vinci's note-books. Duckworth.
McNamara, G. R., & Zanetti, G. (1988). Use of the Boltzmann equation to simulate lattice-gas
automata. Physical Review Letters, 61(20), 2332.
110

Neuzil, C. E. (1994). How permeable are clays and shales? Water resources research, 30(2),
145-150.
Pan, C., Hilpert, M., & Miller, C. T. (2004). Lattice Boltzmann simulation of two phase flow in
porous media. Water Resources Research, 40(1).
Pelesko, J. A., & Bernstein, D. H. (2002). Modeling Mems and Nems. CRC press.
Piri, M. (2003). Pore-scale modelling of three-phase flow (Doctoral dissertation, Imperial
College).
Pomeau, B. H. Y., & Frisch, U. (1986). Lattice-gas automata for the Navier-Stokes
equation. Phys. Rev. Lett, 56(14), 1505.
Ponomarenko, A., Quéré, D., & Clanet, C. (2011). A universal law for capillary rise in
corners. Journal of Fluid Mechanics, 666, 146-154.
Qian, Y. H., d'Humières, D., & Lallemand, P. (1992). Lattice BGK models for Navier-Stokes
equation. EPL (Europhysics Letters), 17(6), 479.
Rahimian, M. H., & Pourshaghaghy, A. (2002). Direct simulation of forced convection flow in a
parallel plate channel filled with porous media. International communications in heat and
mass transfer, 29(6), 867-878.
Riepe, L., Suhaimi, M., Kumar, M., &Knackstedt, M. (2011, January) Application of High
Resolution Micro-CT-Imaging and Pore Network Modeling (PNM) for the Petrophysical
Characterization of Tight Gas Reservoirs-A Case History from a Deep Clastic Tight Gas
Reservoir in Oman. In SPE Middle East Unconventional Gas Conference and Exhibition.
Rothman, D. H., & Zaleski, S. (1994). Lattice-gas models of phase separation: interfaces, phase
transitions, and multiphase flow. Reviews of Modern Physics,66(4), 1417.
Rothman, D. H., & Zaleski, S. (2004). Lattice-gas cellular automata: simple models of complex
111

hydrodynamics (Vol. 5). Cambridge University Press.
Roychaudhuri, B., Tsotsis, T. T., & Jessen, K. (2013). An experimental investigation of
spontaneous imbibition in gas shales. Journal of Petroleum Science and Engineering, 111,
87-97.
Sahimi, M. (2012). Flow and transport in porous media and fractured rock. Wiley. com.
Scheidegger, A. E. (1958). The physics of flow through porous media. Soil Science, 86(6), 355.
Shafer, J., & Neasham, J. (2000). Mercury Porosimetry Protocol for Rapid Determination of
Petrophysical and Reservoir Quality Properties. In International Symposium of the Society
of Core Analysts (pp. 18-22).
Shah, R. K. "A correlation for laminar hydrodynamic entry length solutions for circular and
noncircular ducts." Journal of Fluids Engineering 100.2 (1978): 177-179.
Shan, X., & Chen, H. (1993). Lattice Boltzmann model for simulating flows with multiple
phases and components. Physical Review E, 47(3), 1815.
Shouxiang, M., Morrow, N. R., & Zhang, X. (1997). Generalized scaling of spontaneous
imbibition data for strongly water-wet systems. Journal of Petroleum Science and
Engineering, 18(3), 165-178.
Sing, K. (2001). The use of nitrogen adsorption for the characterisation of porous
materials. Colloids and Surfaces A: Physicochemical and Engineering Aspects, 187, 3-9.
Succi, S. (2001). The lattice Boltzmann equation: for fluid dynamics and beyond. Oxford
university press.
Succi, S., Foti, E., & Higuera, F. (1989). Three-dimensional flows in complex geometries with
the lattice Boltzmann method. EPL (Europhysics Letters),10(5), 433.
Sukop, M. C., & Thorne Jr, D. T. (2007). Lattice Boltzmann modeling: an introduction for
geoscientists and engineers. Springer Publishing Company, Incorporated.
112

Swanson, B. F. (1981). A simple correlation between permeabilities and mercury capillary
pressures. Journal of Petroleum Technology, 33(12), 2498-2504.
Swift, M. R., Osborn, W. R., & Yeomans, J. M. (1995). Lattice Boltzmann simulation of
nonideal fluids. Physical Review Letters, 75(5), 830.
Takahashi, S., & Kovscek, A. R. (2010). Spontaneous countercurrent imbibition and forced
displacement characteristics of low-permeability, siliceous shale rocks. Journal of
Petroleum Science and Engineering, 71(1), 47-55.
Valluri, S. R., Jeffrey, D. J., & Corless, R. M. (2000). Some applications of the Lambert W
function to physics. Canadian Journal of Physics, 78(9), 823-831.
Wagner, A. J. (2008). A practical introduction to the lattice Boltzmann method.Adt. notes for
Statistical Mechanics, 463, 663.
Warren, J. E., and P. J Root. "The behavior of naturally fractured reservoirs." SPE Journal, 1963.
Washburn, Edward W. "The dynamics of capillary flow." Physical review 17.3 (1921): 273.
Weislogel, Mark M. "Compound capillary rise." Journal of Fluid Mechanics 709 (2012): 622-
647.
World Energy Outlook 2014 (EIA, 2014).
Yiotis, A. G., Psihogios, J., Kainourgiakis, M. E., Papaioannou, A., & Stubos, A. K. (2007). A
lattice Boltzmann study of viscous coupling effects in immiscible two-phase flow in porous
media. Colloids and Surfaces A: Physicochemical and Engineering Aspects, 300(1), 35-49.
Zhmud, B. V., Tiberg, F., & Hallstensson, K. (2000). Dynamics of capillary rise. Journal of
Colloid and Interface Science, 228(2), 263-269.
Zou, Q., & He, X. (1996). On pressure and velocity flow boundary conditions and bounceback
for the lattice Boltzmann BGK model. arXiv preprint comp-gas/9611001.
113

8. Appendix-A: Maxwell-Boltzmann Equilibrium Distribution Function

Going back to the setup shown in Figure 3.1 at equilibrium, the number of collisions on one face
(the positive x-plane) of the cube is:

,
A- 1

where

is the number density of the particles, is the area and is time. The particles
with positive x-velocity will have a chance to collide with the positive x wall. By hitting the
wall, particles bounce back and move in the opposite direction. The change in momentum is:

A- 2

Pressure is defined as the normal force divided by area which is related to momentum as:

A- 3

Since we are dealing with a large number of particles and no net forces are applied, then the
number of particles with positive velocity is half of the total number of particles.

A- 4

Assuming that the gas is ideal and replacing :

A- 5

The sum is replaced with a statistical notation:
114

A- 6

Replacing for number density:

A- 7

Cancelling volume:

A- 8

Writing the number of molecules as the product of number of moles and the Avogadro number:

A- 9

Cancelling number of moles:

A- 10

Without loss of generality the equation can be written in three dimensions (

)

,
A- 11

where

is the root-mean-square speed of molecules. Accordingly, the kinetic energy will be:

,
A- 12

which states that the mean kinetic energy of an ideal gas is a function of its temperature only.
This will be used during the next parts of derivation.
Returning to the main problem:
115

A- 13

Taking the derivative of

w.r.t.

A- 14

Taking the derivative of the r.h.s. gives:

A- 15

or

A- 16

Substituting for

:

A- 17

Taking derivative:

A- 18

By analogy:

A- 19

A- 20

This suggests that the left hand side of equations 3-28, 3-29 and 3-30 is constant. For example
from equation 3-30:
116

A- 21

Dividing by

:

A- 22

Separating the first order ordinary differential equation:

A- 23

Integrating:

A- 24

Here is the integration constant. The integral of probability over the domain is equal to one:

A- 25

The value of this integral can be numerically estimated as:

A- 26

Rearranging:

A- 27

Replacing in equation 3-34:
117

A- 28

Determining the constant :

A- 29

A- 30

A- 31

A- 32

A- 33

A- 34

A- 35

Replacing for in equation 3-38:

,
A- 36

Where is the molecular mass.
118

And finally the equilibrium distribution function for a dimensional system is the multiplication
of distribution function for different dimensions:

,
A- 37

Where is the velocity of molecules relative to the container and can be expanded as .
is the microscopic velocity and is the macroscopic velocity.

A- 38

119

9. Appendix-B: Derivation of D
2
Q
9
Scheme

To derive D
2
Q
9
, is defined as (He and Lou, 1997):

B - 1

where

and

are the and components of molecular velocity ( ).
Rewriting the general hydrodynamic integral:

B - 2

The system is two dimensional (D=2) and

:

B - 3

Changing variables as:

B - 4

B - 5

120

B - 6

Separating the integral

where:

B - 7

B - 8

Solving for

:
B - 9

B - 10

Solving for

:

121

B - 11

Solving for

:

122

B - 12

B - 13

Where:

B - 14

Using Gauss-Hermite quadrature, this integral can be estimated as:

B - 15

The weights are:

B - 16

123

is the n
th
Hermite polynomial (Grad, 1949). Here and the Hermite polynomial, in the
weight function, is

Substituting

, the different values of will be

and the weights are

.
The integral of moments (Equation B-13) reduces to:

B - 17

Comparing with Equation B-1, one can see that:

and
the weights corresponding to each direction are

Replacing

and adopting the new notation:

B - 18

Equation

along with equation B-18, weights and
velocities in LBM can be solved numerically to get the equilibrium distribution functions and
consequently macroscopic properties.

Abstract (if available)

Abstract Flow and transport in porous rocks is controlled by the pore structures and surface chemistry that, in part, defines a porous material. Pores are commonly in the micron range for conventional oil and gas formations while often in the nano-meter range for unconventional resources. Traditionally, core experiments are performed for characterization of reservoir rocks. Recently, with advances in imaging and computational power, we are able to image the pores in most relevant rocks and to construct a digital rendering of the rock that can subsequently be used for modeling of flow and transport processes. ❧ In this manuscript, we introduce a workflow that takes the rock images as input. In the first stage of the workflow, a 3D image is analyzed to find effective and total porosities. Single-phase flow simulation, using the Lattice Boltzmann Method (LBM), is then performed to evaluate directional (e.g. x, y and z) permeability of the rock sample. Finally, two-phase flow simulations are performed on the pore geometries dictated by the sample image to determine multiphase characteristics, including relative permeability of the rock sample. ❧ To place imaging technology in context, we compare and discuss the insight and understanding about the pore structure characteristics of rocks, specifically for Monterey shale samples that can be obtained from various characterization techniques. The characterization techniques considered here include: (i) Mercury injection capillary pressure experiments (MICP), (ii) Brunauer, Emmett and Teller (BET) – Nitrogen adsorption experiments, (iii) High-Resolution X-ray CT (HRXCT), and (iv) Focused Ion Beam – Scanning Electron Microscopy (FIB-SEM). The HRXCT and FIB-SEM characterization data are utilized to reconstruct (visualize) the 3D structures of the rocks

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

PDF

Modeling and simulation of complex recovery processes

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

PDF

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

PDF

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

PDF

Mass transfer during enhanced hydrocarbon recovery by gas injection processes

PDF

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

PDF

Multiscale and multiresolution approach to characterization and modeling of porous media: From pore to field scale

PDF

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

PDF

Asphaltene deposition during co₂ injection in enhanced oil recovery applications

PDF

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

PDF

Chemical and mechanical deformation of porous media and materials during adsorption and fluid 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

PDF

Investigation of adsorption and phase transition phenomena in porous media

PDF

Molecular- and continuum-scale simulation of single- and two-phase flow of CO₂ and H₂O in mixed-layer clays and a heterogeneous sandstone

PDF

Flow and thermal transport at porous interfaces

PDF

Investigation of gas transport and sorption in shales

PDF

The projection immersed boundary method for compressible flow and its application in fluid-structure interaction simulations of parachute systems

PDF

Machine-learning approaches for modeling of complex materials and media

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

PDF

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

PDF

The effect of lattice structure and porosity on thermal conductivity of additively-manufactured porous materials

Asset Metadata

Creator Farhadi Nia, Shahram (author)

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

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Petroleum Engineering

Publication Date 11/25/2014

Defense Date 10/20/2014

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

Tag digital rock physics,fluid flow,fractured,image analysis,imbibition,Lattice Boltzmann method,Monterey formation,OAI-PMH Harvest,pore scale,relative permeability,Rock,simulation,unconventional reservoirs

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Jessen, Kristian (committee chair)

Creator Email farhadin@usc.edu,shahram.farhadi@gmail.com

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

Unique identifier UC11297373

Identifier etd-FarhadiNia-3100.pdf (filename),usctheses-c3-518299 (legacy record id)

Legacy Identifier etd-FarhadiNia-3100.pdf

Dmrecord 518299

Document Type Dissertation

Format application/pdf (imt)

Rights Farhadi Nia, Shahram

Type texts

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

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

Repository Name University of Southern California Digital Library

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