EXPERIMENTAL CHARACTERIZATION OF PORE-SCALE CAPILLARY PRESSURE
AND CORNER FILM FLOW IN 2D POROUS MICROMODELS
by
Razin Sazzad Molla
A thesis submitted in partial fulfillment
of the requirements for the degree
of
Master of Science
in
Mechanical Engineering
MONTANA STATE UNIVERSITY
Bozeman, Montana
May 2023
©COPYRIGHT
by
Razin Sazzad Molla
2023
All Rights Reserved
ii
ACKNOWLEDGEMENTS
I would like to thank my advisor Dr. Yaofa Li for the guidance and support throughout
my time at Montana state university, without which this work would not have been possible.
I would also like to thank my thesis reading committee members, Prof. Sarah Codd and
Prof. Eric Johnson for their effort in reviewing this work and their valuable advice and
remarks.
I would like to thank my current and previous lab members in the microfluidics lab
for energy and fluid transport (M-LEFT), Mr. Nishagar Raventhiran, Mr. Rafid Musabbir
Rahman, Mr. Md Ahsan Habib, Ms. Deanna Reynolds and Mr. Elliott Niemur for sharing
the invaluable ideas and helping with setting up the experimental apparatus with me. I
would like to extend my appreciation to Montana Microfabrication Facility (MMF), and
staff members Dr. Andrew Lingley and Dr. Joshua Heinemann for their help with my
sample fabrication.
Finally, I want to thank my family and friends for their support.
iii
TABLE OF CONTENTS
1. CHAPTER 1 INTRODUCTION.............................................................................. 1
Motivation ............................................................................................................. 1
Objectives ............................................................................................................ 14
Outline ................................................................................................................ 14
2. CHAPTER 2 BACKGROUND AND LITERATURE REVIEW ............................... 16
Definitions and Theoretical Background ................................................................. 16
Scaling and Representative Elementary Volume............................................... 16
Porosity ........................................................................................................ 19
Saturation..................................................................................................... 19
Phase and Component: .................................................................................. 20
Interfacial Tension and the Young-Laplace Equation........................................ 21
Contact Angle and Young’s Equation ............................................................. 23
Capillary Pressure ......................................................................................... 24
Drainage and Imbibition ................................................................................ 27
Darcy’s equation ........................................................................................... 28
Literature Review ................................................................................................. 30
Capillary Pressure Characterization in Porous Media....................................... 30
Two-Phase Immiscible Displacement and Film Flow in Porous Media ............... 39
3. CHAPTER 3 EXPERIMENTAL METHOD ........................................................... 52
Micromodels......................................................................................................... 52
Micromodel Design........................................................................................ 53
Micromodel Fabrication ................................................................................. 60
Experimental Methods for Capillary Pressure Quantification ................................... 68
Flow Apparatus and Experimental Procedures ................................................ 68
Image Acquisition ......................................................................................... 70
Contact Angle Measurement .......................................................................... 71
Interfacial Tension Measurement .................................................................... 72
Magnification Calibration .............................................................................. 73
Image Processing........................................................................................... 73
Circle Fit Algorithm...................................................................................... 76
Experimental Method for Film Flow Characterization............................................. 78
Flow Apparatus and Experimental Procedure ................................................. 78
Image Acquisition ......................................................................................... 79
Contact Angle Measurement: ......................................................................... 79
Interfacial Tension Measurement: ................................................................... 79
iv
TABLE OF CONTENTS – CONTINUED
Calibration ................................................................................................... 80
Image Processing........................................................................................... 80
4. CHAPTER 4 PORE-SCALE CAPILLARY PRESSURE QUANTIFICA-
TION................................................................................................................... 83
Capillary Pressure in Homogeneous Micromodels .................................................... 83
Effect of Spatial Resolution............................................................................ 88
Effect of Heterogeneity .................................................................................. 92
Effect of Phase Connectivity .......................................................................... 97
5. CHAPTER 5 FILM FLOW CHARACTERIZATION .............................................104
Single-Phase Flow Validation................................................................................104
Comparison Between Calculated and Applied Flow Rate ................................104
Single Phase Velocity field ............................................................................105
Drainage and Imbibition in Homogeneous Micromodel ...........................................109
Drainage......................................................................................................109
Flow rate vs. time ........................................................................................109
Velocity Fields .............................................................................................111
Film Flow....................................................................................................115
Imbibition....................................................................................................117
Flow rate vs. time ........................................................................................117
Velocity Fields .............................................................................................119
Film Flow....................................................................................................121
Drainage and Imbibition in Heterogeneous Micromodel ..........................................123
Drainage......................................................................................................123
Flow rate vs. time ........................................................................................123
Velocity plots ...............................................................................................125
Film flow .....................................................................................................127
Imbibition....................................................................................................128
6. CONCLUSION....................................................................................................130
REFERENCES CITED.............................................................................................131
v
LIST OF TABLES
Table Page
3.1 Description of the micromodels used........................................................... 67
3.2 Details of the objectives ............................................................................ 72
vi
LIST OF FIGURES
Figure Page
1.1 Application examples of multiphase flow in porous media.
Schematics showing from top left clockwise: Generalized
microfluidic devices such as Lab on a Chip (LOC), repro-
duced from Ref. [137]; Enhanced Oil recovery in operation
such as steam injection; Contamination from industrial
sources into the subsurface, reproduced from Ref. [47];
Hydrology such as the water cycle involves multiphase flow
in porous media; Microbial fuel cell (MFC) involves porous
filament thus flow through porous materials, reproduced
from Ref. [171]; CO2 sequestrations from major sources;
Blood and bile flow through liver lobules can be simulated
as a flow through porous media, reproduced from Ref. [126];
Heat sinks with improved performances applying porous
media, reproduced from Ref. [146] ................................................................ 2
1.2 Different oil entrapment at the pore scale. Here the rock
grains and oil are in brown and black, respectively. (a) oil
film adhered to pore walls; (b) oil trapped in dead ends;
(c) oil trapped due to heterogeneous porous media; and
(d) oil ganglia in pore throats. Figure was reproduced based
on work in Ref. [43]. .................................................................................... 4
1.3 (a) two-phase flow with trapped water (wetting phase), here
blue is connected water, red is trapped water, black is solid
grain and green is connected n-heptane (nonwetting phase);
(b) zoomed in view of the interfaces, connected interfaces in
yellow and disconnected ones in purple. (c) disconnected n-
heptane as yellow and blue and black represents connected
water and solid grain respectively as before................................................... 8
vii
LIST OF FIGURES – CONTINUED
Figure Page
1.4 Images illustrating scenarios where film flow plays an role:
(a) 3D phase distribution obtained by numerical simulation
with brine and super-critical CO2 under strong imbibition
(wetting phase contact angle close to zero). Precursor
wetting film ahead of the main meniscus dominates the flow
(Figure was reproduced based on work in Ref. [78]); (b) µ-
PIV images of the film flow in a homogeneous micromodel.
Water is displacing air from left to right, where water
and air appear bright and dark respectively. (c) Primary
drainage in a heterogeneous micromodel at a fixed boundary
pressure. It is visible that initially trapped water drains
away over time via corner film flow that is connected to the
outlet. ...................................................................................................... 12
2.1 Different scales in porous media with relevant applications.
Figure was reproduced based on work in Ref. [76]. ...................................... 17
2.2 Volume-averaged porosity averaged over cubes of different
length l, for the statistical realization of Fontainebleau
sandstone. This shows the effect of REV. As the averaging
volume is increased the porosity value converges. Figure
was reproduced based on work in Ref. [18]. ................................................. 18
2.3 Two phase system with an ideal interface. .................................................. 22
2.4 Two fluid phases on a solid surface with different wettabil-
ities. Figure was reproduced based on work in Ref. [18]. .............................. 24
2.5 Water and oil in contact with a solid, contact angle is
measured from the denser fluid. Figure was reproduced
based on work in Ref. [18].......................................................................... 25
2.6 hysteretic capillary pressure, P c as a function of wetting
phase saturation. Figure was reproduced based on work in Ref. [56]. ............ 29
2.7 A nonhysteretic model applied and predicted nonwetting
phase volume fraction with the actual volume fraction for
two different models (a and b). (cf. first and second line of
table 5 in [120]).Figure was reproduced based on work in Ref. [120].............. 35
viii
LIST OF FIGURES – CONTINUED
Figure Page
2.8 Equilibrium data of saturation, capillary pressure, and
specific interfacial area for drainage and imbibition. (a)
actual data points (b) fitted data. Figure was reproduced
based on work in Ref. [92].......................................................................... 36
2.9 Curvature-based capillary pressure vs. transducer-based
measured capillary pressure. Figure was reproduced based
on work in Ref. [141]. ................................................................................ 37
2.10 Curvatures measured from micro-CT images of the con-
nected nonwetting phase (air), and curvatures estimated
from transducer-based capillary pressure are compared.
Figure was reproduced based on work in Ref. [71]. ...................................... 39
2.11 Phase diagram for drainage (viscosity ratio M = µnw/µw;
nw and w refers to the nonwetting/advancing and wetting
defending phases respectively, Capillary number Ca =
µnwunw/γ(nw−w)cosθ; unw is the velocity of the advancing
phase, γ(nw−w) is the interfacial tension between the two
phase and θ is the contact angle. Three stability areas are
bounded by dashed lines and the locations of the PEG200
and water (two different wetting phase used) displacement
experiments are represented by pink squares and blue
circles respectively. The gray zones represent the stability
areas originally indicated by Lenormand et al. This works
extends the original work of Lenormand et al. Figure was
reproduced based on work in Ref. [176]....................................................... 41
2.12 Quasi steady state CO2 -water phase distribution for
drainage experiments at various capillary number with the
same viscosity ratio which exhibits the gradual transition
from capillary fingering to viscous fingering as the flow rate
increases. The invading phase CO2, the resident phase
water, and the solid grains are shown in green, blue, and
black, respectively. Flow is from left to right. Figure was
reproduced based on work in Ref. [107]....................................................... 42
ix
LIST OF FIGURES – CONTINUED
Figure Page
2.13 (a) weak imbibition (contact angle 60 degree) where water
displacing silicone oil, cooperative pore filling is shown, (b)
strong imbibition (contact angle 7 degree) shows corner
flow (c) with decreasing contact angle interfacial curvature
forced into the corners (d) strong imbibition with film
swelling over time. Figure was reproduced based on work in Ref. [178] ......... 44
2.14 Micro-PIV measurement in a 2D micromodel experiment
of supercritical CO2 and water. Quantitative velocity
information can tell a lot more about flow unsteadiness
and preferential pathways. Figure was reproduced based
on work in Ref. [94] ................................................................................... 45
2.15 Meniscus displacement mechanism during imbibition (orig-
inal figure from Lenormand et al. 1984). Figure was
reproduced based on work in Ref. [104]....................................................... 46
2.16 Imbibition considering film flow and micro roughness using
a pore network simulator (pore wall microroughness and
wetting film thickness are exaggerated for illustrative pur-
poses). Figure was reproduced based on work in Ref. [35]. ........................... 48
2.17 Film thickness swelling and eventually snap off in a triangu-
lar pore, black is the liquid and white is the vapor. Figure
was reproduced based on work in Ref. [164]. ............................................... 49
2.18 Final trapping nonwetting phase distribution and satura-
tion levels for different injection rates (capillary numbers).
(a),(b) and (c) are imbibition experiments done by Tsao
et al. [163] in a 2D 4-layer micromodel. (d),(e) and (f)
are the ICB-DPNM simulation by Zhao et. al (2021) with
the same geometry and conditions. (g) is the comparison
of trapped nonwetting phase (air) saturation vs capillary
number for both the experiments and the dynamic pore
network modeling. Figure was reproduced based on work in Ref. [180]. ........ 50
x
LIST OF FIGURES – CONTINUED
Figure Page
2.19 Variation of liquid distribution in the pore network during
evaporation obtained from experiment, (Direct pore net-
work modeling with corner films) DPNM wC, and (Direct
pore network modeling without corner films) DPNM woC.
St is the total liquid saturation in the pore network. Figure
was reproduced based on work in Ref. [168]. ............................................... 51
3.1 A workflow chart illustrating the fabrication process used
to make the micromodels in this study........................................................ 53
3.2 (a) The heterogeneous micromodel mask for the capillary
pressure experiment (Type II), (b) single micromodel de-
sign, (c) dimensions of the porous section with the air barrier. ..................... 54
3.3 (a) Design of the idealized micromodel mask for the cap-
illary pressure experiment (Type I), (b) dimensions of the
porous section with the air barrier.............................................................. 55
3.4 Micromodels used for the film flow characterization ex-
periment: (a) and (c) single micromodel design for the
homogeneous and heterogeneous micromodel respectively
(b) and (d) dimensions of the porous section with the air
barrier for the homogeneous and heterogeneous micromodel
respectively. .............................................................................................. 55
3.5 Nonwetting phase barrier outlined in red at the end of the
porous section.(a) large gap for PIV experiment micromodel
(b)smaller barrier gap in heterogeneous design for capillary
pressure experiment................................................................................... 56
xi
LIST OF FIGURES – CONTINUED
Figure Page
3.6 (a) Pore body (radius of pore body) and (b) pore throat
(throat width) distribution of the heterogeneous porous
medium for the capillary pressure quantification experiment
(c) Pore body (radius of pore body) and (d) pore throat
(throat width) distribution of the idealized porous medium
for the capillary pressure quantification experiment (e) Pore
body (radius of pore body) and (f) pore throat (throat
width) distribution of the homogeneous porous medium for
the film flow characterization experiment (g) Pore body
(radius of pore body) and (h) pore throat (throat width)
distribution of the heterogeneous porous medium for the
film flow characterization experiment.......................................................... 57
3.7 (a) Cityblock distance method and watershed segmentation
separates pore body and throats. Smallest distance between
grains gives the throat width and equivalent pore body
radius is calculate√d from approximating each pore body
area as a circle ( (πr2/π)). For a detailed discussion
of pore size distribution from digital images readers are
referred to the paper in Ref. [143] .............................................................. 58
3.8 (a) Bare silicon wafer (b) dehydration bake (c) HMDS vapor
prime oven to coat HMDS layer (d) Laurel spin coater (e)
spin coating PR (f) soft bake. .................................................................... 60
3.9 (a) Photomask (b) ABM contact aligner (c) aligning the
substrate and photomask (d) UV exposure (e) developing
(f) cleaning with DI water. ........................................................................ 61
3.10 Features under bright field microscope after completing
photolithography:(a) and (b) homogeneous design;(c) and
(d) heterogeneous design. ......................................................................... 63
3.11 Dry etching process: (a) Oxford PlasmaLab 100; (b) load-
ing the wafer; (c) plasma striking; (d) etching completed.
................................................................................................................ 64
3.12 SEM image of the idealized micromodel after etching and
photoresist removal. ................................................................................. 65
xii
LIST OF FIGURES – CONTINUED
Figure Page
3.13 SEM image of the heterogeneous micromodel after etching
and photoresist removal. .......................................................................... 65
3.14 Post lithography and etching operations: (a) drilling holes;
(b) piranha cleaning; (c) anodic bonding (d) dicing and (e)
final micromodel after nano-port attachment. ............................................ 66
3.15 Schematic diagram showing the setup to apply a fixed
pressure at the inlet by controlling the tank height. ................................... 69
3.16 Schematic diagram showing the n-heptane injection through
the cut-off valve. ...................................................................................... 70
3.17 Photo illustrating the optical setup used in the experiment. ......................... 71
3.18 A closeup view of the optical path inside the microscope. ............................ 72
3.19 Static contact angle between DI water n-heptane measured
on a silicon wafer of the same surface properties as that of
the micromodel. For the micromodels used in this study,
the contact angle was measured to be ∼30 ◦. ............................................... 73
3.20 (a) mask image (b) RAW image (c) segmented binarized
(green is n-heptane, black is grain and blue is water). .................................. 74
3.21 (a) n-Heptane-water interfaces as white (b) zoomed in view
of image a................................................................................................. 75
3.22 (a) Fitted circle on every interface (b) zoomed-in view of
image a (fitted circles are in green)............................................................. 76
3.23 Static contact angle between DI-water and air measured on
a silicon wafer of the same surface properties as that of the
micromodel. The static contact angle was determined to be 5 ◦. ................... 80
3.24 Interfacial tension measurement using the pendant drop
method between air and DI water. ............................................................ 81
3.25 Pixel to mm scaling for the PIV and PTV experiments. ............................. 81
xiii
LIST OF FIGURES – CONTINUED
Figure Page
4.1 Equilibrium phase distributions during drainage at different
applied boundary pressures. Here, n-heptane displaces
water from right to left. Images shown here were captured
at 5x magnification, where, green, blue and black represent
n-heptane, water and solid grains, respectively. ........................................... 84
4.2 Equilibrium phase distributions during imbibition at differ-
ent applied boundary pressures. Similar to Figure 4.1, but
for the imbibition cases.............................................................................. 85
4.3 (a) Major steps followed to get the radius of curvature;
(b) the curvature in the depth direction, R2 was approx-
imated based on the contact angle in the horizontal plane.
Image is reproduced based on the work in Ref. [90]. .................................... 86
4.4 Comparison between the curvature-based capillary pressure
and the applied bulk pressure differences at various applied
pressures for the drainage case imaged at 5x. .............................................. 87
4.5 Similar to Figure 4.4, but for the imbibition case. The
two data points correspond to the two images shown in
Figure 4.2. .............................................................................................. 87
4.6 Images captured at four different magnifications for the
same equilibrium phase configuration during drainage. While
the FOVs cover the entire porous section at 5x and 8x, the
FOVs at 20x and 40x only cover a portion of the porous
section. For 20x and 40x, it was made sure that the leading
front of the invading phase is always in the FOV......................................... 89
4.7 Interfacial length between wetting and nonwetting phases
for the drainage case as a function of applied pressure at
two different magnifications. ...................................................................... 90
4.8 Applied pressure vs. water saturation. The applied
pressure is calculated based on the water tank height,
whereas the water saturation is calculated from images
captured at 5x and 8x. .............................................................................. 91
xiv
LIST OF FIGURES – CONTINUED
Figure Page
4.9 Comparison between the curvature-based capillary pressure
and the applied bulk pressure differences at various applied
pressures for the drainage case imaged at all four magnifications. ................. 92
4.10 Equilibrium phase distributions during drainage at different
applied boundary pressure in the heterogeneous micro-
model. N-heptane displaces water from right to left. These
particular set of images are captured at 5x magnification.
Green represents n-heptane, blue is water and black is solid grains. .............. 93
4.11 Equilibrium phase distributions during imbibition at differ-
ent applied boundary pressures in the heterogeneous micro-
model. Water displaces n-heptane from left to right. These
particular set of images are captured at 5x magnification.
Green, blue and black represent n-heptane, water and solid
grains, respectively. ................................................................................... 94
4.12 Sample images captured at all for different magnifications
for the same equilibrium state during drainage. Again for 5x
and 8x, the FOVs cover the entire porous section, whereas
FOVs in 20x and 40x only cover portion of the porous section. .................... 95
4.13 Curvature-based macroscale capillary pressure vs. the
applied boundary pressure for the heterogeneous case mea-
sured with four different magnifications. .................................................... 96
4.14 Applied pressure vs. saturation calculated for the hetero-
geneous case with both 5x and 8x magnifications. ....................................... 97
4.15 Connected and disconnected water is shown separately for
the case of applied pressure of 5.6 kPa and 8x magnification.
Disconnected water and connected water regions shown
in pink and green, respectively. It is visible that the
interfacial curvature (blue) associated with disconnected
water appear smaller (i.e., radius is larger) than (yellow)
that associated with the connected water region, indicating
that the pressure within the disconnected water region is
slightly higher than that in the connected water region,
given that the pressure in the nonwetting phase is uniform across. ............... 98
xv
LIST OF FIGURES – CONTINUED
Figure Page
4.16 Histograms of the radii occurring for the drainage case at
5.6 kPa and 8x. (a) all the interfaces (yellow and blue
circles); (b) only the interfaces created between the con-
nected water and n-heptane (yellow circles only); (c) only
the interfaces created between the disconnected water and
n-heptane (blue circles only). ..................................................................... 99
4.17 Corrected curvature-based macroscale capillary pressure
vs. the applied boundary pressure for the heterogeneous
micromodel for the drainage case imaged at 8x. It is seen
that this corrected capillary pressure shows much better
agreement for higher applied pressure during drainage. .............................101
4.18 Imbibition case where the applied pressure is 1355 Pa. This
case is for 5x. Connected and trapped water is shown in
green and purple respectively. Blue circles are associated
with the disconnected water and yellow circles are the ones
between connected water and n-heptane. ...................................................102
5.1 Match between calculated flow rate and applied flow rate in
the homogeneous micromodel for a single phase flow. Errors
bars represent the combined uncertainities while calculating
flow rate from average velocity at the outlet. Maximum
relative difference is 14% that occurs for the 0.4 µL/minute
applied flowrate case. We attribute this difference to the
syringe pump instability. The syringe pump instability has
been demonstrated in Figure 5.2. ..............................................................105
5.2 It is visible that calculated flow rate for 0.6 µL/minute
applied flow rate fluctuates in the higher range. Possible
reasons for this is syringe pump’s inability to stabilize the
flow rate from one value to another, poor calibration and
built-in inaccuracy....................................................................................106
5.3 Steady state velocity field in the homogeneous micromodel.
DI water is flowing from right to left. The applied flow
rate is 0.2 µL/minute. The maximum velocity occurs at
the throats which is 0.5024 mm/s. ...........................................................107
xvi
LIST OF FIGURES – CONTINUED
Figure Page
5.4 Steady state velocity field in the heterogeneous micromodel.
DI water is flowing from right to left. The applied flow
rate is 0.05 µL/minute. The maximum velocity occurs at
the throats which is 0.6333 mm/s. The location of the
maximum velocity is marked with a red star. .............................................108
5.5 The instantaneous flow rate is plotted against time on a
semi-log plot. Time zero corresponds to the first pair where
the PIV processing began and the air phase just enters
the porous medium. This is frame 891 shown in Figure
5.6a. The last frame in this recording is 70000 which is
shown in Figure 5.6b. The whole invasion takes about 760
seconds. During this displacement process, initially, the
air interface velocity is high and maximum velocity (thus
flow rate) occurs at around 14.5 seconds. The sharp peaks
represent the Haines jumps (i.e., pore-scale event when the
nonwetting phase passes a narrow pore and very rapidly fills
the next wider pore space) [5, 6, 14]. The applied pressure
is 8.95 kPa and the flow rate is calculated at the outlet. .............................110
5.6 (a) The first frame in the series where PIV processing
began. (b) last frame in the series. ...........................................................111
5.7 Instantaneous velocity vector field between the frame pairs
902-903. The applied pressure is 8.95 kPa. Maximum
velocity is 0.8324 mm/s. White and gray represent air and
grain respectively. ...................................................................................112
5.8 Instantaneous velocity vector field between the frame pairs
1116-1117. The applied pressure is 8.95 kPa. Maximum
velocity is 0.5222 mm/s. White and gray represent air and
grain respectively. ...................................................................................113
5.9 Instantaneous velocity vector field between the frame pairs
14695-14706. The applied pressure is 8.95 kPa. Maximum
velocity is 0.0384 mm/s. White and gray represent air and
grain respectively. ...................................................................................114
xvii
LIST OF FIGURES – CONTINUED
Figure Page
5.10 The region marked in red is one of the regions of active film
flow. The two images show that air has displaced most of
the regions. The high velocity at the top in outlet (barrier)
means that the flow contribution of all the film flow. ..................................115
5.11 Film flow around the solid posts are tracked via PTV. The
average velocity magnitude is shown by the color map.
For better understanding of the flow pathway, film flow
directions are indicated by the three superimposed lines
around the grains (red, green and black). Also, it is worth
mentioning that only the particles around the grains walls
are tracked and larger water pockets are not included in the
PTV analysis. ..........................................................................................116
5.12 The instantaneous flow-rate for imbibition is plotted against
time on a semi-log plot. Time zero corresponds to the first
pair from where the PIV processing began after the valve
has just been opened. The whole invasion takes about 265
seconds. The applied pressure is 2.55 kPa and the flow rate
is measured at the outlet as usual. ............................................................118
5.13 Instantaneous velocity vector field between the frame pair
2681-2682. The applied pressure is 2.55 kPa. Maximum
velocity is 0.6695 mm/s. White and gray represent air and
grain respectively. Water displaces air from left to right
during imbibition. ...................................................................................119
5.14 Instantaneous velocity vector field between the frame pair
3169-3170. The applied pressure is 2.55 kPa. Maximum
velocity is 0.6691 mm/s. White and gray represent air and
grain respectively. Water displaces air from left to right
during imbibition. ...................................................................................120
5.15 PIV raw images showing film flow and snap-off. Water
films around grains results in swelling of films. The red
rectangle in the two images show one of the regions where
this happens. Time interval between these two frames is
3.982 seconds. .........................................................................................121
xviii
LIST OF FIGURES – CONTINUED
Figure Page
5.16 PTV image averaged over 20 frames. Water films around
grains results in swelling of films. The direction of the films
are shown by the 4 lines with arrows around the grains.
Water is displacing air from left to right. ...................................................122
5.17 The instantaneous flow-rate is plotted against time on a
semi-log plot. Time zero corresponds to the first frame of
the first pair where the PIV processing began and the air
phase just enters the porous medium. This is frame 666
shown on top left. The last framae in this recording is 65000
which is on top right. The whole invasion takes about 715
seconds. During this displacement process, initially the air
interface velocity is high and slows down eventually. The
applied pressure is 9.12 kPa. .....................................................................124
5.18 Instantaneous velocity vector field between the frame pars
727-728. The applied pressure is 9.12 kPa. Maximum
velocity is 0.5500 mm/s. White and gray represent air and
grain respectively. Air displaces water from right to left
during drainage........................................................................................125
5.19 Instantaneous velocity vector field between the frame pars
16315-16515. The applied pressure is 9.12 kPa. Maximum
velocity is 0.0038 mm/s. White and gray represent air and
grain respectively. Air displaces water from right to left
during drainage........................................................................................126
5.20 Film flow connecting the outlet with trapped water down-
stream. The superimposed red and green lines (with
arrows) show the direction of two films and the veloc-
ity magnitude is averaged over 60 frames (12551-12600)
displayed with the color map around grains. The larger
trapped water phase is not included in the PTV analysis,
only the films around the grains are considered. ........................................127
5.21 Instantaneous velocity vector field between the frame pars
5286-5287. The applied pressure is 3.91 kPa. Maximum
velocity is 0.7233 mm/s. White and gray represent air and
grain respectively. Water displaces air from left to right
during imbibition. ....................................................................................128
xix
NOMENCLATURE
Greek Letters
µnw Dynamic viscosity of the nonwetting phase
µw Dynamic viscosity of the wetting phase
∇′ Two-dimensional surface divergence operator
γ Interfacial tension
γwn Interfacial tension between the wetting-nonwetting phase
γws Interfacial tension between the wetting-solid phase
γnws Interfacial tension between the nonwetting-solid phase
θ Static contact angle
ϕ Porosity
µ′ Chemical potential
ρ Density
ρw Density of the wetting phase
ρwn Density of the nonwetting phase
Symbols
J Twice the mean curvature
U Internal Energy
G Gibbs free energy
F Helmholtz free energy
V Volume
N Number of component of a species
T Thermodynamic temperature
P Thermodynamic pressure
S ′ Entropy
A Area
R1 First principal radius of curvature
R2 Second principal radius of curvature
r Pore radius
pc Micro/pore-scale capillary pressure
Pc Macroscale capillary pressure
P c Capillary pressure
Sw Wetting phase saturation
Snw Nonwetting phase saturation
d Micromodel depth
w Micromodel width
g Gravitational acceleration
xx
NOMENCLATURE – CONTINUED
Q Volumetric flow rate
q Darcy flux/Darcy velocity
nw Unit normal vectors at the interface positive outwards from the wetting
phase
nn Unit normal vectors at the interface positive outwards from the nonwetting
phase
vα Fluid velocity inside a porous medium/interstitial velocity
awn Specific interfacial area between the wetting-nonwetting phase
ans Specific interfacial area between the nonwetting-solid phase
aws Specific interfacial area between the wetting-solid
Krα Relative permeability
K Intrinsic permeability
M Viscosity ratio
Ca Capillary number
Bo Bond number
Abbreviations
CFD Computational Fluid Dynamics
CMOS Complementary Metal-oxide Semiconductor
CCS Carbon Capture and Sequestration
DI Water Deionized Water
DPNM Dynamic Pore Network Model
EOR Enhanced Oil Recovery
FOV Field of View
FFT Fast Fourier Transform
ICAL Imaging and Chemical Analysis Laboratory
LED Light-emitting Diode
LBM Lattice Boltzmann Methods
Micro-CT Micro-Computed Tomography
M-LEFT Microfluidics Lab for Energy & Fluid Transport
MMF Montana Microfabrication Facility
NMR Nuclear Magnetic Resonance
NAPLs Non-aqueous Phase Liquids
ppm Parts per Million
PIV Particle Image Velocimetry
PTV Particle Tracking Velocimetry
PEMFC Protonexchange Membrane Fuel Cells
PDMS Polydimethylsiloxane
xxi
NOMENCLATURE – CONTINUED
PMMA Poly(methyl methacrylate)
REV Representative Elementary Volume
SEM Scanning Electron Microscopy
TCAT Thermodynamically Constrained Averaging Theory
xxii
ABSTRACT
Multiphase flow in porous media is ubiquitous in natural and engineering processes. A
better understanding of the underlying pore-scale physics is crucial to effectively guiding,
predicting and improving these applications. Traditional models describe multiphase flows
in porous media based on empirical constitutive relations (e.g., capillary pressure vs.
saturation), which, however, are known to be hysteretic. It has been theoretically shown that
the hysteresis can be mitigated by adding new variables in the functional form. However,
experiments are still needed to validate and further develop the theories. In particular, our
understanding of capillary pressure characterization and numerous pore-scale mechanisms
is still limited. For instance, during capillary pressure measurement, fluid phases become
disconnected, making the bulk pressure an inaccurate measure for the actual capillary
pressure. In a strongly wetting medium, wetting phase always remains connected by corner
films, through which trapped water continues to drain until a capillary equilibrium is reached,
but the effects of corner film flow are minimally characterized.
In this thesis, two different experiments are presented. In the first experiment, we
focused on the capillary pressure characterization and the effect of measurement resolution.
Microscopic capillary pressure along with other geometric measures are characterized during
drainage and imbibition. By strategically varying the pressure at the boundary, different
equilibrium states were achieved and imaged at four different magnifications (i.e., 2, 1.25, 0.5,
0.25µm/pixel). In the second experiment, we for the first time characterized the corner film
flow again during drainage and imbibition condition employing particle image velocimetry.
Overall, our results suggest that the calculated macroscale pressure Pc and the bulk
pressure drop agree reasonably well when only interfaces associated with the connected
phases are considered. A spatial resolution of 2µm/pixel seems to sufficiently resolve the
interface, and further increasing the resolution does not have a significant impact on the
results. Additionally, corner film flow was found to be an active transport mechanism.
During drainage, trapped water is continuously drained over time via thin film, whereas
during imbibition snap-off events are enhanced by wetting films. These observations call for
future studies to carefully treat corner film flows when developing new predictive models.
1
CHAPTER 1 INTRODUCTION
Motivation
Multiphase flow in porous media in general refers to the simultaneous movements of
more than one fluid phase through a solid porous medium [12]. In such a system, a porous
medium is a permeable solid structure with connected pores and passages through which
the moving phases can flow, whereas the moving phases themselves are separated from
each another by identifiable boundaries (i.e., the interfaces), where sharp discontinuities
in composition and physical and chemical properties exist. On a broader perspective,
multiphase flow in porous media can have numerous variations. For instance, the solid
porous matrix can be rigid or deformable and the porous structures can be homogeneous
or heterogeneous. The flowing phases can be fluids (e.g., water) as well as solids (e.g.,
solid suspensions). When more than two fluid phases are present, they can be miscible or
immiscible, and there can even be fluid-fluid and fluid-solid reactions. While the multiphase
flow system can be extremely complex, in this thesis we restrict our discussion to the flow of
two immiscible fluid phases (e.g., water and air) through a rigid and nonreactive solid porous
matrix. Some naturally occurring porous media of this kind include soil, quartz sandstone
and shale gas formation [18]. Infiltration of rainwater into a soil and movement of subsurface
water through a geologic formation are examples of immiscible two-phase flow of this kind.
Multiphase flow in porous media is central to a broad range of natural and engineering
processes, such as enhanced oil recovery (EOR) [9, 13, 34], shale gas production [154], fluid
separation in fuel cells [100], groundwater contamination and remediation [70, 152], and
geological sequestration of carbon dioxide [181] (Figure 1.1). For instance, global warming
has become a serious challenge facing us in the century, primarily due to the continuously
2
Generalized Microfluids Enhanced Oil Recovery Environment
Pattanayak, P. et al. 2021 https://melzerconsulting.com/co2-enhanced-oil-recovery-explained/ Fagerlund, Fritjof. (2022)
Industry Hydrology
Applications of
  Multi-phase Flow 
 in Porous Media
Saman Rashidi et al. 2020 https://www.weather.gov/lot/hydrology_education
Biomedical CO2 Sequestration Energy, Fuel Cell
Mehdi Mosharaf-Dehkordi (2019) Image Credit: VectorMine/Shutterstock.com Jiseon You et al. 2017
Figure 1.1: Application examples of multiphase flow in porous media. Schematics showing
from top left clockwise: Generalized microfluidic devices such as Lab on a Chip (LOC),
reproduced from Ref. [137]; Enhanced Oil recovery in operation such as steam injection;
Contamination from industrial sources into the subsurface, reproduced from Ref. [47];
Hydrology such as the water cycle involves multiphase flow in porous media; Microbial fuel
cell (MFC) involves porous filament thus flow through porous materials, reproduced from
Ref. [171]; CO2 sequestrations from major sources; Blood and bile flow through liver lobules
can be simulated as a flow through porous media, reproduced from Ref. [126]; Heat sinks
with improved performances applying porous media, reproduced from Ref. [146]
3
increasing CO2 level in the atmosphere. The global average level of atmospheric CO2 has
reached from about 280 parts per million (ppm) at the beginning of the industrial age in the
mid-18th century [49] to 414.72 ppm in 2021 [110]. In order to mitigate the green house gas
effect and eventually achieve the 1.5 ◦C goal set by the Intergovernmental Panel on Climate
Change (IPCC), negative emission technologies are critical [81, 87, 136]. Among others,
geological sequestration of CO2 into saline aquifers represents a viable method. In this
so-called carbon capture and sequestration (CCS) technology, CO2 is captured from major
stationary sources, transported by pipelines and injected into suitable deep rock formations
[10, 29, 63, 98]. In this operation, the multiphase flow of resident brine and injected CO2
through porous rock (e.g., sandstone) at the microscopic scale plays a defining role in CO2
storage safety, capacity and security at the macroscopic scale. For example, under reservoir
conditions, the injected CO2 often has a lower viscosity compared with the resident brine
with a viscosity ratio of ∼ 0.01, which leads to instability at the fluid-fluid interface, causing
CO2 to form so-called fingering through resident brine. The fingering phenomenon in turn
causes CO2 to bypass portions of the pore space, resulting in low CO2 trapping and storage
[118]. Therefore, the detailed pore-scale physics underlying the displacement patterns in
porous media is crucial to the prediction of CO2 migration within the porous formations of
geologic storage reservoirs [30, 94].
As another example, EOR is an innovative technology to increase oil production.
Following the conventional primary and secondary recovery processes, EOR is applied by
means of injecting high pressure and high temperature water, chemicals, or gases into the
reservoir to displace residual oil, effectively forming a multiphase flow in porous media.
Our ability to successfully, safely, and economically apply EOR is directly determined by
our understanding of pore-scale physics and flow. For instance, the presence of thin oil
films on rock walls can limit the overall recovery efficiency by causing the non-wetting
phase entrapment as illustrated in Figure 1.2 [43, 124, 142]. Additionally, contamination
4
Figure 1.2: Different oil entrapment at the pore scale. Here the rock grains and oil are in
brown and black, respectively. (a) oil film adhered to pore walls; (b) oil trapped in dead
ends; (c) oil trapped due to heterogeneous porous media; and (d) oil ganglia in pore throats.
Figure was reproduced based on work in Ref. [43].
by hazardous organic chemicals of soil and groundwater poses serious risks to human health
and the environment. In particular, the presence of non-aqueous phase liquids (NAPLs)
(i.e., petroleum hydrocarbon fuels and solvents) in the subsurface serves as a long-term
source for groundwater contamination [2, 3, 41, 138]. In the sector of clean energy, fuel
cell technology exhibits a huge economic and environmental potential in the future power
markets, from small portable cells to large residential power plants. In particular, proton-
exchange membrane fuel cells (PEMFC) consist of porous materials through which fuel and
oxidant are delivered to the reaction sites, where electrochemical reactions take place to
generate power [20, 112]. In this context, reactant transport and water management in
the porous electrodes and membranes are major challenges in fuel cell design. Finally, in
biological systems and tissue engineering, water, minerals, and nutrient transport are often
associated with transport through natural and artificial porous media [38]. It has been shown
5
that tumor growth can be modeled as a multiphase system composed of an extracellular
matrix [153].
All the applications described above directly involve multiphase flow in porous media
and an accurate understanding of fundamental flow physics at the pore scale is directly
linked to our ability to develop engineering systems and predicting larger scale phenomena
[142]. In porous medium studies, two spatial scales are typically identified: the micro/pore
scale (e.g., discrete pores), at which fundamental flow and transport take place, and
the macro/continuum scale, at which predictions are desired for practical applications.
Macroscopic phenomena are however governed collectively by various pore-scale processes.
An accurate continuum description of these processes always requires a rigorous description
of the pore-scale physics, which is often incorporated through constitutive relations. It is
increasingly accepted that developing subgrid-scale models capable of accurately representing
pore-scale processes is critical to improving the accuracy of continuum-scale predictive
tools. Even though it is not always practical to translate pore-scale results directly into
reliable large-scale predictions, the understanding and knowledge learned from pore-scale
investigations do improve our ability to advance such efforts. Enabled by improved imaging
ability, innovative fabrication techniques and better numerical simulation capabilities, our
understanding of pore-scale multiphase flow has been greatly improved over the past few
decades [18]. Techniques such as x-ray computed micro-tomography (micro-CT) and nuclear
magnetic resonance (NMR) have made it possible to visualize and quantify microscale
characteristics associated with interfacial properties and fluid distributions of relevance
to multi-phase fluid systems [33, 106, 155, 156, 167]. Consequently, image analysis of
the resulting three-dimensional microtomographic image data has become an essential
component of pore-scale investigations [106]. Additionally, microfluidic model porous media
or so-called micromodels have become a very powerful tool over the past two decades and
allowed direct real-time observation of pore-scale processes with very high temporal and
6
spatial resolutions [24, 91, 93, 94, 103, 121, 147, 147, 175, 176]. By building surrogate
porous media using transparent materials such as glass, these micromodels overcome
the visualization challenge in the real rock experiments. With the aid of micromodels,
several important pore-scale mechanisms were discovered and investigated, including Haines
jumps and snap-off [5, 14, 54, 70, 157]. Meanwhile, numerical modeling such as direct
numerical simulation and pore network modeling have enabled us to resolve very small flow
patterns [85, 167, 169]. In particular, computational fluid dynamics (CFD) and Lattice-
Boltzmann Method (LBM) have gained increasing popularity with an advent of microscale
characterization methods and advanced computational resources [169].
Despite the exciting research advances, a complete description of multiphase flow in
porous media is still not available, due to the inherently complex pore geometries, fluid-fluid
interfaces and the associated pore-scale physics [139]. For a long time, flow in porous media
has been modeled by the empirical Darcy’s equation which was formulated based on averaged
properties such as bulk pressure drop and permeability,
−kAc
Q = ∆p (1.1)
µL
where Q is the volumetric flow rate, Ac, L and k are the cross-section area, length and
permeability of the porous medium, respectively, and ∆p is the total pressure drop across
the porous medium. For multiphase flows, Darcy’s law has been used in various extended
forms, for which constitutive relations such as capillary pressure (Pc) vs. saturation (S) are
used to close the equation system [12, 44]. While the Pc − S relation curves are easy to
determine experimentally or numerically for various conditions [56], they are, however, well
known to be hysteretic and non-unique. That is one saturation value can be corresponding
to many different values of capillary pressure, which in turn depend on flow history. The
validity of this modeling approach has been questioned over the years and there have been
7
increasing efforts to modify and improve such models. Hassanizadeh et al. employed the
thermodynamically constrained averaging theory (TCAT) and rational thermodynamics
approach and proposed that the apparent hysteresis in Pc − S relations is a result of
incomplete parametrization [61, 66, 67, 68, 134].
Therefore, a long-standing problem in two-phase flow in porous media has been to find
out the complete set of Darcy-scale state variables to model macroscopic fluid displacement in
a theoretically sound fashion [134, 150, 152]. Following the theoretical work of Hassanizadeh
and Gray, numerous experimental, numerial and theoretical efforts have been dedicated to
the study of capillary hysteresis. In addition to the fluid phase saturation (S), these studies
have been trying to correlate macroscale capillary pressure with other variables such as fluid-
fluid interfacial areas and Euler characteristic, which not only account for the amount of each
fluid phase in the porous media, but also the distribution and arrangement of each phase.
For example, Hassanizadeh and Gray first hypothesized that fluid-fluid interfacial area is
the missing state variable since capillary pressure is an interfacial property and therefore
specific interfacial area (i.e, interfacial area between wetting and nonwetting phases per unit
volume) should be included in the macroscopic description. They also posited that capillary
pressure as a function of saturation and specific interfacial area should fall on a unique
surface [57, 58, 59, 66, 67]. This hypothesis was tested and it was observed that hysteresis
was not completely removed for either static or transient conditions [84, 92]. Attempts have
also been made to include integral geometrical measure as a state variable such as Euler
characteristic based on the understanding that this variable provides information on the
distribution of fluid phases especially non-wetting phase. Euler characteristic is considered
as an appropriate candidate to be included in the non-hysteretic functional form [116, 120].
Moreover there has been more general approaches to find a link between capillary pressure
and disconnected phases which again suggests that some variable that distinguishes between
connected and disconnected phases in a two-phase flow system should be included in the
8
flow description [72, 73].

Figure 1.3: (a) two-phase flow with trapped water (wetting phase), here blue is connected
water, red is trapped water, black is solid grain and green is connected n-heptane (nonwetting
phase); (b) zoomed in view of the interfaces, connected interfaces in yellow and disconnected
ones in purple. (c) disconnected n-heptane as yellow and blue and black represents connected
water and solid grain respectively as before.
To obtain a Pc − S curve in the traditional way, a core sample of the size of a few
millimeters to a few centimeter is often used. The core sample is initially saturated with a
9
wetting phase (i.e., water) at a low pressure, following which a nonwetting phase is injected to
the core sample by executing a step change in the boundary pressure values and monitoring
the resultant saturation change over time until an equilibrium value is achieved. Often the
saturation values are averaged over the domain of the system, and the bulk measured pressure
values (i.e., the pressure drop between the inlet and outlet of the core sample measured with
a differential pressure transducer, traditionally called the macroscopic capillary pressure) are
assumed to represent the pressure of all of a given fluid phase present within the system.
However, enabled by the aforementioned theoretical advancement and recent progresses in
advanced imaging techniques, such as micro-CT, NMR and microfluidics [7, 71], a few issues
have been pointed out by recent studies [7, 7, 56, 71, 120].
The first issue concerns the bulk pressure. Essentially, it is not clearly how well
the bulk pressure measured between the inlet and outlet boundaries represents the real
pressure difference between different phases. In fact, thanks to the advanced image
analysis techniques, it has been confirmed that often times, fluid phases, especially, the
nonwetting phase can be disconnected from the boundaries a shown in Figure 1.3. When we
consider a strongly wetting system, it is often acceptable to assume that a singular pressure
measurement is representative of the pressure state of all the wetting fluid because the
wetting phase is connected via thin films through which, the pressure may be transmitted
throughout the domain. However, for nonwetting phase, the assumption cannot be applied
by default, because during flow, nonwetting fluid becomes disconnected and isolated within
pore bodies of the solid matrix, so pressure changes in the bulk nonwetting phase are not
transmitted to the isolated nonwetting clusters [71]. To alleviate this issue, a few researchers
attempted to obtain the macroscopic pressure by directly averaging the microscopic capillary
pressure, which is in turn calculated from the interfacial curvatures within individual pores
based on Young-Laplace equation,
10
1 1
pc = γ( + ) (1.2)
R1 R2
where pc is the microscopic capillary pressure across a interface, γ is the interfacial tension
between the nonwetting and wetting phases, and R1 and R2 are the two principal radii of
curvature of a surface. This approach is believed to yield a more accurate measurement
of capillary pressure as the connected regions and disconnected regions can be evaluated
separately. For example, Armstrong et al. calculated the curvature-based capillary pressure
based on micro-CT images, and achieved reasonable agreement with the transducer-based
measurement, but only when only the connected fluid phases are considered [7]. In a separate
study, Herring et al. demonstrated similar results [71].
However, as the curvature-based approach becomes a potential solution, a new issue
has arisen: most micro-CT measurement suffers from low spatial resolution (and often low
temporal resolution) resulting in inaccuracy for low interfacial area and high curvature
measurements, which in turn affects curvature-based capillary pressure measurements [71].
For instance, Herring et al. [71] reported that variation in curvatures seen at 2x was not
seen at 4x for their glass bead samples. Yet, for their sandstone rock samples, there was
variation in curvatures for disconnected phase at both resolutions [106]. In the measurement
by Armstrong et al., although reasonable agreement was achieved between curvature-based
and transducer-based capillary pressure measurements when disconnected fluid interfaces
are excluded, there is still a discrepancy as high as 10%, And it is not yet clear whether
this difference is due to some unidentified physics or simply due to inaccurate curvature
measurement which in turn is caused by insufficient spatial resolution. Given that typical
micro-CT is limited at a spatial resolution of 5µm/voxel, the question we must answer is
that: Is the current spatial resolution sufficient to fully resolve the interface curvature? If
not, what is the required resolution to achieve so for a given pore size distribution? Can we
11
take advantage of optical based methods, which can offer much higher spatial resolution to
perform a resolution dependence test, similar to the mesh dependence test that is routinely
performed in numerical simulation.
The second issue concerns the establishment of equilibrium state during capillary
pressure. Essentially, it is unclear how much waiting time is needed to achieve true
equilibrium following a step change of inlet pressure. Some studies waited as long as one
day to collect one data point, whereas some other studies waited as short as 15 min, both
of which were selected on a random basis [7, 12]. So far, very little experimental evidence
has been reported to address the relaxation and/or rearrangement of interfaces during the
transition from dynamic to static conditions. However, it has been pointed out that the
relaxation time strongly depends on the film flow that runs along the wetting corners and
crevices of the porous medium. As illustrated in Figure 1.4, as the inlet pressure is increased,
part of the wetting phase (e.g., water) is displaced downstream, whereas the other part of
the wetting phase is left behind forming residual water ganglion. Although these water
ganglion appear to be trapped and disconnected from the outlet in the bulk, they are indeed
connected through thin corner film, through which the water ganglion can be contiguously
drained, until a true equilibrium is achieved. Therefore, to further our understanding of the
relaxation time and capillary pressure, a careful quantification of thin film flow is critical.
In addition to its relevance to capillary pressure measurement, studying the relative
importance of pore-scale flow mechanisms such as film flow is crucial to furthering our
fundamental understanding of multiphase flow in porous media and allowing us to develop
more accurate predictive models [128]. It is well accepted that film flow phenomena can
enhance the connectivity of porous networks, thus introducing new pathways for fluid
transport. It is thus important to consider the relative contribution of such effects when
designing a simulation protocol for porous media systems. Ignoring such effects can lead to
large errors in the estimation of quantities such as the phase saturation and other related
12

  
 
 
Figure 1.4: Images illustrating scenarios where film flow plays an role: (a) 3D phase
distribution obtained by numerical simulation with brine and super-critical CO2 under strong
imbibition (wetting phase contact angle close to zero). Precursor wetting film ahead of the
main meniscus dominates the flow (Figure was reproduced based on work in Ref. [78]); (b) µ-
PIV images of the film flow in a homogeneous micromodel. Water is displacing air from left
to right, where water and air appear bright and dark respectively. (c) Primary drainage in
a heterogeneous micromodel at a fixed boundary pressure. It is visible that initially trapped
water drains away over time via corner film flow that is connected to the outlet.
13
measures that depend on it (e.g., relative permeability curves). Since real porous media is
never smooth, presence of roughness significantly gives rise to film and corner flow as well
[35, 160]. Presence of wetting films also contributes to the total wetting-solid interfacial area
(i.e, interfacial area between the solid and wetting phases) present in a two-phase system
and must be accounted for accurately in the description. It is usually assumed for a perfectly
wetting medium that the wetting phase completely coats the solid grains [67, 83]. Despite
the important of film flow, very little studies have focused on its quantification. Although
a few studies have been reported in the context of drying, evaporation and imbibition (i.e.,
displacement of a nonwetting phase by a wetting phase in porous media) in simplified pore
geometry [142, 168]. However, there has been very little, if any, studies specifically focusing
on the quantification of film flow during drainage processes (i.e., displacement of a wetting
phase by a nonwetting phase in porous media).
To that end, two different experimental studies are presented in the thesis, both
performed in 2D porous micromodel employing epi-fluorescence microscopy. The first
experiment explores the spatial resolution dependency of the capillary pressure derived
from averaged Young-Laplace equation (i.e., equation 1.2) and how that relates to the bulk
pressure difference (boundary pressure) in a quasistatic drainage experiment. The connected
and disconnected phases are separated employing advanced image processing techniques, and
the resulting capillary pressure is compared against the bulk pressure. The experiments are
done for micromodels with both homogeneous and heterogeneous structures. While the
porous structure in the homogeneous micromodel is composed of regularly arranged circular
pillars, then heterogeneous porous section represents realistic structures of sandstone.
In the second experiment, film flow during drainage and imbibition are studied in
homogeneous and heterogeneous micromodels using microscopic particle image velocimetry
(micro-PIV) and particle tracking velocimetry (PTV). Pore velocity profiles and film flow are
quantified using micro-PIV and PTV, respectively. And the effect of film flow on final fluid-
14
fluid distribution are quantified during both drainage and imbibition with different applied
pressure boundary conditions. These studies will potentially enrich our understanding of
pore scale process of two-phase flow in porous media.
Objectives
The overall objectives of this thesis are to:
• Characterize the pore-scale capillary pressure and compare the boundary pressure
termed as capillary pressure in traditional model with the macroscale capillary pressure
derived from Young-Laplace equation and studying the effect of spatial resolution on
the pixel-based curvature calculation in 2D porous micromodels.
• Characterize film flow for both drainage and imbibition conditions employing micro-
PIV and PTV techniques in both homogeneous and heterogeneous micromodels.
Outline
The rest of the thesis is organized as follows:
• Chapter 2 presents necessary concepts, theories and definitions relevant to multiphase
flow in porous media followed by a short literature review on two-phase flow in porous
media, particularly on characterization of capillary pressure and film/corner flow.
• Chapter 3 presents a detailed description of the experimental methodology for both the
capillary pressure measurement and film flow experiment. The micromodel fabrication
procedures are discussed together followed by two separate subsections which discusses
flow apparatus, experimental procedure, image acquisition, image processing and
calibration procedures.
15
• Chapter 4 presents the results of capillary pressure measurements. Results include
segmented images for different magnifications (5x, 8x, 20x, 40x), and comparison
between the calculated macroscale capillary pressure and the applied pressure at
different magnifications for both homogeneous and heterogeneous micromodels. This
chapter ends with a discussion on the results and their significance.
• Chapter 5 first presents the micro-PIV results at different boundary condition. Single-
phase flow measurement is shown as a validation of our method. Instantaneous
flow rates are measured for the drainage and imbibition experiments for both the
homogeneous and heterogeneous micromodels applying different boundary pressure
and compared. Vector flow fields and saturation are calculated and presented for all
the cases. Film flow during both drainage and imbibition are calculated using PTV at
different pores of active film flow. Their velocity magnitude, relative effects in relation
to the final equilibrium distribution is discussed. Several pore-scale mechanisms
including recirculating flows and liquid bridges are also identified and discussed.
• Chapter 6 concludes the thesis with major conclusions and recommendations for future
work.
16
CHAPTER 2 BACKGROUND AND LITERATURE REVIEW
Definitions and Theoretical Background
Scaling and Representative Elementary Volume
Concept of scaling is very important in the description of fluid dynamics in porous
media. Over the decades a great challenge in the modeling of flow in porous media has
been how to establish a theoretically sound linking across the scales (e.g., sub-pore to pore-
scale, and nanoseconds to days/years). Generally, processes from smaller scales must be
properly accounted for in larger-scale models. Upscaling theories exists for porous media
to achieve that goal [60, 64]. Different spatial and temporal scales exist in porous media
to describe the flow phenomena. Different spatial scales are the molecular scale, the pore
scale, the representative elementary volume (REV) scale and the field scale as shown in
figure 2.1 [12]. The molecular scale emerges from the inter-molecular processes, such as
the motion and collisions of molecules in the fluid phases and of special importance for the
determination of the continuum-scale chemical, physical and transport properties such as
mass density, Young’s modulus, diffusion coefficient, viscosity etc., which are continuously
defined in space or macroscopic observables such as pressure and temperature.
On the pore scale the phase distributions and the phase interfaces can be spatially
resolved with 3D measurement technology like x-ray tomography or neutron tomography
[155, 156]. With fast advancements of these technologies in the last decades allowing more
detailed scans with a higher spatial and temporal resolution, they are providing important
insights about the pore space [167]. Besides artificially generated geometries, these 3D scans
can be used as pore geometry information for pore-scale models in order to compute, for
example, the flow field and transport phenomena in the pore space. In real life, however,
most applications in a larger extent require models with quantities averaged over the REV.
As defined by Blunt [18], in any description of the pore space (either from an image or
17
Figure 2.1: Different scales in porous media with relevant applications. Figure was
reproduced based on work in Ref. [76].
a statistical reconstruction) it is implicit that the averaged properties (such as porosity,
permeability etc.) measured or computed are representative of any similar-sized or larger
region of the same porous medium. This averaging volume is the REV [18].
The size of the REV must be within a range where it would still result the same value
for the averaged quantity even if the averaging volume is increased or decreased as shown in
figure 2.2 [12]. Lower bound for the averaging is given by the microscopic variability, which
can cause strong fluctuations of the averaged quantities, if the averaging volume chosen is
too small. On the other hand the upper bound is mostly relevant for heterogeneous porous
media where an increase in the REV size would result in a low resolution of the medium
and, therefore, the macroscopic heterogeneity would strongly influence the resulting averaged
quantities [127]. In the REV scale, averaged or effective quantities are employed, such as
18
15.0
14.5
14.0
13.5
13.0
12.5
0.0 2.5 5.0 7.5 10.0 12.5 15.0
[mm]
Figure 2.2: Volume-averaged porosity averaged over cubes of different length l, for the
statistical realization of Fontainebleau sandstone. This shows the effect of REV. As the
averaging volume is increased the porosity value converges. Figure was reproduced based on
work in Ref. [18].
porosity, tortuosity, permeability, Darcy velocity etc., which come usually along with an
information loss about the distribution within the averaging volume [127].
REV-scale models with varying complexity are widely used for engineering applications
in porous-medium modeling, such as the simulation of CO2 sequestration, oil production,
subsurface remediation, and, groundwater modeling to name a few [135, 177]. REV scale can
again be classified as laboratory scale and field scale. Most engineering applications require a
continuum approach on a scale in the order of tens of meters to kilometers, which is referred
to as field scale or continuum scale (i.e., at a scale much larger than a characteristic pore-
size). The choice of scale depends on the spatial variability of the system and the nature of
φi( ) [%]
19
the phenomenon being studied [127].
Porosity
Porosity (ϕ) is an important averaged property of any porous medium. It is the ratio
of the void space to the total volume of the medium in the averaging domain. Connected
porosity or effective porosity quantifies the total connected pore space. The averaging volume
must be greater than the REV as mentioned before [12, 133]. Mathematically, porosity is
defined as
void volume
ϕ = (2.1)
total bulk volume
Or, more explicitly in terms of the void volume distribution function [166],
∫
1
ϕ = α(r)dV (2.2)
V V
where the void volume distribution function α(r) is defined as,
1 if r lies in the fluid region
α(r) =  (2.3)
0 if r lies in the solid region
Here, V is the averaging volume which contains the fluid volume and the solid volume.
Generally, porosity is treated as a constant and not a function of space and time in the type
of flow considered here.
Saturation
Saturation (S) or volume fraction of an individual phase is another macroscopic variable
of interest in any description of flow in porous media. In a two-phase system there are wetting
and the non-wetting phase along with solid phase. As an example, the saturation (volume
20
fraction, used interchangeably) of the wetting phase is given by,
∫
1 Vw
Sw = dV = (2.4)
Vf V V
w f
where Vw is the total wetting phase volume within the averaging domain and Vf is the void
volume (i.e., the total fluid volume, and that is the total volume of the domain V multiplied
by the porosity ϕ). According to the definition in Equation 2.4, the sum of all the fluid phase
saturation is unity.
∑
Sα = 1 (2.5)
α=w,nw
where w and nw denote wetting and nonwetting phases, respectively, for a two-phase flow
system.
Phase and Component:
A phase is one of the several different homogeneous materials present in the portion
of matter under study. A phase may be either a pure substance in a particular state or a
solution in a particular state (i.e., solid, liquid, or gaseous). Also, the portion of matter
under consideration may consist of several phases of the same substance or several phases of
different substances. In a more precise thermodynamic sense, phase refers to a homogeneous
region of matter in which there is no spatial variation in average density, energy, composition,
or other macroscopic properties. Phases can also be distinct in their molecular structure.
For example, water has multiple ice phases that differ in their crystallographic structure.
A phase can be considered a distinct “system” with boundaries that interface with either
container walls or other phases. On the other hand, a component is a chemically independent
constituent of a system. The number of components in a system is the minimum number
of independent species necessary to define the composition of all the phases present in the
21
system [8, 25, 162].
Interfacial Tension and the Young-Laplace Equation
Interfacial tension (γ) or surface tension originates from the molecular interactions
between fluid-fluid or fluid-solid interfaces. For an interface between two different fluids,
γ is often referred to as the interfacial tension, whereas for an interface between a liquid
and its own vapor, γ is often called the surface tension. In many instances, these terms are
used less precisely and interchangeably [139]. According to Gibbs treatment of interfaces, in
classical thermodynamics the energy (U) representation of the state of an ideal system with n
components is U = U(S ′, V,N1, N2, . . . , Nn), with S ′, V and N denoting entropy, volume and
the number of species of each component, respectively. The system can also be represented
in other potentials such as Gibbs free energy G = G(T, P,N1, N2, . . . , Nn) in terms of
temperature (T ) and pressure (P ) or the Helmholtz free energy F = F (T, V,N1, N2, . . . , Nn)
in terms of temperature (T ) and volume of the system (V ). Fundamental thermodynamic
equation in Gibbs’ coordinate for energy,
∑n
dU = TdS ′ − PdV + µidNi (2.6)
i=1
For a non-simple system inclusion of non-moving boundary work (such as electric
∑charge transport or surface deformation work) gives, dU = TdS ′ − PdV + Fdx +
n ′
i=1 µi∑dNi.Similarly, the total differential of other potentials gives dG = −S ′dT∑+ V dP +
Fdx + n ′
i=1 µidNi for the Gibbs free energy and dF = −S ′dT − PdV + Fdx + n ′
i=1 µidNi
for the Helmholtz free energy where F is any generalized force and x is the conjugate
variable [23, 162]. According to Gibbs’ convention the two phases α and β are thought
to be separated by an infinitesimal thin boundary layer, the Gibbs dividing plane which
has zero volume but surface area. All the extensive variables sum over the three different
phases: α phase, β phase and interface phase σ. The change of total Helmholtz free energy:
22


 
 

	

Figure 2.3: Two phase system with an ideal interface.
dF = −S ′ ∑ ∑
dT−PαdV−(P β−Pα)dV β+ n µ′αdNα+ n µ′ ∑
β
i=1 i i i=1 i dNβ+ n ′σ σ
i i=1 µi dNi +γdA. At
equilibrium, with constant volume, temperature, and constant amounts of material, the free
energy is minimal. At a minimum the derivatives with respect to all independent variables
must be zero. After simplification:
∑n
dF = −(P β − Pα)dV β + γdA + µidNi
i=1
The total differential for the Helmholtz free energy can also be written as the derivatives of
the natural variables and equating terms gives:
( )
∂F
γ =
∂A T,V,V β ,Ni
Which says that interfacial tension between the phases α and β is the change of Helmholtz
23
free energy of the system with respect to the change of surface area of the interface while
keeping the temperature, the total volume of the system, the volume of phase β and the total
β
numbers of all components constant. At equilibrium:∂F = ∂F ∂A + ∂F ∂V + ∂F ∂Ni = 0.
∂A ∂A ∂A ∂V β ∂A ∂Ni ∂A
Simplifying yields (last term being zero and substituting for γ):
β
γ − (P β − Pα ∂V
) = 0
∂A
β
From differential geometry it is known that ∂V = ( 1 + 1 )−1 where R1 and R
∂A R R 2 are
1 2
the two principal radii of curvatures of a surface in 3D which is also equal to twice the mean
curvature [23]. Which leads to:
P β − Pα 1 1
= γ( + )
R1 R2
The pressure difference between the two phases at equilibrium is the capillary pressure which
is the product of the interfacial tension and the two principal radii of curvatures. This is
the Young Laplace equation for capillary pressure which can also be derived from energy
balance of a static interface. Interfacial tension can also be expressed in terms of Gibbs free
energy: ( )
∂G
γ =
∂A T,P,Ni
Where the natural variables are temperature, pressure and number of all components [23].
Contact Angle and Young’s Equation
In a system with two fluid phases, the solid surface generally exhibits a stronger affinity
for one of the fluids, which is termed the “wetting” (W) phase, whereas the other phase
is termed the “nonwetting” (NW) phases. The strength of the wetting properties of the
solid can be quantified in terms of contact angle (θ) as illustrated in Figures 2.4 and 2.5.
24
γ
wetting non-wetting
θ
γws γnws
Figure 2.4: Two fluid phases on a solid surface with different wettabilities. Figure was
reproduced based on work in Ref. [18].
By applying a force balance in the horizontal direction along the solid surface, the Young’s
equation can be recovered,
γnws = γws + γcosθ (2.7)
where γnws, γws and γ are interfacial tension between the nonwetting-solid, wetting-solid and
nonwetting-wetting phases, respectively. It is worth noting that the aforementioned contact
angle assumes a smooth surface and a static condition (i.e., the contact line is stationary).
Contact angle under this static condition is the intrinsic/static contact angle which differs
from advancing and receding contact angles occurring during dynamic events [18].
Capillary Pressure
Capillary pressure is a direct outcome of having more than one fluid phase and key
variable in the description of immiscible two-phase flow in porous media. A detailed
description of capillary pressure from the first principle is long and the mathematics involved
is cumbersome. Readers are referred to the paper by Boruvka et. al [19] which presents
fundamental equations of bulk phases, dividing surfaces, lines and points for a multiphase
multicomponent system and progressively derives simplified forms of those equations for
25
θ = 0 ° θ =45 ° θ = 90 °
Completely water-wet Water-wet Neutrally wet
oil oil
water oil
θ θ
solid surface
θ = 135 ° θ = 180 °
water Oil-wet Completely oil-wet
Oil
θ oil
solid surface
Figure 2.5: Water and oil in contact with a solid, contact angle is measured from the denser
fluid. Figure was reproduced based on work in Ref. [18].
specific cases [19, 23]. Here we only present some derived results.
Capillary pressure at the microscale can be defined at every point on the interface by,
pc := −γnwJw = γnwJn (2.8)
where, γnw is the interfacial tension between the wetting-nonwetting phase and Jw and Jn
are twice the mean curvature of the interface (sum of the two principal curvatures) measured
from the wetting and nonwetting side respectively.
Jw = −Jn = ∇′.nw = −∇′.nn (2.9)
26
where ∇′ is the two-dimensional surface divergence operator and nw and nn are the unit
normal vectors at the interface positive outwards from the wetting and nonwetting phases
respectively. The definitions of unit normal and surface divergence for any general surface in
three-dimensional (3D) space can be found in the literature [19]. This definition of capillary
pressure applies to both equilibrium and dynamic conditions and doesn’t require solid phase
and contact angle to be defined. This is basically a normal stress balance (if ∇′.nw or
∇′.nn ≠ 0, there must be a normal stress jump at the interface, which generally involves both
pressure and viscous terms). Now, only at equilibrium this Pc equals the pressure difference
between the two phases on the two sides of the interfaces (pn and pw for nonwetting and
wetting phases, respectively) [19, 56],
pc = pn − pw = −γnwJw = γnwJn (2.10)
So far, all these definitions are for micro/pore scale. The macroscale counterpart comes
from averaging microscale Pc over the fluid-fluid interfaces within an averaging volume V
that contains both the fluid phases, the solid phase and the interfaces,
∫
1
Pc = Pc dA; where Ω ∈ V (2.11)
Awn
Ω
where Ω is the interfacial area domain within the total volume of V [56]. The last equation
applies to both equilibrium and dynamic condition. If we know the instantaneous microscale
pc for every interface over a domain, we can find macroscale Pc for the whole domain for
that moment.
However, traditionally macroscale capillary pressure is defined as, Pc = Pn −Pw ,where
Pn and Pw are the pressures at the inlet connected with a nonwetting phase reservoir and
outlet connected with a wetting phase reservoir, respectively, in an experimental cell. For a
measurement in a 2D porous medium with a constant depth (d), a simplified 2D versions of
27
the Young-Laplace equation can be written as [5, 90],
1 2cosθ
P c = γnw( + ) (2.12)
R1 d
where R1 is the principal radius of curvature in the horizontal direction, d is the micromodel
depth and θ is the contact angle.
Drainage and Imbibition
Drainage is the displacement process where a nonwetting phase displaces the wetting
phase. Drainage process is restricted by the smallest passage/constriction in a porous
medium. If we consider drainage in a capillary tube with a circular constriction, non wetting
phase can enter a pore throat only if the applied pressure at the inlet is greater then the
capillary pressure of that pore throat (Pc = 2γcosθ/r, where r is the radius of circular
throat). While a pore throat is defined as the narrowest point of the constriction, pore body
is the widest region in a neighborhood connected with one or multiple pore throats, and it
allows entry of the nonwetting phase spontaneously once it passes a throat connecting the
pore body during drainage.
Imbibition is essentially the opposite of drainage process, where the wetting phase
invades the porous medium to displace the nonwetting phase. During imbibition, whether
or not an interface can advance forward is determined primarily by the size of the pore
body (i.e., the widest point of the porous structure), whereas in contrary, whether or not
an interface can advance forward during drainage is determined primarily by the size of the
throat (i.e., the narrowest point). Imbibition is impeded by pore body or the wide regions
since wetting phase tends to reside in the narrower region/throats throat filling is favorable in
imbibition. Drainage and imbibition are central to porous media flow and their mechanistic
description is different [18].
As mentioned above, capillary pressure saturation relation (Pc − S) use in traditional
28
models is hysteretic (Figure 2.6). Depending on the displacement history, the same saturation
can be corresponding to different capillary pressures. During drainage, the nonwetting phase
first invades the porous media, which is originally saturated with the wetting phase. As
the applied pressure (i.e., macroscopic capillary pressure) increases, the wetting saturation
decreases up to a point, after which further increasing pressure does not displacement more
wetting phase (so-called the irreducible wetting phase saturation). This first drainage process
is called the primary drainage. Following that, upon reducing the applied pressure, wetting
phase returns to the porous medium by displacing nonwetting phase, going through an
imbibition process. This is the main imbibition. As shown in Figure 2.6, the two curves are
quite different due to hysteresis. This processes can be repeated many more times, and more
scanning loops can be obtained all with different shapes [18, 129].
Darcy’s equation
From the single-phase flow Darcy’s equation (Equation 1.1), the extended two-phase
Darcy description can be written by introducing additional variables [37, 129],
−KrαK
qα = (∇Pα − ραg); α = wetting or nonwetting (2.13)
µα
where qα = Q/A is the Darcy flux/velocity (m3/s/m2), which is not the actual velocity of
the fluid, but instead the fluid velocity vα is related to the Darcy flux qα = Q/A by the
porosity ϕ.
qα
vα = (2.14)
ϕ
When Darcy’s law is extended from single-phase flow to multiphase flow, the actual
permeability of a porous medium is replaced by the effective permeability of the respective
phase. Relative permeability is a dimensionless measure of the conductance of a porous
29
primary
drainage
secondary
drainage
main
imbibition
s w
ir sw
0.0 1.0
Figure 2.6: hysteretic capillary pressure, P c as a function of wetting phase saturation. Figure
was reproduced based on work in Ref. [56].
medium for the fluid phase, when the medium is saturated with more than one fluid.
effective permeability Kα
Relative permeability, Krα = (2.15)
intrinsic permeability K
Intrinsic permeability K is a 2nd order tensor (matrix) in 3D coordinates for any general
anisotropic porous medium. For an isotropic medium a single scalar value can represent the
intrinsic permeability. Darcy’s equation along with continuity and a constitutive relation
(Pc − S) serves the traditional two-phase flow modeling in porous media. A large set of
experimental data and parameters is available for a broad range of materials and applications
for the constitutive relation. Although Darcy approach has shortcomings, it is still the most
used approach today and implemented in most numerical modeling software [134].
pn — pw
30
Literature Review
Capillary Pressure Characterization in Porous Media
There have been numerous studies over the past 150 years regarding two-phase flow in
porous media in all the areas: theory, experiment, and numerical computation. Henry Darcy
was the first to propose his empirical linear relation for single phase flow in porous media
derived through experimentation [37, 79, 129]. Darcy’s empirical law, Equation 1.1, has since
been the foundation for mathematical analyses of porous media flow [129]. Although physical
description of multiphase flow in porous media ideally should be based on conservation
principles, in practice, however, Darcy’s law is employed as the foundation of multiphase
flow studies. Darcy’s law is an empirical substitute for the momentum conservation obtained
from the data of experimental study of one-dimensional single-phase flow [129]. Extension
to multiple phases from the Darcy’s law was done by Muskat et al.[130, 131]. In 1940,
Leverett pointed out that in order to include capillary pressure effects in the flow equation,
the pressure must be phase dependent [105].
There have been many modifications of Darcy’s equation: Darcy-Brinkman for non-
negligible velocity gradients correction [21, 99], extensions for larger Reynolds numbers
like the Forchheimer equation [65] and many more to account for other effects (e.g.,
electromagnetic effects [80]. Use of the Darcy’s law (any of its versions such as 3D vector
form) has been analyzed, debated, and scrutinized in numerous papers over the last 50 years.
It has always been well known about the validity of Darcy’s law and its limitations since its
emergence [129]. For instance, it has been established experimentally that in high-velocity
groundwater flows, deviations from Darcy’s law occur [15, 16, 40, 46, 48, 62, 101, 113]. Bear
experimentally concluded that Darcy’s law holds correct only for Reynolds number, based
on grain size, below the range 1-10 [12].
To model multiphase flow, extended Darcy’s law introduces the saturation-dependent
31
parameter such as relative permeability. Since the systems of equations for a two-
phase flow has 10 equations (2 mass balance equations, 3 vector velocity equation for
each phase, equation of state for density for both fluid phases) but with 11 variables
(ρw, ρnw, vw, vnw, pw, pnw, sw) for two phase isothermal flow, the introduction of an
additional constitutive relation becomes necessary for the closure of the system. Capillary
pressure as a function of saturation serves this purpose [56, 119]. This is called the traditional
model for two-phase flow which has been used extensively in hydrology, petroleum industry
and many other areas and still being used [53, 140, 158, 159, 161, 174]. However, there are
a number of well-known problems associated with this model [56, 67, 68]. First, constitutive
relationships (capillary pressure vs. saturation, and also relative permeability–saturation
relations) are found to be hysteretic, suggesting that the set of independent parameters
might be incomplete [97]. Also, in the traditional model capillary pressure is not correctly
defined [56, 120]. Bulk pressure difference (pnw − pw) is termed as the capillary pressure and
equilibration is assumed without justification. Since P c − Sw curve derived experimentally
should be based on equilibrium points only, any dynamic effects would impose a lot of
uncertainty on the plots [56]. Usually, to derive the P c −Sw plot experimentally, a drainage
or imbibition curve for a specific porous medium is fitted to a set of data points [22, 165].
Also, since any pore-scale event is completely missing in the simple P c − Sw relation, any
phase configuration or connectivity description is therefore lost. It is also known that when
nonwetting phase becomes disconnected in a two-phase displacement process, the boundary
pressures of the system are not a reliable measure of the averaged pressures in the system
for each of the fluid phases [56].
Because of all these aforementioned deficiencies in the traditional model, for the last
quarter century researchers have been attempting to find a new, a more theoretically sound,
properly upscaled two-phase flow modeling approach [120]. Hassanizadeh and Gray are
among the first ones to develop a new modeling approach for two-phase flow at the macroscale
32
from rational thermodynamics and volume averaging [66, 67, 68]. For example, Hassanizadeh
and Gray proposed an extended continuum theory for two-phase flow considering all the
phases and interfaces and proposed mass, momentum and energy equations for all the phases
and interfaces and derived a revised constitutive relation. To replace the traditional Darcy’s
equation for velocity of individual phase they proposed momentum equation for each phase
that is rigorously derived and supplemented by balance equations (based on Helmholtz free
energy) of interfaces and entropy inequality. The total system of equations has 14 equations
and 14 unknowns providing a firm foundation for a theoretically sound model of two-phase
flow in porous media (see table 2 in [68]). In this system of equations, macroscopic capillary
pressure appear in two following equations besides the two momentum equations (Equations
27b and 38 in [68]).
∂sw
L = pc − (pn − pw) (2.16)
∂t
− γwn∇awn + ϵ(Ωn − Ωw − pc)∇sw = 0 (2.17)
The macroscopic capillary pressure is related to material coefficient L, rate of change of
saturation and the phase pressures in the two equations. The macroscale capillary pressure
is related to other thermodynamic variables by[67, 68],
w n ∑ αβ αβ αβ
pc = − w ∂A ∂A a Γ ∂A
ρ sw( ) − ρnsn( ) − (2.18)
∂sw ∂sn ϵ ∂sw
αβ
where pc, ρn, ρw, sn, sw, Aw, An, Aαβ, Γαβ, Aαβ, ϵ are macroscale capillary pressure,
nonwetting phase density, wetting phase density, nonwetting saturation, wetting phase
saturation, macroscopic Helmholtz free energy of the wetting phase, macroscopic Helmholtz
free energy of the nonwetting phase, macroscopic Helmholtz free energy of the αβ interface
33
per unit mass of interface, the average mass of the αβ interface (i.e., the mass of all
αβ interfaces in an averaging volume per unit area), specific interfacial area between αβ
phase (interfacial area per unit volume) and porosity, respectively. This thermodynamic
representation gives the most general functional dependence of macroscale capillary pressure
as the Equation 41 in [67]. After proper simplification, functional dependence of macroscopic
capillary pressure becomes (Equation 46 in [67]):
P c = P c(awn, ans, aws, sw) (2.19)
Upon further simplification,
P c = P c(awn, sw) (2.20)
where awn, ans, aws, sw are interfacial area per unit volume (specific interfacial area)
between the wetting-nonwetting, nonwetting-solid, and wetting-solid phases and wetting
phase saturation, respectively [67, 68]. This revised macroscopic two-phase flow description
was purely theoretical at the time and is non-hysteretic and doesn’t have relative permeability
in it. assanizadeh and Gray acknowledged that the rigorous proposed theory and range of
applicability of its simplified versions have to be verified experimentally [68].
There also have been other approaches to defining a revised model such as phe-
nomenologically derived constitutive relationships of Marle [115] and Kalaydjian [88]. In
recent years, the two phase Darcy’s equation (also called the extended Darcy’s equation)
has been compared with the thermodynamically consistent Darcy velocity originally given
by Hassnizadeh and Gray in [66] by Niessneret et al. [134], which includes saturation
gradient and rate of change of Helmholtz free energy with respect to saturation into the
equation of Darcy velocity [134]. After careful comparison they concluded that during fluid
flow, by introducing either interfacial area and/or a more complex relationship for relative
permeability the traditional Darcy’s equation and the thermodynamically consistent flow
34
theory by Hassnizadeh and Gray can be comparable [66, 134]. It was suggested in all these
approaches, that capillary function is not a unique function of saturation and that incomplete
parametrization results in the apparent hysteresis. While saturation is a volume dependent
property, capillary pressure is inherently a pore-scale phenomenon directly linked to fluid-
fluid interfaces and curvatures, and thus require more than one independent variable to fully
describe it. Hysteresis is therfore a result of projecting a 3D surface (P c − Sw − anw) on a
2D plane (P c − Sw).
All these new modeling approach along with new constitutive relations have been
tested both experimentally and numerically in a number of studies [29, 42, 52, 74, 75, 85,
86, 92, 182, 183]. Most of these works included additional variables to model capillary
pressure such as wetting-nonwetting interfacial area, interfacial curvature and/or Euler
characteristic. Significant improvement was seen in the hysteretic behavior for quasi static
process. It was shown that capillary pressure, interfacial area and saturation do fall on
a unique plane mostly for quasi-static displacement process. For example Karadimitriou
et al. [90] presented 2D micromodel study and quantified effect of spatial interfacial area
in a transient flow experiment for the first time. They defined their macroscale capillary
pressure as averaging microscale capillary pressure over interfacial area and they fitted a
unique P c−Sw−anw surface but noticed discrepancy between the surfaces for transient and
quasi-static cases, suggesting that one interfacial area surface is unable to sufficiently describe
two-phase flow under both transient and quasi-static conditions. They conjectured that one
possible reason can be the phase connectivity since the original theory by Hassanizadeh
and Gray [66] assumed continuous phases while they had discontinuous nonwetting phase.
Their experimental work is also in agreement with numerical work in the same conditions
by Joekar-Niasar et al. and Hassanizadeh et al. by pore network modeling [66, 84, 90].
In another paper, Karadimitriou et al. [92] presented a similar P c − Sw − anw plane
obtained from a quasi 2D micromodel drainage and imbibition experiment which shows
35
Figure 2.7: A nonhysteretic model applied and predicted nonwetting phase volume fraction
with the actual volume fraction for two different models (a and b). (cf. first and second line
of table 5 in [120]).Figure was reproduced based on work in Ref. [120].
unique plane for the fitted equilibrium data as shown in figure 2.8. There have been similar
approaches in the search for a nonhysteretic functional form of capillary pressure and an
extensive review of the research conducted to verify Hassnizadeh and Gray’s hypothesis can
be found in Tables 1 and 2 in Ref. [120].
Central to all the experimental work in capillary hysteresis is the description of capillary
pressure (e.g., how they are defined). As mentioned before and systematically outlined in
[56, 120], the traditional definition of macroscale capillary pressure as the difference between
the boundary pressure in two phases is theoretically inexact and heuristic. First of all, it
does not apply for dynamic condition but might apply at equilibrium depending on how the
phase pressure are measured. Usually, boundary pressure difference is some sort of volume
averaging but macroscale capillary pressure should be derived from microsclae capillary
pressure averaged over all the areas of the interfaces since capillary pressrue is an interfacial
property [56]. This is where the phase connectivity comes into play. Since pressure in a
fluid phase is only representative of the boundary pressure if the fluid phase continuously
36
Figure 2.8: Equilibrium data of saturation, capillary pressure, and specific interfacial area
for drainage and imbibition. (a) actual data points (b) fitted data. Figure was reproduced
based on work in Ref. [92].
connected to the boundary, any presence of trapped ganglia would deem the boundary
pressure difference unreliable as the true macroscale capillary pressure. Also, if the system
is not at true equilibrium, dynamic effects may cause the equilibrium capillary pressure to
deviate.
37
Since transducer-based pressure difference or column-based pressure (pressure difference
between the inlet and outlet) has long been called the capillary pressure in many research in
porous media flows, research has also been done to determine the relationship between the
transducer based capillary pressure and the averaged Young-Laplace pressure in equilibrium
drainage and imbibition experiments [7, 71]. Work has been done to capture interfacial
curvature optically or using x-ray tomography in 3D micromodels to link pore scale capillary
pressure with transducer-based capillary pressure [7, 141]. In this context Pyrak-Nolte et
al. [141] presented comparison between curvature-based capillary pressure and transducer-
based capillary pressure utilizing optical methods in two different micromodels as presented
in Figure 2.9. In this work it is acknowledged that for each equilibrium state, the images
did not exhibit a single curvature value, but instead a range of curvatures were visible.
By averaging all the curvatures, a capillary pressure was calculated from individual image.
The equilibrium time was around 3-5 minutes and phase connectivity was not considered.
Moreover, it is not known to which extent the accuracy of capturing the curvature depends
on the spatial resolution of the optical system.
Figure 2.9: Curvature-based capillary pressure vs. transducer-based measured capillary
pressure. Figure was reproduced based on work in Ref. [141].
Similar studies in the literature to investigate in detail the relationships between
38
curvature-based and transducer-based pressure is present mostly by applying micro-CT in
their experiments [7, 71, 106]. Aside from applying new image processing methods for micro-
CT images, these studies distinguished between interfaces associated with connected and
disconnected wetting-nonwetting phases and achieved better agreement with the transducer-
based pressure when only connected interfaces between wetting-nonwetting phases were
considered. For example, Herring et al. [71] employed micro-CT and a novel 3D curvature
algorithm in Bentheimer sandstone cores, and performed a detailed study of two phase
drainage and imbibition experiments with brine and air [71]. Overall, it was seen that the
image-based capillary pressure (from interfacial curvatures) agrees moderately well when
only connected air-brine interfaces are taken into account (Fig 2.10). Despite this overall
agreement there were some discrepancies which they attribute to low voxel resolution for
small curvatures, low interfacial areas, which did not yield accurate capillary pressure, and
a lack of equilibrium when brine saturation is low. It is worth noting that, disconnected
nonwetting phase has received attention in many studies in the literature but it is uncertain
how trapped wetting phases can influence the true macroscale capillary pressure.
To summarize, we state that, capillary pressure should always be measured at the
microscale and then properly upscaled using the right averaging scheme to apply to
macroscale. Transducer-based measurement is a volume averaging approach and it is still
not clear how that corresponds to the true macroscale capillary pressure. As our first goal
in this work, we aim to compare the two approaches taken into consideration the phase
connectivity. One of the reasons for comparison is to be able to use image based capillary
pressure calculation while evaluating two-phase flow and know how it would correspond to
transducer-based measurements. Also, we want to investigate if spatial resolutions (i.e.,
optical magnification) play a role in the curvature calculations.
39
Figure 2.10: Curvatures measured from micro-CT images of the connected nonwetting phase
(air), and curvatures estimated from transducer-based capillary pressure are compared.
Figure was reproduced based on work in Ref. [71].
Two-Phase Immiscible Displacement and Film Flow in Porous Media
Extensive studies have been dedicated to the visualization and characterization of
immiscible displacement processes at the pore scale experimentally, theoretically and
numerically employing various methods over the past few decades [36, 70, 91, 176]. Starting
with simplified etched networks, now real-time 3D visualization is possible with the
advancement of micro-PIV, synchrotron based microtomography and high-speed cameras
[5, 14, 107]. It is neither possible nor the scope of this short review section to summarize
even a small portion of the literature present in their entirety, but some relevant points will
be drawn.
40
Important aspects of immiscible two-phase flow displacement process have been to
understand the flow physics and predictions of displacement stability, fluid flow pathways,
pore-saturation levels, the spatial and temporal distributions of fluid phases and the relevant
parameters such as fluid viscosity and density, interfacial tension, wetting properties,
geometrical heterogeneity, fluid flow rates, and the considered length scales during both
drainage and imbibition [36, 70, 176]. For example, in a water-wet porous media the wetting
fluid will preferentially coat the walls of the solid grain, and occupy small crevices and
pore throats, while the nonwetting fluid will tend to invade and fill the larger pore spaces
and reside in the center of pore bodies [70]. As mentioned before, drainage is the process
where a nonwetting fluid invades the pore spaces of the medium by displacing the wetting
fluid, whereas in imbibition the opposite happens. Fluid displacement mechanisms are
fundamentally different for drainage and imbibition, and both processes have been studied
in an array of different micromodels under numerous conditions [102, 103].
The complexity of pore space geometry and the interplay among gravity (characterized
by Bond number Bo = ∆ρga2/γ, where a is the pore size [117]), capillarity (characterized
by capillary number Ca = µV/γcosθ), viscous (characterized by viscosity ratio M =
µnw/µw) and inertial forces limit determination of the dynamics and states of drainage
[6, 111, 117, 122]. Our fundamental understanding or modeling of drainage has followed
mostly from experiments performed in 2D micromodels. Lenormand et al. have done many
2D micromodel experiments for both drainage and imbibition and documented physical
mechanisms of meniscus displacements [102, 103]. Depending the relative strengths of
different forces, three flow regimes were discovered, including capillary fingering, viscous
fingering and stable displacement (Figures 2.11 and 2.12). Their equivalent statistical
description (invasion percolation for capillary fingering, diffusion limited aggregation (DLA)
for viscous fingering and anti-DLA for stable displacement) was first explored by Lenormand
et al. and have been reproduced by many different authors in drainage using different fluid
41
Figure 2.11: Phase diagram for drainage (viscosity ratio M = µnw/µw; nw and w refers to the
nonwetting/advancing and wetting defending phases respectively, Capillary number Ca =
µnwunw/γ(nw−w)cosθ; unw is the velocity of the advancing phase, γ(nw−w) is the interfacial
tension between the two phase and θ is the contact angle. Three stability areas are bounded
by dashed lines and the locations of the PEG200 and water (two different wetting phase used)
displacement experiments are represented by pink squares and blue circles respectively. The
gray zones represent the stability areas originally indicated by Lenormand et al. This works
extends the original work of Lenormand et al. Figure was reproduced based on work in
Ref. [176]
42
Figure 2.12: Quasi steady state CO2 -water phase distribution for drainage experiments at
various capillary number with the same viscosity ratio which exhibits the gradual transition
from capillary fingering to viscous fingering as the flow rate increases. The invading phase
CO2, the resident phase water, and the solid grains are shown in green, blue, and black,
respectively. Flow is from left to right. Figure was reproduced based on work in Ref. [107]
pairs and boundary conditions [30, 176]. Considering the microscale nature of these flows in
2D micromodels, gravity effect is negligible, and hence capillary number and viscosity ratio
were the controlling factors. It is worth noting that there have been experiments considering
gravity force as well such as gravity assisted drainage, which however is not relevant to
our experiments [96, 111, 117]. Many attempts have been made in regard to drainage to
indicate which flow regime will dominate a given displacement process (by comparing the
relative importance of gravity, viscosity, and capillarity) both experimentally and numerically
43
(Dynamic pore network model, DNS, LBM) and from a force-balance perspective by many
previous studies.
Imbibition dynamics are much more complicated because of the effect of the pore
geometry of the solid and the flow by film at a very low flow rate [102]. During imbibition
wetting phase travels preferentially through the small crevices and pore throats in porous
media and thickens films on the surface of grains (Figure 2.13). These films swell from
the grain surface into the pore space. In some cases the films merge and completely
emptynonwetting phase fluid out of pore bodies (termed “piston-like” displacement), and in
other cases, they enclose and trapnonwetting phase within pore bodies as thickening menisci
merge at constrictions in the porous media and “snap-off” ganglia (Roof, 1970)[149]. Snap-off
is enhanced during imbibition by lower flow rates, high aspect ratios (the ratio of pore radius
to throat radius), and low contact angles [26, 82, 114, 114, 125, 132, 172, 172, 173]. To better
understand the imbibition mechanisms in porous media, different experiments, mathematical
models, and numerical simulations have been developed using different geometries such as
square single capillary, capillary bundle, square etched network, patterned porous media, etc.
Detailed pore scale mechanism of meniscus displacement was documented by Lenormand et
el. in some of his early works [102, 103]. He published one of the first papers documenting
the piston type displacement and snap-off (called choke-off by Mohanty et al. [123]) during
imbibition and the effect of varying capillary number and pore throat pore body aspect ratio
on the displacement.
Micro-PIV as an experimental technique measures flow field directly, and can provide
a better insight into instantaneous dynamics in micromodel experiments during pore scale
displacements. This optical method of flow visualization relies on the detection of tracer
particles in the flow field. The motion of the seeded particles is used to calculate the velocity
field. In the context of microfluidics, Santiago et al. in 1998 introduced micro-PIV that
used a microscope, micron-sized seeding particles, and a charge couple device (CCD) camera
44
Figure 2.13: (a) weak imbibition (contact angle 60 degree) where water displacing silicone
oil, cooperative pore filling is shown, (b) strong imbibition (contact angle 7 degree) shows
corner flow (c) with decreasing contact angle interfacial curvature forced into the corners
(d) strong imbibition with film swelling over time. Figure was reproduced based on work in
Ref. [178]
to record high-resolution particle-image fields, and the velocity field is calculated by cross-
correlating the acquired particle images [109, 151]. More recently, there have been a number
of studies reporting pore-scale measurements in micromodels using micro-PIV to capture
45
spatially and temporally resolved velocity data of immiscible multiphase flow [17, 93, 94, 148].
Inertial effects and small scales events can be quantitatively evaluated using micro-PIV and
have been done in the context of geologic carbon dioxide sequestration. By using idealized
porous structures or cylindrical posts, Kazemifar et al. visualized and quantified the drainage
process of water by CO2 at the reservoir conditions (80 bar, 24 ◦C and 40 ◦C) and identified
three distinct flow stages (Figure 2.14) [94].
Bulk ow direction 
400 µm
CO2
Figure 2.14: Micro-PIV measurement in a 2D micromodel experiment of supercritical CO2
and water. Quantitative velocity information can tell a lot more about flow unsteadiness
and preferential pathways. Figure was reproduced based on work in Ref. [94]
The literature on films on surfaces is vast. Herein we only focus on larger-scale
phenomena relevant to porous media or the so called macroscopic films or “thick” films
(relative to the nanoscale films), since the hydraulic conductivity of very thin films (at the
46
Figure 2.15: Meniscus displacement mechanism during imbibition (original figure from
Lenormand et al. 1984). Figure was reproduced based on work in Ref. [104]
nanoscale) is very low [39, 104]. Most research on film flow in porous media have been done
in the context of imbibition since the competition between main meniscus and corner film
flow makes imbibition in strongly wetting porous media more complex. Experiments show
that during imbibition, corner film flow in strongly wetting porous media becomes even more
significant with slow displacement rates [178, 180]. Flow through macroscopic films has been
studied for imbibition process by many authors, including Lenormand and Zarcone [104],
Lenormand et al. [103], Dullien et al. [45], and Constandinides and Payatakes [35], among
many others. Lenormand described in detail the expected mechanisms attributed to film
47
flow in imbibition along with drainage displacement mechanism qualitatively in a micromodel
experiment. In addition to elaborating displacement mechanism in drainage based on relative
effect of capillary force, viscous force and injection rate (capillary fingering, stable and
unstable displacement), the particular study also showed how film-flow during imbibition
along the pore wall leads to pinch-off (snap-off), a pore filling mode, where wetting film flow
along grain and filling up a pore eventually [102]. Lenormand also noted that wetting film
flow during imbibition occurs only for strong wettability and very low flow rate [102]. In
another paper Lenormand et al. showed detailed meniscus displacement mechanism during
imbibition in a glass etched micromodel [104] (Figure 2.15). Most of these studies mentioned
above were in single capillary or 2D network of capillaries.
Constantinides et al. developed a pore network simulator for imbibition using 3D
network of pores. In this simulation, water displaced oil in primary imbibition. They take
into consideration the effect of pore wall roughness and wettability. Among their results
they showed that wetting films cause extensive disconnection and entrapment of nonwetting
fluid, especially for small capillary number (Ca ≤ 10−6). Wetting films cause a substantial
increase of the residual oil saturation and the mean volume of residual ganglia increases as
Ca decreases. Swelling of wetting films is illustrated in Figure 2.16 [35].
Tuller et al. developed an analytical model to calculate sample scale saturated and
unsaturated hydraulic conductivity and liquid saturation upscaled from pore-scale based on
chemical potentials. They considered angular pore and wetting film flow in the corners.
They compared analytical results for hydraulic conductivity and saturation with published
soil data and the van Genuchten-Mualem model for drainage and imbibition. Their work
elucidated the important contribution of film flow to the overall conductivity, especially at
lower chemical potential (Figure 2.17) [164].
Yiotis et al. investigated the effect of capillarity-driven viscous flow through macro-
scopic liquid films during the isothermal drying of porous materials, effectively a drainage
48
Figure 2.16: Imbibition considering film flow and micro roughness using a pore network
simulator (pore wall microroughness and wetting film thickness are exaggerated for
illustrative purposes). Figure was reproduced based on work in Ref. [35].
process. Using theoretical and numerical approach they quantitatively characterize film
radius and drying rates among other variables. They showed that film flow is a major
transport mechanism which has dominant effect in a capillarity-controlled displacement
process, which is the case in typical porous media applications [170].
In recent years there have been a few numerical studies to incorporate film flow into
the displacement modeling during imbibition. Zhao et. al demonstrated a single-pressure
49
Figure 2.17: Film thickness swelling and eventually snap off in a triangular pore, black is
the liquid and white is the vapor. Figure was reproduced based on work in Ref. [164].
dynamic pore network model (DPNM) to simulate imbibition in strongly wetting porous
media with corner film flow. Similar as before, they detailed the snap-off mechanism caused
by the thickening of wetting corner film in imbibition. The numerical results were compared
with real microfluidic experiments done by Tsao et al. as depicted in Figure 2.18 [163, 180].
In a more fundamental work Zhao et al. developed the modified interacting capillary
bundle model to describe the imbibition dynamics of a liquid-gas system in a square tube
with corner films. Although it only applies to gas-liquid systems as the viscous dissipation
in the nonwetting gas phase is neglected. The model assumes the pores to be straight tubes
and cannot directly be used in more complex porous media [179].
Hoogland et al. pointed out the importance of film flow in drainage process in their
experiment to measure drainage dynamics from sand columns and assessed the amount of
water retained behind the drainage front, which is expected to drain according to the foam
50
Figure 2.18: Final trapping nonwetting phase distribution and saturation levels for different
injection rates (capillary numbers). (a),(b) and (c) are imbibition experiments done by Tsao
et al. [163] in a 2D 4-layer micromodel. (d),(e) and (f) are the ICB-DPNM simulation by
Zhao et. al (2021) with the same geometry and conditions. (g) is the comparison of trapped
nonwetting phase (air) saturation vs capillary number for both the experiments and the
dynamic pore network modeling. Figure was reproduced based on work in Ref. [180].
drainage equation (FDE). Their results suggested that dynamics of corner flow in porous
media can be represented by a similar governing equation as FDE [77]. In another paper,
they used sand and glass beads and micro-CT to study drainage transition from the rapid
piston-like regime to corner flow dominated regime. For this experiment, a vertical sample
initially saturated with water was used, and air displaced water from top to bottom [77].
Film flow has also been studied in the context of evaporation in a single capillary
[27, 28, 69, 95]. Wu et al. simulated evaporation using dynamic pore network model with
corner flow and obtained good agreement with the evaporation experiment using a quasi-2D
micromodel. They distinguished continuous and discontinuous corner flows in a gas injection
process. Their experiment and dynamic pore network with corner films (called DPNM wC in
their paper) model show good agreement for evaporation rate results. However, the network
model that does not consider corner flow predicts a low evaporation rate, thus confirming
the importance of corner flow (Figure 2.19) [168].
Although corner film flow plays a crucial role in both drainage and imbibition
51
Figure 2.19: Variation of liquid distribution in the pore network during evaporation obtained
from experiment, (Direct pore network modeling with corner films) DPNM wC, and (Direct
pore network modeling without corner films) DPNM woC. St is the total liquid saturation
in the pore network. Figure was reproduced based on work in Ref. [168].
processes, to the best of our knowledge, there has been very few, if any, studies that
focus on the visualization and velocity estimation of corner flow using micro-PIV and PTV
techniques. Additionally, most previous studies on film flow during imbibition and drainage
are performed in extremely simple pore geometry. Given the importance of corner flow
in pore scale displacement process of a strongly wetting medium, a detailed pore-scale
micromodel study of the corner flow dynamics is warranted. Understanding the direct pore-
scale displacement using micro-PIV and PTV will allow us to understand pore-scale physics
and possibly reveal unidentified mechanisms.
52
CHAPTER 3 EXPERIMENTAL METHOD
Micromodels
Micromodel is a microfluidic model porous system. It is an artificial representation
of a natural porous media fabricated from optically transparent materials such as glass
or transparent polymers using either additive or nonadditive manufacturing techniques [4].
Micromodel has become a popular tool to study two-phase flow in the past two decades,
as they effectively overcome the challenges of visualizing the complex flow structures and
fluid dynamics in real rock porous media. Since the emergence of the first micromodel,
dating back to 1950s [4, 89], micromodels of various types including 2D, 2.5D and 3D ones
made from different materials (glass, silicon, PDMS, PMMA, etc.) with different techniques
are widespread. A 2D micromodel consists of single-layer microfluidic channels with an
arbitrary porous structure. The complex flow behavior in 2D micromodels can be visualized
using optical microscopes such as bright-field or epi-fluorescence microscopy [4].
The micromodels used in this study were fabricated using microfabrication techniques,
developed primarily for semi-conductor industries. The fabrication consists of two steps:
i) lithography, which transfers a pre-designed pattern from the photomask to the silicon
wafer; and ii) etching, which generates the desired shape on the micromodel material (i.e.,
silicon) [4, 55, 89]. Typically, the lithography and etching steps are completed in a cleanroom
using silicon wafers. The rest of the process includes photoresist stripping, inspection and
depth measurement using a profilometer, drilling the inlet and outlet holes on the silicon
wafer, piranha cleaning, anodic bonding, dicing, and finally attaching nanoports to serve as
fluid delivery ports. Below we provide more details pertinent to the fabrication process.
53
  

 


		
	Ž		
ƒ	

‹

	
 ƒ	
 
	
Œ
Š




ˆ 

	
	
 	
‰
……†‡€
	 


	
	
 


­ €‚ ƒ€„


Figure 3.1: A workflow chart illustrating the fabrication process used to make the
micromodels in this study.
Micromodel Design
The micromodel fabrication procedure begins with the design of the photomask, which
carries the porous patterns to be etched on a silicon wafer. The masks are deigned with
Adobe Illustrator (2021) and printed on a 4-inch × 4-inch soda lime glass substrate using
direct laser writing method (Front Range Photomask Inc.). A total of four different designs
of micromodels were used in this study. For the capillary pressure experiments (i.e., results
presented in Chapter 4), two different micromodels were designed with one being the idealized
pattern (Type I) featuring a regularly arranged array of circles (Figure 3.3) and the other
ˆ

 	

54
one is heterogeneous (Type II) featuring a porous pattern inspired by real sandstone geology
(Figure 3.2). These two different designs were printed onto two different photomasks with
each mask containing ten identical units as shown in Figure 3.2. For the film flow experiments
(i.e., Chapter 5), one mask that contains both homogeneous (Type III) and heterogeneous
(Type IV) patterns was used as shown in Figure 3.4.
    
Capillary Pressure Micromodel -  V3 - #1 1 mm M-LEFT Capillary Pressure Micromodel -  V3 - #2 1 mm M-LEFT
 


Capillary Pressure Micromodel -  V3 - #3 1 mm M-LEFT Capillary Pressure Micromodel -  V3 - #4 1 mm M-LEFT
Capillary Pressure Micromodel -  V3 - #2 1 mm M-LEFT
Capillary Pressure Micromodel -  V3 - #5                     1 mm M-LEFT Capillary Pressure Micromodel -  V3 - #6 1 mm M-LEFT

Capillary Pressure Micromodel -  V3 - #7 1 mm M-LEFT Capillary Pressure Micromodel -  V3 - #8 1 mm M-LEFT
Capillary Pressure Micromodel -  V3 - #9 1 mm M-LEFT Capillary Pressure Micromodel -  V3 - #10 1 mm M-LEFT
Capillary Pressure Micromodel - V3
Designed by Razin
M-LEFT
08/15/2021 
	
Figure 3.2: (a) The heterogeneous micromodel mask for the capillary pressure experiment
(Type II), (b) single micromodel design, (c) dimensions of the porous section with the air
barrier.
As shown in Figures 3.2, 3.3 and 3.4, a typical micromodel consists of an inlet, an
outlet, and a porous section. While the inlet and outlet are connected to the upstream and
downstream plumbing to serve as fluid delivery, the porous section is the region of interest,
where interesting fluid flow phenomena occur. While the porous sections in designs Type I
and Type III are composed of regularly arranged cylinders, that porous sections in designs
Type II and Type IV are constructed based on the micro-CT scan of a thin sandstone slice.
Quantitative information about the pore body and pore throat size distribution of all the
micromodels is illustrated in Figure 3.6. The algorithm utilized to separate the pore body
and pore throat from an image of a porous medium is also briefly explained in Figure 3.7
	
55

 



	
 




 



Capillary Pressure Micromodel -  V6 - #2 1 mm M-LEFT


Figure 3.3: (a) Design of the idealized micromodel mask for the capillary pressure experiment
(Type I), (b) dimensions of the porous section with the air barrier.
 

 
	


 
		 

 
 
	


 
		 





 


		
	 		
	
Capillary Pressure Micromodel -  V12 - #1 1 mm M-LEFT Capillary Pressure Micromodel -  V13 - #1 1 mm M-LEFT
1 mm
 
 
Figure 3.4: Micromodels used for the film flow characterization experiment: (a) and (c)
single micromodel design for the homogeneous and heterogeneous micromodel respectively
(b) and (d) dimensions of the porous section with the air barrier for the homogeneous and
heterogeneous micromodel respectively.



56
1 mm




Figure 3.5: Nonwetting phase barrier outlined in red at the end of the porous section.(a)
large gap for PIV experiment micromodel (b)smaller barrier gap in heterogeneous design for
capillary pressure experiment.
and also in the image procession section.
57
 200  200
150 150
100 100
50 50
0 0
0 20 40 60 80 0 20 40 60 80 100
Pore radius [ m] Throat width [ m]
 100  100
80 80
60 60
40 40
20 20
0 0
0 10 20 30 40 50 0 20 40 60 80 100
 Pore radius [ m]  Throat width [ m]
100 60
80 50
40
60
30
40
20
20 10
0 0
0 20 40 60 0 10 20 30 40 50 60
 Pore radius [ m]  Throat width [ m]
30 35
25 30
25
20
20
15
15
10
10
5 5
0 0
0 20 40 60 80 100 0 20 40 60 80 100 120
Pore radius [ m] Throat width [ m]
Figure 3.6: (a) Pore body (radius of pore body) and (b) pore throat (throat width) distribu-
tion of the heterogeneous porous medium for the capillary pressure quantification experiment
(c) Pore body (radius of pore body) and (d) pore throat (throat width) distribution of the
idealized porous medium for the capillary pressure quantification experiment (e) Pore body
(radius of pore body) and (f) pore throat (throat width) distribution of the homogeneous
porous medium for the film flow characterization experiment (g) Pore body (radius of pore
body) and (h) pore throat (throat width) distribution of the heterogeneous porous medium
for the film flow characterization experiment.
Count Count Count Count
Count Count Count Count
58
 


 

 
  
Figure 3.7: (a) Cityblock distance method and watershed segmentation separates pore body
and throats. Smallest distance between grains gives the throat width and equi√valent pore
body radius is calculated from approximating each pore body area as a circle ( (πr2/π)).
For a detailed discussion of pore size distribution from digital images readers are referred to
the paper in Ref. [143]
59
In addition to the standard features, a special design of all the micromodels used in this
study is the addition of the non-wetting phase barrier which provides capillary constrictions
narrower than the smallest pore throat in the porous section as shown in Figure 3.5. The
non-wetting phase barrier allows the wetting phase to pass, but prevents the nonwetting
phase (i.e., n-heptane or air) from breaking through the porous section by means of capillary
pressure, which is crucial to running drainage and imbibition in the same experiment in a
pressure-controlled way. The gap of the barrier is within the range of 15 – 4µm in the
capillary pressure micromodels, whereas the the minimum gap was kept at 15µm in the film
flow micromodels, to prevent particle clogging at the barrier.
Revisiting the concept of REV size of any porous medium we want to emphasize that
the REV size is insignificant for the objectives of this thesis. The REV size is especially
significant for the determination of averaged properties (e.g., Pc vs. S) as outlined in Ref.
[92]. In their work the goal was to test a new two-phase flow modeling approach (i.e.,
by inclusion of the interfacial area as a state variable) on elongated micromodel which is
larger than just one REV in order for the model to be valid on macroscale as well. They
determined the REV size by obtaining through simulation capillary pressure vs. saturation
plot for gradually changing subdomains and the smallest domains beyond which the plots
remained unchanged gave the REV size. Their micromodel was four times the size of REV.
Had their micromodel been smaller than one REV their P c−Sw−anw plot might have been
different. One REV and larger gives almost identical results as they have acknowledged as
well.
Now, for the purpose of resolution dependence and phase connectivity study REV size is
not so important since capillary pressure is always defined at the pore scale and comparison
between average capillary pressure based on pore scale curvature vs. the bulk pressure
difference would be insensitive to micromodel domain. Moreover, the type of film flow
studied here always occur at the pore scale and qualitative nature of study deem the REV
60
size insignificant for our micromodels. In general, the micromodels used in the experiments
presented in this thesis are considered as one REV [92].
Micromodel Fabrication
 
 
 
Figure 3.8: (a) Bare silicon wafer (b) dehydration bake (c) HMDS vapor prime oven to coat
HMDS layer (d) Laurel spin coater (e) spin coating PR (f) soft bake.
61
 
 
 
Figure 3.9: (a) Photomask (b) ABM contact aligner (c) aligning the substrate and photomask
(d) UV exposure (e) developing (f) cleaning with DI water.
62
The photolithography process was carried out in the Montana Microfabrication Facility
(MMF) at the Montana State University as illustrated in Figures 3.8, 3.9 and 3.10.
The silicon wafers used in the process are 4-inch standard one-side polished wafer of
500µm (±25µm) thick (University wafer). First of all, a wafer was dehydrated by baking
it at at 110 ◦C for 10 minutes on a hotplate (EMS 1000-3). Following that, the wafer
was coated with a thin layer of Hexamethyldisilazane (HMDS) adhesion promoter using a
vapor prime oven (Yield Engineering System III Vapor Prime Oven) to improve photoresist
adhesion. A positive photoresist (Microchemicals GmbH AZ1512) was then coated on the
wafer using a spin coater (Laurell EDC-650-8B Spin Processor) running at 3000 rpm (3000
rpm/s acceleration) for 30 second, and the final depth of the photoresist is ∼1.5µm according
to the manufacturer data.
Following spin coating, the wafer was soft baked for 1 minute at 110◦C, which removes
the solvent from the photoresist and improves the adhesion of the resist to the wafer. The
wafer was then exposed to UV light for 3 seconds at a dose of 15 mW/cm2 using a contact
aligner (ABM-USA Inc.). The photoresist used herein (AZ1512) is mercury g-line/broadband
(436 nm) sensitive materials which can produce features as small as 1µm. After exposure the
wafer is dipped into the developer (Integrated Micro Materials, AZ300 MIF) and stirred for
50 seconds with an orbital shaker. The wafer was then rinsed and cleaned with DI water and
blown dry with N2, and the features on the wafer were examined under an optical microscope
(Figure 3.10), following which it was post-baked for 3 minutes at 110 ◦C. Post-baking was
the final step of the photolithography process, which is necessary to stabilize and harden the
photoresist prior to etching.
Once the lithography was completed, the wafer was etched in a plasma etcher (Oxford
ICP PlasmaLab 100) as shown in Figure 3.11. A standard recipe was used with about 40
etching cycles to achieve a depth of 20 – 28µm. Briefly, C4F8 and SF6 are used as the
etchant gases, which were alternatively supplied at a rate of 100 standard cubic centimeters
63
 
 
Figure 3.10: Features under bright field microscope after completing photolithography:(a)
and (b) homogeneous design;(c) and (d) heterogeneous design.
per minute (sccm) to the etching chamber. Following the etching step, the photoresist was
stripped off with acetone and cleaned with isopropyl alcohol (IPA). The wafer was then
inspected under profilometer (Filmetrics Profilm 3D) to measure the etching depth of all
the micromodels. Inlet and outlet holes were then drilled on the wafer using diamond drill
bit (Dremel 3000) on all the micromodels, following which the wafer was cleaned in piranha
solution for 15 minutes. Additionally, a borofloat glass wafer (University wafer) of the same
size was also piranha cleaned for bonding with the silicon wafer.
To perform piranha cleaning, a mixture of sulfuric acid (H2SO4) and hydrogen peroxide
(H2O2) was prepared at a ratio of 3:1, which helps to remove organic residues from substrates
64
 
 
Figure 3.11: Dry etching process: (a) Oxford PlasmaLab 100; (b) loading the wafer;
(c) plasma striking; (d) etching completed.
(e.g., photoresist and etchant residuals). Piranha solution also alters the surface property
and makes it hydrophilic by coating a SiO2 layer on top of the silicon wafer. For this study,
60 ml of 98% H2SO4 (VWR Chemicals) and 20 ml of 30% H2O2 (Fisher Chemical) were
used and the cleaning was carried out at 140 ◦C. Once both the glass and silicon wafers
were cleaned, they were rinsed with excess amount of DI water and blown dry to get ready
for anodic bonding. To perform anodic bonding, the wafers were stacked up, placed in
between of two graphite plates, and heated up to 430 ◦C in a furnace (Lindberg/Blue).
Once the temperatures stabilized, a 900 Volt DC was applied between the graphite plates
serving as the anode and cathode for 30 minutes. When bonded, the wafers were diced
65
Figure 3.12: SEM image of the idealized micromodel after etching and photoresist removal.
Figure 3.13: SEM image of the heterogeneous micromodel after etching and photoresist
removal.
into individual micromodels using a dicing machine (DISCO DAD 3350). Finally, nanoports
(IDEX NanoPort UX-02013-65) were attached using epoxy glue to the inlet and outlet for
fluid delivery. Figure 3.13 show the sample scanning electron microscopy (SEM) images after
66
  
 
Figure 3.14: Post lithography and etching operations: (a) drilling holes; (b) piranha cleaning;
(c) anodic bonding (d) dicing and (e) final micromodel after nano-port attachment.
silicon etching for the Type I and Type II designs, respectively, confirming the high precision
of the fabrication methods used herein and Figure 3.14 illustrates the major steps for anodic
bonding and nanoport attachment. The key parameters of all four types of micromodels are
summarized in Table 3.1.
67
Table 3.1: Summary of the micromodels’ properties used in the experiments
Fluid
Average Average pairs
Static
pore pore used Interfacial
Micro- Depth Contact
Porosity throat body (wetting- tension
model (µm) angle
width radius non- (mN/m)
(degree)
(µm) (µm) wet-
ting)
Capillary Water
pressure 0.55 41.29 31.70 23 and n- 30 34
Idealized Heptane
Capillary
Water
pressure
0.49 27.65 19.66 23 and n- 30 34
hetero-
Heptane
geneous
Film
flow Water
0.47 39.1173 38.3210 28 5 72
homoge- and air
neous
Film
flow Water
0.49 40.24 26.5124 28 5 72
hetero- and air
geneous
68
Experimental Methods for Capillary Pressure Quantification
Flow Apparatus and Experimental Procedures
The experiment was done at the Microfluidics Lab for Energy and Fluid Transport
(M-LEFT) at Montana state university. After completing the fabrication process individual
micromodel was prepared before each experiment and a new micromodel was used for each
experiment. The wetting fluid is DI water and the nonwetting fluid is n-heptane (Acros
Organics, 99.5% CAS: 142-82-5). The n-heptane phase was dyed with the nile red dye (Acros
Organics, 99%, CAS: 7385-67-3) at a concentration of 5 mg nile red in 8 ml n-Heptane. A
new work fluid solution was made before every experiment. The entire setup was built on
an optical table (Newport I-2000 series). The micromodel was first affixed on 2 microscope
slides on the two ends using double sided tape and then placed on an inverted microscope
stage (Olympus IX-71) securely taped with electrical tape to make sure no further movement
can occur during experiments. A hydrostatic pressure setup was used to apply known fixed
pressure at the inlet. The system can conveniently vary the applied pressure in a range of 0
– 12 kPa by adjusting the water tank height as shown in Figure 3.15.
Prior to each experiment, all the tubing and plumbing system were cleaned with acetone,
methanol and DI water and blown dry with nitrogen. The micromodel was first saturated
with DI water through the cutoff valve using a syringe and a PTFE tubing (Zeus Industrial
Products, Inc. 1/32 inch ID). The micromodel was made sure to be fully saturated with DI
water and no air bubble was present in the porous section or anywhere upstream (tubing) by
injecting enough water through the tubing and the micromodel. After that, dyed n-heptane
was injected using a 5 ml syringe with a syringe filter attached again through the cut-off
valve (IDEX Health and Science LLC) side. When n-heptane filled about halfway through
the tubing, the injection was stopped and the tubing from the hydrostatic pressure apparatus
was attached to the cutoff valve with the valve closed as shown in Figure 3.16. Then the water
69
Figure 3.15: Schematic diagram showing the setup to apply a fixed pressure at the inlet by
controlling the tank height.
tank was set to the desired height to apply a constant pressure to the inlet. The cutoff valve
was then opened and n-heptane-water interface started to move in the upstream tubing,
eventually entering the porous media to displace the resident water partially out of the
micromodel. We then waited for about 20 minutes for the system to achieve quasi-equilibrium
state (i.e., water and n-heptane configurations remained unchanged) at the applied pressure
before taking images.
For each case the fluid-fluid pattern were imaged four times using four different
objectives corresponding to four different magnifications (i.e., 5x, 8x, 20x, and 40x). While
the field of view (FOV) under 5x and 8x magnifications covers the entire porous section,
the FOVs under 20x and 40x only cover part of the porous section, so for the higher
70
	 





 	 








Figure 3.16: Schematic diagram showing the n-heptane injection through the cut-off valve.
magnifications, only the part next to the outlet was imaged. Once imaging with all four
objectives was done for one pressure, the cut-off valve was closed, the LED turned off to
reduce photobleaching effect and the tank was raised by around 5 cm (i.e., corresponding to
∼500 Pa). The cutoff valve was then opened again and the same process was repeated.
Image Acquisition
As shown in Figure 3.17, all images were acquired using the inverted microscope
(Olympus IX71). The excitation and emission filters in the filter cube are 465(±60) nm
and 500 nm long pass, respectively, which were chosen and optimized for the nile red dye
with excitation and emission peaks of around 480 nm and 550 nm, respectively (Figure 3.18).
A 2x coupler was used to mount the camera to the side port of the microscope, and the
properties of different magnification objectives are listed in Table 3.2. A Phantom VEO
440 camera (Scientific high-speed CMOS camera) and a Thorlab blue LED (SOLIS-470C,
dominant wavelength 470 nm, full width half maximum 34 nm) equipped with an LED
controller (Thorlab DC2200) were used for image capture and illumination, respectively.
The camera has a CMOS sensor of 2560×1600 pixels with a pixel size of 10µm resulting
71






	

		


	




  

	
Figure 3.17: Photo illustrating the optical setup used in the experiment.
in a physical sensor size of 25.6×16 mm2. The camera and LED were synchronized with
a function generator (BNC 565 Pulse/Delay Generator). The exposure time was set to 1
millisecond to help to minimize photobleaching. For each magnification setting under each
applied pressure, 100 frames were captured at a frame rate of 25 fps. The Phantom Camera
Control (PCC) software which was supplied by the camera manufacturer was used to save
RAW cine image files at 12 bit.
Contact Angle Measurement
Contact angle was measured with a goniometer employing the sessile drop method.
The “DropAnalysis” plugin in ImageJ was used to find static contact angle as shown in
Figure 3.19.
72



 
	
 LED
 









Figure 3.18: A closeup view of the optical path inside the microscope.
Table 3.2: Summary of the Objectives’ list and corresponding properties
Spatial
Numerical Effective mag-
Objectives FOV (mm2) resolution
Aperture nification
(µm/pixel)
2.5x 0.08 5x 5.12×3.2 2
4x 0.13 8x 3.2×2 1.25
10x 0.30 20x 1.28×0.8 0.5
20x 0.45 40x 0.64×0.4 0.25
Interfacial Tension Measurement
The interfacial tension between DI water and dyed n-heptane solution was measured
using the pendant drop method (VCA2500XE). With multiple data points measured, an
average interfacial tension between DI-water and n-Heptane was determined to be 34.5 mN/m
73
Figure 3.19: Static contact angle between DI water n-heptane measured on a silicon wafer
of the same surface properties as that of the micromodel. For the micromodels used in this
study, the contact angle was measured to be ∼30 ◦.
at the experimental temperature (i.e., 25◦C).
Magnification Calibration
To accurately convert from the image coordinate to physical coordinate, the magnifi-
cation of each objective needs to be precisely determined, which often slightly deviates from
the manufacturer-specified values. To perform such calibrations, a standard resolution test
target with known division and grids was imaged to back out the magnification. Once the
magnification is determined, the following relationship was used to convert between pixel
number and physical size.
sensor pixel size[µm/pixel]
[µm] = × [pixels] (3.1)
effective magnification
Image Processing
The images were processed in MATLAB R2021a using an in-house code. The code was
based on the image processing toolbox and a number of built-in MATLAB functions. To
74



Figure 3.20: (a) mask image (b) RAW image (c) segmented binarized (green is n-heptane,
black is grain and blue is water).
process the images, the raw cine files were first aligned with the pore mask image, which is a
binary image as shown in Figure 3.20a using image transformation functions. The purpose
of this step is to determine precisely the locations of solid phase (i.e., solid grains) in the
75
 
Figure 3.21: (a) n-Heptane-water interfaces as white (b) zoomed in view of image a.
images, to prepare for the wetting and nonwetting phase segmentation in the next step.
Thanks to the use of fluorescence, the nonwetting phase (i.e., n-heptane), which is dye with
the fluorescent nile red, appears bright, whereas the wetting phase (i.e., water), which is free
of fluorescent dye appears dark. By leveraging the intensity difference, a thresholding scheme,
which automatically and dynamically selects a threshold based on the local intensity and
noise level, was applied to separate the nonwetting phase from the wetting phase as shown in
Figure 3.20b, c. Once segmented, the interfaces between the nonwetting and wetting phases
were then determined, each of which was in turn fitted with a polynomial curve as shown in
Figure 3.21. Sum of the lengths of all interfaces directly yielded the total wetting-nonwetting
interfacial length.
To evaluate the curvature of each interface, here also known as a meniscus, a circle was
fitted to every meniscus with the x and y locations (center of the fitted circle) and radius
R calculated, as shown in Figure 3.22. The microscopic capillary pressure Pc of a given
meniscus was then calculated based on the calculated radius of the meniscus and the Young-
Laplace equation (i.e., Equation 2.12). The microscopic capillary pressure values were then
averaged using the interface length as a weight to get the macroscopic capillary pressure. It
is worth noting that for the averaging, the smallest 20% and the largest 20% values were
76
 
Figure 3.22: (a) Fitted circle on every interface (b) zoomed-in view of image a (fitted circles
are in green).
eliminated to minimize the effect of spurious curvatures, if present.
Circle Fit Algorithm
In order to fit the circle and find out the curvature of each meniscus direct least-squared
fitting was used as demonstrated in Ref. [32]. Any least-squared fitting method is essentially
minimizing the squared differences between the given points and the fitted shape which is
a circle in our case with parameters a, b, and R (x and y coordinates of the center and
the radius respectively). There have been many least-squared circle fit methods present in
the literature such as Kasa fit, Pratt fit or Taubin fit to name a few [32]. If the geometric
distances between the fitted circle and the data points are minimized it is usually called
geometric fit and algebraic fit if algebraic distances are minimized. Here in this thesis, we
utilize a MATLAB open-source code that is based on Pratt’s approximation to Gradient-
Weighted Algebraic Fit (GRAF) [31].
Gradient-Weighted Algebraic Fit (GRAF) is an extension of simple algebraic fit. Simple
algebraic fit is the minimization of the sum of squared differences which is the following
function: ∑n
F 2 2 2 2
1 = [(xi − a) + (yi − b) −R ] (3.2)
i=1
77
This can be rewritten as: ∑n
F1 = [zi + Bxi + Cy 2
i + D] (3.3)
i=1
Where, zi = x2 + y2i i , B = −2a, C = −2b and D = a2 + b2 + R2. Here a, b, R, xi, and yi
are the x coordinate of the center of the fitted circle, y coordinate of the center of the fitted
circle, the radius of the fitted circle, x and y coordinates of the data points respectively.
Differentiating F1 with respect to B, C, and D gives a system of equations and the a, b
and R can be retrieved by solving them. This simple algebraic fit has certain limitations.
A better method is called Gradient-Weighted Algebraic Fit (GRAF). GRAF minimizes the
following function: ∑n [P (x , y )]2i i
F2 = (3.4)
∥[∇P (xi, y 2
i=1 i)] ∥
Where P (xi, yi) = A(x2
i + y2i ) + Bxi + Cyi + D is the polynomial function representing the
circle and gradient of the function ∇P (xi, yi) = (2Axi +B)̂i+ (2Ayi +C )̂j. After calculating
the magnitude of the squared gradient vector and substituting GRAF becomes a problem of
minimizing the following:
∑n [Az 2
i + Bxi + Cyi + D]
F3 = (3.5)
[4A(Az + Bx + Cy + D) + B2 + C2 − 4AD]2
i=1 i i i
Where zi = x2 + y2i i as usual and A = 1. The minimization of F3 is equivalent to the
minimization of the following according to Pratt’s approximation for cases where the data
points (xi, yi) lie close to the circle (Azi + Bxi + Cyi + D ≈ 0):
∑n [Azi + Bxi + Cyi + D]2
F4 = (3.6)
[B2 + C2 − 4AD]2
i=1
This objective function was proposed by Pratt and he achieved the minimization using the
matrix method with the constraint being B2 + C2 − 4AD = 1. Introducing a Lagrange
multiplier η one eventually minimizes the objective function written in matrix form, F5 =
78
ATMA− η(ATBA− 1) where A, B and M are the matrices of the coefficient, constant
and the moments respectively. Minimization of F5 is the problem at hand. It follows (by
differentiating with respect to A) that η is actually a generalized eigenvalue. This problem
leads to an eigenvalue problem which is solved by finding the roots of a quartic equation
(determinant of the matrix M − ηB ) by Newton’s method through an iterative process.
Once the smallest nonnegative root is found one eventually solves for a, b, and R after the
matrix A is recovered where A = (A,B,C,D)T [31, 32].
Experimental Method for Film Flow Characterization
The experimental method for film flow characterization is quite similar to that used for
capillary pressure measurement. Herein, in this section we will only focus on the things that
are done differently in the film flow characterization experiment.
Flow Apparatus and Experimental Procedure
The flow apparatus is the same as the capillary pressure experiment and pressure was
applied at the inlet using the hydrostatic pressure apparatus. Wetting and nonwetting
fluids used in the film flow experiments are DI water and air (instead of n-heptane for the
capillary pressure experiments), respectively. The DI water phase was seeded with fluorescent
tracer particles of 0.5µm in diameter (Thermo Fisher Scientific, F8813) at a concentration
of 150µl stock solution (contains 2 vol./vol.% polystyrene solids) in every 5 ml DI water.
The fluorescent dye contained in the particles works well with the blue (SOLIS 470C) LED
and fluorescent filters used before. As mentioned previously, a new micromodel and a new
particle solution were prepared for every experiment.
To prepare the tracer particle solution, the stock solution was first mixed in a ultrasound
cleaner, following which 150µl particle solution was transferred to a vial containing 5 ml Di
water using a calibrated pipette. The solution in the vial was sonicated for 10 minutes to
79
ensure the tracer particles were well dispersed in water. In the meantime, the micromodel was
mounted onto the microscope as described previously. Then the micromodel was saturated
with the particle solution through the cut-off valve (Figure 3.16) using a 5 mL syringe
equipped with a syringe filter (5µm, Milipore Sigma-Aldrich). Extra care was taken to make
sure no air bubble got into the micromodel or tubing. Once the micromodel is saturated
and all other systems are in perfect order (e.g., camera in focus and acquisition system set
as intended), the hydrostatic pressure apparatus is attached with the cut-off valve at the
desired height of the tank. The valve was then opened and recording was initiated as the
interface approached to the inlet of the micromodel.
Image Acquisition
A scientific CMOS camera (Andor Zyla 5.5mp) with a sensor size of 2560 × 2160 pixels
(5.5 Megapixels and 6.5µm×6.5µm per pixel) was used for this experiment. Compared with
the Phantom camera used in the capillary pressure characterization, this camera offers better
sensitivity and longer imaging period, which are both necessary to image the tracer particles
for as long as 15 minutes. A 10x objective was used, and again the camera and the LED
(Blue LED, SOLIS 470C) were triggered externally using the function generator. Images
were captured at 90 fps with an effective exposure time of 400µs.
Contact Angle Measurement:
The contact angle of DI water seeded with tracer particles on a silicon substrate was
measured using the same goniometer as shown in Figure 3.23. The static contact angle was
determined to be 5 ◦.
Interfacial Tension Measurement:
Interfacial tension between air and DI-water was measured again using the pendant
drop method (VCA2500XE) at room temperature. An average interfacial tension between
80
Figure 3.23: Static contact angle between DI-water and air measured on a silicon wafer of the
same surface properties as that of the micromodel. The static contact angle was determined
to be 5 ◦.
DI-water (seeded with particles at the same concentration used in the experiment) and air
was determined to be ∼73 mN/m at ∼ 25◦C as shown in Figure 3.24, in good agreement
with values in the literature [11].
Calibration
The pixel to millimeter conversion was was done in the same way. An image of the test
target was first taken, where 1540 pixels on the image represent 1 mm in physical dimension,
as shown in Figure 3.25.
Image Processing
The particle images were processed in a commercial software, LaVision Davis 2022.
For a more detailed working principle of micro-PIV and PTV, readers are referred to the
books by [1, 145]. Briefly, 16-bit grayscale images were exported from Andor Solis software
81
Figure 3.24: Interfacial tension measurement using the pendant drop method between air
and DI water.
1540 pixel
Figure 3.25: Pixel to mm scaling for the PIV and PTV experiments.
and imported to Davis project as IM7 files. Using the “Time resolved 2D -PIV (2D2C)”
method in the FlowMaster module, the images were first divided into small regions, called
82
interrogation windows. Cross correlation between the interrogation window in the first
image and the corresponding window in the next image was then calculated employing
a fast Fourier transform (FFT) algorithm. For this study, a final interrogation size of
32×32 pixels was used with 50% overlap, which yielded spatial spacing of the velocity
vectors of 16 pixels. To ensure high-quality calculation, every cross-correlation peak with
a primary-secondary-ratio less than 1.3 was removed as an outlier detection method. The
peak location of the cross-correlation map essentially provides the displacement information
between the two interrogation windows, which combined with the known time interval yields
the instantaneous velocity field. It is worth noting that, the nonwetting phase, which is
not seeded with tracer particles, and the solid grains, which remained stationary, were both
masked out by a thresholding step [108]. For quantification of the film flow, a PTV algorithm
implemented in an in-house MATLAB code was utilized which is motivated by the routines
in Ref. [50].
PTV processes of the image sequences consist of particle detection in an individual
frame and tracking particles across frames in order to find the average velocity magnitude.
A band pass filter is first applied to all the images in the sequence to filter out background
intensity and subsequently increase the signal to noise ratio (SNR). Since our region of
interest is around the posts where the film flow occurs any region outside of the edge is
masked out during particle tracking using the same thresholding step as PIV processing.
The particle locations are detected from the brightest pixel locally and subsequently refined
from the intensity weighted centroid surrounding the brightest pixel to achieve sub-pixel
accuracy [51]. Subsequently all the particle position lists from all the images in the sequence
are generated. By choosing a suitable value for maximum displacement of particles between
two frames particles are tracked and average velocity is found for the whole sequence since
the time interval is known between frames.
83
CHAPTER 4 PORE-SCALE CAPILLARY PRESSURE QUANTIFICATION
This chapter presents the results of the capillary pressure characterization experiments
in both homogeneous and heterogeneous micromodels using n-heptane and DI water as
working fluids. Herein we focus on the effects of measurement spatial resolution and the
effects of phase connectivity on the curvature-based capillary pressure measurement.
Capillary Pressure in Homogeneous Micromodels
In the drainage case, the nonwetting phase (i.e., n-heptane, displaces resident water
as the applied pressure at inlet was increased incrementally. Figure 4.1 shows the phase
configurations at equilibrium under a few selected applied pressures. In this figure, n-heptane
displaces water from right to left, and the images were captured at 5x magnification, where
green, blue and black represent n-heptane, water and solid grains, respectively. It can be
seen that, as the applied pressure increases, the nonwetting phase (i.e., green) front invades
more and more into the porous section in a nearly uniform way by displacing the wetting
phase (i.e., blue) out of the porous section. The phase configurations are expected and
essentially results of the porous structures. While the pore throat size is uniform in the
vertical direction, it gradually decreases from upstream to downstream (i.e., right to left).
According to the Young-Laplace equation, a smaller pore throat corresponds to a large entry
capillary pressure, which requires a higher applied pressure to break through. As a result,
when the applied pressure is increased, the nonwetting-wetting interfaces proceed to smaller
pore throat, where the capillary entry pressure is just high enough to sustain the applied
pressure, and a equilibrium state is achieved therein. It is worth mentioning that during
a drainage process like this, where the invading phase (n-heptane) is less viscous than the
defending phase (water), instability in the form of fingers (i.e., uneven fronts) are expected
to occur during the dynamics process [148, 178]. However, this is not seen in Figure 4.1
84
for two reasons: i) only the final equilibrium stage is shown, which provides sufficient time
for the trapped water to drain through the connected films around the solid grains; ii) the
flow herein is extremely slow, potentially helping to mitigate the dynamic instability. This
observation is in agreement with experimental results by Rabbani et al. in Ref. [144],
where they showed that below a capillary number (flow rate) with certain pore gradient
range (reduced pore size in the direction of flow), stable and uniform displacement can be
achieved.
   
 

 	
 


   

 
 

 

 

 

   


 
 	
 		



  

 	 

 
 

 


 
 
 
Figure 4.1: Equilibrium phase distributions during drainage at different applied boundary
pressures. Here, n-heptane displaces water from right to left. Images shown here were
captured at 5x magnification, where, green, blue and black represent n-heptane, water and
solid grains, respectively.
Figure 4.2 presents similar sample images in the same micromodel for the imbibition
case. The imbibition experiment was performed as soon as the drainage experiment ended
with a sufficiently high pressure (7.45 kPa here). Start from this value, the applied pressure

85
was gradually reduced in steps, causing the wetting phase to imbibe back into the porous
structures. It can be seen from Figure 4.2 that unlike the drainage case, the interfacial
front during imbibition is quite uneven. This is because the wetting phase, instead moving
with the front uniformly, preferentially flows through the corners and boundaries around
the nonwetting phase, leaving a large amount of nonwetting phase trapped in the porous
structure at the end of the imbibition process. The trapped nonwetting phase is expected
to be permanently trapped, as it is completely disconnected from the inlet and outlet
boundaries, and a longer waiting time doesn’t help to mobilize it. In fact, the dynamics
and mobilization of trapped nonwetting phase ganglion is a very interesting research topic,
which is, however, out of the scope of the thesis.
 
 


 
	
   
Figure 4.2: Equilibrium phase distributions during imbibition at different applied boundary
pressures. Similar to Figure 4.1, but for the imbibition cases.
Based on the segmented images as shown in Figures 4.1 and 4.2, curvature-based
capillary pressure can be calculated as illustrated in Figure 4.3. To do this, the pore scale
capillary pressure can be calculated for each individual meniscus based on the Young-Laplace

	
86
   

	 	
R1

Figure 4.3: (a) Major steps followed to get the radius of curvature; (b) the curvature in the
depth direction, R2 was approximated based on the contact angle in the horizontal plane.
Image is reproduced based on the work in Ref. [90].
equation at every pressure step in the entire domain. These individual pressures represent the
pore-scale equilibrium capillary pressure for each pressure step. While the in-plane curvature
was calculated based on the shape of the meniscus in the 2D image, the curvature in the
wall-normal direction is assumed constant due to the constant depth of the micromodel, as
shown in Figure 4.3b. In this regard, the principal radius of curvature R1 in the horizontal
direction for every interface was calculated based on the image, whereas the curvature in the
wall-normal direction was approximated with 2 cos θ/d, where d is the micromodel depth.
These pore-scale capillary pressures for every interface were then averaged over the domain
to give the true macroscale capillary pressure.
Figure 4.4 presents the curvature-based capillary pressure as a function of the applied
pressure for drainage in the homogeneous micromodel captured at 5x. It is seen that all
87
10
8  5x
6
4
2
2 4 6 8 10
Applied Pressure (kPa)
Figure 4.4: Comparison between the curvature-based capillary pressure and the applied bulk
pressure differences at various applied pressures for the drainage case imaged at 5x.
5
4
3 all interfaces
2 only connected
interfaces
1 only
disconnected
interfaces
0
0 1 2 3 4 5
Applied Pressure (kPa)
Figure 4.5: Similar to Figure 4.4, but for the imbibition case. The two data points correspond
to the two images shown in Figure 4.2.
Calculated Pc (kPa) Calculated Pc (kPa)
88
the data points line up reasonably well with a 1:1 ratio, suggesting that the curvature-
based pressure agrees well with the applied pressure for this particular case. Similarly,
Figure 4.5 shows the curvature-based capillary pressure as a function the applied pressure
for the imbibition case. Only two data points were obtained for this case, due the extremely
low pressure required for imbibition to happen. Nevertheless, this plot demonstrates an
interesting behavior. While the calculated pressure and the applied pressure agree reasonably
well when the applied pressure is at 3.08 kPa, a huge deviation is observed when the applied
pressure is reduced to 0.1, suggesting that the applied pressure is not representative of the
capillary pressure at all for the low pressure case. A closer look at Figure 4.2 reveals that
the nonwetting phase (green) is completely disconnected from the inlet and out boundaries
when the applied pressure is reduced to 0.1. In this case, the curvature-based pressure, which
is essentially the pressure jump across the interface between the wetting and nonwetting
phases, is not directly linked with the applied pressure difference between the inlet and
outlet anymore. This simple yet important observation suggest that phase connectivity play
an important role in curvature-based capillary pressure calculation, which will be discussed
in detail in a later section.
Effect of Spatial Resolution
As mentioned previously, one of the focuses of the thesis is to understand the effect
of imaging resolution on curvature-based pressure measurement. Essentially, what is the
necessary spatial resolution to sufficiently resolve a given interfacial meniscus? To answer
that question, the same flow configuration and menisci were imaged at four different
magnifications, yielding four different spatial resolutions. The curvature-based capillary
pressure and other relevant flow properties such as interfacial length were calculated and
compared between all four magnifications. Figure 4.6 presents several sample images
captured at four different magnifications for the same equilibrium phase configuration during
89
drainage. While the FOVs cover the entire porous section at 5x and 8x, the FOVs at 20x
and 40x only cover a portion of the porous section. For 20x and 40x, it was made sure that
the leading front of the invading phase is always in the FOV.
Figure 4.7 presents the total interfacial length between water and n-heptane at various
applied pressures for two different magnifications, 5x and 8x. As shown by the figure, the
interfacial length monotonically decreases as n-heptane progresses to narrower regions with
smaller pore throats and thus short interfaces. The interfacial lengths calculated at 5x and 8x
magnifications agree well, again suggesting that spatial resolution does not have a significant


		

 
 
Figure 4.6: Images captured at four different magnifications for the same equilibrium phase
configuration during drainage. While the FOVs cover the entire porous section at 5x and
8x, the FOVs at 20x and 40x only cover a portion of the porous section. For 20x and 40x,
it was made sure that the leading front of the invading phase is always in the FOV.
90
1
5x
0.8
8x
0.6
0.4
0.2
0
3 4 5 6 7 8
Applied Pressure (kPa)
Figure 4.7: Interfacial length between wetting and nonwetting phases for the drainage case
as a function of applied pressure at two different magnifications.
effect on interfacial length measurement in the homogeneous porous media with fairly stable
displacement patterns. It is worthing noting that only 5x and 8x data are presented here,
as the FOVs in 20x and 40x only covers partially the porous section, which prevents us
from getting the total interfacial length for the entire domain. Nevertheless, we believe the
experimental evidence presented is sufficient to support our conclusion herein.
Figure 4.8 shows the Pc − S curve plotted in the traditional way for both 5x and 8x
cases. The micromodel is initially saturated with water (S = 1) at low pressures. As the
applied pressure increases, more and more n-heptane displaces resident water and the water
saturation decreases, reaching to the irreducible wetting saturation. Upon decreasing the
pressure at the inlet, main imbibition occurs, where water reconnects with the inlet leaving
large oil ganglia behind (called the residual nonwetting phase as shown in Figure 4.2). It is
Interfacial length [mm]
91
8
7 5x
6 Drainag 8x
e
5
4
3
2
1
0
0 0.2 0.4 0.6 0.8 1
Water Satura�on [-]
Figure 4.8: Applied pressure vs. water saturation. The applied pressure is calculated based
on the water tank height, whereas the water saturation is calculated from images captured
at 5x and 8x.
apparent that the curves during drainage and imbibition are quite different, resembling the
classic hysteretic behavior in the PC − S relations. Again very little difference is observed
between the 5x and 8x data, suggesting that spatial resolution does not have a significant
effect on the saturation calculation. Or more precisely speaking, this data shows that 5x
(i.e., 2µm) is sufficient to resolve the phase patterns and menisci in a porous structure
with a pore throat size as small as 8µm. Increasing the measurement resolution does not
significantly improving the measurement accuracy.
Finally, curvature-based macroscale capillary pressure are calculated and plotted for all
four magnifications as shown in Figure 4.9. Again, the spatial resolution of measurement
does not appear to have a significant effect of the results. We note the good agreement and
Applied Pressure (kPa)
i�o
n
bib
Im
92
10
 5x
 8x
8
 20x
 40x
6
4
2
2 4 6 8 10
Applied Pressure (kPa)
Figure 4.9: Comparison between the curvature-based capillary pressure and the applied
bulk pressure differences at various applied pressures for the drainage case imaged at all four
magnifications.
weak dependence on the spatial resolution are in part due to the idealized structure of the
porous section, where the nonwetting-wetting interfacial menisci are uniform and perfectly
circularly, posing relatively less challenge to imaging and curvature calculation.
Effect of Heterogeneity
Similar to the cases in the homogeneous micromodel, drainage and imbibition experi-
ments were performed in the heterogeneous micromodel and the final phase configurations
at equilibrium states were captured for each preset pressure as depicted in Figures 4.10
(drainage) and 4.11 (imbibition). For the heterogeneous micromodel, the entry pressure
was determined to be ∼4.2 kPa, which is slightly higher than the homogeneous case. As
the applied pressure on the inlet is increased, more water is displaced out of the porous
section, forming highly complex displacement patterns due to pore-scale heterogeneity. The
Calculated Pc (kPa)
93
displacement patterns are partially due to the variable sizes of pore throats: at a given applied
pressure the larger pore throats can be invaded due to lower capillary entry pressure, whereas
the smaller pore throats cannot be invaded. The effect is further enhanced by capillary
fingering during the dynamics pore invasion process, where nonwetting phase grows in all
directions and even backwards due to flow instability. It is also apparent in Figure 4.10 that
this process leaves a large amount of the wetting phase (water) surrounded and trapped by
the nonwetting phase, and the wetting phase appear to be disconnected from the inlet and
outlet boundaries, which, as discussed later, has a significant effect on the curvature-based
capillary pressure.
For the imbibition case shown in Figure 4.11, again the wetting phase imbibes back
into the porous section increasing the water saturation therein, as the applied pressure is
   
 
  


 
    
  


  


 


	 
	
 
	  	
  	


   
Figure 4.10: Equilibrium phase distributions during drainage at different applied boundary
pressure in the heterogeneous micromodel. N-heptane displaces water from right to left.
These particular set of images are captured at 5x magnification. Green represents n-heptane,
blue is water and black is solid grains.

94
reduced at the inlet. It can be seen that the wetting phase trapped during the drainage
process remains trapped for nearly the entire imbibition process due to the heterogeneous
pore structure. Towards the end of the imbibition process (i.e., 100 Pa), while the bulk
wetting phase get connected with the inlet boundary, a large portion of the nonwetting
phase get trapped and disconnected from the boundaries, similar to the homogeneous case.
Moreover, at this final stage of imbibition, the fluids within the porous section can be
classified into four different groups: connected wetting phase, disconnected wetting phase,
connected nonwetting phase and disconnected nonwetting phase. As this is not observed
in the homogeneous case, it is presumably due to the pore-scale heterogeneity of porous
structures. Again as discussed later, the connectivity of the phases turns out to play an
  

 
 
	
 

  

  

­
  
	
 	
 
 


  

­­
  

		



 
 
 


Figure 4.11: Equilibrium phase distributions during imbibition at different applied boundary
pressures in the heterogeneous micromodel. Water displaces n-heptane from left to right.
These particular set of images are captured at 5x magnification. Green, blue and black
represent n-heptane, water and solid grains, respectively.

95
important role in curvature-based capillary pressure quantification.
Similar to the homogeneous cases, each case was imaged with four different magnifica-
tions (Figure 4.12) and the pore-scale capillary pressure was calculated for every meniscus
and averaged over the measurement domain to get the macroscale curvature-based capillary
pressure for every pressure steps. Figure 4.13 shows the curvature-based capillary pressure
at various applied pressures with all menisci (interfacial curvatures) considered for all four
magnifications. While the pressures line up relatively well for low pressure applied (i.e., high
water saturation), large deviations are observed for high applied pressures. This is true for all
four magnifications, again suggesting negligible effect of measurement resolution. We believe
the large discrepancy at high applied pressures is due the disconnected phases that are vastly
present in the heterogeneous porous structures under both drainage and imbibition cases.
Briefly, when a blob of nonwetting or wetting phase get disconnected from the inlet and
outlet boundaries, the interfacial curvature formed at the blob surface become uncorrelated


		

 
 
Figure 4.12: Sample images captured at all for different magnifications for the same
equilibrium state during drainage. Again for 5x and 8x, the FOVs cover the entire porous
section, whereas FOVs in 20x and 40x only cover portion of the porous section.
96
with the applied boundary pressure. Therefore, the microscale capillary pressure at the
particular blob is not effectively represented by the boundary pressure difference, causing
discrepancies in the two sets of data. This point will be revisited in detail in the next section.
For the sake of completeness, the Pc − S curve is plotted in Figure 4.14 for 5x
and 8x cases, which again resembles the classic hysteretic capillary pressure elation. As
before, the pore-scale capillary pressure was calculated for each individual meniscus from
the Young-Laplace equation at every pressure step in the entire domain and then averaged
over the domain to give the true macroscale capillary pressure. Again this curvature-based
macroscale capillary pressure represent the true pore-scale equilibrium capillary pressure,
which, however should be distinguished from the applied bulk pressure difference between
the inlet and outlet. The spatial resolution again demonstrates minimal effect on capillary
10
8
6
4
2
2 4 6 8 10
Applied Pressure (kPa)
Figure 4.13: Curvature-based macroscale capillary pressure vs. the applied boundary
pressure for the heterogeneous case measured with four different magnifications.
Calculated Pc (kPa)
97
8
6
4
8x
2
 5x
0
0 0.2 0.4 0.6 0.8 1
Water satura�on [-]
Figure 4.14: Applied pressure vs. saturation calculated for the heterogeneous case with both
5x and 8x magnifications.
pressure measurement.
Effect of Phase Connectivity
Based the results shown previously, we note the following: i) the lowest measurement
resolution (i.e., 2µm) we are using is sufficient to resolve the curvature, and further improving
the accuracy do not significantly affect the measurement; ii) the curvature-based capillary
pressure agrees relatively well with the applied bulk pressure difference when all fluid phases
are sufficiently connected with the boundaries. Following the discussion, it is apparent
that one must be very careful when comparing bulk pressure and curvature-based capillary
pressure, as the pressure at the disconnected fluids is not necessarily the same as that
maintained at the inlet or outlet due to a lack of connectivity. To that end, this section
Applied Pressure (kPa)
98




Figure 4.15: Connected and disconnected water is shown separately for the case of applied
pressure of 5.6 kPa and 8x magnification. Disconnected water and connected water regions
shown in pink and green, respectively. It is visible that the interfacial curvature (blue)
associated with disconnected water appear smaller (i.e., radius is larger) than (yellow)
that associated with the connected water region, indicating that the pressure within the
disconnected water region is slightly higher than that in the connected water region, given
that the pressure in the nonwetting phase is uniform across.
is dedicated to a detailed discussion on the role of phase connectivity.
Figure 4.15 shows a sample image captured at 8x for a drainage case at 5.6, where
the connected wetting phase (water) region and disconnected water regions are treated
separately. In this particular image, the nonwetting phase (gray) is all connected with the
inlet, and the connected wetting phase and disconnected wetting phase are demarcated in
green and pink, respectively. It is visible that the curvature of the interfaces (blue) associated
with disconnected water appears smaller (i.e., radius is larger) than that associated with the
connected water region (yellow), indicating that the pressure within the disconnected water
region is slightly higher than that in the connected water region, given that the pressure in
the nonwetting phase is uniform across.
To test this hypothesis, the histograms of the radii of all fitted circles for this particular
99
 
all Interfaces connected water-oil
interfaces

disconnected 
water-oil interfaces
Figure 4.16: Histograms of the radii occurring for the drainage case at 5.6 kPa and 8x.
(a) all the interfaces (yellow and blue circles); (b) only the interfaces created between the
connected water and n-heptane (yellow circles only); (c) only the interfaces created between
the disconnected water and n-heptane (blue circles only).
case (i.e., drainage at 5.6 kPa) are calculated and presented in Figure 4.16. It can be seen
that when only the interfaces associated with the connected wetting phase are considered
(Figure 4.16b), the radius values are populated around the value of 9 pixels with minimal
values greater than 15 pixels. However, when the interfaces associated with the disconnected
phase are considered (Figure 4.16c), the radius values are populated around 10 pixels, with
many more values greater than 15 pixels. It is apparent that the systematic shift of radius
100
values has caused the underestimation of the interfacial curvature, and thus the lower values
of the curvature-based capillary pressure previously as shown in Figure 4.13. With that
in mind, the curvature-based capillary pressure is recalculated by ignoring the interfaces
associated with the disconnected water regions for the drainage case and the results are
plotted in Figure 4.17. It is seen that this corrected curvature-based capillary pressure
shows much better agreement for higher applied pressure during drainage, confirming our
hypothesis that only the pressure in the connected phase is related to the applied boundary
pressure difference. However, it is worth noting that, although a better agreement is
achieved between the correct capillary pressure and the applied pressure, it is not guaranteed
that the corrected capillary pressure (compared with the uncorrected value) offers a better
representation of the state of the pore-scale flow. In fact, a good macroscale pressure should
be able to account for both connected phases and disconnected phases and be used in
predictive models.
The same behavior is also observed in the imbibition case in the heterogeneous
micromodel. Figure 4.18 shows a sample image during the imbibition case, where again the
interfaces associated with disconnected regions and connected regions are treated separately.
In this case, the interfaces associated with the connected wetting phase display larger radii of
curvatures (yellow), whereas the interfaces associated with the disconnected wetting phase
display smaller radii of curvatures (blue), again suggesting different pressures within the
connected and disconnected wetting phases. This also suggests that if there are wetting
water films linking between the connected wetting phase regions (green) and the disconnected
wetting phase regions (purple), there must be a flow through the film driven by this pressure
difference, which we will focus on in Chapter 5. It turns out that this discovery about phase
connectivity can similarly explain the “outlier” data point in Figure 4.5. Essentially, in
traditional capillary pressure measurement, a big assumption is that the nonwetting phase
(the invading phase) is always connected to the inlet, whereas the wetting phase is always
101
10
8x Connected
Interfaces Only
8
8x All Interfaces
6
4
2
2 4 6 8 10
Applied Pressure (kPa)
Figure 4.17: Corrected curvature-based macroscale capillary pressure vs. the applied
boundary pressure for the heterogeneous micromodel for the drainage case imaged at 8x.
It is seen that this corrected capillary pressure shows much better agreement for higher
applied pressure during drainage.
connected to the outlet. Therefore at equilibrium (i.e., when there is minimal flow within
the porous structure), one can say the pressure difference between the inlet and the outlet
equals to the pressure jump across the interface (i.e., capillary pressure). However, when one
of more phases get disconnected from the inlet or outlet boundaries, the assumption breaks
down, causing the two quantities to be independent, and thus no agreement is observed in
Figure 4.5.
To summarize, the first point we want to make is that in all these cases whether drainage
or imbibition and in any micromodel, the true macroscale capillary pressure is always the
ones calculated based on the pore-scale curvature of all the interfaces whether connected
Modified Macroscale Pc (kPa)
102


 

 


	





		
Figure 4.18: Imbibition case where the applied pressure is 1355 Pa. This case is for 5x.
Connected and trapped water is shown in green and purple respectively. Blue circles are
associated with the disconnected water and yellow circles are the ones between connected
water and n-heptane.
or disconnected. Our goal was to compare this measurement with the applied boundary
pressures and systematically study the effect measurement resolution and phase connectivity.
Our results have shown that measurement resolution (i.e., 2µm or smaller/pixel) does not
seem to have a significant effect on the calculated capillary pressure, saturation and interfacial
length for both homogeneous and heterogeneous cases. The results have also shown that
103
boundary pressure and macroscale capillary pressure calculated from averaging pore-scale
capillary pressure over the interfaces are two separate measures in their own domain. One
must be careful when using them interchangeably. As explicitly outlines in chapter 2 and 3,
pressure with individual phase can only have the same pressure as that of a reservoir if they
are connected. As we have seen in drainage for the heterogeneous case and imbibition for
both micromodels, any presence of disconnected phases significantly give rise to pore-scale
pressure variation within disconnected phases. More importantly, when nonwetting phase
is disconnected we have a significantly higher macroscale capillary pressure that cannot be
estimated based on pressure at the boundary. We also note that the curvature-based capillary
pressure is always defined precisely both at equilibrium and transient conditions, and the
right way to go from pore scale to macroscale is to average the right pore-scale capillary
pressure (depending on equilibrium or dynamic) over fluid-fluid interfacial area. Transducer-
based pressure measured at flow boundaries is theoretically inexact and can cause deviations.
Finally, understanding the corner film flows that run between the connected and disconnected
wetting phases could be a key to furthering our understanding of the equilibration process
and curvature-based capillary pressure.
104
CHAPTER 5 FILM FLOW CHARACTERIZATION
As demonstrated by our results in Chapter 4, the interfaces (menisci) associated with
connected and disconnected phases have very different curvatures for both drainage and
imbibition, suggesting that a true equilibrium may not be achieved within the porous media.
Given that the pressure in the nonwetting phase is nearly a constant with minimal flow, the
pressures within the connected wetting phase and the disconnected wetting phase must be
different. With that in mind, if there are wetting films linking the connected wetting phase
regions and the disconnected wetting phase regions, there must be flows through the films
driven by this pressure difference. Proper characterization of this film flow can provide insight
into the flow equilibration process, further our understanding of pore-scale capillary pressure,
and help us identify new pore-scale physics and flow mechanisms. Inspired by this, in this
Chapter, we attempt to present a characterization of corner film flows that occur during the
equilibration process under drainage and imbibition conditions in both homogeneous and
heterogeneous micromodels. Although this current characterization is relatively qualitative
and preliminary, to the best of our knowledge, this study represents the first one to quantify
corner film flows under both drainage and imbibition conditions.
Single-Phase Flow Validation
Comparison Between Calculated and Applied Flow Rate
Single phase flow of DI water in both homogeneous and heterogeneous micromodel were
performed at different constant flow rates and compared against the calculated flow rates
based on the average velocity at the outlet. As it is shown in Figure 5.1 reasonable agreement
between calculated and applied flow rates has been achieved. This serves as a validation of
PIV algorithm used in this study.
We want to mention that, the deviations are mostly due to the syringe pump fluctuation.
105
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.1 0.2 0.3 0.4 0.5 0.6 0.7
Applied Flowrate (µL/min)
Figure 5.1: Match between calculated flow rate and applied flow rate in the homogeneous
micromodel for a single phase flow. Errors bars represent the combined uncertainities while
calculating flow rate from average velocity at the outlet. Maximum relative difference is 14%
that occurs for the 0.4 µL/minute applied flowrate case. We attribute this difference to the
syringe pump instability. The syringe pump instability has been demonstrated in Figure 5.2.
As seen in Figure 5.2, for applied flow rate of 0.6 µL/minute, the calculated steady-state
flow rate appear to me pulsating indicating syringe pump instability.
Single Phase Velocity field
The velocity field of the steady-state flow of seeded DI water was calculated as shown
in Figures 5.3 and 5.4, for the homogeneous and heterogeneous micromodels respectively.
In the homogeneous porous section, due to regularly arranged pore geometry the velocity
distribution is periodic whereas in heterogeneous porous section the velocity magnitudes
vary randomly depending on pore size. The maximum velocities are 0.5024 mm/s and
0.6333 mm/s for the homogeneous and heterogeneous micromodels respectively. Although
Calculated Flowrate (µL/min)
106
Fluctua�on in the calculated flowrate 
for 0.6 µL/min applied flowrate 
0.68
0.67
0.66
0.65
0.64
0.63
0.62
0.61
0.6
0 1 2 3 4 5 6
Time (Second)
Figure 5.2: It is visible that calculated flow rate for 0.6 µL/minute applied flow rate fluctuates
in the higher range. Possible reasons for this is syringe pump’s inability to stabilize the flow
rate from one value to another, poor calibration and built-in inaccuracy.
the applied flow-rate is less for the heterogeneous case the maximum velocity is higher than
the homogeneous case because of the pore heterogeneity and smallest pore size being smaller
than the constant pore size of homogeneous design (18 µm as opposed to 38 µm). The
location where the maximum velocity occurs is marked as red star in Figure 5.4.
Calculated Flowrate (µL/min)
107
mm/s
Figure 5.3: Steady state velocity field in the homogeneous micromodel. DI water is flowing
from right to left. The applied flow rate is 0.2 µL/minute. The maximum velocity occurs at
the throats which is 0.5024 mm/s.
108
mm/s
Figure 5.4: Steady state velocity field in the heterogeneous micromodel. DI water is flowing
from right to left. The applied flow rate is 0.05 µL/minute. The maximum velocity occurs
at the throats which is 0.6333 mm/s. The location of the maximum velocity is marked with
a red star.
To summarize, from the single-phase flow experiments, we have learned that our PIV
algorithm used in this study is robust and reasonable as the velocity fields and calculated
flow rates explain. In the subsequent discussions, we present the results of the two-phase
flow experiments aimed at characterizing the corner film flow.
109
Drainage and Imbibition in Homogeneous Micromodel
This section describes in detail the results of the experiment performed in the homo-
geneous micromodel. Results of both the drainage and imbibition process is systematically
presented with the primary goal of film flow study in mind.
Drainage
In order to perform the drainage experiment, a fixed boundary pressure was applied
similar to the experiments for capillary pressure characterization in Chapter 4. The only
difference is that instead of varying the pressure, this time one fixed pressure was applied
which was 8.95 kPa. As the applied pressure was set to 8.95 kPa, the nonwetting phase (i.e.,
air here) invaded the porous section and displaced resident water. During this transient
process, the air-water interface velocity is initially high and gradually slows down as more
pores are invaded. What’s noticeable and different from the idealized geometry is that
instead of stable displacement (as in Chapter 4) the air entraps water downstream initially
as these trapped water regions are connected to the outlet by the thin films, which helps
drain these water regions over time. The total recording time for this drainage was 12.83
minutes. Next, we present the results of the volumetric flow rate vs. time measured at the
outlet, a few representative velocity vector plots, and film flow plots obtained via PTV.
Flow rate vs. time
The instantaneous flow rate is calculated for the entire recording time and plotted as
a function of time as shown in Figure 5.5. Time zero corresponds to the first pair from
where the PIV processing began. This starts from frame 891 and the last frame is 70000.
These frames are shown in Figure 5.6 a and b respectively. The initial and final distribution
can be seen on these figure. Most of the pore invasion occurred during the first 30 seconds
and multiple trapped water regions were left behind most of which were eventually drained
110
out by film flow over time. The flow rates at the latter part correspond to these film flow.
Maximum flow rate is 9.6 µL/minute.
Frame number: 891 Frame number: 70000
10 Maximum Flow Rate
Applied Pressure: 8.95 kPa
1
0.1
0.01
0.001
0.0001
0.00001
0.000001
0.0000001
1E-08
0 30 80 130 180 230 280 330 380 430 480 530 580 630 680 730 780
Time (second)
Figure 5.5: The instantaneous flow rate is plotted against time on a semi-log plot. Time zero
corresponds to the first pair where the PIV processing began and the air phase just enters
the porous medium. This is frame 891 shown in Figure 5.6a. The last frame in this recording
is 70000 which is shown in Figure 5.6b. The whole invasion takes about 760 seconds. During
this displacement process, initially, the air interface velocity is high and maximum velocity
(thus flow rate) occurs at around 14.5 seconds. The sharp peaks represent the Haines jumps
(i.e., pore-scale event when the nonwetting phase passes a narrow pore and very rapidly fills
the next wider pore space) [5, 6, 14]. The applied pressure is 8.95 kPa and the flow rate is
calculated at the outlet.
Flow Rate (µL/minute)
111
Frame 891 t=0s Frame 70000 t=760.199s
(a) (b)
Figure 5.6: (a) The first frame in the series where PIV processing began. (b) last frame in
the series.
Velocity Fields
In this section the instantaneous velocity vector plots at certain time during the drainage
process are presented. Since water is the seeded phase cross correlating 2 frames directly
gives the velocity vectors in the water phase from the known time step between those frames.
The vector plots presented herein also have been segmented to show the air and grain phase
separately. The velocity of the water after opening the valve (before air first appears in the
FOV) is very high. Similarly, the first time air invades the pores the displacement rate is
high exceeding the sustainable frame rates needed to resolve the flow. Also during first burst
events (called Haines jumps) velocity in the vicinity of pores exceeds by order of magnitude.
Therefore, We only present a few representative cases here. As shown in Figures 5.7, 5.8
and 5.9, with time more and more air (white) displaces water. Air is displacing water from
right to left in all the drainage cases. The velocity vectors corresponds to the instantaneous
velocity in the water suggesting nonuniform flow. Whenever a new pore is invaded at a
location, or flow is taking place via multiply connected thin films, velocity of the water
phase is high in the vicinity of that location. As we will see in the next section, a continuous
112
film flow is actually present and connected to the outlet which is essentially responsible for
the trapped water to be displaced as shown in Figure 5.11.
Frame 902-903 t=9.922s
mm/s
Figure 5.7: Instantaneous velocity vector field between the frame pairs 902-903. The applied
pressure is 8.95 kPa. Maximum velocity is 0.8324 mm/s. White and gray represent air and
grain respectively.
113
Frame 1116-1117 t=12.276s
mm/s
Figure 5.8: Instantaneous velocity vector field between the frame pairs 1116-1117. The
applied pressure is 8.95 kPa. Maximum velocity is 0.5222 mm/s. White and gray represent
air and grain respectively.
114
Frame 14695-14706 t=161.645s
mm/s
Figure 5.9: Instantaneous velocity vector field between the frame pairs 14695-14706. The
applied pressure is 8.95 kPa. Maximum velocity is 0.0384 mm/s. White and gray represent
air and grain respectively.
115
Frame 1116-1117 t=12.276s Frame 14695-14706 t=161.645s
Figure 5.10: The region marked in red is one of the regions of active film flow. The two
images show that air has displaced most of the regions. The high velocity at the top in outlet
(barrier) means that the flow contribution of all the film flow.
Film Flow
As one of the goals of this study is to verify and characterize the film flow during
drainage we present some qualitative results in this section. Overall, film flow was seen to
be present and be an active transport mechanism. For example, if we look at Figures 5.8
and 5.9, we see the air almost fills the pores that were not initially invaded and trapped
water was present. This comparison is outlined in Figure 5.10. This occurred because of the
presence of multiple films that are connected by each other eventually leading to the outlet.
It is worth mentioning that it is possible for completely isolated water ganglia to occur that
are disconnect from the outlet and remain trapped for the entire recording period. Our focus
here is the trapped water that are connected to the outlet. To demonstrate this, We have
processed frames 14650-14700 via PTV and averaged the velocity magnitude over time as
shown in Figure 5.11. Is is evident that after around 2.7 minutes (first frame being t=0) the
film flow is the dominant mechanism of drainage which suggests that pressure in the trapped
116
water phase is higher than the bulk water phase at the outlet driving the flow. This is in
agreement with the findings in Chapter 4 that the curvature between trapped water and
the nonwetting phase is smaller (larger radius). Moreover, it is observed that water films
that contribute to flow are all connected although disconnected films were seen around just
one post (or multiple) as mentioned before. The velocity magnitude is also varied as can be
seen in the Figure 5.11 as s simple result of conservation of mass. The film that is near the
outlet is relatively faster (represented by the single black line) than the branch films that
are connected downstream (represented by the green and red lines).
6 5 mm/s
4
3
1
2
Frame 14650-14700
Figure 5.11: Film flow around the solid posts are tracked via PTV. The average velocity
magnitude is shown by the color map. For better understanding of the flow pathway, film
flow directions are indicated by the three superimposed lines around the grains (red, green
and black). Also, it is worth mentioning that only the particles around the grains walls are
tracked and larger water pockets are not included in the PTV analysis.
Outlet
117
Imbibition
After drainage, we perform similar imbibition experiment as before. This time only
one pressure step is applied. The pressure at the inlet was decreased from 8.95 kPa to 2.55
kPa. After a very short quiescent period the water starts to imbibe back displacing the air.
The total recording time was approximately 4.42 minutes. Next, similar to the drainage
experiment we present the results of the flow rate vs. time, velocity plots and the film flow
study via PTV.
Flow rate vs. time
Similar to drainage, instantaneous flow-rate has been calculated and plotted as a
function of time as shown in Figure 5.12. There are a few observations here. Firstly, the
first frame in this series is the last frame from the drainage image set (frame number 1 as
seen on the top left of Figure 5.12). Secondly, the last frame of the series (top right of Figure
5.12) indicates that almost complete displacement of air by water is achieved at the end
of imbibition. Also, the initial quiescent period after the valve is open is noticeable on the
plot. The maximum flow-rate is 2.77 µL/minute which unlike drainage is low but flow is
nearly continuous after the initial pause until water reconnects to the inlet. Moreover, the
irreducible air saturation is also low (almost no trapped air at the end of imbibition).
118
Frame number: 1 Frame number: 24111
10
Maximum Flow Rate
1 Applied Pressure: 2.55 kPa
0.1
0.01
0.001
0.0001
0.00001
0.000001
0.0000001
1E-08
0 30 80 130 180 230 280
Time (second)
Figure 5.12: The instantaneous flow-rate for imbibition is plotted against time on a semi-log
plot. Time zero corresponds to the first pair from where the PIV processing began after
the valve has just been opened. The whole invasion takes about 265 seconds. The applied
pressure is 2.55 kPa and the flow rate is measured at the outlet as usual.
Flow Rate (µL/minute)
119
Velocity Fields
The instantaneous velocity vector plots are again calculated during imbibition for some
representative frame pairs. Water is displacing air from left to right in all the figures. As we
have seen in the flow rate plot that after the quiescent period the flow is almost continuous
until water almost completely fills the porous section displacing air (and reconnects to the
inlet) is also understood from the velocity plots in Figures 5.13 and 5.14. We can see within
5 seconds water saturation almost doubled if we compare the two figures indicating that rate
of change of saturation is high during imbibition.
Frame 2681-2682 t=29.491s
mm/s
Figure 5.13: Instantaneous velocity vector field between the frame pair 2681-2682. The
applied pressure is 2.55 kPa. Maximum velocity is 0.6695 mm/s. White and gray represent
air and grain respectively. Water displaces air from left to right during imbibition.
120
Frame 3169-3170 t=34.859s
mm/s
Figure 5.14: Instantaneous velocity vector field between the frame pair 3169-3170. The
applied pressure is 2.55 kPa. Maximum velocity is 0.6691 mm/s. White and gray represent
air and grain respectively. Water displaces air from left to right during imbibition.
121
Frame 2683 Frame 3045
(a) (b)
Figure 5.15: PIV raw images showing film flow and snap-off. Water films around grains
results in swelling of films. The red rectangle in the two images show one of the regions
where this happens. Time interval between these two frames is 3.982 seconds.
Film Flow
Similar to drainage film flow is seen to be equally important and an active mode of
transport during imbibition. We more specifically focus on one mechanism of invasion which
is snap-off. For example snap-off (i.e, an event where spontaneous filling of the throat by
wetting phase occurs [157]) is enhanced by film flow. As the pressure in the air decreases
water starts to fill up the narrowest regions and eventually completely fills the pore-space.
This event is caused by the film flow. As can be seen in Figure 5.15 within 3.982 seconds
water filled regions that are apparently disconnected but are actually connected by water
films around the grains which the PTV image (Figure 5.16) shows with the direction of
the films by the four superimposed lines (red, green, orange, and blue lines with arrows)
around the grains. The magnitude of the average velocity is also shown by the color map.
What we can understand from film flow during imbibition is that water films ahead of the
main meniscus aids in filling the throats which can also disconnect the nonwetting phase.
Moreover, we saw in Chapter 4 that during imbibition radii of curvature of the interfaces
122
between connected water-nonwetting phase and the disconnected water-nonwetting phase
were different. The interfaces associated with the connected wetting phase displayed larger
radii of curvatures, whereas the interfaces associated with the disconnected wetting phase
displayed smaller radii of curvatures suggesting different pressures within the connected and
disconnected wetting phases. This phenomenon can again be explained by this film flow
assisting in snap-off during imbibition as displayed by the film flow PTV plot (Figure 5.16).
mm/s
Frame 3016-3035
Figure 5.16: PTV image averaged over 20 frames. Water films around grains results in
swelling of films. The direction of the films are shown by the 4 lines with arrows around the
grains. Water is displacing air from left to right.
Outlet
Inlet
123
Drainage and Imbibition in Heterogeneous Micromodel
This section describes very briefly the results of the experiment performed in the
heterogeneous micromodel. Similar to the homogeneous micromodel a fixed applied pressure
at the inlet was applied and the displacement process was recorded during both drainage
and imbibition. This experiment supplements the results already presented in the earlier
section.
Drainage
Similar to the homogeneous drainage experiment, a fixed boundary pressure is applied
again. As the applied pressure is set to 9.12 kPa the nonwetting phase (i.e., air here) invades
the porous section and displaces resident water. During this transient process initially the
interface velocity is high and gradually slows down as more pores are invaded. What’s
noticeable again is that the air entraps water downstream initially as these trapped water
regions are connected to the outlet by the thin films, which helps draining these water regions
over time. The total recording time for this drainage is 11.92 minutes. Next, we present
the results of the volumetric flow rate vs. time measured at the outlet, a few representative
velocity vector plots, and film flow plots obtained via PTV.
Flow rate vs. time
The instantaneous flow rate is calculated as before for the entire displacement sequence
and plotted as a function of time as shown in Figure 5.17. The maximum flow-rate is 5.45
µL/minute which corresponds to the largest jump (air displacing water) occurring within the
very first few seconds from where the PIV processing begins. After the initial displacement
the rest of the pore filling takes place over few minutes. The trapped water drains away via
film flow similar to homogeneous micromodel. Moreover, we can get an idea of the change of
water saturation at the end of drainage process from the top left and right frames as shown
124
Frame number: 666 Frame number: 65000
10
Maximum Flow Rate Applied Pressure: 9.12 kPa
1
0.1
0.01
0.001
0.0001
-20 30 80 130 180 230 280 330 380 430 480 530 580 630 680 730
Time (second)
Figure 5.17: The instantaneous flow-rate is plotted against time on a semi-log plot. Time
zero corresponds to the first frame of the first pair where the PIV processing began and the
air phase just enters the porous medium. This is frame 666 shown on top left. The last
framae in this recording is 65000 which is on top right. The whole invasion takes about 715
seconds. During this displacement process, initially the air interface velocity is high and
slows down eventually. The applied pressure is 9.12 kPa.
in Figure 5.17 which are the first and the last frame of the series respectively.
Flow Rate (µL/minute)
125
Velocity plots
Similar to the homogeneous case we present a few representative velocity vector plots
here since our frame rate was not sufficient to resolve the smallest time scale events.
Two velocity plots are presented here as shown in Figures 5.18 and 5.19. Air saturation
significantly increases from the first pair (727-728) to the second pair (16315-16515) which
are 172 seconds apart. It is noticeable how the saturation changed and trapped water is
already visible in Figure 5.18. Similarly, in Figure 5.19 some of the initially trapped water
was displaced indicating again the presence of the corner film flow as the homogeneous case.
Frame 727-728 t=7.997s
mm/s
Figure 5.18: Instantaneous velocity vector field between the frame pars 727-728. The applied
pressure is 9.12 kPa. Maximum velocity is 0.5500 mm/s. White and gray represent air and
grain respectively. Air displaces water from right to left during drainage.
126
Frame 16315-16515 t=179.465s
mm/s
Figure 5.19: Instantaneous velocity vector field between the frame pars 16315-16515. The
applied pressure is 9.12 kPa. Maximum velocity is 0.0038 mm/s. White and gray represent
air and grain respectively. Air displaces water from right to left during drainage.
127
Film flow
mm/s
Frame 12551-12600
Figure 5.20: Film flow connecting the outlet with trapped water downstream. The
superimposed red and green lines (with arrows) show the direction of two films and the
velocity magnitude is averaged over 60 frames (12551-12600) displayed with the color map
around grains. The larger trapped water phase is not included in the PTV analysis, only
the films around the grains are considered.
Wetting films are seen to be present and be more dominant for heterogeneous pore-
structure as well. Although we had some seemingly disconnected films in the homogeneous
micromodel, for a heterogeneous micromodel we noticed that all the films are connected
Outlet
Inlet
128
along with all the trapped water phase. As usual we have processed via PTV frame pairs
12551-12600 as shown in Figure 5.20.
Imbibition
For completeness we present imbibition results again for the heterogeneous case. After
the drainage, imbibition was similarly brought upon by reducing the pressure to 3.91 kPa.
For this imbibition we only present the velocity plot as shown in Figure 5.21.
Frame 5286-5287 t=58.146s
mm/s
Figure 5.21: Instantaneous velocity vector field between the frame pars 5286-5287. The
applied pressure is 3.91 kPa. Maximum velocity is 0.7233 mm/s. White and gray represent
air and grain respectively. Water displaces air from left to right during imbibition.
129
To summarize, film flow was seen to be present and play an important role in pore-
scale displacement in both drainage and imbibition in both the micromodels. Revisiting the
question about phase connectivity and the nature of equilibration, we state that as seen in
the drainage cases, a new configuration of the meniscus is possible by displacing trapped
wetting phases leading to the equilibrium condition. Although we have presented film flow
within the initial 5 minutes of the beginning of the displacement process for the drainage
cases in homogeneous and heterogeneous micromodels respectively, film flow was seen to
be present even after 10 minutes, although relatively slower and almost stationary at some
point. We also want to mention the fact that some trapped water was not displaced even at
the end of the recording for both the homogeneous and heterogeneous drainage experiment.
Also, we verify the fact that wetting layers is always present and trapped wetting phase is
indeed connected to the outlet and can be displaced until final equilibrium state is reached.
Therefore we reinstate the fact that capillary equilibrium and film flow are linked and
any localized imbalance would give rise to film flow. We also verify the fact that pressure
imbalance does exist between the connected and trapped wetting phase driving the film flow.
Besides the importance in regard to capillary pressure and phase connectivity, these films
are important in a general sense as they can significantly contribute to the mass transport
process.
130
CONCLUSION
In this thesis, two different experiments have been presented, both performed in 2D
porous micromodels employing epi-fluorescence microscopy coupled with high-speed imaging.
Our first goal was to find the relation between the traditionally defined capillary pressure and
macroscale capillary pressure averaged from the pore-scale curvature and find the effect of
spatial resolution on interfacial curvature and resulting macroscale capillary pressure along
with other variables of interest such as interfacial area between wetting and nonnwetting
phases and saturation. Effect of phase connectivity was also in mind to study the role it
plays in resulting capillary pressure. The second objective was to study and characterize film
flow which exist in a strongly wetting medium and serves as a pathway for trapped phases
to flow until true equilibrium state is reached.
Our results have shown that measurement resolution (i.e., 2µm and smaller/pixel)
does not seem to have a significant effect on the calculated capillary pressure, saturation
and interfacial length for both homogeneous and heterogeneous cases. The results have also
shown that boundary pressure and macroscale capillary pressure calculated from averaging
pore-scale capillary pressure over the interfaces are two separate measures in their own
domain. One must be careful when using them interchangeably.
Film flow was seen to be an active mode of pore-scale displacement process. During
drainage they seemed to be giving way to trapped wetting phase at a slower rate which is
result of the pressure gradient in trapped phase connected by the films. During imbibition
they were also dominant filling up the narrowest throats.
Additional studies are needed to understand capillary equilibrium and interfacial
relaxation to enrich our understanding. Film flow should be studied more carefully as well
such in 3D porous media.
131
REFERENCES CITED
[1] Lara Adrian, Ronald J Adrian, and Jerry Westerweel. Particle image velocimetry.
Number 30. Cambridge university press, 2011.
[2] Abdurrahman Akyol, Orhan Taner Can, Erhan Demirbas, and Mehmet Kobya. A
comparative study of electrocoagulation and electro-fenton for treatment of wastewater
from liquid organic fertilizer plant. Separation and Purification Technology, 112:11–19,
2013.
[3] Nihat Hakan Akyol and Sevgi Turkkan. Effect of cyclodextrin-enhanced dissolution on
mass removal and mass-flux reduction relationships for non-uniformly organic liquid
distribution in heterogeneous porous media. Water, Air, & Soil Pollution, 229(2):1–11,
2018.
[4] Alimohammad Anbari, Hung-Ta Chien, Sujit S Datta, Wen Deng, David A Weitz,
and Jing Fan. Microfluidic model porous media: Fabrication and applications. Small,
14(18):1703575, 2018.
[5] Ryan T Armstrong and Steffen Berg. Interfacial velocities and capillary pressure
gradients during haines jumps. Physical Review E, 88(4):043010, 2013.
[6] Ryan T Armstrong, Nikolay Evseev, Dmitry Koroteev, and Steffen Berg. Modeling
the velocity field during haines jumps in porous media. Advances in Water Resources,
77:57–68, 2015.
[7] Ryan T Armstrong, Mark L Porter, and Dorthe Wildenschild. Linking pore-scale
interfacial curvature to column-scale capillary pressure. Advances in Water resources,
46:55–62, 2012.
[8] Peter Atkins, Peter William Atkins, and Julio de Paula. Atkins’ physical chemistry.
Oxford university press, 2014.
[9] Rimsha Aziz, Vahid Joekar-Niasar, Pedro J Mart́ınez-Ferrer, Omar E Godinez-
Brizuela, Constantinos Theodoropoulos, and Hassan Mahani. Novel insights into pore-
scale dynamics of wettability alteration during low salinity waterflooding. Scientific
reports, 9(1):1–13, 2019.
[10] Stefan Bachu. Sequestration of co2 in geological media: criteria and approach for
site selection in response to climate change. Energy conversion and management,
41(9):953–970, 2000.
[11] Shazia Bashir, Muhammad Bashir, Julia Margaret Rees, William Bauer Jay Zimmer-
man, et al. Dynamic wetting in microfluidic droplet formation. BioChip Journal,
8(2):122–128, 2014.
132
[12] Jacob Bear. Dynamics of fluids in porous media. Courier Corporation, 1988.
[13] Jacob Bear and Alexander H-D Cheng. Modeling groundwater flow and contaminant
transport, volume 23. Springer, 2010.
[14] Steffen Berg, Holger Ott, Stephan A Klapp, Alex Schwing, Rob Neiteler, Niels Brussee,
Axel Makurat, Leon Leu, Frieder Enzmann, Jens-Oliver Schwarz, et al. Real-time 3d
imaging of haines jumps in porous media flow. Proceedings of the National Academy
of Sciences, 110(10):3755–3759, 2013.
[15] Brian Berkowitz. Boundary conditions along permeable fracture walls: Influence on
flow and conductivity. Water Resources Research, 25(8):1919–1922, 1989.
[16] R Byron Bird. Transport phenomena. Appl. Mech. Rev., 55(1):R1–R4, 2002.
[17] Gianluca Blois, Julio M Barros, and Kenneth T Christensen. A microscopic particle
image velocimetry method for studying the dynamics of immiscible liquid–liquid
interactions in a porous micromodel. Microfluidics and Nanofluidics, 18(5):1391–1406,
2015.
[18] Martin J Blunt. Multiphase flow in permeable media: A pore-scale perspective.
Cambridge university press, 2017.
[19] L Boruvka and AW Neumann. Generalization of the classical theory of capillarity. The
Journal of Chemical Physics, 66(12):5464–5476, 1977.
[20] NP Brandon and DJ Brett. Engineering porous materials for fuel cell applications.
Philosophical Transactions of the Royal Society A: Mathematical, Physical and
Engineering Sciences, 364(1838):147–159, 2006.
[21] Hendrik C Brinkman. A calculation of the viscous force exerted by a flowing fluid on
a dense swarm of particles. Flow, Turbulence and Combustion, 1(1):27–34, 1949.
[22] Royal Harvard Brooks. Hydraulic properties of porous media. Colorado State
University, 1965.
[23] Hans-Jürgen Butt, Karlheinz Graf, and Michael Kappl. Physics and chemistry of
interfaces. John Wiley & Sons, 2013.
[24] C Chalbaud, M Robin, JM Lombard, F Martin, P Egermann, and H Bertin. Interfacial
tension measurements and wettability evaluation for geological co2 storage. Advances
in water resources, 32(1):98–109, 2009.
[25] David Chandler. Introduction to modern statistical. Mechanics. Oxford University
Press, Oxford, UK, 5:449, 1987.
[26] Ioannis Chatzis and Norman R Morrow. Correlation of capillary number relationships
for sandstone. Society of Petroleum Engineers Journal, 24(05):555–562, 1984.
133
[27] Fabien Chauvet, Paul Duru, Sandrine Geoffroy, and Marc Prat. Three periods of
drying of a single square capillary tube. Physical review letters, 103(12):124502, 2009.
[28] Fabien Chauvet, Paul Duru, and Marc Prat. Depinning of evaporating liquid films
in square capillary tubes: Influence of corners’ roundedness. Physics of Fluids,
22(11):112113, 2010.
[29] Daiquan Chen, Laura J Pyrak-Nolte, Jessica Griffin, and Nicholas J Giordano.
Measurement of interfacial area per volume for drainage and imbibition. Water
Resources Research, 43(12), 2007.
[30] Yu Chen, Yaofa Li, Albert J Valocchi, and Kenneth T Christensen. Lattice boltzmann
simulations of liquid co2 displacing water in a 2d heterogeneous micromodel at reservoir
pressure conditions. Journal of contaminant hydrology, 212:14–27, 2018.
[31] Nikolai Chernov. Circle fit (pratt method). Available at https://www.mathworks.
com/matlabcentral/fileexchange/22643-circle-fit-pratt-method.
[32] Nikolai Chernov and Claire Lesort. Least squares fitting of circles. Journal of
Mathematical Imaging and Vision, 23:239–252, 2005.
[33] Daniel A Clarke, Fabian Dolamore, Conan J Fee, Petrik Galvosas, and Daniel J
Holland. Investigation of flow through triply periodic minimal surface-structured
porous media using mri and cfd. Chemical Engineering Science, 231:116264, 2021.
[34] Torsten Clemens, Kostas Tsikouris, Markus Buchgraber, Louis Castanier, and Anthony
Kovscek. Pore-scale evaluation of polymers displacing viscous oil—computational-
fluid-dynamics simulation of micromodel experiments. SPE Reservoir Evaluation &
Engineering, 16(02):144–154, 2013.
[35] George N Constantinides and Alkiviades C Payatakes. Effects of precursor wetting
films in immiscible displacement through porous media. Transport in porous media,
38(3):291–317, 2000.
[36] Christophe Cottin, Hugues Bodiguel, and Annie Colin. Drainage in two-dimensional
porous media: From capillary fingering to viscous flow. Physical Review E,
82(4):046315, 2010.
[37] Henry Darcy. Les fontaines publiques de la ville de Dijon: Exposition et application des
principes à suivre et des formules à employer dans les questions de distribution d’eau:
Ouvrage terminé par un appendice relatif aux fournitures d’eau de plusieurs villes, au
filtrage des eaux et à la fabrication des tuyaux de fonte, de plomb, de tôle et de bitume,
volume 2. V. Dalmont, 1856.
[38] Malay K Das, Partha P Mukherjee, and K Muralidhar. Porous media applications:
biological systems. In Modeling Transport Phenomena in Porous Media with Applica-
tions, pages 123–154. Springer, 2018.
134
[39] Howard Ted Davis. Statistical mechanics of phases, interfaces, and thin films. Wiley-
VCH, 1996.
[40] Roger JM DeWiest. History of the dupuit-forchheimer assumptions on groundwater
hydraulics. Transactions of the ASAE, 8(4):508–0509, 1965.
[41] Erica L DiFilippo and Robert P Eganhouse. Assessment of pdms-water partition
coefficients: implications for passive environmental sampling of hydrophobic organic
compounds. Environmental science & technology, 44(18):6917–6925, 2010.
[42] F Doster and R Hilfer. A comparison between simulation and experiment for hysteretic
phenomena during two-phase immiscible displacement. Water Resources Research,
50(1):681–686, 2014.
[43] P Druetta, P Raffa, and F Picchioni. Chemical enhanced oil recovery and the role of
chemical product design. Applied energy, 252:113480, 2019.
[44] Francis AL Dullien. Porous media: fluid transport and pore structure. Academic press,
2012.
[45] Francis AL Dullien, Cesar Zarcone, Ian F Macdonald, Ann Collins, and Ron DE
Bochard. The effects of surface roughness on the capillary pressure curves and the
heights of capillary rise in glass bead packs. Journal of Colloid and Interface Science,
127(2):362–372, 1989.
[46] Sabri Ergun. Fluid flow through packed columns. Chem. Eng. Prog., 48:89–94, 1952.
[47] Fritjof Fagerlund. Experimental and modelling studies on the spreading of non-aqueous
phase liquids in heterogeneous media. 11 2022.
[48] Philipp Forchheimer. Wasserbewegung durch boden. Z. Ver. Deutsch, Ing., 45:1782–
1788, 1901.
[49] H Friedli, H Lötscher, Hj Oeschger, Urs Siegenthaler, and Bernhard Stauffer. Ice core
record of the 13c/12c ratio of atmospheric co2 in the past two centuries. Nature,
324(6094):237–238, 1986.
[50] Eric M. Furst. Particle tracking with matlab. Available at https://lem.che.udel.
edu/wiki/index.php?n=Main.Microrheology.
[51] Eric M. Furst. Particle tracking with matlab. Available at https://lem.che.udel.
edu/wiki/uploads/Main/Microrheology/Particle_tracking_handout.pdf.
[52] Huhao Gao, Alexandru Bogdan Tatomir, Nikolaos K Karadimitriou, Holger Steeb,
and Martin Sauter. A two-phase, pore-scale reactive transport model for the kinetic
interface-sensitive tracer. Water Resources Research, 57(6):e2020WR028572, 2021.
135
[53] L Gascón, JA Garćıa, F LeBel, E Ruiz, and F Trochu. Numerical prediction
of saturation in dual scale fibrous reinforcements during liquid composite molding.
Composites Part A: Applied Science and Manufacturing, 77:275–284, 2015.
[54] Helmut Geistlinger, Iman Ataei-Dadavi, and Hans-Jörg Vogel. Impact of surface rough-
ness on capillary trapping using 2d-micromodel visualization experiments. Transport
in Porous Media, 112(1):207–227, 2016.
[55] Alireza Gerami, Yara Alzahid, Peyman Mostaghimi, Navid Kashaninejad, Farzan
Kazemifar, Tammy Amirian, Nader Mosavat, Majid Ebrahimi Warkiani, and Ryan T
Armstrong. Microfluidics for porous systems: fabrication, microscopy and applications.
Transport in Porous Media, 130(1):277–304, 2019.
[56] William G Gray, Kelsey Bruning, and Cass T Miller. Non-hysteretic functional form
of capillary pressure in porous media. Journal of Hydraulic Research, 57(6):747–759,
2019.
[57] William G Gray and S Majid Hassanizadeh. Paradoxes and realities in unsaturated
flow theory. Water resources research, 27(8):1847–1854, 1991.
[58] William G Gray and S Majid Hassanizadeh. Unsaturated flow theory including
interfacial phenomena. Water Resources Research, 27(8):1855–1863, 1991.
[59] William G Gray and Cass T Miller. Tcat analysis of capillary pressure in non-
equilibrium, two-fluid-phase, porous medium systems. Advances in water resources,
34(6):770–778, 2011.
[60] William G Gray and Cass T Miller. A generalization of averaging theorems for porous
medium analysis. Advances in water resources, 62:227–237, 2013.
[61] William G Gray and Cass T Miller. Introduction to the thermodynamically constrained
averaging theory for porous medium systems, volume 696. Springer, 2014.
[62] Leon Green Jr and Pol Duwez. Fluid flow through porous metals. 1951.
[63] WD Gunter, B Wiwehar, and EH Perkins. Aquifer disposal of co2-rich greenhouse
gases: extension of the time scale of experiment for co2-sequestering reactions by
geochemical modelling. Mineralogy and petrology, 59(1):121–140, 1997.
[64] Majid Hassanizadeh and William G Gray. General conservation equations for multi-
phase systems: 1. averaging procedure. Advances in water resources, 2:131–144, 1979.
[65] S Majid Hassanizadeh and William G Gray. High velocity flow in porous media.
Transport in porous media, 2(6):521–531, 1987.
[66] S Majid Hassanizadeh and William G Gray. Mechanics and thermodynamics of
multiphase flow in porous media including interphase boundaries. Advances in water
resources, 13(4):169–186, 1990.
136
[67] S Majid Hassanizadeh and William G Gray. Thermodynamic basis of capillary pressure
in porous media. Water resources research, 29(10):3389–3405, 1993.
[68] S Majid Hassanizadeh and William G Gray. Toward an improved description of the
physics of two-phase flow. Advances in Water Resources, 16(1):53–67, 1993.
[69] Feng He, Zhijun Wang, Lilin Wang, Junjie Li, and Jincheng Wang. Effects of surfactant
on capillary evaporation process with thick films. International Journal of Heat and
Mass Transfer, 88:406–410, 2015.
[70] Anna L Herring, FJ Gilby, Zhe Li, JE McClure, Michael Turner, JP Veldkamp,
L Beeching, and AP Sheppard. Observations of nonwetting phase snap-off during
drainage. Advances in Water Resources, 121:32–43, 2018.
[71] Anna L Herring, Jill Middleton, Rick Walsh, Andrew Kingston, and Adrian Sheppard.
Flow rate impacts on capillary pressure and interface curvature of connected and
disconnected fluid phases during multiphase flow in sandstone. Advances in Water
Resources, 107:460–469, 2017.
[72] R Hilfer. Capillary pressure, hysteresis and residual saturation in porous media.
Physica A: Statistical Mechanics and its Applications, 359:119–128, 2006.
[73] R Hilfer. Macroscopic capillarity and hysteresis for flow in porous media. Physical
Review E, 73(1):016307, 2006.
[74] R Hilfer, RT Armstrong, S Berg, A Georgiadis, and H Ott. Capillary saturation and
desaturation. Physical Review E, 92(6):063023, 2015.
[75] R Hilfer and F Doster. Percolation as a basic concept for macroscopic capillarity.
Transport in porous media, 82(3):507–519, 2010.
[76] Johannes Hommel, Edward Coltman, and Holger Class. Porosity–permeability
relations for evolving pore space: a review with a focus on (bio-) geochemically altered
porous media. Transport in Porous Media, 124(2):589–629, 2018.
[77] Frouke Hoogland, Peter Lehmann, and Dani Or. Drainage dynamics controlled by
corner flow: Application of the foam drainage equation. Water Resources Research,
52(11):8402–8412, 2016.
[78] Ran Hu, Jiamin Wan, Zhibing Yang, Yi-Feng Chen, and Tetsu Tokunaga. Wettability
and flow rate impacts on immiscible displacement: A theoretical model. Geophysical
Research Letters, 45(7):3077–3086, 2018.
[79] M King Hubbert. Darcy’s law and the field equations of the flow of underground fluids.
Transactions of the AIME, 207(01):222–239, 1956.
137
[80] Kh Kh Imomnazarov. Modified darcy laws for conducting porous media. Mathematical
and computer modelling, 40(1-2):5–10, 2004.
[81] IPCC. Global warming of 1.5 ºc. Available at https://www.ipcc.ch/sr15/, July
2022.
[82] GR Jerauld and SJ Salter. The effect of pore-structure on hysteresis in relative
permeability and capillary pressure: pore-level modeling. Transport in porous media,
5(2):103–151, 1990.
[83] Hao Jiang, Bo Guo, and Mark L Brusseau. Pore-scale modeling of fluid-fluid interfacial
area in variably saturated porous media containing microscale surface roughness.
Water resources research, 56(1):e2019WR025876, 2020.
[84] V Joekar-Niasar and SM Hassanizadeh. Specific interfacial area: The missing state
variable in two-phase flow equations? Water Resources Research, 47(5), 2011.
[85] V Joekar-Niasar, SM Hassanizadeh, and A Leijnse. Insights into the relationships
among capillary pressure, saturation, interfacial area and relative permeability using
pore-network modeling. Transport in porous media, 74(2):201–219, 2008.
[86] Vahid Joekar-Niasar and S Majid Hassanizadeh. Effect of fluids properties on non-
equilibrium capillarity effects: Dynamic pore-network modeling. International Journal
of Multiphase Flow, 37(2):198–214, 2011.
[87] Fortunat Joos and Renato Spahni. Rates of change in natural and anthropogenic
radiative forcing over the past 20,000 years. Proceedings of the National Academy of
Sciences, 105(5):1425–1430, 2008.
[88] François Kalaydjian. A macroscopic description of multiphase flow in porous media
involving spacetime evolution of fluid/fluid interface. Transport in Porous Media,
2(6):537–552, 1987.
[89] NK Karadimitriou and SM Hassanizadeh. A review of micromodels and their use in
two-phase flow studies. Vadose Zone Journal, 11(3):vzj2011–0072, 2012.
[90] NK Karadimitriou, SM Hassanizadeh, V Joekar-Niasar, and PJ Kleingeld. Micromodel
study of two-phase flow under transient conditions: Quantifying effects of specific
interfacial area. Water Resources Research, 50(10):8125–8140, 2014.
[91] NK Karadimitriou, V Joekar-Niasar, SM Hassanizadeh, PJ Kleingeld, and LJ Pyrak-
Nolte. A novel deep reactive ion etched (drie) glass micro-model for two-phase flow
experiments. Lab on a Chip, 12(18):3413–3418, 2012.
[92] NK Karadimitriou, Michiel Musterd, PJ Kleingeld, MT Kreutzer, SM Hassanizadeh,
and V Joekar-Niasar. On the fabrication of pdms micromodels by rapid prototyping,
and their use in two-phase flow studies. Water Resources Research, 49(4):2056–2067,
2013.
138
[93] Farzan Kazemifar, Gianluca Blois, Dimitrios C Kyritsis, and Kenneth T Christensen.
A methodology for velocity field measurement in multiphase high-pressure flow of co2
and water in micromodels. Water Resources Research, 51(4):3017–3029, 2015.
[94] Farzan Kazemifar, Gianluca Blois, Dimitrios C Kyritsis, and Kenneth T Christensen.
Quantifying the flow dynamics of supercritical co2–water displacement in a 2d porous
micromodel using fluorescent microscopy and microscopic piv. Advances in Water
Resources, 95:352–368, 2016.
[95] Emmanuel Keita, Stephan A Koehler, Paméla Faure, David A Weitz, and Philippe
Coussot. Drying kinetics driven by the shape of the air/water interface in a capillary
channel. The European Physical Journal E, 39(2):1–10, 2016.
[96] H Khorshidian, LA James, and SD Butt. The role of film flow and wettability in
immiscible gas assisted gravity drainage. In Proceedings of International Symposium
of the Society of Core Analysts held in Snowmass, CO, USA, pages 1–16, 2016.
[97] JE Killough. Reservoir simulation with history-dependent saturation functions. Society
of Petroleum Engineers Journal, 16(01):37–48, 1976.
[98] H Koide, Y Tazaki, Y Noguchi, S Nakayama, M Iijima, K Ito, and Y Shindo.
Subterranean containment and long-term storage of carbon dioxide in unused aquifers
and in depleted natural gas reservoirs. Energy Conversion and management, 33(5-
8):619–626, 1992.
[99] Joel Koplik, Herbert Levine, and A Zee. Viscosity renormalization in the brinkman
equation. The Physics of fluids, 26(10):2864–2870, 1983.
[100] Jacob M LaManna, James V Bothe Jr, Feng Yuan Zhang, and Matthew M Mench.
Measurement of capillary pressure in fuel cell diffusion media, micro-porous layers,
catalyst layers, and interfaces. Journal of Power Sources, 271:180–186, 2014.
[101] SL Lee and JH Yang. Modeling of darcy-forchheimer drag for fluid flow across a bank of
circular cylinders. International journal of heat and mass transfer, 40(13):3149–3155,
1997.
[102] Roland Lenormand. Liquids in porous media. Journal of Physics: Condensed Matter,
2(S):SA79, 1990.
[103] Roland Lenormand, Eric Touboul, and Cesar Zarcone. Numerical models and
experiments on immiscible displacements in porous media. Journal of fluid mechanics,
189:165–187, 1988.
[104] Roland Lenormand and Cesar Zarcone. Role of roughness and edges during imbibition
in square capillaries. In SPE annual technical conference and exhibition. OnePetro,
1984.
139
[105] MoC Leverett. Capillary behavior in porous solids. Transactions of the AIME,
142(01):152–169, 1941.
[106] Tianyi Li, Steffen Schlüter, Maria Ines Dragila, and Dorthe Wildenschild. An improved
method for estimating capillary pressure from 3d microtomography images and its
application to the study of disconnected nonwetting phase. Advances in Water
Resources, 114:249–260, 2018.
[107] Yaofa Li, Gianluca Blois, Farzan Kazemifar, and Kenneth T Christensen. High-
speed quantification of pore-scale multiphase flow of water and supercritical co2 in
2-d heterogeneous porous micromodels: flow regimes and interface dynamics. Water
Resources Research, 55(5):3758–3779, 2019.
[108] Yaofa Li, Gianluca Blois, Farzan Kazemifar, and Kenneth T Christensen. A particle-
based image segmentation method for phase separation and interface detection in
piv images of immiscible multiphase flow. Measurement Science and Technology,
32(9):095208, 2021.
[109] Ralph Lindken, Massimiliano Rossi, Sebastian Große, and Jerry Westerweel. Micro-
particle image velocimetry (µpiv): recent developments, applications, and guidelines.
Lab on a Chip, 9(17):2551–2567, 2009.
[110] Rebecca Lindsey. Climate change: Atmospheric carbon dioxide. Available
at https://www.climate.gov/news-features/understanding-climate/
climate-change-atmospheric-carbon-dioxide.
[111] Grunde Løvoll, Yves Méheust, Knut Jørgen Måløy, Eyvind Aker, and Jean Schmit-
tbuhl. Competition of gravity, capillary and viscous forces during drainage in a two-
dimensional porous medium, a pore scale study. Energy, 30(6):861–872, 2005.
[112] L Ma, DB Ingham, and MC Pourkashanian. Application of fluid flows through porous
media in fuel cells. In Transport phenomena in porous media III, pages 418–440.
Elsevier, 2005.
[113] IF Macdonald, MS El-Sayed, K Mow, and FAL Dullien. Flow through porous media-
the ergun equation revisited. Industrial & Engineering Chemistry Fundamentals,
18(3):199–208, 1979.
[114] Walid Mohamed Mahmud and Viet Hoai Nguyen. Effects of snap-off in imbibition in
porous media with different spatial correlations. Transport in porous media, 64(3):279–
300, 2006.
[115] CM Marle. On macroscopic equations governing multiphase flow with diffusion and
chemical reactions in porous media. International Journal of Engineering Science,
20(5):643–662, 1982.
140
[116] James E McClure, Ryan T Armstrong, Mark A Berrill, Steffen Schlüter, Steffen Berg,
William G Gray, and Cass T Miller. Geometric state function for two-fluid flow in
porous media. Physical Review Fluids, 3(8):084306, 2018.
[117] Yves Méheust, Grunde Løvoll, Knut Jørgen Måløy, and Jean Schmittbuhl. Interface
scaling in a two-dimensional porous medium under combined viscous, gravity, and
capillary effects. Physical Review E, 66(5):051603, 2002.
[118] Bert Metz, Ogunlade Davidson, HC De Coninck, Manuela Loos, and Leo Meyer. IPCC
special report on carbon dioxide capture and storage. Cambridge: Cambridge University
Press, 2005.
[119] Cass T Miller, George Christakos, Paul T Imhoff, John F McBride, Joseph A Pedit,
and John A Trangenstein. Multiphase flow and transport modeling in heterogeneous
porous media: challenges and approaches. Advances in Water Resources, 21(2):77–120,
1998.
[120] CT Miller, K Bruning, CL Talbot, JE McClure, and WG Gray. Nonhysteretic capillary
pressure in two-fluid porous medium systems: Definition, evaluation, validation, and
dynamics. Water Resources Research, 55(8):6825–6849, 2019.
[121] Franziska Moebius and Dani Or. Interfacial jumps and pressure bursts during fluid
displacement in interacting irregular capillaries. Journal of colloid and interface
science, 377(1):406–415, 2012.
[122] Franziska Moebius and Dani Or. Inertial forces affect fluid front displacement dynamics
in a pore-throat network model. Physical Review E, 90(2):023019, 2014.
[123] Kishore K Mohanty, H Ted Davis, and LE Scriven. Physics of oil entrapment in
water-wet rock. SPE Reservoir Engineering, 2(01):113–128, 1987.
[124] Norman R Morrow. Wettability and its effect on oil recovery. Journal of petroleum
technology, 42(12):1476–1484, 1990.
[125] Norman R Morrow, I Chatzis, and JJ Taber. Entrapment and mobilization of residual
oil in bead packs. SPE Reservoir Engineering, 3(03):927–934, 1988.
[126] Mehdi Mosharaf-Dehkordi. A fully coupled porous media and channels flow approach
for simulation of blood and bile flow through the liver lobules. Computer Methods in
Biomechanics and Biomedical Engineering, 22(9):901–915, 2019.
[127] Klaus Mosthaf. Modeling and analysis of coupled porous-medium and free flow with
application to evaporation processes. 2014.
[128] Marcel Moura, Eirik Grude Flekkøy, Knut Jørgen Måløy, Gerhard Schäfer, and Renaud
Toussaint. Connectivity enhancement due to film flow in porous media. Physical Review
Fluids, 4(9):094102, 2019.
141
[129] Julia C Muccino, William G Gray, and Lin A Ferrand. Toward an improved
understanding of multiphase flow in porous media. Reviews of Geophysics, 36(3):401–
422, 1998.
[130] M Muskat, RD Wyckoff, HG Botset, and MW Meres. Flow of gas-liquid mixtures
through sands. Transactions of the AIME, 123(01):69–96, 1937.
[131] Morris Muskat and Milan W Meres. The flow of heterogeneous fluids through porous
media. Physics, 7(9):346–363, 1936.
[132] Viet Hoai Nguyen, Adrian P Sheppard, Mark A Knackstedt, and W Val Pinczewski.
The effect of displacement rate on imbibition relative permeability and residual
saturation. Journal of Petroleum Science and Engineering, 52(1-4):54–70, 2006.
[133] Donald A Nield, Adrian Bejan, et al. Convection in porous media, volume 3. Springer,
2006.
[134] Jennifer Niessner, Steffen Berg, and S Majid Hassanizadeh. Comparison of two-phase
darcy’s law with a thermodynamically consistent approach. Transport in porous media,
88(1):133–148, 2011.
[135] SO Ochs, H Class, A Färber, and R Helmig. Methods for predicting the spreading of
steam below the water table during subsurface remediation. Water Resources Research,
46(5), 2010.
[136] Dlugokencky E Trans P. Trends in atmospheric carbon dioxide (noaa/esrl). Available
at https://gml.noaa.gov/ccgg/trends/, July 2022.
[137] Prapti Pattanayak, Sachin Kumar Singh, Monica Gulati, Sukriti Vishwas, Bhupinder
Kapoor, Dinesh Kumar Chellappan, Krishnan Anand, Gaurav Gupta, Niraj Kumar
Jha, Piyush Kumar Gupta, et al. Microfluidic chips: recent advances, critical strategies
in design, applications and future perspectives. Microfluidics and nanofluidics,
25(12):1–28, 2021.
[138] Kurt D Pennell, Gary A Pope, and Linda M Abriola. Influence of viscous and buoyancy
forces on the mobilization of residual tetrachloroethylene during surfactant flushing.
Environmental Science & Technology, 30(4):1328–1335, 1996.
[139] George F Pinder and William G Gray. Essentials of multiphase flow and transport in
porous media. John Wiley & Sons, 2008.
[140] GA Pope and RC Nelson. A chemical flooding compositional simulator. Society of
Petroleum Engineers Journal, 18(05):339–354, 1978.
[141] Laura J Pyrak-Nolte, David D Nolte, Daiquan Chen, and Nicholas J Giordano.
Relating capillary pressure to interfacial areas. Water Resources Research, 44(6), 2008.
142
[142] Zhipeng Qin, Soheil Esmaeilzadeh, Amir Riaz, and Hamdi A Tchelepi. Two-phase
multiscale numerical framework for modeling thin films on curved solid surfaces in
porous media. Journal of Computational Physics, 413:109464, 2020.
[143] Arash Rabbani, Saeid Jamshidi, and Saeed Salehi. An automated simple algorithm for
realistic pore network extraction from micro-tomography images. Journal of Petroleum
Science and Engineering, 123:164–171, 2014.
[144] Harris Sajjad Rabbani, Dani Or, Ying Liu, Ching-Yao Lai, Nancy B Lu, Sujit S Datta,
Howard A Stone, and Nima Shokri. Suppressing viscous fingering in structured porous
media. Proceedings of the National Academy of Sciences, 115(19):4833–4838, 2018.
[145] Markus Raffel, Christian E Willert, Jürgen Kompenhans, et al. Particle image
velocimetry: a practical guide, volume 2. Springer, 1998.
[146] Saman Rashidi, Faramarz Hormozi, and Mohammad Hossein Doranehgard. Abilities
of porous materials for energy saving in advanced thermal systems. Journal of thermal
analysis and calorimetry, 143(3):2437–2452, 2021.
[147] Masoud Riazi, Mehran Sohrabi, Christian Bernstone, Mahmoud Jamiolahmady, and
Shaun Ireland. Visualisation of mechanisms involved in co2 injection and storage in
hydrocarbon reservoirsand water-bearing aquifers. Chemical Engineering Research and
Design, 89(9):1827–1840, 2011.
[148] Sophie Roman, Cyprien Soulaine, Moataz Abu AlSaud, Anthony Kovscek, and Hamdi
Tchelepi. Particle velocimetry analysis of immiscible two-phase flow in micromodels.
Advances in Water Resources, 95:199–211, 2016.
[149] JG Roof. Snap-off of oil droplets in water-wet pores. Society of Petroleum Engineers
Journal, 10(01):85–90, 1970.
[150] Subhadeep Roy, Santanu Sinha, and Alex Hansen. Flow-area relations in immiscible
two-phase flow in porous media. Frontiers in Physics, 8:4, 2020.
[151] Juan G Santiago, Steve T Wereley, Carl D Meinhart, DJ Beebe, and Ronald J Adrian.
A particle image velocimetry system for microfluidics. Experiments in fluids, 25(4):316–
319, 1998.
[152] Steffen Schlüter, S Berg, M Rücker, RT Armstrong, H-J Vogel, R Hilfer, and
D Wildenschild. Pore-scale displacement mechanisms as a source of hysteresis for
two-phase flow in porous media. Water Resources Research, 52(3):2194–2205, 2016.
[153] G Sciumè, SE Shelton, WG Gray, CT Miller, Fazle Hussain, Mauro Ferrari, P Decuzzi,
and BA Schrefler. Tumor growth modeling from the perspective of multiphase porous
media mechanics. Molecular & cellular biomechanics: MCB, 9(3):193, 2012.
143
[154] Vahid Shabro, Carlos Torres-Verdin, and Farzam Javadpour. Numerical simulation
of shale-gas production: from pore-scale modeling of slip-flow, knudsen diffusion, and
langmuir desorption to reservoir modeling of compressible fluid. In North American
Unconventional Gas Conference and Exhibition. OnePetro, 2011.
[155] Ebrahim Shahraeeni and Dani Or. Pore-scale evaporation-condensation dynamics
resolved by synchrotron x-ray tomography. Physical Review E, 85(1):016317, 2012.
[156] Nima Shokri, Peter Lehmann, P Vontobel, and Dani Or. Drying front and water
content dynamics during evaporation from sand delineated by neutron radiography.
Water resources research, 44(6), 2008.
[157] Kamaljit Singh, Hannah Menke, Matthew Andrew, Qingyang Lin, Christoph Rau,
Martin J Blunt, and Branko Bijeljic. Dynamics of snap-off and pore-filling events
during two-phase fluid flow in permeable media. Scientific reports, 7(1):1–13, 2017.
[158] M Soltani and Pu Chen. Numerical modeling of fluid flow in solid tumors. PloS one,
6(6):e20344, 2011.
[159] HL Stone. Probability model for estimating three-phase relative permeability. Journal
of petroleum technology, 22(02):214–218, 1970.
[160] Zhonghao Sun, Ayaz Mehmani, and Carlos Torres-Verd́ın. Subpore-scale trapping
mechanisms following imbibition: A microfluidics investigation of surface roughness
effects. Water Resources Research, 57(2):e2020WR028324, 2021.
[161] Ankur Taneja and Jonathan Higdon. A fully-coupled discontinuous galerkin spectral
element method for two-phase flow in petroleum reservoirs. Journal of Computational
Physics, 352:341–372, 2018.
[162] Jefferson W Tester, Michael Modell, et al. Thermodynamics and its Applications.
Prentice Hall PTR, 1997.
[163] Chia-Wen Tsao, Qun-Zhan Huang, Chang-Ye You, Markus Hilpert, Shao-Yiu Hsu,
Krzysztof Lamorski, Liang-Cheng Chang, and Cezary S lawiński. The effect of channel
aspect ratio on air entrapment during imbibition in soil-on-a-chip micromodels with
2d and 2.5 d pore structures. Lab on a Chip, 21(2):385–396, 2021.
[164] Markus Tuller and Dani Or. Hydraulic conductivity of variably saturated porous
media: Film and corner flow in angular pore space. Water resources research,
37(5):1257–1276, 2001.
[165] M Th Van Genuchten. A closed-form equation for predicting the hydraulic conductivity
of unsaturated soils. Soil science society of America journal, 44(5):892–898, 1980.
[166] Stephen Whitaker. Advances in theory of fluid motion in porous media. Industrial &
engineering chemistry, 61(12):14–28, 1969.
144
[167] Dorthe Wildenschild and Adrian P Sheppard. X-ray imaging and analysis techniques
for quantifying pore-scale structure and processes in subsurface porous medium
systems. Advances in Water resources, 51:217–246, 2013.
[168] Rui Wu, Tao Zhang, Chao Ye, CY Zhao, Evangelos Tsotsas, and Abdolreza
Kharaghani. Pore network model of evaporation in porous media with continuous
and discontinuous corner films. Physical Review Fluids, 5(1):014307, 2020.
[169] Xiaofan Yang, Yashar Mehmani, William A Perkins, Andrea Pasquali, Martin
Schönherr, Kyungjoo Kim, Mauro Perego, Michael L Parks, Nathaniel Trask,
Matthew T Balhoff, et al. Intercomparison of 3d pore-scale flow and solute transport
simulation methods. Advances in water resources, 95:176–189, 2016.
[170] AG Yiotis, AG Boudouvis, AK Stubos, IN Tsimpanogiannis, and YC Yortsos. Effect
of liquid films on the drying of porous media. AIChE Journal, 50(11):2721–2737, 2004.
[171] Jiseon You, Richard J Preen, Larry Bull, John Greenman, and Ioannis Ieropoulos.
3d printed components of microbial fuel cells: Towards monolithic microbial fuel cell
fabrication using additive layer manufacturing. Sustainable Energy Technologies and
Assessments, 19:94–101, 2017.
[172] Li Yu and Norman Claude Wardlaw. The influence of wettability and critical pore-
throat size ratio on snap—off. Journal of Colloid and Interface Science, 109(2):461–
472, 1986.
[173] Li Yu and Norman Claude Wardlaw. Mechanisms of nonwetting phase trapping during
imbibition at slow rates. Journal of colloid and interface science, 109(2):473–486, 1986.
[174] Xiaoxiong Zha, Min Yu, Jianqiao Ye, and Ganlin Feng. Numerical modeling of
supercritical carbonation process in cement-based materials. Cement and concrete
research, 72:10–20, 2015.
[175] Changyong Zhang, Mart Oostrom, Jay W Grate, Thomas W Wietsma, and Marvin G
Warner. Liquid co2 displacement of water in a dual-permeability pore network
micromodel. Environmental science & technology, 45(17):7581–7588, 2011.
[176] Changyong Zhang, Mart Oostrom, Thomas W Wietsma, Jay W Grate, and Marvin G
Warner. Influence of viscous and capillary forces on immiscible fluid displacement:
Pore-scale experimental study in a water-wet micromodel demonstrating viscous and
capillary fingering. Energy & Fuels, 25(8):3493–3505, 2011.
[177] Dongxiao Zhang, Raoyang Zhang, Shiyi Chen, and Wendy E Soll. Pore scale study of
flow in porous media: Scale dependency, rev, and statistical rev. Geophysical research
letters, 27(8):1195–1198, 2000.
145
[178] Benzhong Zhao, Christopher W MacMinn, and Ruben Juanes. Wettability control on
multiphase flow in patterned microfluidics. Proceedings of the National Academy of
Sciences, 113(37):10251–10256, 2016.
[179] Jianlin Zhao, Feifei Qin, Robert Fischer, Qinjun Kang, Dominique Derome, and Jan
Carmeliet. Spontaneous imbibition in a square tube with corner films: theoretical
model and numerical simulation. Water Resources Research, 57(2):e2020WR029190,
2021.
[180] Jianlin Zhao, Feifei Qin, Qinjun Kang, Chaozhong Qin, Dominique Derome, and Jan
Carmeliet. A dynamic pore network model for imbibition simulation considering corner
film flow. Water Resources Research, page e2022WR032332, 2021.
[181] Xiang Zhou, Qingwang Yuan, Yizhong Zhang, Hanyi Wang, Fanhua Zeng, and Liehui
Zhang. Performance evaluation of co2 flooding process in tight oil reservoir via
experimental and numerical simulation studies. Fuel, 236:730–746, 2019.
[182] L Zhuang, CR Bezerra Coelho, SM Hassanizadeh, and M Th van Genuchten. Analysis
of the hysteretic hydraulic properties of unsaturated soil. Vadose Zone Journal,
16(5):1–9, 2017.
[183] Luwen Zhuang, S Majid Hassanizadeh, Chao-Zhong Qin, and Arjen de Waal.
Experimental investigation of hysteretic dynamic capillarity effect in unsaturated flow.
Water Resources Research, 53(11):9078–9088, 2017.