University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

Efficient Simulation of Flow and Transport in Complex Images of Porous
Materials and Media Using Curvelet Transformation
by
Abdullah Aljasmi

A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(CHEMICAL ENGINEERING)
August 2021

Copyright 2021 Abdullah Aljasmi

ii

Dedication

To my Mother, Abdulaziz, Reem, Rana and Muhammad.

iii

Acknowledgements

“In the name of God, most Gracious, most Merciful ”

I am very grateful to several people for helping and supporting me during
my graduate study period in the Chemical Engineering Department at the
University of Southern California. Those people made my journey towards
completing my PhD dissertation a precious experience.
First, I would like to express my appreciation to my advisor, Professor
Muhammad Sahimi, to whom I owe most of my accomplishments and
learning in the field of transport phenomena in porous materials. Since the
first day of my journey in graduate studies, he has kindly provided me with
valuable guidance, advice and support. He taught me how to become a
patient and hardworking person. I am really privileged to have an inspiring
advisor and a great friend.
I would also like to thank my dissertation committee members, Professor
Kristian Jessen from Chemical Engineering Department and Professor
Felipe de Barros from Civil Engineering Department at USC for their time
and valuable suggestions.

iv
I am very grateful to the Chemical Engineering Department at the
College of Technological Studies of the Public Authority for Applied
Education and Training of Kuwait for their sponsorship and support during
my graduate studies.
I would like to express my special thanks and gratitude to my family,
starting with my great mother, thank you for everything you did for me
throughout my life. This dissertation is dedicated to you, please accept it as
a simple reward.
My dear children, Abdulaziz, Reem, Rana and Muhammad, thank you for
giving my life a powerful sense of purpose. You are the most important piece
of my life, and I love you all.
My deepest appreciation and special thanks go to Dr. Bader Albusairi of
Chemical Engineering Department at Kuwait University, to whom I am
immensely grateful for inspiring me to complete my graduate studies. Since
the beginning of my journey in this field, he provided me with guidance,
insights and endless support.

v

Table of Contents
Dedication ii
Acknowledgements iii
List of Tables vii
List of Figures viii
Abstract xii
Chapter 1: Introduction 1
1.1 Background and Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Curvelet Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.1 Continuous-Time Curvelet Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2.2 Three-Dimensional Curvelet Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.2.3 Curvelet Transform Features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.2.4 Recent Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.2.5 CurveLab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.3 Organization of This Dissertation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

Chapter 2: Efficient Image-Based Simulation of Flow and Transport in Hetero-
geneous Materials and Media: Application of Curvelet Transforms 21
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2 Thresholding an Image of a Porous Media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3 Simulation of the Transport Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4 Computational Efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.5 Extensions of the Methodology. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.6 Summary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

Chapter 3: Speeding-up Simulation of Multiphase Flow in Digital Images of Het-
erogeneous Porous Media by Curvelet Transformation 36
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2 Spatial Correlations and Removal of Redundancy from Images of Porous Media 39
3.3 The Governing Equations and Details of Simulation of Two-Phase Flow . . . . . . 41
3.3.1 The porous medium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

vi
3.3.2 The governing equations and the numerical simulation procedure . . . . . . . 42
3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.5.1 Accuracy and Speed-up of the Computations with the Denoised Images . . 53
3.5.2 Computations in the Curvelet Space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.5.3 Simulation in Three-Dimensional Images . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

Chapter 4: Fast simulation of two-phase flow in three dimensional digital images
of heterogeneous porous media using multiresolution curvelet transforma-
tion 59
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.2 Denoising of 3D Images of Porous Media. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.2.1 Two-dimensional images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.2.2 Three-dimensional images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.3 Simulations in the curvelet or real space? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.4 Details of the Simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.4.1 The porous medium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.4.2 The governing equations and the numerical simulation procedure . . . . . . . 67
4.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.5.1 Single phase flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.5.2 Two phase flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.6 Speed-up of the computations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

Chapter 5: Highly efficient image-based simulation of two-phase flow in porous
media with lattice-Boltzmann method using three-dimensional curvelet
transforms 80
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.2 The porous medium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
5.3 Lattice-Boltzmann Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.4 Processing of the Three-Dimensional Image . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
5.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

Bibliography 105
Appendix 119

vii

List of Tables

2.1 Comparison of the number of grid blocks N and the computation times
for various images. De is the effective diffusivity . . . . . . . . . . . . . . . . .

32
2.2 Comparison of the number of grid blocks N and the computation times
for various images. Ke is the effective diffusivity . . . . . . . . . . . . . . . . .

33
3.1 Comparison of the number of grid blocks N, the computation times,
and the oil residual saturation Sor, computed with the original and
curvelet transform-thresholded images. . . . . . . . . . . . . . . . . . . . . . . . .

54
4.1 Comparison of the number of grid blocks N, the computation times,
and the brine residual saturation Sbr, computed with the original and
curvelet transform-thresholded images. The CPU times represent
those for simulating the two-phase flow. . . . . . . . . . . . . . . . . . . . . . . .

77
5.1 Comparison of the physical properties of three porous media before
and after they are processed by curvelet transformation and their
dependence on the threshold ε. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92
5.2 Comparison of the number of the lattice nodes M, the computation
times (in CPU seconds), the speed-up factor S of the computations,
and the brine residual saturation Sr, computed with the original and
thresholded images. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

102

viii

List of Figures
1.1 Illustration of Curvelet thresholding process of an image. . . . . . . . . 3
1.2 Wavelet coefficients (a), and curvelets on various scale, directions
and positions (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6
1.3 Curvelets at increasingly fine scales. Spatial domain (a-c).
Frequency domain (d-f). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
1.4 Curve representation comparison between (a) wavelet and (b)
curvelet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15
1.5 Original image (a), Noisy image (b) and curvelet denoising image
(c). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
2.1 Comparison of the images of two porous media with their curvelet-
transformed versions. (a), (b), and (c) show, respectively, the
sandstone, a small zoomed-in portion of it, and its curvelet-
transformed with the threshold ε = 0.8. (d), (e), and (f) show the
same, but for the carbonate medium. . . . . . . . . . . . . . . . . . . . . . . . . .

24

2.2 Transient concentration profiles in the curvelet-transformed
images of the sandstone with two thresholds ε, and their comparison
with those in the original image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26
2.3 Transient concentration profiles in the curvelet-transformed
images of the carbonate porous medium with two thresholds ε, and
their comparison with those in the original image. . . . . . . . . . . . . . .

27
2.4 Steady-state concentration profiles in the curvelet-transformed
images of the sandstone with two thresholds ε, and their comparison
with those in the original image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

ix
2.5 Steady-state concentration profiles in the curvelet-transformed
images of the carbonate porous medium with two thresholds ε, and
their comparison with those in the original image. . . . . . . . . . . . . . .

29
3.1 Sandstone image (a), and Carbonate image (b). . . . . . . . . . . . . . . . . . 42
3.2 Absolute permeability vs N (number of grid points), (a) sandstone
and (b) carbonate image, in order to identify the size of the REV. . . . 46
3.3 Time evolution of wetting phase distribution for the (a-c) Sandstone
original image, (d-f) curvelet threshold ε = 0.5 image and (g-i)
curvelet threshold ε = 0.9 image at time 30 (ms), 60 (ms) and at
breakthrough time, respectively. The invading fluid (water) and the
defending fluid (oil) shown in red and blue, respectively. . . . . . . . . . .

47
3.4 Time evolution of wetting phase distribution for the (a-c) Carbonate
original image, (d-f) curvelet threshold ε = 0.5 image and (g-i)
curvelet threshold ε = 0.9 image at time 30 (ms), 60 (ms) and at
breakthrough time, respectively. The invading fluid (water) and the
defending fluid (oil) shown in red and blue, respectively. . . . . . . . . . .

48
3.5 Saturation profile of the wetting fluid along the flow direction (x
direction) at breakthrough time in the curvelet-transformed images
of the (a) Sandstone and (b) carbonate porous medium with two
thresholds ε, and their comparison with those in the original image.
X is the normalized front distance from the inlet (X = x/L, where L
is the length of the sample in the x direction . . . . . . . . . . . . . . . . . . . .

49
3.6 Saturation of the wetting fluid as a function of dimensionless time
,tD (normalized by the breakthrough time) in the curvelet-
transformed images of the (a) sandstone and (b) Carbonate porous
medium with two thresholds ε, and their comparison with those in
the original image. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50
3.7 Probability density function (PDF) of the normalized velocity
(u/Umax, where Umax is the maximum flow velocity) in the curvelet-
transformed images of the Sandstone porous medium with two
thresholds ε, and their comparison with those in the original image.

51

x
3.8 Probability density function (PDF) of the normalized velocity
(u/Umax, where Umax is the maximum flow velocity) in the curvelet-
transformed images of the Carbonate porous medium with two
thresholds ε, and their comparison with those in the original image.

51
3.9 Water-oil relative permeability curves obtained from Sandstone
original image (continuous curves) and curvelet thresholded images
(dashed curves). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52
3.10 Water-oil relative permeability curves obtained from Carbonate
original image (continuous curves) and curvelet thresholded images
(dashed curves). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52
3.11 Small portion of the (a) fine grid, (b) the coarsened grid using ε = 0.5
and (c) the coarsened grid using ε = 0.9. . . . . . . . . . . . . . . . . . . . . . . . 55
4.1 A 200
3
Berea sandstone sample, (a) original and (b) curvelet
transformed sample with ε = 0.8 used for single and two-phase flow
simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66
4.2 Sample-size dependence of the absolute permeability of the Berea
sandstone image in order to identify the size of the REV. . . . . . . . . . 71
4.3 The pressure distribution in the pore volume at the end of the
simulation for (a) original, (b) curvelet threshold image ε = 0.5 and
(c) curvelet threshold image ε = 0.9. . . . . . . . . . . . . . . . . . . . . . . . . . . .

71
4.4 Time evolution of non-wetting phase distribution for the Berea
sandstone (a-b) original image, (c-d) curvelet threshold ε = 0.5 image
and (e-f) curvelet threshold ε = 0.9 image at time 60 (ms) and at
breakthrough time, respectively. CO2 shown in red. . . . . . . . . . . . . . .

73
4.5 Saturation profile of the non-wetting fluid along the flow direction (x
direction) at breakthrough time in the curvelet-transformed images
of the Berea sandstone porous medium with two thresholds ε, and
their comparison with those in the original image. X is the
normalized front distance from the inlet (X = x/L, where L is the
length of the sample in the x direction). . . . . . . . . . . . . . . . . . . . . . . . .

74

xi
4.6 Saturation of the non-wetting fluid as a function of dimensionless
time ,tD (normalized by the breakthrough time) in the curvelet-
transformed images of the Berea sandstone porous medium with two
thresholds ε, and their comparison with those in the original image.

74
4.7 Drainage relative permeabilities of the CO2-water system in a Berea
Sandstone sample, original image (continuous curves) and curvelet
thresholded images (dashed curves). . . . . . . . . . . . . . . . . . . . . . . . . . .

75
4.8 Probability density function (PDF) of the normalized velocity
(u/Umax, where Umax is the maximum flow velocity) in the curvelet-
transformed images of the Sandstone porous medium with two
thresholds ε, and their comparison with those in the original image.

76
5.1 The image of the Berea sandstone used in the simulations. . . . . . . . . 84
5.2 The image of the Regenerated sandstone utilized in the simulations.

84
5.3 The image of the carbonate porous medium analyzed by curvelet
transformation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
85
5.4 Comparison of the original rock sample (a) with its thresholded
version (b), obtained with a threshold ε = 0.8. Observe that all the
essential features are preserved under thresholding. . . . . . . . . . . . . .

92
5.5 The connectivity function 𝑝 (ℎ;𝑚 ) computed for the Berea sandstone,
and its two coarsened versions, with m = 100. . . . . . . . . . . . . . . . . . .

94
5.6 Dependence of the permeability of the Berea sandstone on the
number of the lattice nodes N in the LB simulation. . . . . . . . . . . . . . 95
5.7 𝐶𝑂
2
invasion patterns in the Berea sandstone. (a) The original
image; (b) in the thresholded image with ε = 0.5 image, (c) in
thresholded image with ε = 0.9 with brine saturation 𝑆 𝑏 = 0.50. 𝐶𝑂
2

is shown in blue and brine in red. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

96
5.8 Saturation profile of 𝐶𝑂
2
along the direction x of macroscopic flow at
the breakthrough time in the original images of Berea sandstone, as
well as two thresholded ones, where 𝑋 = 𝑥 /𝐿 is the normalized
distance of the interface from the inlet, with L being the sample's
length in the x direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

96

xii
5.9 Comparison of time-dependence of saturation of 𝐶𝑂
2
, where 𝑡 𝐷 =
𝑡 /𝑡 𝐵 , with 𝑡 𝑠 being the time at which the system reaches steady state,
in the original image of the Berea sandstone and its curvelet-
transformed images. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97

5.10 Relative permeabilities of 𝐶𝑂
2
and brine in the Berea sandstone as
functions of the brine saturation 𝑆 𝑏 in the sandstone's original
image, as well as its curvelet-thresholded images. . . . . . . . . . . . . . . .

98
5.11 The probability distribution function (PDF) of the local velocity of
𝐶𝑂
2
in the Berea sandstone, normalized by its maximum, in the
curvelet-transformed images of the sandstone with two thresholds ε,
and its comparison with that of the original image. . . . . . . . . . . . . . .

99
5.12 Comparison of time-dependence of saturation of oil, where 𝑡 𝐷 =
𝑡 /𝑡 𝑠 , with 𝑡 𝑠 being the time at which the system reaches steady
state, in the original image of the Regenerated sandstone and its
curvelet-transformed images. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

101
5.13 Relative permeabilities of oil and water in the Regenerated
sandstone as functions of the brine saturation Sb in the sandstone's
original image, as well as its curvelet-thresholded images. . . . . . . . .

101

xiii

Abstract
Flow, transport, reaction, adsorption and deformation (FTRAD) are
phenomena that occur in a wide variety of heterogeneous materials and
media, ranging from biological membranes, to composite solids, and porous
formations. They have been studied for decades based on models of
heterogeneous media that involved approximations, simplifications, and
even unjustified assumptions. With advances in instrumentation, obtaining
highly resolved images of such media has become feasible, opening the way
to direct image-based simulation of the FTRAD. Due to complexity of such
images, however, the simulations still require long intensive computations,
limiting the size of the image that can be used.
We propose a novel methodology based on curvelet transformation (CT), a
multiscale and multidirectional transformation that provides an optimal
sparse representation of objects -images in our case-. Our methodology is
described as following, the image (i) is taken to the curvelet space; (ii) all the
curvelet coefficients that are smaller than a threshold are set to zero, and
(iii) the image is reconstructed. We then demonstrate that simulation of the
FTRAD in the reconstructed image yields as accurate results as those in the
original one, while speeding up the computation by two orders of magnitude
or better, hence taking a very significant step toward achieving the goal of
making direct image-based simulation of the FTRAD a standard practice.
The research presented in this dissertation demonstrates the virtues of
the method, by testing our proposed method on several different
applications. First, to demonstrate the increase in the computational

xiv
efficiency of simulation of the flow and transport processes, we solve for
unsteady-state diffusion problem and a single-phase fluid flow in complex
images of porous media using our method and compare the results with the
original image solution.
We then extend our proposed method for speeding-up simulation of
multiphase fluid flow in images of porous media by curvelet transformations,
which are specifically designed for processing of images that contain complex
curved surfaces. Most porous media contain correlations in their morphology
and, therefore, their images carry redundant information that, in the
curvelet transform space, can be removed efficiently and accurately in order
to obtain a coarser image with which the computations are far less intensive.
We utilize the methodology to simulate two-phase flow of oil and water in
two-dimensional digital images of sandstone and carbonate samples, and
demonstrate that while the results with the curvelet-processed images are
as accurate as those with the original ones, the computations are speeded up
by a factor of 110-150. Thus, the methodology opens the way toward
achieving the ultimate goal of simulation of multiphase flow in porous media.
Next, we use the same approach that utilizes curvelet transformations
(CTs) for speeding-up simulation of multiphase flow on 3D images of porous
media. As described previously, the CTs are designed for denoising of images
that contain complex curved surfaces, such as those of heterogeneous porous
media. This is possible because the morphology of most porous media contain
extended correlations, implying that their images carry redundant
information that can be eliminated by a suitable CT. As a result, simulation
of multiphase flow in the coarser images are sped up very significantly. The
new approach is used to simulate two-phase flow of CO2 and brine in a 3D
image of Berea sandstone. We demonstrate that the results with the CT-
processed image are as accurate as those with the original one, while

xv
computations are significantly faster by a speed-up factor between one order
of magnitude and more than two orders of magnitude.
We also use the same methodology on a prime method for carrying out
multiphase fluid flow simulations, which is the color-fluid lattice Boltzmann
method with multirelaxation time (CFLB-MRT) collision operator. The
simulations are, however, time consuming and intensive. We propose a
method to accelerate image-based computations with the CFLB-MRT
method. In the proposed method, the 3D image is preprocessed by curvelet
transforming it and eliminating those details that do not contribute
significantly to multiphase flow, which is done by thresholding the image.
After inverting the coarser image back to the real space, it is utilized in the
simulation of multiphase flow by the CFLB-MRT approach. As the test of the
method, we carry out simulation of a two-phase flow problem in which the
porous media are initially saturated by brine or water, which is then
displaced by CO2 or oil, injected into the pore space. The simulations are
carried out with two types of sandstone. We show that the method
accelerates the computations significantly by a speed-up factor of up to 35.
The results in all applications indicate that excellent agreement between
the solutions in the original images, and the less noisy and coarser images
by two thresholds, while the speed of the computations increases
significantly.

1

Chapter 1

Introduction

1.1 Background and Motivation
Fluid Flow and transport in complex environments, such as porous media
and materials, constitute a set of important problems, in view of their
relevance to a wide variety of phenomena in natural and industrial
processes. Examples of porous media and material vary from relatively small
length scales which are usually called laboratory scale, such as catalysts, cell
membranes, skin, and biological tissues, to large scale porous media
including oil and gas reservoirs and ground water aquifers.
Understanding the physics of flow, transport, reaction, adsorption and
deformation (FTRAD) in rocks and large-scale porous media, provides a
valuable insight into their properties, structure and fluid content. Predicting
the properties of FTRAD processes are vital to address many important
problems, such as energy production, global warming and managing water
resources. Studying such complex phenomena experimentally is time
consuming and expensive. In addition, experimental studies of FTRAD in
complex porous media and materials to measure the physical properties is a
very difficult task. On the other hand, various theoretical approaches have
been developed to study the FTRAD in porous media, due to the
heterogeneities one must ultimately resort to numerical simulations, since

2

exact solutions of such problems are always very difficult, if not impossible
to derive. Compared to experimental approaches, the numerical methods
have the advantages in terms of flexibility in changing parameters and
studying different flow mechanisms, also they are less expensive than
experiments.
Several numerical approaches for solving the equations governing the
FTRAD, such as the Lattice-Boltzmann's method (Huang et al., 2015; Krüger
et al., 2017) have been developed. Another approach is the pore network
model (Fatt, 1956; Blunt and King, 1991) in which the pore space of the
porous rock is simplified to a network of connected throats. Pore network
models have the advantages of being inexpensive and computationally more
efficient. However, pore network models are based on unjustified
assumptions and approximations in geometry as well as in the physical
details of the problem, which increases the uncertainty of the results.
A classical well-known method is based on modelling the FTRAD by
Eulerian grid-based methods such as, finite difference technique (Davis,
2011; Rahman et al., 2015; Kelly et al., 1976), which estimate the derivative
in a location by selecting values on grid points. It is fast and simple but may
cause high errors, especially when dealing with very complicated geometries,
such as porous media. A more accurate method is the finite volume method
and the finite element method (Lysmer et al., 1972; Meyer, 1992) have been
utilized to solve problems with complexed boundary conditions.
Advancements in instrumentation and enormous increase in
computational power, it has become increasingly possible to not only extract
much more information about the morphology of heterogeneous porous
media, but also to carry out the simulations of FTRAD directly in the image
of the complex media. For example, X-ray computed tomography is
increasingly being used (Arns et al., 2001, 2002) to obtain high resolution

3

images of porous media. Thus, image-based simulation of the FTRAD has
become more common (Piller et al., 2009; Tahmasebi and Sahimi, 2012;
Raeini et al., 2015; Andrew et al., 2015; Tahmasebi et al., 2015; Cooper et al.,
2017; Tahmasebi et al., 2017; Bakhshian et al., 2018). But, although such
simulations eliminate the need for developing models of porous media that
involve assumption and approximations, due to their heterogeneities, the
problem is still very difficult, because the structure of the complex image of
porous media and materials implies that a high resolution computational
grid must be used in order to compute the properties of interest, such grids
may contain several hundreds of thousands of grid points, or millions in some
cases. Thus, developing efficient numerical methods for solving FTRAD in
2D or 3D images of complex porous media has been widely considered in
research in the past few decades.
In this thesis we propose a novel methodology for speeding up image-based
numerical simulation of the flow and transport in porous media. The
methodology is based on processing of the image with curvelet
transformations (CTs) (Cand´es and Donoho, 1999, 2005) in order to smooth
out those details of the morphology that do not contribute significantly to the
phenomenon to be simulated. The basic process is as follows, first, the image
is taken to the curvelet space. Then, apply a threshold to the image and set
all the curvelet coefficients that are smaller than a threshold to 0. Finally,
apply the inverse CT and obtain the reconstructed image.

Figure 1.1: Illustration of Curvelet thresholding process of an image.

4

Finite-volume method (FVM) is used in this study to solve the flow and
transport problems in the original and curvelet transformed images. The
FVM is the most suitable technique for numerical computations in which the
grid blocks, although regularly shaped, do not have the same size, which is
the case with the coarsened grid that we obtain using CT. We demonstrate
that the approach speeds up the computations by about two orders of
magnitude or better, while producing results that are virtually identical with
those obtained with the original image.
On the other hand, another popular numerical method to simulate
multiphase fluid flow through porous media is the Lattice Boltzmann
method. However, LB method is also computationally intensive, which limits
the size of the image used for simulation. We propose a method for image-
based numerical simulation of two-phase flow in porous media using the LB
method, whereby the image is processed with 3D curvelet transforms (CTs)
in order to develop a coarse, yet accurate image of the pore space. We show
that the results are as accurate as those obtained by the original image, while
the computations' time is reduced by more than one order of magnitude.
The basic idea of our proposed methodology is, instead of attempting to
develop an even more efficient numerical algorithm, than what is currently
available for solving even larger sets of equations we focus on reducing the
number of equations to be solved by answering the following key question:
Given a complex image of a porous medium that contains a large number of
grid points, how much of the image must be retained in the model and how
much of them can be discarded, without distorting the image and the true
value of the flow and transport properties? (Daubechies, 1989).
Enhancing sparsity simplifies the numerical model by reducing number of
computations required to achieve high accuracy estimation of the flow or the
transport property of interest. We demonstrate that one can drastically

5

reduce the number of equations to be solved without impacting the important
information about the medium and its image and the true value of the
property to be calculated, hence, improving the computational efficiency.

1.2 Curvelet Transformation
Segmentation and processing of an image in order to obtain its sparse but
accurate representation has been an active research field, going back many
decades. But it was the seminal works of Daubechies (1988,1992) and Mallat
(1989a,b) advanced the field significantly by introducing multiresolution
wavelet transforms (WTs), which decompose an image into a series of high-
and low-pass filter bands, and extract the details of its structure in the
horizontal, vertical, and diagonal directions. Heterogeneous porous media,
particularly 3D ones, are, however, not only highly noisy, but also contain
curved pore surfaces. Thus, the three directions that the WTs employ are not
enough for extracting adequate amount of information from 3D images. Note,
however, that the WTs have been used very fruitfully in upscaling of
geological models of heterogeneous porous media (see, for example, Mehrabi
and Sahimi, 1997; Beylkin et al, 1991Ebrahimi and Sahimi, 2002,2004;
Pancaldi et al. 2007; Rasaei and Sahimi, 2008, 2009, Babaei and King 2011;
Rezapour et al. 2019). in geophysics (for reviews see Kumar and Foufoula-
Georgiou,1997; Grinsted et al., 2004), and many other applications.
Although wavelets have become widely popular in many applications, but
they still have several drawbacks, one major drawback is that wavelets are
suitable only for representing point singularities, since they are non-
directional and ignore the geometric properties of structures, which led to
the poor representation of curved edges.

6

Most natural signals or images contains edges, for example,
discontinuities along curves (line or curve singularities). Therefore, wavelets
and related ideas are computationally inefficient for natural geometric
features with curve singularities. In order to overcome the directional
sensitivity of wavelets, a multiscale multidirectional geometric analysis,
named curvelet transformation (CT) was developed (see, for example,
Donoho, 1999; Cand´es et al., 2005). Figure 1.2 compares between the
wavelet and curvelet coefficients at different scales, positions and directions.

(a) (b)
Figure 1.2: Wavelet coefficients (a), and curvelets on various scale,
directions and positions (b) (CurveLab software by, Cand´es et al., 2005).

The CTs are geometric multiscale transforms that allow optimal sparse
representation of objects - images in the present thesis - with second-order
continuous differentiable singularities by precise representation of the edges
along curves with very high directional sensitivity. They have already been

7

used in geophysics for denoising of seismic data (Herrmann et al., 2007;
Neelamani et al., 2008; Dashtian and Sahimi, 2014) and analysis of
nonstationary time series (Olhede and Walden, 2004); see Ma and Plonka
(2009) for other applications.

(a) (b) (c)

(d) (e) (f)
Figure 1.3: Curvelets at increasingly fine scales. Spatial domain (a-c).
Frequency domain (d-f) (CurveLab software by, Cand´es et al., 2005).

8

In Figure 1.3, the upper row is the spatial representation of curvelets at
scales 1,3 and 5. It can be seen that in 11(a) curvelets are nondirectional at
coarse scale. And the lower row is the frequency domain representation of
curvelets.
The CT has gone through two major stages. First, it was proposed in
(Donoho and Cand´es, 1999). It is based on complex series of steps involving
ridgelet transformation as preprocessing step of curvelets. Its limitation is
that it is accurate mostly for detecting linear radial structures, which are not
the type that one has in many types of images, ranging from those in medical
science (Al-Zubi et al. 2011), to heterogeneous porous media. Thus, the same
researchers developed the so-called curvelet second generation (Cand´es et
al. 2005) to overcome the first-generation drawbacks. Curvelet second
generation which have been shown to be highly effective at identifying an
image's important features along curves, an aspect that is highly important
to the rough pores' surface of a porous medium. It is also simple, easy to
understand and implement and faster, since the use of ridgelet
transformation as a preprocessing step of CT was discarded, which reduces
the amount of redundancy in the transformation.
Curvelet second generation is known as Fast Discrete Curvelet Transform
(FDCT). According to (Cand´es et al., 2005), there are two proposed
implementations of FDCT. First implementation is unequally spaced fast
Fourier transform (USFFT), and the second implementation is wrapping
function. They differ mainly by the choice of the spatial grid used for
transformation. Both implementations return a set of digital curvelet
coefficients in a table indexed by a scale, orientation and spatial location
parameters. Wrapping based curvelet transformation is faster and easier to
implement and understand.

9

1.2.1 Continuous-Time Curvelet Transform

Suppose that f(x) represents a 2D image of a porous medium, with ω being
frequency, and r and 𝜃 the polar coordinates. Consider a pair of functions
W(r) and V (𝜃 ), called the radial and angular windows, respectively, which
are smooth, nonnegative and real-valued functions, with W(r) taking on
positive real arguments and supported in r ϵ (1/2, 2), and V (𝜃 ) taking on real
arguments in 𝜃 ϵ [ − 1, 1]. The two windows satisfy the admissibility
conditions,

∑ 𝑊 2
(2
𝑗 𝑟 )=1
∞
𝑗 =−∞
,with 𝑟 ∈(
3
4
,
3
2
) (1.1)
and,
∑ 𝑉 2
(𝜃 −1)=1
∞
𝑗 =−∞
,with 𝜃 ∈(−
1
2
,
1
2
) (1.2)

For each 𝑗 a frequency window 𝑈 𝑗 is defined in the Fourier space by,

𝑈 𝑗 (𝑟 ,𝜃 )= 2
−3𝑗 4 ⁄
𝑊 (2
−𝑗 𝑟 ) 𝑉 (2
[
𝑗 2
⁄ ]
𝜃 2𝜋 ⁄ ) (1.3)

Where [·] represents the number’s integer part. By definition, the support
of 𝑈 𝑗 is a polar wedge defined by W(r) and V(𝜃 ) and applied with scale-
dependent window widths in each direction. To obtain real-valued curvelets,
we use the symmetric version, namely, 𝑈 𝑗 (r,𝜃 ) + 𝑈 𝑗 (r,𝜃 + π).
The mother curvelet φ
𝑖 (𝑥 ) is defined by its Fourier transform, 𝜑 𝑖 ̂ (ω)=
𝑈 𝑗 (ω), where 𝑈 𝑗 (ω
1
,ω
2
)= 𝑈 𝑗 (ω) is the window in polar coordinates. All the
curvelets at scale 2
−j
are obtained by rotation and translation of φ
𝑖 (𝑥 ). An

10

equispaced, but scale-dependent, sequence of rotation angles, 𝜃 𝑙 =2𝜋𝑙 ×
2
−[
𝑗 2
]
is then introduced with, 𝑙 = 0, 1, . . ., such that 𝜃 𝑙 =0<2, as is a
sequence of translation parameters with 𝑘 =(𝑘 1
,𝑘 2
)∈ 𝑍 2
. Then, the
curvelets are defined as a function of 𝑥 =(𝑥 1
,𝑥 2
) at scale (2
−𝑗 ), orientation 𝜃 𝑙 ,
and position 𝑥 𝑘 (𝑗 ,𝑙 )
= 𝑅 𝜃 𝑙 −1
(𝑘 1
×2
−𝑗 ,𝑘 2
×2
−𝑗 2 ⁄
) by (Donoho, 1999; Cand´s et
al., 2004),

𝝋 𝑗 ,𝑙 ,𝑘 (𝑥 )=𝜑 𝑗 [𝑅 𝜃𝑙
(𝑥 − 𝑥 𝑘 (𝑗 ,𝑙 )
)] (1.4)

where 𝑅 𝜃𝑙
represents a rotation by 𝜃 𝑙 radians, given by

𝑅 𝜃 𝐼 =(
𝑐𝑜𝑠 𝜃 𝑙 𝑠𝑖𝑛 𝜃 𝑙 −𝑠𝑖𝑛 𝜃 𝑙 𝑐𝑜𝑠 𝜃 𝑙 ) (1.5)

With 𝑅 𝜃 𝐼 −1
= 𝑅 𝜃 𝐼 𝑇 = 𝑅 −𝜃 𝐼 , and T representing the transpose operation. The
curvelet coefficients (CCs), 𝑐 𝑗 ,𝑙 ,𝑘 are defined as the inner product of the image
𝑓 (x) and the curvelet 𝝋 𝑗 ,𝑙 ,𝑘 ,

𝑐 𝑗 ,𝐼 ,𝑘 = <𝑓 ,𝝋 𝑗 ,𝐼 ,𝑘 > = ∫𝑓 (𝑥 )𝜑 𝑗 ,𝑙 ,𝑘 (𝑥 )
̅̅̅̅̅̅̅̅̅̅
𝑑𝑥 (1.6)

where the overline represents a complex conjugate. Since the CTs operate in
the frequency domain, we use the Plancherel’s theorem to express the inner
product as an integral over the frequency plane,

𝑐 𝑗 ,𝐼 ,𝑘 =
1
(2𝜋 )
2
∫𝑓 ̂
(ω)𝜑̂
𝑗 ,𝑙 ,𝑘 (ω)
̅̅̅̅̅̅̅̅̅̅̅
𝑑 ω=
1
(2𝜋 )
2
∫𝑓 ̂
(ω)𝑈 𝑗 (𝑅 𝜃 𝑙 ω)𝑒 (𝑖 <𝑥 𝑘 (𝑗 ,𝑙 )
,𝜔 >)
𝑑 ω (1.7)

11

As in the application of the wavelet transforms (Mehrabi and Sahimi,
1997; Ebrahimi and Sahimi, 2002,2004; Rasaei and Sahimi, 2008, 2009) to
simulation and upscaling of flow and transport in heterogeneous porous
media, we also have coarse-scale elements when working with the CTs. We
introduce a low-pass window 𝑊 0
by,

|𝑊 0
(𝑟 )|
2
+ ∑ |𝑊 (2
−𝑗 𝑟 )|
2
=1
𝑗 =0
(1.8)

and define for (𝑘 1
,𝑘 2
) ∈ Z the coarse scale curvelets by,

𝜑 𝑗 0,𝑘 (𝑥 )= 𝜑 𝑗 0
(𝑥 − 2
−𝑗 0
𝑘 ) (1.9)
and,
𝜑̂
𝑗𝑜
(ω)=2
−𝑗 0
𝑊 0
( 2
−𝑗 0
|ω|) (1.10)
Note that the coarse-scale curvelets are nondirectional, and that the
complete CT consists of the fine-scale directional curvelets (𝜑 𝑗 ,𝑙 ,𝑘 )
𝑗 =𝑗𝑜 ,𝑙 ,𝑘 and
the coarse-scale isotropic father wavelets (Φ
𝑗𝑜 ,𝑘 )
𝑘 .

1.2.2 Three-Dimensional Curvelet Transforms
Although the original CTs that were developed for image processing were
two dimensional (Donoho, 1999; Cand´s et al., 2005), there are many
scientific problems and engineering applications, such as processing of 3D
seismic and medical images, that motivated the development of 3D CTs,
particularly because such CTs preserve the important features of complex
3D systems, which is why their use is becoming more common.

12

Three-dimensional CTs are constructed in a manner similar to their 2D
counterparts (Ying et al., 2005; Esmaeili et al., 2017). Thus, in continuous
CTs a frequency window 𝑈 𝑗 ̂
(ω) is constructed by introducing two other
windows. One is a radial window 𝑊 (𝑟 ), defined by

∑ 𝑊 2
(2
𝑗 𝑟 )=1
∞
𝑗 =−∞
(1.11)

which smoothly extracts the frequency near the dyadic corona, {2
𝑗 −1
≤𝑟 ≤
2
𝑗 +1
}, with r being the polar coordinate. The second window 𝑉 (𝜃 ) is defined
by

∑ 𝑉 2
(𝜃 −1)=1
∞
𝑗 =−∞
,with 𝜃 ∈(−
1
2
,
1
2
) (1.12)

Implying that for each scale 𝑗 a unit sphere 𝑆 2
, representing all the
orientations in 𝑅 3
is considered, and is partitioned into 𝑂 (2
𝑗 /2
· 2
𝑗 /2
) =
𝑂 (2
𝑗 ) smooth angular windows 𝑉 (𝜃 ) with a circular support of radius
𝑂 (2
−𝑗 /2
), the squares of which form a partition of unity on 𝑆 2
. The window
𝑈 𝑗 (𝑟 ,𝜃 ) is then defined by

𝑈 𝑗 (𝑟 ,𝜃 )= 2
−3𝑗 4 ⁄
𝑊 (2
−𝑗 𝑟 ) 𝑉 (2
[
𝑗 2
⁄ ]
𝜃 2𝜋 ⁄ ) (1.13)

Here [·] denotes the integer part of the number. Then 𝑈 𝑗 ̂
(ω) is simply the
Fourier transform of 𝑈 𝑗 (𝑟 ,𝜃 ).
The mother curvelet 𝜑 𝑖 (𝑥 ), based on which all other curvelets at scale 2
−j

are constructed by rotation and translation of 𝜑 𝑖 (𝑥 ), is defined by 𝜑 𝑖 ̂ (ω)=
𝑈 𝑗 (ω). The 3D curvelets are defined as a function of 𝑥 =(𝑥 1
,𝑥 2
,𝑥 3
) at scale
( 2
−𝑗 ), orientation 𝜃 𝑗 ,𝑙 and position 𝑥 𝑘 (𝑗 ,𝑙 )
= R
𝜃 𝑗 ,𝑙 −1
(𝑘 1
×2
−𝑗 ,𝑘 2
×2
−𝑗 2 ⁄
,𝑘 3
×

13

2
−𝑗 2 ⁄
), i.e., 𝝋 𝑗 ,𝑙 ,𝑘 (𝑥 )=𝜑 𝑗 [𝑅 𝜃𝑙
(𝑥 − 𝑥 𝑘 (𝑗 ,𝑙 )
)] , where (𝑘 1
,𝑘 2
,𝑘 3
) are translation
parameters, and 𝑅 𝜃𝑙
is the rotation by the angle 𝜃 𝑗 ,𝑙

𝑅 𝜃 𝐼 =(
𝑐𝑜𝑠 𝜃 𝑙 𝑠𝑖𝑛 𝜃 𝑙 −𝑠𝑖𝑛 𝜃 𝑙 𝑐𝑜𝑠 𝜃 𝑙 ) (1.14)

Suppose that the 3D image of a porous medium, a 𝑛 ×𝑛 × 𝑛 array of pixels,
is represented by I(x). Then, its CT, usually referred to as the curvelet
coefficients (CCs) 𝐶 𝑗 ,𝑙 ,𝑘 , is defined by

𝑐 𝑗 ,𝐼 ,𝑘 == ∫𝐼 (𝑥 )𝜑 𝑗 ,𝑙 ,𝑘 (𝑥 )
̅̅̅̅̅̅̅̅̅̅
𝑑𝑥 (1.15)

which after using the Plancherels theorem, is rewritten in the frequency
space:

𝑐 𝑗 ,𝐼 ,𝑘 =
1
(2𝜋 )
3
∫𝐼̂
(ω)𝜑̂
𝑗 ,𝑙 ,𝑘 (ω)
̅̅̅̅̅̅̅̅̅̅̅
𝑑 ω=
1
(2𝜋 )
3
∫𝐼̂
(ω)𝑈 𝑗 (𝑅 𝜃 𝑙 ω)𝑒 (𝑖 <𝑥 𝑘 (𝑗 ,𝑙 )
,𝜔 >)
𝑑 ω (1.16)

where the overline denotes a complex conjugate.
Since the digitized image of the porous medium is discrete, we need to
rewrite the above formulation in discrete form. The 3D discrete CT operates
on a 3D Cartesian grid - the 3D image 𝐼 (𝑛 1
,𝑛 2
,𝑛 3
) with ,0 ≤𝑛 1
,𝑛 2
,𝑛 3
< 𝑛 of
the porous medium - and generates a set of the CCs 𝐶 𝐷 (𝑗 ; 𝑙 ; 𝑘 )
that are defined
by

𝐶 𝐷 (𝑗 ,𝑙 ,𝑘 )= ∑ ∑ ∑ 𝐼 (𝑛 1
,𝑛 2
,𝑛 3
) 𝜑 𝑗 ,𝑙 ,𝑘 𝐷 (𝑛 1
,𝑛 2
,𝑛 3
)
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
𝑛 3
𝑛 2
𝑛 1
(1.17)

with (𝑗 ,𝑙 ) 𝜖 𝑍 𝑎𝑛𝑑 𝑘 = (𝑘 1
,𝑘 2
,𝑘 3
). Since, as described below, by curve-
transforming of the image we obtain a new coarse-scale image, we also need
to consider such scales. Thus (Ying et al., 2005), suppose that the frequency
window 𝑈̃
𝑠 0
at scale 𝑗 = 𝑠 0
is defined in a rectangular box of integer size,

14

𝐿 1,𝑠 0
×𝐿 2,𝑠 0
×𝐿 3,𝑠 0
is defined by the usual relation, 𝑈̃
𝑠 𝑜 ,
0
(𝜔 )=𝑊̃
𝑠 𝑜 (𝜔 ), Then,
the 3D discrete curvelets at coarsest scale are defined by their Fourier
transform:

𝜑 𝑠 0
,0,𝑘 𝐷 (𝜔 )=
1
𝑁 𝑈̃
𝑠 𝑜 ,
0
(𝜔 ) ·𝑒𝑥𝑝 [−2𝜋 𝑖 (𝑘 1
𝜔 1
𝐿 1,𝑠 𝑜 ⁄ +𝑘 2
𝜔 2
𝐿 2,𝑠 𝑜 +𝑘 3
𝜔 3
𝐿 3,𝑠 𝑜 ⁄ ⁄ )] (1.18)

Where N is given by

𝑁 =(𝐿 1,𝑠 0
𝐿 2,𝑠 0
𝐿 3,𝑠 0
)
1/2
, (1.19)

for 0 ≤ 𝑘 1
< 𝐿 1,𝑠 𝑜 ,0 ≤ 𝑘 2
< 𝐿 2,𝑠 𝑜 𝑎𝑛𝑑 0 ≤ 𝑘 3
<𝐿 3,𝑠 𝑜 .
On the other hand, for the fine scales 𝑠 0
< 𝑗 every Cartesian corona has
six components, one for each face of the unit cube, and is partitioned into
2
𝑗 /2
·2
𝑗 /2
=2
𝑗 same-volume wedges. For example, suppose that for the 𝑙 th
wedge in the first component, (1, 𝛼 𝑙 , 𝛽 𝑙 ) indicates the direction of the center
line of the wedge. Then, an angular window is defined by

𝑉̃
𝑗 ,
𝑙 (𝜔 )= 𝑉̃
(2
𝑗 /2
·
𝜔 2
− 𝛼 𝑙 ·𝜔 1
𝜔 1
)·𝑉̃
(2
𝑗 /2
·
𝜔 3
− 𝛽 𝑙 ·𝜔 1
𝜔 1
) , (1.20)

with 𝑉̃
given by Eq. (1.12). Similar definitions hold for the other five
components by exchanging the roles of 𝜔 1
with 𝜔 2
or 𝜔 3
. Then, the frequency
window 𝑈̃
𝑗 ,
𝑙 is defined by, 𝑈̃
𝑗 ,
𝑙 (𝜔 )= 𝑊̃
𝑗 (𝜔 )·𝑉̃
𝑗 ,
𝑙 (𝜔 ) which isolates the
frequencies near the wedge {(𝜔 1
,𝜔 2
,𝜔 3
) ∶ 2
𝑗 −1
≤ 𝜔 1
≤2
𝑗 +1
,−2
−𝑗 /2
≤
𝜔 2
𝜔 1
⁄ −𝛼 𝑙 ≤2
−𝑗 /2
,−2
−𝑗 /2
≤ 𝜔 3
𝜔 1
⁄ −𝛽 𝑙 ≤2
−𝑗 /2
} . Then 𝜑 𝑗 ,𝑙 ,𝑘 𝐷 (𝜔 ) , the
discrete curvelet at the fine scales 𝑠 0
< 𝑗 with index 𝑘 at scale 𝑗 and angle 𝑙 ,
are defined in a manner similar to Eq. (1.18), but with 𝐿 𝑖 ,𝑠 0
replaced by 𝐿 𝑖 ,𝑗 ,𝑙
with 𝑖 = 1,2,𝑎𝑛𝑑 3.

15

1.2.3 Curvelet Transformation Features
Curvelet transformation is a special member of the multiscale geometric
transforms (Donoho, 1999; Cand´s et al., 2004) which was basically
developed to overcome the limitations of traditional multiscale
transformations such as wavelets. What makes curvelets special and
different from other multiscale representations is that curvelets contain a
very high degree of directional sensitivity, as it is shown in Figure 1.2. In
principle, the main objective of curvelets is to deal with interesting
phenomena that occur along curved edges in an image. As illustrated in
Figure 1.4, curvelets require fewer elements for representations, and it
produces much smoother edges than wavelets (Boubchir and Fadili, 2005).

(a) (b)
Figure 1.4: Curve representation comparison between (a) wavelet and (b)
curvelet.

Curvelets are superior to other multiscale representations as they
efficiently address very important limitations where wavelets are far from
ideal. As in the following:

16

(a) As mentioned earlier, curvelets allow optimal sparse representation of
objects (images) with second-order continuous differentiable singularities.
Optimal in the sense that there is no other multiscale analysis tool can yield
a smaller error with the same number of elements. They represent edges and
singularities along curves precisely by elements with very high directional
sensitivity. Moreover, it is sparse in the sense that most of the elements have
very small magnitudes, except for a few large ones that represent high
frequency features such as curved edges.

(b) Curvelets provide optimal image reconstruction. They have a special
geometric feature which can be utilized to certain reconstruction problems
with missing or noisy data with very high accuracy. For instance, in many
medical applications one wishes to reconstruct an image from noisy and
incomplete tomographic data. In fact, this is considered one of the
challenging problems especially when we have incomplete data. A great
application on this is seismic data exploration and denoising which will be
described in the next section.

(c) Optimal sparse representation of wave propagation. Curvelets are very
powerful tools for analyzing and computing partial differential equations.
Since, curvelets greatly model the geometry of wave propagation. In other
words, curvelets can be viewed as waveforms with enough frequency
localization which makes them behave like waves, and at the same time they
have spatial localization which makes them simultaneously behave like
particles. As mentioned in (a), curvelets provide sparse representation, the
sparsity simplifies the mathematical analysis which can be utilized for
speeding up the simulations by reducing the number of computations
required to achieve a certain accuracy (Cand´s et al., 2005).

17

1.2.4 Recent Applications
Curvelet transforms has been widely applied in many fields of science and
engineering, especially in the field of image processing, seismic exploration
and scientific computing. In image processing, CT may be used to denoise
images and detect features embedded in them, extracting patterns from
large digital images and compress large images. Figure 1.5 shows an
example of image denoising using curvelets.

(a) (b) (c)
Figure 1.5: Original image (a), Noisy image (b) and curvelet denoising image
(c).

In seismic data exploration and denoising, a critical problem in this field
is to preserve the smoothness of the wave fronts when denoising the seismic
data. Curvelets can be viewed as local plane waves. They provide optimal
sparse representation of local seismic patterns and can effectively preserve
the wave fronts in seismic processing. Curvelets have been successfully
applied in seismic data exploration and denoising (Herrmann et al., 2007;
Neelamani et al., 2008; Dashtian and Sahimi, 2014).

18

The third application is using curvelets in scientific computing, the CT
may be used for speeding up numerical simulations, such as the numerical
propagation of waves in heterogenous media and material of special interest.
Some examples of using curvelets for scientific computing is turbulence
analysis and solving partial differential equations. A crucial problem in
turbulence simulation is to efficiently compress the flow field with minimum
loss of the geometric flow structures. Recently, curvelets have been used in
extraction of the coherent vortex field in turbulent flows (Moreno and Pullin,
2008; Ma and Hussaini, 2007,2009). Moreover, curvelets may be used in
solving PDEs numerically, as they provide an optimal sparse representation
of the solution operator (Cand´s et al., 2005). The key idea is to transform
the initial field (images in our case) to curvelet space, then compute the
significant curvelet coefficients through solving the associated governing
equation. Finally, reconstruct the coefficients at all scales by an inverse CT
and obtain the approximated solution (see, for example, Ma et al., 2009) . By
using a suitable threshold in curvelet domain, one can select the significant
curvelet coefficients and discard unnecessary coefficients that do not affect
the results. Hence, speeding up the simulations with high accuracy.

1.2.5 CurveLab
A software package of curvelets called curvelab made by (Cand´es et al.,
2005) was used to curvelet transform the images. Curvelab includes several
Matlab and C++ programs for the fast discrete curvelet transform to
demonstrate how to use this software, curvelab is available at
http://www.curvelet.org. There are two distinct implementations included in
the software package for 2D and 3D curvelet transforms: the wrapping based
transform and the unequally spaced fast Fourier transform (USFFT). In our

19

research, we used the wrapping based transform, which is simpler, faster
and more accurate.

1.3 Organization of This Dissertation
The research presented in this dissertation demonstrates the virtues of
our proposed method. First, to demonstrate the accuracy and the increase in
the computational efficiency of simulation of flow and transport processes,
we solve the diffusion problem and the single-phase fluid flow problem in two
types of complex images of porous media, Mt. Simon sandstone and
carbonate samples. To simulate both problems we have two options: (i)
determining the inverse CT of the new image, and carry out the simulation
of the FTRAD in it, or (ii) take the CT of the equations governing the FTRAD,
carry out the numerical simulations in the CT space, and then determine the
inverse CT of the quantities of the results. We have used both approaches,
as described in Chapter 2.
Next, solving more complex phenomena in the images of all types of
heterogenous materials and media, such that two-phase flow of immiscible
fluids. Understanding and modelling two phase immiscible flow in porous
media is one of the well-known challenging and important problem, due to
its wide range of applications including enhancing oil and gas recovery. In
this research, we are focusing on simulating two phase immiscible flow in a
complex sandstone and carbonate images. Due to the complexity of the
geometry and flow patterns happening on different spatial and temporal
scales make it difficult to model and simulate, that requires the use of high-
resolution computational grid to achieve high accuracy. Here, we will
compare the original image-based simulation of the problem with our
proposed method which is based on processing the image with curvelet

20

transformation, in order to speed-up the simulation of the same problem. To
do so, we first CT the original image, apply the threshold value. Then, set
the curvelet coefficients below the threshold to zero, finally apply the inverse
CT to obtain the reconstructed image. Then, we use it for simulation
(Chapter 3).
We then describe our approach for 3D image-based simulation of
multiphase flow in porous media. The main idea is to focus on the image,
process it in order to denoise it and remove the complexities that do not
contribute significantly to flow and transport in the pore space, and then
carry out numerical simulation of multiphase flow in the coarse image. In
this research, we divide a 3D image into multiple 2D slices, apply the
methodology to each 2D image, and then stack together the curvelet-
transformed 2D denoised images to reconstruct a 3D image and use it for
simulation, more details in chapter 4.
Lattice-Boltzmann simulation of multiphase flow in heterogeneous porous
media, represented by their high-resolution images, has become practical,
but the computations are still intensive and take a long time. Thus, one must
address the issue of high computational cost of such simulations. We applied
our proposed method to speed up simulations of two-phase flow with color-
fluid Lattice-Boltzmann with multirelaxation time collision operator (CFLB-
MRT) method. To do so, the 3D image is preprocessed by curvelet
transforming it to eliminate those details of the morphology that do not
contribute significantly to multiphase flow simulation. Then, inverting the
coarser image and utilize it in the simulation of multiphase flow by CFLB-
MRT method. As shown in chapter 5, the methods accelerate the
computations very significantly by one order of magnitude.

21

Chapter 2
Efficient Image-Based Simulation of Flow and
Transport in Heterogeneous Materials and Media:
Application of Curvelet Transforms

2.1 Introduction
Flow, transport, reaction, adsorption, and deformation (FTRAD) in
heterogeneous porous media are encountered in a vast number of systems of
fundamental scientific interest and industrial applications and in particular
in large-scale porous formations, such as oil, gas, and geothermal reservoirs
and groundwater aquifers (for reviews see, Blunt, 2017; Hunt & Sahimi,
2017; Sahimi, 2011). The heterogeneity is manifested by the variations in the
size, shape, and connectivity of the pores, as well as their physical properties.
Although various theoretical approaches have been developed to study the
FTRAD in porous media, due to the heterogeneities one must ultimately
resort to numerical simulations. Thus, various models of heterogeneous
porous media and several numerical approaches for solving the equations
governing the FTRAD, such as the Lattice-Boltzmann's method (Huang et
al., 2015; Krüger et al., 2017), have been developed.
Considerable advancements in the development of nondestructive two-
and three-dimensional (3D) imaging (Kak and Slaney, 2001; Herman, 2009)

22

have made it possible to not only extract much more information about the
morphology of heterogeneous porous media, but also use the images directly
in the numerical simulation of the FTRAD. In particular, X-ray imagining
eliminates the need for destructive sectioning; high-resolution focus ion-
beam scanning electron microscopy has become an essential part of
characterization of porous materials’ microstructure, and X-ray computed
tomography is increasingly being used (Arns et al., 2001, 2002). Thus, image-
based simulation of the FTRAD has become more common (Piller et al., 2009;
Tahmasebi and Sahimi, 2012; Raeini et al., 2015; Andrew et al., 2015;
Tahmasebi et al., 2015; Cooper et al., 2017; Tahmasebi et al., 2017;
Bakhshian et al., 2018). But, although such simulations eliminate the need
for developing models of porous media that involve assumption and
approximations, due to their heterogeneities, image-based simulations entail
use of a high-resolution computational grid with hundreds of thousands of
grid blocks, if not larger.
In this chapter we propose a novel methodology for speeding up image-
based numerical simulation of the FTRAD. The methodology is based on
processing of the image with curvelet transformations (CTs) in order to
smooth out those details of the morphology that do not contribute
significantly to the phenomenon to be simulated. We demonstrate that the
approach speeds up the computations by about two orders of magnitude or
better, while producing results that are virtually identical with those
obtained with the original image. In the next section we briefly explain how
the original image is made smooth in the curvelet space. Image-based
simulation of flow and transport processes in two distinct types of porous
media are described in section 3. Section 4 compares the computational
efficiency of the new method with that of the traditional simulations.

23

Extension of the method to 3D images is described in section 5. The chapter
is summarized in the last section.

2.2 Thresholding an Images of a Porous Medium
As is well known in the theory of image processing with the CTs (Starck et
al., 2002), 𝐶 𝑗 ,𝑙 ,𝑘 represents a measure of the local roughness or complexity of
an image’s morphology. Thus, the larger the CCs are, the more significant is
the roughness contribution to f(x) and, therefore, to image-based simulation
of the FTRAD in f(x). Therefore, after computing the CCs 𝐶 𝑗 ,𝑙 ,𝑘 of the image
f(x), we normalize them by the largest CC and introduce a threshold 0 ≤ ε <
1. Each of the CCs is then examined. If each 𝐶 𝑗 ,𝑙 ,𝑘 < 𝜀 , we set it to zero, hence
ignoring the local roughness in that particular neighborhood. After sweeping
the CT of f(x) once, we obtain the new CT of the image in which many of the
CCs are zero. As discussed below, the choice of the threshold depends on the
computational time that one can afford and the computational resources that
one has access to. We then have two options: (i) computing the inverse CT of
the thresholded image and carrying out the image-based simulation in it, or
(ii) taking the CT of the governing equations, carrying out the numerical
simulations in the CT space and determining the inverse CT of the results.
We have used both methods, as described in the following.
Note that this scheme is different from upscaling of computational grids
by wavelet transforms (Rasaei & Sahimi, 2008, 2009), since in the wavelet
upscaling the local permeabilities are upscaled explicitly, whereas in the
present method it is the image itself that is coarsened by the CT. We should
also point out that the entire curvelet thresholding of the images and
bringing them back to real space represents only about 2% of the total
computation time for simulating the FTRAD in the images.

24

(a) (b) (c)

(d) (e) (f)
Figure 2.1: Comparison of the images of two porous media with their
curvelet-transformed versions. (a), (b), and (c) show, respectively, the
sandstone, a small zoomed-in portion of it, and its curvelet-transformed with
the threshold ϵ = 0.8. (d), (e), and (f) show the same, but for the carbonate
medium.

To demonstrate the virtues of the method, we used 2D images of two types
of porous media, a fossiliferous outcrop carbonate, Mt. Gambier limestone in
South Australia with a porosity of 43%. The second example was (Bakhshian
et al., 2018; Tahmasebi et al., 2017), an image of a core plug from Mt. Simon
formation (named for Mount Simon escarpment in Eau Claire County,

25

Wisconsin) with a porosity of 25 %. We utilized one 2D slice of size 512
2
pixels
from each image. Figure 2.1 presents a small portion of each, both before and
after thresholding them with ε = 0.8, showing clearly that the roughness of
the internal pore surfaces, which is a prime reason for long duration of
image-based simulations, has become smooth. Generally speaking, small
CCs are always associated with such noisy features in the imaging data, and
therefore, their removal (or smoothening) does not alter the results
significantly.

2.3 Simulation of Flow and Transport
We simulated transient diffusion, as well as single-phase fluid flow in both
the sandstone and the carbonate samples. Thus, we first determined the size
of the computational grids in the original images and the thresholded ones
that yield macroscopic diffusivity 𝐷 𝑒 and Darcy permeability 𝐾 𝑒 that no
longer change, if the grid size is increased. Diffusion is, of course, governed
by the standard equation

∂C/∂t = D ∇
2
C (2.1)

where 𝐶 (𝑥 ,𝑡 ) is the concentration field at time 𝑡 , and D is the diffusivity. To
simulate the process, concentrations C0 and 0 were imposed between two
opposite faces of the image in one direction, and no-flux boundary condition
was used in the transverse direction, with the initial condition being
𝐶 (𝑥 ,0) = 0. As for the fluid flow, we solved the Stokes’ equation,

μ∇
2
v −∇ P = 0 (2.2)

26

for the fluid velocity v in the same two images, where μ is the fluid viscosity
(assumed to be water), and 𝑃 is the pressure. A constant pressure gradient
was imposed across the images in one direction, and periodic boundary
conditions were utilized in the transverse direction. We used the finite-
volume method with square (but unequal) grid blocks. All the discretized
equations were solved by the conjugate-gradient method.

Figure 2.2: Transient concentration profiles in the curvelet-transformed
images of the sandstone with two thresholds ϵ, and their comparison with
those in the original image.

Figure 2.2 compares the concentration normalized C∕C0, averaged over the
transverse direction, in the sandstone at time, t = 5 s, a short time after the
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
C/C0
X
Sandstone
Original, t = 5 s
= 0.5
= 0.9
Original, t = 100 s
= 0.5
= 0.9
ε
ε
ε
ε

27

image was exposed to the concentration difference, and a second one at an
intermediate time, t = 100 s. Figure 2.3 presents the same, but in the
carbonate porous medium. In both cases X = x/L with X representing the
dimensionless distance from the image’s inlet in the direction of the
macroscopic concentration difference, and L being the images’ length. Both
figures indicate excellent agreement between the solutions in the original
images and two images obtained by the two thresholds. Differences will
emerge only when 𝜀 → 1.

Figure 2.3: Transient concentration profiles in the curvelet-transformed
images of the carbonate porous medium with two thresholds ε, and their
comparison with those in the original image.

0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
C/C0
X
Carbonate
Original, t = 5 s
= 0.5
= 0.9
Original, t = 100 s
= 0.5
= 0.9
ε
ε
ε
ε

28

Figures 2.4 and 2.5 present the steady-state concentration profiles in the two
porous media. The agreement between all the profiles is excellent. In fact,
except for a small region near the inlet of the system at 𝑋 = 0 (which is due
to boundary effects), there is very little difference, if any, between all the
profiles. As for fluid flow in the same porous media, we computed the
effective permeability, with the results described in the next section.

Figure 2.4: steady-state concentration profiles in the curvelet-transformed
images of the sandstone with two thresholds ε, and their comparison with
those in the original image.

0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
C/C0
X
Sandstone
Original
= 0.5
= 0.9
ε
ε

29

Figure 2.5: steady-state concentration profiles in the curvelet-transformed
images of the carbonate porous medium with two thresholds ε, and their
comparison with those in the original image.

2.4 Computational Efficiency
Table 2.1 compares the details of the computations for the diffusion
simulation, where N is the number of grid blocks. Also presented are the
results for the cases in which the images were curvelet transformed, but no
threshold was applied, and numerical simulations of the diffusion and Stokes’
equations were carried out in the curvelet-transformed images. With the
results inversed curvelet-transformed back to the real space. To do so, we took
the CT of the two governing equations. Since digital CTs operate in the
frequency domain, we substituted the Fourier transform 𝑓 ̂
(𝜔 ) and 𝑈 (𝑅 𝜃 𝑙 𝜔 ) in
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
C/C0
X
Carbonate
Original
= 0.5
= 0.9
ε
ε

30

Eq. (1.7). We now show how to estimate the derivatives, by using the pseudo-
spectral method that gives a relation between Fourier transform of a function
and its derivatives. Using this idea will lead us to a new expression of the
diffusion equation in the curvelet domain. Let, 𝜑̂
𝑗 ,𝐼 ,𝑘 (𝜔 ) be the curvelet
function in the frequency domain

𝜑̂
𝑗 ,𝐼 ,𝑘 (𝑤 )=𝑈 𝑗 ̃
(𝑆 𝜃𝑙
𝜔 )𝑒 𝑖 <𝑥 𝑘 (𝑗 ,𝐼 )
,𝜔 >
(2.3)

the curvelet coefficients are defined by

𝐶 𝑗 ,𝐼 ,𝑘 =
1
(2𝜋 )
2
∫𝐶̂
(𝜔 )𝜑̂
𝑗 ,𝑙 ,𝑘 (ω)
̅̅̅̅̅̅̅̅̅̅̅
𝑑𝜔 (2.4)

𝐶 𝑗 ,𝐼 ,𝑘 ∆
=
1
(2𝜋 )
2
∫𝐹𝑇 (∆ 𝐶 ) 𝜑̂
𝑗 ,𝑙 ,𝑘 (ω)
̅̅̅̅̅̅̅̅̅̅̅
𝑑𝑤 (2.5)
from Fourier transform properties we know

𝐹𝑇 (
𝜕 𝑛 𝑓 𝜕 𝑡 𝑛 )=(𝑖𝜔 )
𝑛 𝑓̂
(𝜔 ) (2.6)

Then, the curvelet coefficients becomes

C
j,I,k
∆
=−
1
(2π)
2
∫|ω|
2
𝐶̂
(ω)𝜑̂
𝑗 ,𝑙 ,𝑘 (ω)
̅̅̅̅̅̅̅̅̅̅̅
dω (2.7)

Using eq. (2.3) and perform change of variables from 𝑆 𝜃 𝐼 𝜔 to ω, yields the
following

31

𝐶 𝑗 ,𝐼 ,𝑘 =
1
(2𝜋 )
2
∫𝑆 𝜃𝐼
−
𝐶̂
(𝜔 )𝑈 𝑗 ̃
(𝜔 )𝑒 𝑖 <𝑦 ,𝜔 >
(2.8)
𝐶 𝑗 ,𝐼 ,𝑘 ∆
=−
1
(2𝜋 )
2
∫(
𝜔 1
2
𝑐𝑜𝑠
2
𝜃 1
+𝜔 2
2
−2𝜔 1
𝜔 2
𝑡𝑎𝑛 𝜃 1
) 𝑆 𝜃𝐼
−
𝐶̂
(𝜔 )𝑈 𝑗 ̃
(𝜔 )𝑒 𝑖 <𝑦 ,𝜔 >
(2.9)
From eq. (2.9) 𝐶 𝑗 ,𝐼 ,𝑘 ∆
compose from three parts, 𝐶 𝑘 1
∆
,𝐶 𝑘 2
∆
,𝐶 𝑘 3
∆

𝐶 𝑘 1
∆
=
1
(2𝜋 cos𝜃 1
)
2
∫−𝜔 1
2
𝐶̃
(𝜔 )𝑒 𝑖 <𝑦 ,𝜔 >
𝑑𝜔 (2.10)

𝐶 𝑘 2
∆
=
1
(2𝜋 )
2
∫−𝜔 2
2
𝐶̃
(𝜔 )𝑒 𝑖 <𝑦 ,𝜔 >
𝑑𝜔 (2.11)

𝐶 𝑘 1𝑘 2
∆
=
2𝑡𝑎𝑛 𝜃 1
(2𝜋 )
2
∫𝜔 1
𝜔 2
𝐶̃
(𝜔 )𝑒 𝑖 <𝑦 ,𝜔 >
𝑑𝜔 (2.12)

The integral parts can be seen as an inverse Fourier transform of −𝜔 1
C
𝑎 𝜙 ̃
(𝜔 )
thus,
𝐶 𝑘 1
∆
=
1
𝑐𝑜𝑠
2
𝜃 l
𝜕 2
𝐶 𝑗 ,𝑙𝑘
𝜕 𝑦 2
(2.13)

𝐶 𝑘 2
∆
=
𝜕 2
𝐶 𝑗 ,𝑙𝑘
𝜕 𝑦 2
2
(2.14)

𝐶 𝑘 1𝑘 2
∆
=−2𝑡𝑎𝑛 𝜃 𝑙 𝜕 2
𝐶 𝑗 ,𝑙𝑘
𝜕 𝑦 1
𝜕 𝑦 2
(2.15)
Now, the CT of the diffusion equation is then given by:

𝜕 𝐶 𝑗 ,𝑙𝑘
𝜕𝑡
=𝐷 [(1+𝑡𝑎𝑛 2
𝜃 𝑙 )
𝜕 2
𝐶 𝑗 ,𝑙𝑘
𝜕 𝑦 1
2
+
𝜕 2
𝐶 𝑗 ,𝑙𝑘
𝜕 𝑦 2
2
−2𝑡𝑎𝑛 𝜃 𝑙 𝜕 2
𝐶 𝑗 ,𝑙𝑘
𝜕 𝑦 1
𝜕 𝑦 2
] , (2.16)

32

where 𝐶 𝑗 ,𝑙 ,𝑘
are the CCs associated with 𝐶 (𝑡 ,𝑥 ), and (𝑦 1
, 𝑦 2
) are the coordinates
in the CT space. We then solved for 𝐶 𝑗 ,𝑙 ,𝑘
in the curvelet-transformed image by
discretizing Eq. (2.16) and defining the scale 𝑗 , direction 𝑙 , and translation
parameter 𝑘 =(𝑘 1
, 𝑘 2
) ∈ 𝑍 2
, with (𝑘 1
, 𝑘 2
) = (𝑘 1
· 2
−j
,𝑘 2
·2
−j/2
). The inverse CT
of the results was then computed. Note that the Stokes’ equation has a
Laplacian term similar to the diffusion equation (but with D → μ). The rest of
the computations for the flow problem is similar to those for the diffusion
problem.

Table 2.1. Comparison of the number of grid blocks N and the computation
times for various images.
Sandstone N Threshold ε De × 10
7
(cm
2
/s) Time (CPU sec)
Original image 238941 23.72 602
Original image in CT space 54181 0 23.95 134
Curvelet-transformed image 3659 0.5 23.97 10
Curvelet-transformed image 2182 0.7 24.08 7
Curvelet-transformed image 1896 0.9 24.24 5
Carbonate
Original image 543069 2.02 1556
Original image in CT space 129052 0 2.07 391
Curvelet-transformed image 8641 0.5 2.07 26
Curvelet-transformed image 5174 0.7 2.08 15
Curvelet-transformed image 4366 0.9 2.08 13
Note: De is the effective diffusivity.

33

Table 2.2. Comparison of the number of grid blocks N and the computation
times for various images.
Note: Ke is the effective permeability.

Table 2.2 makes the same comparison as in Table 2.1 for fluid flow in the two
images. There are several noteworthy features of Table 2.1 and 2.2. (i) Even
with a threshold of ε = 0.9, we still obtain highly accurate solutions that differ
from those of the original images by at most 2–3%, well within the possible
errors caused by the images themselves. (ii) The speedup in the computations
in both types of porous media is up to a factor of 120 for ε = 0.9, for both
problems, more than two orders of magnitude improvement in the
computation time. (iii) Even with a relatively low threshold of ε = 0.5 the
Sandstone N Threshold ε Ke (mD) Time (CPU sec)
Original image 243046 32.38 714
Original image in CT space 58395 0 32.71 184
Curvelet-transformed image 3967 0.5 32.93 12
Curvelet-transformed image 2355 0.7 33.47 8
Curvelet-transformed image 1985 0.9 33.57 6
Carbonate
Original image 583383 5.54 1704
Original image in CT space 154337 0 5.71 489
Curvelet-transformed image 9865 0.5 5.74 31
Curvelet-transformed image 5684 0.7 5.77 17
Curvelet-transformed image 4993 0.9 5.78 15

34

speed-up factor is about 60, a very significant improvement in the efficiency
of the computations. (iv) When we carry out the simulations in the CT space
without any thresholding of the image, we still have a speed-up factor of 4 or
better. The reason is that one requires a coarser computational grid in the
curvelet space, as many of the CCs of the image are either very small and,
thus, do not make significant contributions, or that they are close to those of
their neighbors, indicating a relatively smooth local environment, allowing
use of larger grid blocks. (v) As 𝜖 increases, so also do Ke and De, albeit
slightly. This is due to the removal of the noisy pixels in the original image
by the CT thresholding. Although the volume fraction of such pixels is very
small, they serve as obstacles for flow and transport. Therefore, their removal
leads to (slightly) increase in Ke and De. As such noisy pixels may be due to
artifacts during image processing, their removal would actually lead to more
accurate estimates of the flow and transport properties and closer to their
true values in the porous media than what would obtain with the original
image.

2.5 Extensions of the Methodology
To demonstrate the accuracy and efficiency of the computations, we solved
the diffusion equation, the same levels of accuracy and speed-up are expected
for image-based simulation of more complex phenomena in all type of
heterogeneous media, such as deformation induced by multiphase flows,
propagation of waves in highly heterogeneous media, and fracture
propagation in composite solids. Our preliminary simulation of two-phase
flow of two immiscible fluids in the same images confirms this.
Extension of the method to simulation in 3D images is also
straightforward. Two approaches may be taken to extend the method to 3-D

35

images: (i) One may use 3D curvelet transforms that have been developed
recently (Woiselle et al., 2010; Ying et al., 2005), the proposed methodology
is directly applicable. We suspect that the speed-up factor for computations
with 3-D images will be even larger than what we have presented here. (ii)
Alternatively, one can split the 3-D image into a large number of 2-D slices,
apply the 2-D methodology to each, and then stack together the curvelet-
transformed 2-D images to reconstruct a 3-D image. In this case, care should
be taken to make sure that the transition from one slice to the next is
seamless.
More generally, the proposed methodology is applicable to the simulation
of any physical processes in any types of materials and media in which their
complex morphology plays a fundamental role. As such, the proposed
methodology is a universal approach to image-based simulation of physical
phenomena in complex media.

2.6 Summary
As image-based simulation of the FTRAD in heterogeneous porous media
becomes more practical, one must address the issue of its high computational
cost. We have proposed a novel methodology, whereby the image is curvelet
transformed in order to remove those details of the morphology that do not
contribute significantly to the simulation of the FTRAD by setting to 0 all
the CCs that are smaller than a threshold. The simulation in the less noisy
images yields as accurate results as those in the original ones but with the
speed of the computation increasing by 2 orders of magnitude or better. Thus,
the method represents a significant step toward making image-based
simulation of the FTRAD in heterogeneous porous media a standard
practice.

36

Chapter 3

Speeding-up Simulation of Multiphase Flow in Digital
Images of Heterogeneous Porous Media by Curvelet
Transformation

3.1 Introduction
Multiphase flow in heterogeneous porous media is important to many
fundamental scientific problems, as well as industrial applications (Helmig
and Schulz 1997; Sahimi 2011; Blunt 2017). Examples include evaporation
and drying in soil (Shokri et al. 2012), ink imbibition in printing paper
(Ghassemzadeh et al. 2001; Ghassemzadeh and Sahimi 2004; Swerin 2018;
Aslannejad and Hassanizadeh 2017; Aslannejad et al. 2019), extraction of
oil, gas and geothermal energy from porous formations, transport of
contaminants in soils and aquifers, and sequestration of CO2 in depleted oil
reservoirs (Friedlingstein and Solomon 2005; Nordbotten and Celia 2011;
Kohanpur et al. 2020). Aside from the classical continuum models of such
phenomena (Bear 1972), studies of multiphase flows in a porous medium
entail, as the first step, developing models of its pore space. There is a wide
variety of such models in the literature, ranging from percolation models (for
a recent review see Hunt and Sahimi 2017), both random (Larson et al. 1977;
Heiba et al. 1992; Sahimi et al. 1986; Kantzas and Chatzis 1988; Blunt and
Scher 1995) and invasion percolation (Chandler et al. 1982; Wilkinson and

37

Willemsen 1983; Knackstedt et al. 2001; for a simple introduction to invasion
percolation see Ebrahimi 2010), detailed pore-network models (Blunt and
King 1991; Blunt et al. 1992; Blunt 1997; Piri and Blunt 2005a,b), and a
combination of effective-medium approximation and percolation theory
(Ghanbarian et al. 2016). Although such models have provided much insight
into the physics of multiphase flow in porous media, the way the morphology
of the pore space, i.e., its pore-size distribution, pore connectivity, and the
pore shapes, is represented is based on assumptions, approximations, and
sometimes even unreasonable simplifications.
Two- and three-dimensional (3D) imaging of porous media has been
becoming an important tool for their characterization because, depending on
their resolution, such images provide detailed information about the
morphology of the pore space. With advances in instrumentation, acquisition
of high-resolution 2D and 3D images by focused ion-beam scanning electron
microscopes (Lemmens et al. 2011), magnetic resonance imaging (Sheppard
et al. 2003; Sankey et al. 2009) and X-ray computed tomography (Arns et al.
2001, 2002; Porter et al. 2010; Iglauer et al. 2010; Berg et al. 2013;
Wildenschild et al. 2014) that can be used to display the fluid distributions
in 3D have become an essential part of studying various phenomena in
porous media. But, in addition to enabling their precise characterization,
high-resolution images of porous media have another advantage: provided
that their size is larger than the representative elementary volume (REV),
they may be used in direct numerical simulation of fluid flow and transport
in the porous media without any need to develop models of their pore space
that entails making assumptions and approximations. The REV is the
minimum size of a porous medium such that for sizes larger than the REV
the medium’s properties are independent of its dimensions.

38

Numerical simulation of the governing equations of fluid flow and
transport in the images of porous media is computationally very expensive.
Pores in natural and even many synthetic porous media have irregular
shapes. Their surface is rough, often representing a self-affine fractal
structure. The connectivity of the pores is also complex, giving rise to highly
tortuous flow paths. The computational grid that must represent all such
complexities faithfully must have high resolution, particularly near the
pores’ surface. Moreover, a computational grid with equal grid blocks is
highly inefficient, and increases unnecessarily the computational burden.
Use of irregular computational grids with higher resolution segments near
the pores’ surface and coarser ones away from it is a step in the right
direction, but the computations are still highly intensive. Simulation of flow
phenomena with the lattice-Boltzmann approach is even slower. For direct
image-based simulation of fluid flow and transport in porous media to
become the standard approach, particularly when one deals with more than
one fluid phase, the speed of the calculation must increase very significantly
without losing accuracy.
In this chapter we present a novel approach for image-based simulation of
multiphase flow in porous media. Instead of focusing on developing methods
that can solve the governing equations more efficiently, we focus on the
image itself. We propose a method whereby, taking advantage of extended
correlations in the image, we first process it in order to develop a coarser
image, and then carry out numerical simulation of two-phase flow in the
processed image. The processing of the image is done by curvelet
transformations (CTs), a powerful method developed by computer scientists
for image processing. We demonstrate that the new method increases the
speed of computations by two orders of magnitude or more, while preserving
the accuracy of the solution.

39

The rest of this chapter is organized as follows. In the next section we
explain how CTs are used to process the images in order to remove their non-
essential features. Next, the two-phase flow problem is described, followed
by the presentation of the results. We then discuss a few key aspects of the
proposed method. The chapter is summarized in the last section.

3.2 Spatial Correlations and Removal of Redundancy
from Images of Porous Media
A common characteristic of the morphology of porous media is extended
spatial correlations between the pores and grains. The same correlations
exist, of course, in the neighboring pixels (voxels in 3D) of the media’s images
and, therefore, the image contains redundant information that may be
removed. In other words, the foremost task for speeding-up image-based
simulation of fluid flow and transport in porous media is perhaps developing
a compressed representation of the image by removing the redundant
information. This is the main idea that we propose in this chapter for
speeding-up image-based simulation of multiphase flow in porous media: we
denoise the image and remove the redundant information, i.e., those features
of the morphology that are not essential to fluid flow.
The denoising is based on the fact that (Starck et al. 2002) the CCs
represent a measure of the local complexity of an image. In an image of a
porous medium the complexity is caused by the pores’ rough surface, which
is typically a self-affine fractal (see Sahimi 2011 for a comprehensive
discussion) containing the aforementioned correlations, their associated
irregular shapes. In other words, the CCs, particularly those associated with
the interface between the pores and the solid matrix and their neighborhood,
represent a measure of roughness or curviness of the interface. Because of

40

the extended correlations, however, even the pores are affected. Therefore,
the larger the CCs, the more significant are the curviness and the roughness.
Clearly, if the roughness is severe, its contribution to image-based
simulation of two-phase flow is significant. Thus, if we keep only the most
curved portions of the pore space and smoothen the rest, i.e., if we denoise
the image, a coarser computational grid should suffice for the simulations.
Denoising of images by the WTs was considered the ideal approach, until
its aforementioned shortcomings became clear, hence motivating the
development of the second generation of the CTs. While there are several
ways of exploiting the CTs to denoise or compress an image (Starck et al.
2002), we take the simplest approach that, as we demonstrate shortly,
suffices for our purpose. The approach represents what is usually referred to
as thresholding: After curvelet transforming of an image and computing the
CCs 𝐶 𝑗 ,𝑙 ,𝑘
, we normalize the coefficients by their largest value. We then
introduce a threshold 0 ≤ ε < 1, such that if 𝐶 𝑗 ,𝑙 ,𝑘 < 𝜀 , we set it to zero, hence
ignoring the local roughness and complexity in that particular neighborhood.
We, therefore, obtain a sparser representation of the image in the curvelet
space, as many of the CCs of the original image are set to zero. The choice of
the numerical value of ε depends on the computation time and resources that
one can afford and has access to. We shall return to this point shortly. We
then have two options:
(i) We compute the inverse CT of the thresholded image, and carry out
numerical simulation of two-phase flow in it, or
(ii) The governing flow equations are curvelet transformed and solved in
the CT space. Then, the inverse CT of the numerical results are computed
numerically.
The results presented in this chapter were obtained with option (i), but we
will also discuss the numerical procedure with (ii).

41

3.3 The Governing Equations and Details of Simulation
of Two-Phase Flow
We now describe the porous media whose images we have used in the
computations, and explain the computational procedure.

3.3.1 The Porous Media
We utilized 2D images of two distinct types of porous media. One was from
a sample of Mt. Simon sandstone (named for Mount Simon escarpment in
Eau Claire County, Wisconsin) at a depth of 6700 feet (Tahmasebi et al. 2017;
Bakhshian et al. 2018; Kohanpur et al. 2020). The formation is located at
verification well number 2 of a study site in Decatur, Illinois, where Illinois
State Geological Survey carried out a pilot injection study to better
understand the feasibility of full-scale carbon dioxide sequestration by
injecting it into the formation. A core plug from the formation was scanned
by micro-CT imaging to produce a series of grey-scale scans. We utilized one
2D slice of the image with 512
2
pixels and a porosity of 0.25.
The second image, also of size 512
2
pixels, was from a fossiliferous outcrop
carbonate, Mt. Gambier limestone in South Australia with porosity of 0.43.
We showed previously (Kohanpur et al. 2020) that the selected size of both
samples exceeds that of a REV. The physical sizes of the images of the
sandstone and carbonate rock were, respectively, 1.4 mm
2
and 2.76 mm
2
.
Figure 3.1 shows the images.

42

(a) (b)
Figure 3.1: Sandstone image (a), and Carbonate image (b).

3.3.2 The Governing Equations and the Numerical
Simulation Procedure
We simulated two-phase flow of water and oil in the images of the two porous
media. In both cases, water was injected into the (image) of the pore space,
saturated by oil. The contact angle for the water phase was 55◦. To carry out
the numerical simulations, we used the OpenFOAM program (Ubink, 1997),
an open source software that uses finite-volume discretization to solve the
continuity and the Stokes’ equations in the pore space. The computational
grid was unstructured with grid sizes of 2
3
and 1 pixel in the pores and
throats, and corner points between the solid boundaries, respectively. Thus,
the grid blocks in the throats and the corners were half the size in each
direction of those used in the pores. To do so, we used a modified OpenFOAM
mesh generator (Ubink 1999; Raeini 2013).

43

In order to carry out numerical simulations of two-phase flow in porous
media, one must upscale the microscopic dynamic and capillary pressures to
the Darcy-scale pressure, for which there are various methods. We used the
velocity-weighted average of the viscous forces, as well as the velocity-
weighted average of the pore-scale pressure gradient that has been shown
(Raeini et al. 2014) to produce numerical results that are in agreement with
experimental data most closely. The same equations with some modifications
can be used to calculate the macroscopic pressure drop in two-phase flow.
Thus, the method was used to calculate the total macroscopic pressure drop
∆Pα in the fluid phase α:

∆𝑃 𝛼 =−
1
Q
𝛼 ∫v.[𝛁 .(μ𝛁 v)]dV
𝛼 , α=1,2 (3.1)

where the density ρ and viscosity µ are given by

µ = φµ
1
+ (1 − φ) µ
2
, (3.2)
ρ = φ𝜌 1
+ (1 − φ) 𝜌 2
,
(3.3)

Here, Qα is the flow rate of phase α, v is the velocity field, Vα is the portion
of pore volume filled with phase α, µ
1
and µ
2
are viscosities of the two fluids,
and φ is the volume fraction of fluid 1 (water) in each grid block. More details
are given by Raeini et al. (2014).
The computational grid is Eulerian, for which various methods have been
developed in order to include the interfacial tension in it and reduce
numerical diffusion at the interface between the two fluids. Such methods
include the continuous surface stress (CSS) (Gueyffier et al. 1999),
continuous surface force (CSF), sharp surface force (SSF) (Francois et al.

44

2006), and filtered surface force (FSF) (Raeini et al. 2012). The CSF suffers
from spurious velocity in the flow field. The shortcoming is, to some extent,
controlled by the SSF method that controls the sharpness of the capillary
pressure by defining a sharpened indicator function, which is the same as
the volume fraction φ of fluid phase 1. Thus, it takes on a value of one in fluid
1 the wetting fluid - and zero in fluid 2 - the non-wetting fluid - and varies
smoothly between 0 and 1 in the interface region. A slightly modified
indicator function with a capillary force that, relative to the SSF approach,
is smoother is used in the FSF approach for the interface motion, which
compresses the transition area of the capillary pressure and, hence, resolves
the issue with the non-physical velocities that may arise. Then, the capillary
pressure transition area is only one grid block. Therefore, we utilized the
FSF formulation, since the morphology of the two porous media in which we
simulate two-phase flow is complex. The magnitude f of the interfacial force
f is given by

𝑓 =𝑘𝜎
2𝜌 𝜌 𝑤 + 𝜌 𝑛𝑤
𝛻𝜙 (3.4)

where σ is the interfacial tension, and κ is the curvature given by

𝑘 = − 𝛻 ·(
𝛻𝛼
|𝛻𝛼 |
) (3.5)

The governing equation, in addition to the usual mass conservation, is the
momentum equation,

𝜕 𝜕𝑡
(𝜌 v)=− 𝛁 𝑃 +𝛁 ·(𝜇𝛻 v)+𝐟 (3.6)

45

At each time step the volume fraction φ, which indicates the interface
position of the two-phases, is first updated, followed by the calculation of the
interfacial force f. Then, the velocity and pressure fields are computed based
on the Stokes’ and continuity equations. After the steady-state is reached,
the Darcy’s law for two-phase flow,

𝐯 α
=−
k
r,α
(S
α
)k
μ
α
𝛥𝑃
𝛼 𝛥𝑥
, α=1,2 (3.7)

is used to calculate the relative permeability Kr,α of fluid phase α as a
function of its saturation. As for the fluid properties, we used those given by
Oak et al. (1990) who also reported experimental data for the relative
permeabilities. They are, µ
1
= 1.05×10
−3
for water, µ
2
= 1.4×10
−3
for oil,
both in Pa.s, and 𝜌 1
= 𝜌 2
= 10
3
kg/m
3
for the densities. The interfacial tension
was, σ = 30 mN/m. Water was injected into the pore space with a velocity of
1.4 cm/s, while the pressure was atmospheric on the opposite face of the
injection surface. Thus, the capillary number was 𝐶𝑎 = 5 × 10
−5
.

3.4 Results
To compare the results computed with the original images and their denoised
counterparts, we first determined the resolution (number of grid blocks) of
the computational grids in the images that yields macroscopic Darcy
permeability that no longer changes if the resolution is increased further.
Figure 3.2 presents the results. Thus, we carried out a series of simulations,
one with the original image with a computational grid indicated by Figure
3.2, and a series of simulations with coarsened images after they had been
denoised by the CT.

46

(a) (b)
Figure 3.2: Absolute permeability vs N (number of grid points), (a) sandstone
and (b) carbonate image, in order to identify the size of the REV.

Figures 3.3 and 3.4 compare, respectively, the evolution with time of the
water and oil distributions in the original image of the sandstone and
carbonate rock with those obtained with the curvelet-thresholded images
with the thresholds ε = 0.5, 0.7 and 0.9. The essential features of the
distributions are completely similar in all the images. Numerically, the
differences between the saturations in the original images and the
thresholded ones is less than 5 percent.
To make quantitative comparison between the results obtained with the
original images and their thresholded versions, we present in Figure 3.5 the
water saturation profiles along the macroscopic direction of flow at the
breakthrough time 𝑡 𝐵 , where 𝑋 = 𝑥 /𝐿 , with 𝐿 being the linear size of the
images. The agreement in all cases is excellent.

0
10
20
30
40
50
0 30 60 90 120 150
Absolute Permeability (mD)
N (number of grid points) x 𝟏𝟎
𝟒 0
2
4
6
8
0 30 60 90 120 150
Absolute Permeability (mD)
N (number of grid points) x 𝟏𝟎
𝟒

47

(a) (b) (c)

(d) (e) (f)

(g) (h) (i)
Figure 3.3: Time evolution of wetting phase distribution for the (a-c)
Sandstone original image, (d-f) curvelet threshold ε = 0.5 image and (g-i)
curvelet threshold ε = 0.9 image at time 30 (ms), 60 (ms) and at breakthrough
time, respectively. The invading fluid (water) and the defending fluid (oil)
shown in red and blue, respectively.

48

(a) (b) (c)

(d) (e) (f)

(g) (h) (I)
Figure 3.4: Time evolution of wetting phase distribution for the (a-c)
Carbonate original image, (d-f) curvelet threshold ε = 0.5 image and (g-i)
curvelet threshold ε = 0.9 image at time 30 (ms), 60 (ms) and at breakthrough
time, respectively. The invading fluid (water) and the defending fluid (oil)
shown in red and blue, respectively.

49

Figure 3.6 presents the evolution with the time 𝑡 𝐷 of the water saturation
𝑆 𝑤 in the original images of the two porous media, and compares it with those
computed with the thresholded images, where 𝑡 𝐷 = 𝑡 /𝑡 𝐵 , with 𝑡 being the
actual time. Once again, the agreement is excellent in all the cases.

Figure 3.5: Saturation profile of the wetting fluid along the flow direction (x direction)
at breakthrough time in the curvelet-transformed images of the Sandstone and carbonate
porous medium with two thresholds ε, and their comparison with those in the original
image. X is the normalized front distance from the inlet (X = x/L, where L is the length
of the sample in the x direction).

0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
Sw
X
Sandstone
Original
= 0.5
= 0.9
ε
ε
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
Sw
X
Carbonate
Original
= 0.5
= 0.9
ε
ε

50

Figure 3.6: Saturation of the wetting fluid as a function of dimensionless time ,t D
(normalized by the breakthrough time) in the curvelet-transformed images of the (a)
sandstone and (b) Carbonate porous medium with two thresholds ε, and their comparison
with those in the original image.

Figures 3.7 and 3.8 compare the distribution of local flow velocities in the
two original images with those computed with the thresholded images. The
velocities were normalized with respect to their maximum value. For
velocities up to the maximum of the distributions, the agreement between all
the distributions is perfect. Beyond the maximum, the agreement is good,
with all the distributions having the same shape and almost identical tails.
Figures 3.9 and 3.10 present the comparison between the relative
permeabilities Kr of the two fluid phases, computed in the original images of
the two porous media, and those calculated with the thresholded images. The
agreement is excellent in all the cases, indicating the accuracy of the
computations with the thresholded images.

0
0.1
0.2
0.3
0.4
0.5
0.6
0 0.2 0.4 0.6 0.8 1
Sw
t
D
Sandstone
Original
= 0.5
= 0.9
ε
ε
0
0.1
0.2
0.3
0.4
0.5
0 0.2 0.4 0.6 0.8 1
Sw
t
D
Carbonate
Original
= 0.5
= 0.9
ε
ε

51

Figure 3.7: Probability density function (PDF) of the normalized velocity
(u/Umax, where Umax is the maximum flow velocity) in the curvelet-transformed
images of the Sandstone porous medium with two thresholds ε, and their comparison
with those in the original image.

Figure 3.8: Probability density function (PDF) of the normalized velocity
(u/Umax, where Umax is the maximum flow velocity) in the curvelet-transformed
images of the Carbonate porous medium with two thresholds ε, and their comparison
with those in the original image.
0.2
0.4
0.6
0.8
1
1.2
1.4
0 0.2 0.4 0.6 0.8 1
PDF
Normalized Velocity
Sandstone
Original
= 0.5
= 0.9
ε
ε
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
0 0.2 0.4 0.6 0.8 1
PDF
Normalized Velocity
Carbonate
Original
= 0.5
= 0.9
ε
ε

52

Figure 3.9: Water-oil relative permeability curves obtained from Sandstone original
image (continuous curves) and curvelet thresholded images (dashed curves).

Figure 3.10: Water-oil relative permeability curves obtained from Carbonate original
image (continuous curves) and curvelet thresholded images (dashed curves).
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
Kr
Sw
Sandstone
Original, Kro
Original, Krw
= 0.5, Kro
= 0.5, Krw
= 0.9, Kro
= 0.9, Krw
ε
ε
ε
ε
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
Kr
Sw
Carbonate
Original, Kro
Original, Krw
= 0.5, Kro
= 0.5, Krw
= 0.9, Kro
= 0.9, Krw
ε
ε
ε
ε

53

Finally, an important quantity in two-phase flow in porous media is the
oil residual saturation Sor, i.e., its saturation at the water breakthrough
point. In Table 3.1 we compare the computed residual saturations in the
original images and their thresholded counterparts. In both the sandstone
and carbonate porous sample the maximum difference between the residual
saturations in the original image and the thresholded ones with a threshold
as high as ε = 0.9 is 4.5 percent.

3.5 Discussion
A few important aspects of the results should be discussed in detail.

3.5.1 Accuracy and Speed-up of the Computations with the
Denoised Images
We point out that the computational cost of curvelet transformation of a
digitized image of size n × n is O(n
2
logn). But the entire computation time
for curve transforming of an image, thresholding it, and obtaining the
inverse CT of the coarser image represents only a small fraction of the total
computation time for simulating two-phase flow in the image.
In Table 3.1 we compare the details of the computations: the total number N
of grid blocks, the CPU times for all the cases, and the oil residual
saturations.

54

Table 3.1 Comparison of the number of grid blocks N, the computation
times, and the oil residual saturation Sor, computed with the original and
curvelet transform-thresholded images.

Two aspects of Table 3.1 are noteworthy. One is that the computed results
with the denoised images with a threshold as high as ε = 0.9 are still very
accurate, and differ from those of the original images by less than 5 percent,
which is well within the possible inaccuracies that the images themselves
may contain. The second aspect is the increase in the speed of the
computations for both types of porous media. With a threshold ε = 0.5, for
which the difference between the computed results with the original and
denoised images is about 3.5 percent, the speed-up factor is between 105 and
111, better than two orders of magnitude improvement. When the threshold
is increased to 0.9, the speed-up factor is between 144 and 151.

Sandstone N Threshold ε Sor Time (CPU sec)
Original image 859973

0.2758 19241
Curvelet-transformed image 12646 0.5 0.2861 173
Curvelet-transformed image 8187 0.7 0.2879 139
Curvelet-transformed image 6871 0.9 0.2893 127
Carbonate
Original image 879103
0.2855
37827
Curvelet-transformed image 13833 0.5 0.2916 359
Curvelet-transformed image 9116 0.7 0.2962 283
Curvelet-transformed image 7727 0.9 0.2984 261

55

(a) (b) (c)
Figure 3.11: Small portion of the (a) fine grid, (b) the coarsened grid using ε = 0.5 and
(c) the coarsened grid using ε = 0.9.

The reason is removal of the redundant and irrelevant information from
the images that make it possible to use a coarser but still accurate
computational grid. Figure 3.11 shows the structure of the grid in a small
portion of the original image of the sandstone, and compares it with those in
two of its coarsened versions. First, note that the grid blocks in Figure 3.11(a)
next to the solid walls are highly resolved, because the surface is very
complex. Second, after the images are denoised, one no longer has to use a
highly resolved computational grid everywhere, but only in selected areas,
which are determined by the CT, including portions of the pores since they
are correlated with the regions near the solid surface.

3.5.2 Computations in the Curvelet Space
As already mentioned, one may carry out the simulations in the curvelet
space by the following procedure. After the images are curvelet-transformed,
no threshold is applied, and numerical simulations of two-phase flow are

56

carried out in the curvelet-transformed images, with the numerical results
inversed curvelet-transformed back to the real space. To carry out such
simulations, one must take the CT of the governing equations. Since digital
CTs operate in the frequency domain, one substitutes the Fourier transforms
𝑓 ̂
(ω) and U(R
θl
𝜔 ) in Eq. (1.6).
We first rewrite the second-order term of the momentum equation, Eq.
(3.6), as, ∇ · (µ∇v) = ∇µ · ∇v + µ∇
2
v. The curvelet transform of the gradient ∇F
of a function F is given by

1
𝑐𝑜𝑠 𝜃 𝑙 𝜕 𝐹 𝑗 ,𝑙 ,𝑘
𝜕 𝑦 1
𝑒 1
+
𝜕 𝐹 𝑗 ,𝑙 ,𝑘
𝜕 𝑦 2
𝑒 2
, (3.8)

where 𝐹 𝑗 ,𝑙 ,𝑘 is the CCs associated with F, (𝑦 1
, 𝑦 2
) are the coordinates in the
curvelet space, and 𝑒 1
and 𝑒 2
are unit vectors associated with 𝑦 1
and 𝑦 2
. The
CT of the Laplacian ∇
2
G of a function G is given by

(1+𝑡𝑎𝑛 2
𝜃 𝑙 )
𝜕 2
𝐺 𝑗 ,𝑙𝑘
𝜕 𝑦 1
2
+
𝜕 2
𝐺 𝑗 ,𝑙𝑘
𝜕 𝑦 2
2
−2𝑡𝑎𝑛 𝜃 𝑙 𝜕 2
𝐺 𝑗 ,𝑙𝑘
𝜕 𝑦 1
𝜕 𝑦 2
, (3.9)

with 𝐺 𝑗 ,𝑙 ,𝑘 being the CC associated with G. One then solves the curvelet-
transformed governing equation in the curvelet-transformed image by
discretizing it and defining the scale j, direction 𝑙 , and translation parameter
𝑘 = (𝑘 1
,𝑘 2
) ∈ Z
2
, with (𝑦 1
, 𝑦 2
) = (2
−𝑗 𝑘 1
,2
−𝑗 /2
𝑘 1
). The inverse CT of the
numerical solution is then computed.
It has been shown (Aljasmi and Sahimi 2020) that if this procedure is
followed, the computations will be speeded up by a factor of four or better,
even if the image in the curvelet space is not processed. Clearly, then, if the
image is denoised, and the computations are still carried out in the CT of the
denoised image, the speed-up factor will be very large, comparable with what

57

is presented in Table 3.1. The reason for the higher efficiency even without
thresholding is that a coarser computational grid suffices in the curvelet
space, because many of the CCs of the image are either very small and, thus,
do not make significant contributions, or that due to the correlations they
are close to those of their neighbors, indicating a relatively smooth local
environment, hence allowing use of larger grid blocks. As such noisy pixels
may be due to artifacts during image processing, their removal would
actually lead to more accurate simulation of flow and transport processes
than what would obtain with the original image.

3.5.3 Simulation in Three-Dimensional Images
Although the CTs were originally developed for 2D images, their extension
to 3D systems has also been proposed and developed (Ying et al. 2005;
Woiselle et al. 2010). Thus, one may use such 3D CTs in the computations
with 3D images. Alternatively, one may divide a 3D image into a large
number of 2D slices, apply the methodology described here to each 2D image,
and then stack together the curvelet-transformed 2D images to reconstruct
a 3D coarser image. Care must be taken to ensure that after curve-
transforming and thresholding, the transition from one 2D slice to the next
is seamless, and does not generate any artifacts. Simulations of two-phase
flow in 3D images of porous media using the CTs, the results of which will
be reported in Chapter 4 and 5.

58

3.6 Summary
Advances in instrumentations, together with rapid increase in the speed of
computations, have made it possible to begin simulation of multiphase flow
in 2D and 3D images of heterogeneous porous media, hence avoiding
altogether the assumption and approximations that are involved in the
development of models of porous media. However, as image-based simulation
of multiphase flow in heterogeneous porous media becomes more practical,
one must still address the issue of its high computational cost. In this chapter
we described a new methodology to address the issue, which is based on
curvelet transforming the image and denoising it by removing those details
of the image that do not contribute significantly to the simulation. The
computations in the less noisy, and hence more coarsened images yield
results that are as accurate as those in the original ones, but with a speedup
factor of 100-150 or higher. Thus, the methodology represents a significant
step toward achieving the ultimate goal of such simulations, namely, making
image-based simulation of multiphase flow in heterogeneous porous media a
standard practice.

59

Chapter 4

Fast simulation of two-phase flow in three dimensional
digital images of heterogeneous porous media using
multiresolution curvelet transformation

4.1 Introduction
Although the physics of two-phase flow in heterogeneous porous media is
reasonably well understood (Helmig and Schulz, 1997; Sahimi, 2011; Blunt,
2017), its efficient simulation is still the subject of active research. In
particular, due to societal interest in geological sequestration of CO2, there
is increasing incentive for studying two-phase flow of CO2 and brine in such
porous media as depleted oil or gas reservoirs (Friedlingstein and Solomon,
2005; Metz et al., 2005; Nordbotten and Celia, 2011). Such simulations are
either based on the classical continuum formulation of two-phase flows
(Bear, 1972), which necessitates supplying the simulator with various
properties of the pore space, such as its porosity and permeability, or entails,
as the first step of the simulation, developing or assuming a model of the pore
space that represents faithfully its morphology and the complexities that are
associated with it.
Up until relatively recently, the most realistic model of a porous medium
at laboratory scale was the pore-network model (PNM) in which the pore
space, extracted from an image of the medium, was represented by a network

60

of interconnected pore throats and pore bodies, the effective sizes of which
were selected according to some representative statistical distributions that
can, in principle, be measured or obtained from the image. Based on
representation of a pore space by the PNMs, various approaches were
developed in order to simulate two-phase flows in porous media. They
included a variety of percolation models, ranging from random (see, for
example, Sahimi et al., 1986), to correlated (see, for example, Kantzas and
Chatzis, 1988) and invasion percolation (for a simple introduction to invasion
percolation see Ebrahimi, 2010), site-bond invasion percolation (Sahimi et
al., 1998) that considers the effect of both pore bodies and pore throats on the
invasion process, detailed pore-to-pore simulations (Valvatne and Blunt,
2002; Piri and Blunt, 2005), and a combination of effective-medium
approximation and random percolation (Ghanbarian et al. 2016). These
models have been described and discussed comprehensively by Sahimi
(2011), Blunt (2017), and Hunt and Sahimi (2017) and, therefore, need not
be described here.
There have been tremendous advances in instrumentation, which have
made it possible to obtain high-resolution two- or three-dimensional (3D)
images of porous media at laboratory scale. Such images have been becoming
important tools for characterization of porous media, as they provide,
depending on their resolution, detailed information about the morphology of
the pore space. The imaging techniques include focused ion-beam scanning
electron micro-scopes (Lemmens et al., 2011), magnetic resonance imaging
(Sheppard et al., 2003) and X-ray computed tomography (Kantzas et al.,
1992; Arns et al., 2001, 2002; Porter et al., 2010; Iglauer et al., 2010; Berg et
al., 2013; Wildenschild et al., 2014). The images have been used to display
fluid distributions in 3D porous media (Sheppard et al., 2003), even under
reservoir conditions (Andrew et al., 2015).

61

With advances in computational power, high-resolution images of porous
media offer the opportunity for precise simulation of fluid flow and transport
in porous media. If the physical size of the images is larger than the
representative elementary volume (REV) - the minimum volume of a porous
medium for its properties to be independent of the dimensions – they may be
used directly in numerical simulation of fluid flow and transport in the
porous media (Arns et al., 2001, 2002; Piller et al., 2009; Tahmasebi et al.,
2017), as well as adsorption and deformation (Bakhshian et al., 2018).
Despite the great increase in computational power, numerical simulation
of fluid flow and transport in the images of porous media is still
computationally very expensive. Natural and many synthetic porous media
have pores with irregular shapes, with their surface being rough and often a
self-affine fractal. The pores’ connectivity is also a spatially distributed
property, giving rise to highly tortuous flow and transport paths. Thus, one
needs a high-resolution computational grid that represents such
complexities accurately.
Simulation of flow in porous media with the lattice-Boltzmann (LB)
algorithm (Huang et al.,2015) is still slower than direct numerical simulation
based on the Stokes' equation, although progress has recently been made in
speeding it up. An et al. (2020) developed a method that integrates the GPU
- accelerated volumetric LB method (VLBM) with an upscaling scheme to
solve the pore-scale two-phase flow. They utilized the LB simulator to solve
both the level-set equation for image segmentation and the governing
equations for multiphase flow dynamics. The method accelerates the
standard LB simulations significantly. Aside from the LB-based methods, for
image-based direct simulation of multiphase flow and transport in porous
media, based on the Navier-Stokes equation, to become the standard
approach, the speed of the calculations must increase very significantly

62

without losing accuracy (see, for example, Etemad et al., 2017, for a
significant effort in this direction).
In this chapter we describe an approach for image-based simulation of
multiphase flow in 3D porous media. The main idea is to focus on the image,
process it in order to denoise it and remove the complexities that do not
contribute significantly to flow and transport in the pore space, and then
carry out numerical simulation of multiphase flow in the coarse image. The
image is denoised by curvelet transformations (CTs), a powerful method
developed for image processing. As we show in this chapter, the method
increases the speed of computations very significantly, between one and
more than two orders of magnitude, while preserving the accuracy of the
solution.
The rest of this chapter is structured as follows. In Section 2 we explain
the method for denoising the images. The two-phase flow problem is
described in Section 3, followed by the presentation and discussion of the
results in Section 4. The efficiency of the computation is discussed in Section
5, while the chapter is summarized in the last section.

4.2 Denoising of Images of Porous Media
Having set up the CTs in chapter 1, we now describe their use for processing
of 2D and 3D images of porous media.

4.2.1 Two-dimensional images
The idea for processing images of porous media by the CTs (or by the WTS,
for that matter) and denoising them exploits the physical fact that the
morphology of such media contains extended spatial correlations between

63

the pores and grains (see, for example, Knackstedt et al., 2001). Thus, the
same correlations exist between the pixels (voxels in 3D) in the images of
porous media. The existence of the correlations implies that images contain
redundant information, and that their removal will not damage the essential
information that the images carry. Once the redundant information is
removed, simulation of flow and transport in the resulting coarser images is
accelerated very significantly, as we demonstrate shortly.
The denoising exploits the fact that (Starck et al., 2002) the CCs represent
a measure of the local complexity of an image, which in a porous medium is
due to the pores' rough surface, which is typically a self-affine fractal (see
Sahimi 2011 for a comprehensive discussion) and its correlation with the
surrounding pore space. In other words, the CCs represent a measure of
roughness or curviness of the pore-grain interface and, due to the extended
correlations, the pores next to the interface are also affected. Therefore, the
larger the CCs, the more significant is the curviness and, hence, the
roughness. Clearly, if the roughness is severe, its contribution to image-
based simulation of multiphase flow should be significant. Thus, if we keep
only the most curved portions of the pore space, i.e., if we denoise the image,
a coarser computational grid should suffice for the simulations.
There are several ways of denoising or compressing an image (Starck et
al., 2002), but we choose to take the simplest approach that, as we
demonstrate shortly, suffices for our purpose.
The approach, which is referred to as thresholding in image processing,
consists of the following steps: (i) Curvelet transform the image and compute
the CCs 𝐶 𝑗 ,𝑙 ,𝑘 . (ii) Normalize the CCs by their largest value. (iii) Introduce a
threshold 0≤𝜀 <1, such that if 𝐶 𝑗 ,𝑙 ,𝑘 < 𝜀 , it is set to zero. In this way, the
algorithm generates a sparser representation of the image in the curvelet
space, as many of the CCs of the original image are set to zero, while it still

64

contains all the essential information. The numerical value of 𝜀 depends on
the computation time and resources that one can afford and has access to.

4.2.2 Three-dimensional images
The CTs were originally developed for 2D images. Their extension to 3D
initially encountered some difficulties associated with curved surfaces, but
curvelets for some 3D systems have also been developed (Ying et al., 2005;
Woiselle et al., 2010). But, for reasons that will be explained in Section 5, in
the present research we use an alternative approach, namely, we divide a n
× n × n 3D image into n 2D images of size n × n, apply the methodology
described in subsection 4.2.1 to each 2D image, and then stack together the
curvelet-transformed 2D denoised images to reconstruct a 3D image. The
orientations of the 2D images relative to the original 3D image does not play
a significant role, unless the pore space is strongly anisotropic. If strong
anisotropy is present, one may use a CT-based anisotropic feature extraction
method (Shinde et al., 2017). In the present work we take the 2D images in
the planes perpendicular to the direction of macroscopic flow. With the
development of fast 3D CTs, the issue of the orientation of the 2D slices
relative to the 3D images becomes irrelevant.
Care must, however, be taken to ensure that after curvelet-transforming
and thresholding, the transition from one processed 2D slice to the next is
seamless, and does not generate any artifacts. To do this, we used an open-
access image processing package that does the segmentation and the
stacking process. The software is called Fiji (Schindelin et al., 2012), which
generates directly 3D images from a series of multiple 2D images. It is also
designed to automatically address issues that may arise from image

65

thresholding, which is why the software is a very power tool (for more details
see, https://imagej.net/Fiji).

4.3 Simulations in the curvelet or real space?
Once the coarser image is developed, one has two options:
(i) One computes the inverse CT of the thresholded image, and carries out
numerical simulation of two-phase flow in it, or
(ii) one curvelet transforms the governing flow equations, solves them in
the CT space and then, computes the inverse CT of the numerical results.
The results presented in this chapter were obtained with option (i). It can
be shown (Aljasmi and Sahimi 2020) that if procedure (ii) is followed, the
speed-up factor for computations will be a factor of four or better, even if the
images in the curvelet space is not denoised. Clearly, then, if the image is
denoised, and the computations are still carried out in the CT of the denoised
image, the speed-up factor will be quite large.
The reason that one has better efficiency even without thresholding is
that, in the CT space a coarser computational grid suffices, because many of
the CCs of the image are either very small and, thus, do not make significant
contributions, or that due to the correlations they are close to those of their
neighbors, indicating a relatively smooth local environment, hence allowing
use of larger grid blocks.

4.4 Details of the Simulations
We now describe the porous medium whose image we have used in the
computations, and explain the computational procedure.

66

4.4.1 The porous medium
We utilized a 3D image of Berea sandstone, which is considered a benchmark
for testing various computational approaches to simulating multiphase flow
in porous media, and has been used extensively in the past (Valvatne and
Blunt, 2004; Ramstad et al., 2012; Raeini et al., 2014). The digital and
physical sizes of the sample are, respectively, 200
3
voxels and 2.14
3
mm
3
. The
porosity and absolute permeability of the sample are, respectively, 0.196 and
1368 mD. Figure 4.1 presents the image of the sandstone and its 3D coarse
version. The size of the image is larger than its corresponding REV.

(a) (b)
Figure 4.1: A 200
3
Berea sandstone sample, (a) original and (b) curvelet transformed
sample with ε = 0.8 used for single and two-phase flow simulation.

67

4.4.2 The governing equations and the numerical
simulation procedure
We first carried out simulation of single-phase flow in the image of the Berea
sandstone, both in the original one and in the processed images in order to
not only estimate the absolute permeability of the pore space, but also
determine the size of the computational grid that yield reliable estimate (see
below). We then simulated two-phase flow of CO2 and brine in the same
images.
Simulating two-phase flow in porous media requires upscaling the
microscopic dynamic and capillary pressures to the Darcy scale. While there
are several methods for the upscaling, we used the velocity-weighted
averages of both the viscous forces and the pressure gradient that have been
demonstrated (Raeini et al., 2014) to yield results that are in agreement with
experimental data. As for macroscopic pressure drop ΔP in fluid phase of
the two-phase system, we utilized a slightly modified versions of the same
equations, which is given by

∆𝑃 𝛼 =−
1
Q
𝛼 ∫v.[𝛁 .(μ𝛁 v)]dV
𝛼 , α=1,2, (4.1)

where the density ρ and viscosity µ are given by

µ = 𝜑 µ
1
+ (1 − 𝜑 )µ
2
, (4.2)
𝜌 = 𝜑 𝜌 1
+ (1 − 𝜑 )𝜌 2
,
(4.3)

Here, Q is the flow rate of phase , v is the velocity field, V is the portion of
pore volume filled with phase , 𝜇 1
and 𝜇 2
are viscosities of the two fluids,
and 𝜙 is the volume fraction of fluid 1 (brine) in each grid block.

68

Since the computational grid is Eulerian, one needs a method for including
the interfacial tension in grid, so as to reduce numerical diffusion at the
interface between the two fluids. Several methods have been proposed for
this purpose, such as the continuous surface stress (CSS) (Gueyffier et al.,
1999), continuous surface force (CSF), sharp surface force (SSF) (Francois et
al., 2006), and filtered surface force (FSF) (Raeini et al., 2012). But the CSF
method may generate spurious fluid velocity field that, to some extent, is
controlled by the SSF method. To do so, one introduces in the SSF method a
sharpened indicator function 𝜙 of fluid phase 1 that takes on a value of one
in fluid 1 (brine) - the wetting fluid - and zero in fluid 2 (CO2) - the non-
wetting fluid - and varies smoothly between the two in the interface region.
The approach is capable of controlling the sharpness of the capillary
pressure. The FSF approach utilizes a slightly modified indicator function
for the interface motion with a capillary force that, compared with the SSF
method, is smoother. This compresses the transition zone of the capillary
pressure and, hence, resolves the issue with the spurious velocities that may
arise. Then, the transition area for capillary pressure is only one grid block.
Therefore, we utilized the FSF formulation. The magnitude f of the
interfacial force f is given by

𝑓 =𝑘𝜎
2𝜌 𝜌 𝑤 + 𝜌 𝑛𝑤
𝛻𝜙 (4.4)

where σ is the interfacial tension, and κ is the curvature given by

𝑘 = − 𝛻 ·(
𝛻𝛼
|𝛻𝛼 |
) (4.5)

69

The governing equation, in addition to the usual mass conservation, is the
momentum equation,

𝜕 𝜕 𝑡 (𝜌 v)=− 𝛁 𝑃 +𝛁 ·(𝜇𝛻 v)+𝐟 (4.6)

The indicator function 𝜙 , indicating the position of the interface between the
two fluids, is updated at each time step. This is followed by the calculation of
the interfacial force f. Then, the velocity and pressure fields are computed by
solving the Stokes’ and continuity equations. At steady state Darcy’s law for
two-phase flow,

𝐯 α
=−
𝑘 𝑟 ,𝛼 𝑘 𝜇 𝛼 𝛥𝑃
𝛼 𝛥𝑥
, 𝛼 =1,2 (4.7)

is utilized for calculating the relative permeability 𝑘 𝑟 ,𝛼 of fluid phase as a
function of its saturation, where K is the absolute permeability. As for the
fluid properties, we used the typical values (Kohanpur et al., 2020): 𝜈 1
= 1.0
× 10
−6
for brine, 𝜈 2
= 1.0 × 10
−7
for CO2 (ν is the kinematic viscosity), both
in m
2
/s, while the densities were 𝜌 1
= 𝜌 2
= 10
3
kg/m
3
. The interfacial tension
was assumed to be, σ = 30 mN/m. The contact angle for the brine was 0 ◦.
The OpenFOAM program (Ubink, 1997) and its modified version (Raeini,
2013), which are open-source softwares that use finite-volume discretization
to solve the continuity and the momentum equation in the pore space, were
utilized for carrying out the numerical simulations. For the two-phase flow
simulation two grid blocks per pixel were used. CO2 was injected into the
pore medium with a speed of 4.65 ×10
−2
m/s at one face, so that the capillary
number was Ca = 1.55 ×10
−4
, and a fixed atmospheric pressure was applied
at the opposite face. The other four faces were assumed to be impermeable.

70

All the computations were carried out using a HP Spectre x360 computer
with a speed of 2GHz and 16 GB of memory.

4.5 Results
As pointed out earlier, we carried out simulation of both single- and two-
phase flows. In what follows we present and discuss the results.

4.5.1 One-phase flow
We first carried out a series of simulations of single-phase flow, varying the
resolution of the computational grid (the number of grid blocks), in order to
determine the grid size that yields a permeability independent of the grid’s
resolution. Figure 4.2 presents the results, which indicates that a grid with
about 3.34 × 10
6
blocks is necessary to obtain a stable value of the
permeability K, independent of grid resolution. The calculated K turned out
to be K = 1368 mD. We utilized the same grid in the simulation of two-phase
flow in the original image without any processing by the CT; see below.
Two thresholds, ε = 0.5 and 0.9, were utilized to denoise the image and
carry out simulation of single- and two-phase flows. Figure 4.3 compares the
pressure distributions in the original image and in the two coarser ones. All
the essential features of the three distributions are identical. The computed
permeabilities are K = 1443 and 1467 mD for ε = 0.5 and 0.9, respectively.
Even the latter estimate with a very high threshold of 0.9 is only 7.3 percent
higher than the estimate of K computed with the original image.

71

Figure 4.2: Sample-size dependence of the absolute permeability of the Berea sandstone
image in order to identify the size of the REV.

(a) (b) (c)
Figure 4.3: The pressure distribution in the pore volume at the end of the simulation for
(a) original, (b) curvelet threshold image ε = 0.5 and (c) curvelet threshold image ε = 0.9.

1000
1200
1400
1600
1800
0 100 200 300 400 500 600
Absolute Permeability (mD)
N (number of grid points) x 𝟏𝟎
𝟒

72

4.5.2 Two-phase flow
Figure 4.4 compares the evolution with time of the spatial distribution of
CO2 and brine in the original image of the sandstone with those obtained
with the CT-denoised images with the two thresholds. Similar to the
pressure distribution in single-phase flow, all the important features of the
distributions are completely similar in all the images. More quantitatively,
the differences between the saturations in the original images and the
thresholded ones is less than 7 percent.
In Figure 4.5 we compare the saturation profile of CO2 along the
macroscopic direction x of flow at the breakthrough time t
B
, where X = x/L,
with L being the linear size of the image. The saturation was computed for
each plane at position X perpendicular to the direction of macroscopic flow.
The agreement is excellent.
Figure 4.6 presents the evolution with time 𝑡 𝐷 of CO2 saturation SCO2 in
the original image of the porous medium, and compares it with those
computed with the thresholded images, where 𝑡 𝐷 = 𝑡 /𝑡 𝐵 , with 𝑡 being the
actual time. At short times the three profiles are identical. At longer times,
all the way to the breakthrough time, the maximum difference between the
profile with ε = 0.9 and the original image is about 3 percent.

73

(a) (b)

(b) (d)

(e) (f)

Figure 4.4: Time evolution of non-wetting phase distribution for the Berea sandstone
(a-b) original image, (c-d) curvelet threshold ε = 0.5 image and (e-f) curvelet threshold
ε = 0.9 image at time 60 (ms) and at breakthrough time, respectively. CO2 shown in red.

74

Figure 4.5: Saturation profile of the non-wetting fluid along the flow direction (x
direction) at breakthrough time in the curvelet-transformed images of the Berea
sandstone porous medium with two thresholds ε, and their comparison with those in the
original image. X is the normalized front distance from the inlet (X = x/L, where L is
the length of the sample in the x direction).

Figure 4.6: Saturation of the non-wetting fluid as a function of dimensionless time ,tD
(normalized by the breakthrough time) in the curvelet-transformed images of the Berea
sandstone porous medium with two thresholds ε, and their comparison with those in the
original image.
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
Snw
X
Original
= 0.5
= 0.9
ε
ε
0
0.1
0.2
0.3
0.4
0.5
0.6
0 0.2 0.4 0.6 0.8 1
Snw
t
D
Original
= 0.5
= 0.9
ε
ε

75

Figure 4.7 makes a comparison between the relative permeabilities Kr of
the two fluids that were computed in the original image of Berea sandstone
and those calculated with the thresholded images. The agreement is
excellent in all the cases, indicating the accuracy of the computations with
the thresholded images.

Figure 4.7: Drainage relative permeabilities of the CO2-water system in a Berea
Sandstone sample, original image (continuous curves) and curvelet thresholded images
(dashed curves).

To understand the excellence of the accuracy of the computed relative
permeabilities, we present in Figure 4.8 the distribution of local flow
velocities in the original image of the sandstone with those computed with
the thresholded images, where we have normalized the velocities by their
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
Kr
Sw
Original, Krg
Original, Krw
= 0.5, Krg
= 0.5, Krw
= 0.9, Krg
= 0.9, Krw
ε
ε
ε
ε

76

maximum value. Figure 4.8 indicates that for velocities up to the maximum
of the distributions, the agreement between all the distributions is perfect.
Beyond the maximum, the agreement is good, with all the distributions
having the same shape and almost identical tails.

Figure 4.8: Probability density function (PDF) of the normalized velocity
(u/Umax, where Umax is the maximum flow velocity) in the curvelet-transformed
images of the Sandstone porous medium with two thresholds ε, and their comparison
with those in the original image.

An important quantity in two-phase flow in porous media is the residual
saturation 𝑆 𝑏𝑟
of the displaced fluid (brine in the present simulation), its
saturation at the CO2 breakthrough point. In Table 4.1 we compare the
computed residual saturations in the original image and its thresholded
counterparts. The maximum difference, computed with ε = 0.9, between the
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
0 0.2 0.4 0.6 0.8 1
PDF
Normalized Velocity
Original
= 0.5
= 0.9
ε
ε

77

residual saturations in the original image and the thresholded one is 7.5
percent.

4.6 Speed-up of the computations
We used 2D CT in order to denoise the 3D image, for which the
computational cost of a digitized image of size n × n is O(n
2
log n). For the
2D slices of size 200
2
pixels the CPU time is very small, but since we denoise
a large number of them, the computation time is relatively significant (see
below). Stacking the 200 denoised slices took only 82 CPU seconds.
In Table 4.1 we compare the details of the computations: the total number
N of grid blocks, the CPU times for all the cases, and the brine residual
saturations.

Table 4.1: Comparison of the number of grid blocks N, the computation
times, and the brine residual saturation Sbr, computed with the original and
curvelet transform-thresholded images. The CPU times represent those for
simulating the two-phase flow.

We note that even the computed results with the denoised image with a
threshold ε = 0.9 are still very accurate, and differ from those of the original
images by at most 7 percent. Moreover, the speed-up factors for the
computations with the thresholds ε = 0.5 and 0.9 are, respectively, 263 and
Sandstone N Threshold ε Sbr Time (CPU sec)
Original image 3346168

0.3883 144231
Curvelet-transformed image 15005 0.5 0.4086 547
Curvelet-transformed image 11996 0.9 0.4172 489

78

294. These speed-up factors are for simulating the two-phase flow problem.
If we add the CPU time for denoising the 2D slices and reconstructing a
coarser 3D image to the CPU time for simulating the flow problem, then the
speed-up factor is about one order of magnitude, a factor of 7-10, which is
still very significant. Note, however, that the denoising is done only once,
while the denoised image can be used in the simulation of a variety of
processes in the same image.
The reason for the very significant speed-up is removal of the redundant
information from the image, which makes it possible to use a coarser
computational grid. After the image is denoised, one no longer must use a
highly resolved computational grid everywhere, but only in selected areas
that are determined by the CT. This includes parts of the pores, as the
difference between their CCs and those of the neighbors and in the solid
phase is large and, moreover, they are also correlated with the areas next to
the interface between the pores and the matrix.
Can the efficiency of image-based simulation of two-phase flows in porous
media be increased by other methods? It can, of course, be increased, but we
are not aware of any method that uses a single processor and achieves the
order of the speed-up in the computations that we report while preserving
the accuracy, based on a transformation whose mathematical foundation is
rigorous and was developed precisely for curvy complex space and,
specifically, for image processing.
According to Ying et al. (2005), a 256
3
grid takes 250 CPU seconds to
curvelet transform by a 3D curvelet. This does not include its thresholding
(which is insignificant), and inverse transforming its coarser version. Despite
this, direct use of 3D CT for denoising a 3D image result in better speed-up
factor than what we obtained here, even after adding the time of the image

79

processing to the cost of simulating the flow problem. We will report on this
aspect in Chapter 5.

4.7 Summary
With the rapid increase in the speed of computations and significant
advances in instrumentations it has become possible to simulate multiphase
flows in 2D and 3D images of heterogeneous porous media. As image-based
simulation of multiphase flow in heterogeneous porous media becomes more
practical, one must, however, address the issue of high computational cost
for such simulations, while a variety of methods have been proposed, such as
use of graphics processing units (GPUs), and parallel computations using
multiple processors, we described in this chapter a new methodology to
address the issue. The approach is based on curvelet transforming the image
and processing it by removing those details that do not contribute
significantly to physics of fluid flow. The computations in the less noisy and
coarser images yield results that are as accurate as those in the original ones,
but with a very large speed-up factor. This was demonstrated by simulation
of two-phase flow of CO2 and brine in the 3D images of a sandstone. All the
properties of the two-phase flow system, computed with the CT-processed
images were in excellent agreement with those computed with the original
resolved image. Thus, the methodology represents a significant step toward
achieving the ultimate goal of direct image-based simulations of two-phase
flows based on the governing equations, namely, making such simulations a
standard practice.

80

Chapter 5

Highly efficient image-based simulation of two-phase
flow in porous media with lattice-Boltzmann method
using three-dimensional curvelet transforms

5.1 Introduction
Porous media and materials are ubiquitous and present in everyday life.
They range anywhere from the better known porous media, such as soil,
sandpiles, oil, gas, and geothermal reservoirs and groundwater aquifers, to
lesser known ones, such as sponge, textile, printing paper, biological tissues,
wood, pavement, and asphalt. Practically all the natural porous media, and
a large fraction of man-made ones are heterogeneous, with the heterogeneity
manifesting itself over a range of length scales, from the pore to core and
field scale. Of particular interest is understanding of how fluid flow and
transport processes occur in the pore space of porous media and materials.
Aside from experimental studies of such phenomena that have been
undertaken for at least two centuries, characterization of porous media and
materials, and modeling of various phenomena in them have also been
studied for decades. To carry out theoretical and computational studies of
various phenomena in porous media one requires, as a first step, a
representative model of its pore space. The earliest theoretical studies of
fluid flow and transport in porous media were based on the assumption
(Bear, 1972) that they could be represented as continuua to which effective

81

flow and transport properties, such as the permeability and diffusivity were
attributed and, thus, the governing equation describing the phenomena of
interest in such continuua were solved. The effective properties were either
measured, or estimated by simple models, such as bundles of parallel
capillary tubes (Scheidegger, 1974). Since such models proved inadequate,
beginning with the work of Meyer (Meyer, 1953) and Fatt (Fatt, 1956, 1956),
more advanced models of pore space in the form of networks of
interconnected pores began to emerge. Although due to the computational
limitations, early pore-network models (PNMs) were relatively limited and
provided mostly qualitative insights into the phenomena under study, with
advances in computational power the PNMs gradually became quantitative,
and are still used in the study of many phenomena in porous media and
materials (Sahimi, 2011, Blunt, 2017, Hunt and Sahimi, 2017).
But even the most accurate PNMs are constructed based on a number of
assumptions, such as how the pore throats and pore bodies are defined, the
manner by which sizes are attributed to both, and where a pore throat ends
and a pore body begins, which are all approximations. On the other hand,
with significant advances in instrumentation, one can obtain high-resolution
two- or three-dimensional (3D) images of porous media that are important,
not only for their characterization (Lemmens et.al., 2011, Sheppard et.al.,
2003, Arns et.al., 2001, 2002, Porter et.al., 2010, Iglauer et.al., 2010, Andrew
et.al., 2015, Berg et.al., 2013) given the detailed information that they
provide for their morphology, but also for modeling of fluid flow and transport
through them. The only constraint for using the images in the computations
is that their physical size must be larger than the representative elementary
volume (REV), the minimum size of the porous medium for its properties to
be independent of its dimensions. Thus, use of high-resolution images in the
numerical simulations is becoming more common. They have been used for

82

computing the elastic properties of porous media (Arns et.al., 2002), the
permeability (Arns et.al., 2001; Piller et.al., 2009), simulation of two-phase
flow (Tahmasebi et.al., 2017, Kohanpur et.al., 2020), and adsorption and
deformation (Bakhshian et.al., 2018) in porous materials.
If an image is to be used in the numerical simulation of fluid flow (and
other phenomena for that matter), one must select a suitable numerical
method for solving the governing equations. One method for doing so is based
on solving the governing equations - the mass conservation and Stokes'
equations - in the image by a finite element or finite-volume method. A
second popular approach is based on the lattice-Boltzmann (LB) method
(Chen and Doolen, 1998, Succi, 2001, Mohamad, 2011, Huang et.al., 2015)
that solves the discrete Boltzmann equations. The method is based on
streaming, collision and relaxation of a set of fluid particle distribution
functions (PDFs) on a lattice. Provided that the lattice has enough symmetry,
the discrete equations reduce in the continuum limit to the Stokes' or the
Navier-Stokes equations.
Whether one solves the governing equations directly, or utilized the LB
method, the computations are intensive, particularly when the image is 3D
and has high resolution. In this chapter we propose a method for image-based
numerical simulation of two-phase flow in porous media using the LB
method, whereby the image is first processed with 3D curvelet transforms
(CTs) in order to develop a coarse, yet accurate image of the pore space. We
show that the results are as accurate as those obtained by the LB method
and the original high-resolution image, while the computations' time is
reduced by more than one order of magnitude or better.
The rest of this chapter is organized as follows. In Sec. 2 we briefly
describe the 3D porous medium that we use in the simulations. Section 3
explains the LB method used for simulating the two-phase flow problem. In

83

Sec. 4 processing of the 3D image of the porous medium by CTs is described.
The results are presented and discussed in in Sec. 5 The last section presents
a summery of the chapter.

5.2 The Porous Medium
We simulated a two-phase flow in the 3D image of two porous media. One
was Berea sandstone, a benchmark for testing accuracy and efficiency of
various computational approaches to numerical solution of multiphase flows
in porous media (Vavantne and Blunt, 2004; Ramstad et.al., 2012; Raeini
et.al., 2014). The resolution of the scan was 10.7 μm, while the physical size
of the sample was 2.14
3
mm
3
, with the size of its digitized image being 200
3

voxels. Since the porosity of the sample is about 0.2 (see Table 5.1 below),
then, roughly speaking, about 1.6 million voxels present the pore space. The
basic physical properties of the sandstone are listed in Table 5.1, while its
3D image is shown in Fig. 5.1.
The second porous medium in which we simulated two-phase flow was also
a sandstone, originally developed by An et al. (An et.al., 2017, 2020). They
drilled a cylindrical core plug from a sample whose length and diameter
were, respectively, 4 cm and 10 mm. Its image, obtained by x-ray, was then
digitized, and a section of it with 256
3
voxels was used in their simulations.
An et al., 2017, showed that the size of the sample was larger than the
representative elementary volume. Figure 5.2 presents its 3D image, while
its basic properties are also listed in Table 5.1. Following (An et al.
2017,2020), we refer to this porous medium as the Regenerated sandstone.

84

Figure 5.1: The image of the Berea sandstone used in the simulations.

Figure 5.2: The image of the Regenerated sandstone utilized in the
simulations.

The third image with a size of 512
2
pixels was from a fossiliferous outcrop
carbonate, Mt. Gambier limestone in South Australia with porosity of 0.43.
The physical size of the image was 2.76 mm
2
. Figure 5.3 shows the image.

85

Its physical properties are listed in Table 5.1. Since we previously carried out
simulation of two-phase flow in this image (Kohanpur et. al., 2020), we only
analyzed its main properties under the CT and coarse-graining, which will
be described shortly.

Figure 5.3: The image of the carbonate porous medium analyzed by curvelet
transformation.

5.3 Lattice-Boltzmann Method
To simulate the two-phase flow problem, we used the color-fluid lattice
Boltzmann (CFLB) multiphase simulator developed by Chen et al. (2018),
which is a variant of the original work of Gunstensen et al. (1991) and T𝑜 ̈lke
et al. (2006). The method produces a relatively sharp interface between
completely immiscible fluids, as well being capable of simulating two-phase
flows with high ratios of the two fluids’ viscosities that give rise to unstable
flows, due to its independent control of the surface tension and viscosity,
which is why it has been widely adopted; see, for example, Ref. (Ahrenholz,
2008).

86

In the CFLB model one labels the various fluids with color. Thus, let red
(r) and blue (b) represent the two fluids in a two-phase system, which are
represented by their own particle distribution function (PDF) 𝑓 𝑖 𝑟 (x,t) and
𝑓 𝑖 𝑏 (x,t) in the ith direction of the lattice. We used the 3D LB model with 19
velocities, i.e., the D3Q19 model. The overall PDF for the fluid system at
position x at time t is given by, 𝑓 𝑖 (𝑥 ,𝑡 ) = 𝑓 𝑖 𝑟 (x,t) + 𝑓 𝑖 𝑏 (x,t) where

𝑓 𝑖 𝑐 (𝒙 + v
𝑖 𝑐 ∆𝑡 ,𝑡 + ∆𝑡 ) = 𝑓 𝑖 𝑐 (𝒙 ,𝑡 ) +Ω
𝑖 𝑐 (𝒙 ,𝑡 ), i = 0,1,···,18 (5.1)

Here, superscript c denoting the colors (c = r, b), v
𝑖 𝑐 is the discrete velocities,
and Ω
𝑖 𝑐 (𝑥 ,𝑡 ) is the usual collision operator, given by (T𝑜 ̈lke, et. al., 2002)

Ω
𝑖 𝑐 =Ω
𝑖 𝑐 (3)
[Ω
𝑖 𝑐 (1)
+Ω
𝑖 𝑐 (2)
] . (5.2)

Where, Ω
𝑖 𝑐 (1)
is the standard LB single relaxation time collision operator
(Chen et.al., 1992; Qian et.al., 1992),

Ω
𝑖 𝑐 (1)
= −(
1
𝑇 𝑐 )[ 𝑓 𝑖 𝑐 − 𝑓 𝑖 𝑐 ,𝑒𝑞
], (5.3)

With 𝑇 𝑐 being the relaxation parameter. Ω
𝑖 𝑐 (2)
is the perturbation operator
that generates surface tension, and Ω
𝑖 𝑐 (3)
)
represents the “recoloring” that
mimics the separation of the two fluids. The macroscopic variables are given
by the standard expressions, namely, 𝜌 𝑟 = ∑ 𝑓 𝑖 𝑟 𝑖 , 𝜌 𝑏 = ∑ 𝑓 𝑖 𝑏 𝑖 , 𝜌 = 𝜌 𝑟 +
𝜌 𝑏 , and 𝜌 𝒖 ∑ 𝑓 𝑖 𝑒 𝑖 𝑖 , for, respectively, the densities of red and blue fluids, the
overall density, and the momentum, with 𝒗 = (𝑣 𝑥 ,𝑣 𝑦 ,𝑣 𝑧 ) being the fluid
velocity vector. The pressure 𝑝 is then given by, 𝑝 = |𝒗 |
2
𝜌 /3, where 𝒗 =
∆𝒙 /∆𝑡 is the speed.

87

The operator Ω
𝑖 𝑐 (2)
is given by (T𝑜 ̈lke, et. al., 2006)

Ω
𝑖 𝑐 (2)
= 𝐴 𝜎 |𝐶 |[
(
v
𝑖 𝑐 ·𝐶 )
2
𝐶 2
−
5
9
] , (5.4)

With the free parameter 𝐴 𝜎 being proportional to the interfacial tension σ,
and the same collision operator used for both fluids. The color gradient C is
given in the Appendix. The recoloring collision step redistributes 𝑓 𝑖 𝑐 ,(2)
to
achieve the separation of the two fluids (T𝑜 ̈lke, et. al., 2002), where 𝑓 𝑖 𝑐 ,(2)
is
the distribution function after applying the second operator. The step is
represented by the following maximizing problem ( Gunstensen, et. al., 1991)

𝑚𝑎𝑥 𝑓 𝑖 𝑟 ,(3)𝐶 ·∑ v
𝑖 𝑟 𝑖 𝑓 𝑖 𝑟 ,(3)
, (5.5)

In practice this is done numerically (T𝑜 ̈lke, et. al., 2002).
For two-phase fluid systems, particularly in porous media, it is more
efficient and accurate to use a multirelaxation time (MRT) collision operator,
proposed by d’Humieres et al. (2002) and Lallemand and Luo (2000), since it
also offers more stability to the computational scheme. Thus, we utilized the
MRT collision operator, in which the post-collision bulk PDF is given by
(d’Humieres et al., 2002)

𝑓 𝑖 𝑝 (𝒙 ,𝑡 ) =𝑓 𝑖 (𝒙 ,𝑡 ) − 𝐌 −1
𝐒 (𝐌𝐟 − 𝐦 eq
) , i = 0,1,···,18 , (5.6)

where M is the matrix that transforms 𝐟 to the momentum space, and S is
the diagonal matrix of the relaxation rates 𝑆 𝑖𝑖
. The size of both matrices is 19
× 19. The vector 𝐦 eq

denotes the equilibrium state of the moments 𝐦 = 𝐌𝐟
of the PDF, with 𝐌 given by d’Humieres et al. (1992), and presented in the

88

Appendix for completeness. The nineteen entries of 𝐦 eq
, as well as the
diagonal matrix S, which were also given by d’Humieres et al. (2002), are
also quoted in the Appendix for completeness.
Since the two fluids that we simulate are immiscible, the effect of the
surface tension σ, as well as the contact angle θ, must be taken into account.
An order parameter 𝜙 is defined by

𝜙 =
𝜌 𝑟 − 𝜌 𝑏 𝜌 𝑟 + 𝜌 𝑏 , (5.7)

To include the effect of contact angle, value of the order parameter 𝜙 on the
solid boundaries of the porous medium is set to be (Latava-Kokko and
Rothman, 2005), ϕ = cosθ. Then, additional terms 𝐦 s

are added to the stress-
related entries of the equilibrium moments 𝐦 eq
. Finally, the effective
viscosity µe of the mixture is taken to be

1
𝜇 𝑒 =
1+𝜙 2𝜇 𝑟 +
1−𝜙 2𝜇 𝑏 , (5.8)

One of the main characteristic quantities of any two-phase flow is the
relative permeabilities of the two fluids. The relative permeability of a fluid
phase α is the ratio of the permeability of that part of the pore space occupied
by α and the single-phase permeability Ke of the entire pore space. To begin
the LB simulation and to compute the relative permeabilities, we utilized the
steady-state method (Sahimi, 2011; Blunt, 2017) in which a predefined
fractional flow of both fluid phases is injected into the pore space at constant
flow rates, while the pressure in the opposite face of the sample is constant.
Steady state is deemed to have been reached when the downstream and
upstream fractional flows are equal.

89

To specify the predefined fractional flow in the LB method, an initial
distribution of the two fluids in the pore space at a given value of the
saturation is imposed and, then, the LB simulation determines the actual
spatial distribution of the two fluids at that saturation. The input geometry
is mirrored and periodic boundary conditions along the direction of
macroscopic flow is also imposed, in order to allow both fluid phases enter
and exit the model smoothly. In addition, a body force is applied to each fluid
phase to achieve the same pressure drop and avoid capillary end effects. The
volume flow rates of both phases are monitored in order to ensure that they
converge to the true steady state, which are then used to compute the relative
permeability at the target saturation.

5.4 Processing of the Three-Dimensional Image
The morphology of natural porous media is not completely random, but
contains extended spatial correlations even at the pore and grain scale
(Knackstedt et.al., 2001). Thus, the voxels in a 3D image of a porous medium
must also be correlated with the same type of correlation function. If the
correlations are positive, which they usually are, then, the images contain
regions with similar properties. As a result, use of resolved lattice structure
in such regions is not necessary, and use of a coarser lattice would not distort
the accuracy of the simulations. As we demonstrate shortly, this is indeed
the case.
To implement the pre-processing of the image prior to the two-phase
simulation, we exploit the fact that (Starch et.al., 2002) the CCs represent a
measure of the local complexity of an image of a porous medium between the
pores, as well as between them and solid matrix through the pores’ rough

90

surface, such that the larger the CCs, the more significant are the
correlations and the local complexity of the pore space, as well as their
contribution to the image-based simulation of multiphase flow. Thus, by
preserving only the most important parts of the pore space based on the CCs,
a coarser computational grid should suffice for the simulations.
To compress the image, we use a simple approach, which is used widely in
image processing and is referred to as thresholding. It consists of, (i) curvelet-
transforming the image to compute the CCs; (ii) normalizing the CCs by their
largest value, and (iii) introducing a threshold 0 ≤ 𝜀 < 1, such that if the
CCs 𝐶 (𝑗 ,𝑙 ,𝑘 )
< 𝜀 , they are set to zero. Thus, the pre-processing generates a
sparser representation of the image in the curvelet space, because many of
the CCs of the original image are set to zero. Other methods of compression
to develop a coarser image with the CTs have also been developed (Rabbani
et.al., 2007).
Before presenting the results, we point out that there are very good open-
source toolbox for image processing with CTs (Verma, A., 2021) which can be
used in the type of problems that we discuss in the present chapter.

5.6 Results and Discussion
We simulated a two-phase flow problem in which the porous medium was
initially saturated by brine or water, and CO2 or oil was injected into the pore
space at one face with constant injection speed 𝑣 to displace the brine or
water, while the pressure was specified on the opposite face. No-flow
boundary condition was used in the remaining four faces of the 3D images.
The capillary number was 𝐶𝑎 = 10
−4
, corresponding to slow, capillary-
dominated displacement, with 𝜇 and σ being, respectively, the viscosity and

91

interfacial tension. The ratio of viscosities of the two fluids was 10,
corresponding to an unstable displacement. For the CO2-brine system the
interfacial tension between the two fluids was assumed to be 30 mN/m; see
below for the second two-phase flow system. The two sandstones in which
the two-phase flow problem was simulated are brine-wet or water-wet and,
therefore, the contact angle for CO2, and oil, was assumed to be 180
◦
,
although, as described earlier, any other contact angle can be implemented.
All the computations were carried out on a HP Spectre x360 laptop with
speed of 2 GHz and 16 GB of memory.

A. Properties of Porous Media under Curvelet Transformation
To see how the essential features of porous media are preserved after being
processed by the CTs, consider, first, the image shown in Fig. 5.4(a), which
is that of a tight porous medium dominated by fractures (Liang, Y., 2016).
The image was processed by the CT with a threshold, 𝜀 = 0.8. The resulting
image is shown in Fig. 5.4(b). All the essential features of the medium,
including the connectivity of the fractures that plays the most important role
in the flow properties of the porous medium, are preserved, even though the
threshold is very large.
Although Fig. 5.4 indicates clearly that the essential features of the image
are preserved, we also address the question of whether any important
physical property of the porous media that we process, such as the porosity,
pore surface area, and others, are also preserved under curvelet
thresholding. To address this issue, we computed the porosity, single-phase
permeability, and pore surface area of the original two sandstones, as well
as those of the carbonate sample. After thresholding their images by the
method described earlier, we recomputed the same properties. The results
are listed in Table 5.1.

92

Figure 5.4: Comparison of the original rock sample (a) with its thresholded
version (b), obtained with a threshold ε = 0.8. Observe that all the essential
features are preserved under thresholding.

Table 5.1: Comparison of the physical properties of three porous media before
and after they are processed by curvelet transformation and their
dependence on the threshold ε.

Berea
Sandstone
Threshold ε Ke (mD) Porosity (%) Pore Surface Area (mm^2)
Original image 1273 19.63 189.74
CT image 0.5 1356 20.54 196.31
CT image 0.9 1383 20.98 199.22
Regenerated
Sandstone

Original image 125 18.36 44.11
CT image 0.5 131 19.12 46.45
CT image 0.9 134 19.77 47.23
Carbonate
Original image 5.54 43.00 1.94
CT image 0.5 5.74 44.13 2.01
CT image 0.9 5.78 44.54 2.06

93

In the case of the Berea sandstone, even with a threshold as high as ε = 0.9
(recall that the maximum value of ε is 1) the effective permeability Ke differs
from that of the original porous medium by only 8 percent (see the next
section below for the details of computing Ke).
The corresponding differences for the porosity and pore surface area are,
respectively, 6.9 and 5 percent. As we describe below, such a high threshold
results in tremendous reduction in the overall computation time. The same
type of observations may be made regarding the other two porous media that
we consider. For example, with a threshold, 𝜀 = 0.9, the permeabilities of
the curvelet-processed images of the regenerated sandstone and the
carbonate sample differ from their original values by, respectively, only 6.4
and 4.3 percent.
We also computed a multiple-point connectivity function 𝑝 (ℎ;𝑚 ) that
quantifies the long-range connectivity of a porous material. 𝑝 (ℎ;𝑚 ) is the
probability of having a sequence of m connected points in the pore space of a
porous material in a specific direction h, and is defined by

𝑝 (ℎ;𝑚 )=𝑃𝑟𝑜𝑏 [𝐼 (𝑥 )=1,𝐼 (𝑥 +ℎ( =1, ···𝐼 (𝑥 +𝑚𝐻 )=1] , (5.9)

where 𝐼 (𝑥 ) is the indicator function, i.e., 𝐼 (𝑥 )=1 if x is in the pore space, and
𝐼 (𝑥 )=0 , otherwise. Note that 𝑝 (ℎ;𝑚 ) accounts for curvilinearity and
complexity of a microstructure, since it represents the probability of finding
multiple connected points in such complex media. As such, it provides a
stringent test of preserving the essential features of the image under CTs.
We computed 𝑝 (ℎ;𝑚 ) for the two sandstones. Figure 5.5 shows the results
for the Berea sandstone, and its curvelet-processed images. As can be seen
there is hardly any difference between the original image and the processed
versions. Similar results were obtained for the Regenerated sandstone. We

94

conclude, therefore, that all the essential properties of the pore space are
preserved under curvelet transformation.

Figure 5.5: The connectivity function 𝑝 (ℎ;𝑚 ) computed for the Berea
sandstone, and its two coarsened versions, with m = 100.

B. Berea Sandstone
We first determined the minimum lattice size for the CFLB simulations that
would yield an effective single-phase permeability Ke that will not change, if
the resolution of the lattice was increased and a larger number of lattice
nodes was used. To do so, we computed the dependence of Ke on the lattice
size. Figure 5.6 presents the results, indicating that the Q19 lattice with
about 1.69 × 10
6
lattice points produces effective permeability that is
independent of the size. The same lattice was then used in the simulation of
two-phase flow. To estimate the speed-up in the computations, we carried
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0 50 100 150 200 250 300
Connectivity Probability
Pixel
Original
= 0.5
= 0.9
ε
ε

95

out the same detailed simulations in the original image and in two coarser
ones obtained by thresholding the image with ε = 0.5 and 0.9.

Figure 5.6: Dependence of the permeability of the Berea sandstone on the
number of the lattice nodes N in the LB simulation.

In Fig. 5.7 we compare the spatial distribution of CO2 and brine in the
original image of the sandstone with those obtained with the CT-processed
images with the two thresholds. The brine saturation is 𝑆 𝑏 = 0.5. All the
important features of three distributions are completely similar.
To make a quantitative comparison between the results obtained with the
three images, we compare in Fig. 5.8 the saturation profile of CO2 along the
macroscopic direction 𝑥 of the flow of CO2 at its breakthrough time 𝑡 𝐵 , i.e.,
the time at which CO2–filled pores form a sample–spanning cluster across
the image. In Fig. 5.8 𝑋 = 𝑥 /𝐿 , where L is the linear size of the image, and
the saturation SCO2(X) was computed in the planes perpendicular to the
macroscopic direction at 𝑋 . The agreement between the three sets of results
is excellent.
1000
1200
1400
1600
1800
0 100 200 300 400 500 600
Absolute Permeability (mD)
N (number of lattice points) x 10
-4

96

(a) (b) (c)
Figure 5.7: CO2 invasion patterns in the Berea sandstone. (a) The original
image; (b) in the thresholded image with ε = 0.5 image, (c) in thresholded
image with ε = 0.9 with brine saturation 𝑆 𝑏 = 0.5. CO2 is shown in blue and
brine in red.

Figure 5.8: Saturation profile of CO2 along the direction x of macroscopic flow
at the breakthrough-time in the original images of Berea sandstone, as well
as two thresholded ones, where 𝑋 = 𝑥 /𝐿 is the normalized distance of the
interface from the inlet, with L being the sample's length in the x direction.
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
S
CO2
X
Original
= 0.5
= 0.9
ε
ε

97

We make another quantitative comparison between the three sets of
results in Fig. 5.9, where we show time–dependence of CO2 saturation SCO2
in the original image of the porous medium and those computed with the
thresholded images. In this figure, 𝑡 𝐷 = 𝑡 /𝑡 𝑠 , where 𝑡 is the actual time, and
𝑡 𝑠 is the time at which the flow system has reached steady state. Saturation
increases linearly with the time as more CO2 is injected into the pore space,
and after the system reaches steady state, it reaches a constant value. More
importantly, not only are the three profiles completely similar, but also the
maximum difference between the profile computed in the coarser image with
the threshold 𝜀 = 0.9 and the original image is practically negligible.

Figure 5.9: Comparison of time-dependence of saturation of CO2, where 𝑡 𝐷 =
𝑡 /𝑡 𝑠 , with 𝑡 𝑠 being the time at which the system reaches steady state, in the
original image of the Berea Sandstone and its curvelet-transformed images.

0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0 0.2 0.4 0.6 0.8 1
S
CO2
t
D
Original
= 0.5
= 0.9
ε
ε

98

Next, we compare the relative permeabilities, the most important
quantity in two-phase flow in porous media. The relative permeability Kr,α of
fluid α at saturation Sα is given by the generalized Darcy’s law:

𝑣 𝛼 =−
𝑘 𝑟 ,𝛼 (𝑆 𝛼 )𝑘 𝑒 𝜇 𝛼 𝛥𝑃
𝛼 𝛥𝑥
, (5.10)

where 𝑣 𝛼 is the Darcy velocity of fluid phase 𝛼 , ∆Pα is the pressure drop in
phase 𝛼 over a distance ∆x, and 𝜇 𝛼 is the viscosity of phase 𝛼 . Figure 5.10
compares the three sets of results for the relative permeabilities of the two
fluids. The agreement between the three sets of results is excellent.

Figure 5.10: Relative permeabilities of CO2 and brine in the Berea sandstone
as functions of the brine saturation 𝑆 𝑏 in the sandstone's original image, as
well as its curvelet-thresholded images.

0
0.2
0.4
0.6
0.8
1
0.0 0.2 0.4 0.6 0.8 1.0
Kr
Sw
Original,
Original,
= 0.5,
= 0.5,
= 0.9,
= 0.9,
ε
ε
ε
ε
Krb
Krb
Krb
KrCO2
KrCO2
KrCO2

99

The high accuracy of the computed relative permeabilities for the
thresholded results provided the motivation for making another quantitative
comparison between the three sets of results. We show in Fig. 5.11 the
distribution of local CO2 velocities in the original image and compare it with
those computed with the thresholded images. All the velocities were
normalized by their maximum value. For flow velocities up to the maximum
of the distributions, the agreement between the three sets of results is
perfect. For flow velocities larger than the maximum, the agreement is good,
with all the distributions having the same shape, and eventually the same
tail.

Figure 5.11: The probability density function (PDF) of the local velocity of
CO2, in the Berea sandstone normalized by its maximum, in the curvelet-
transformed images of the sandstone with two thresholds ε, and its
comparison with that of the original image.

0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
0 0.2 0.4 0.6 0.8 1
PDF
Normalized Velocity
Original
= 0.5
= 0.9
ε
ε

100

From a practical view point, a most important quantity in two-phase flow
in porous media, in addition to the relative permeability, is the residual
saturation 𝑆 𝑟 of the displaced fluid (here, brine) at CO2 breakthrough point.
We compare in Table 5.2 the three computed residual saturations. The
maximum difference, computed with 𝜀 = 0.9, is about 7 percent.

C. The Regenerated Sandstone

We simulated two-phase flow of water and oil in the Regenerated sandstone.
Once again, water was the wetting fluid. The interfacial tension between the
two fluids was 20 mN/m. Oil was injected into the pore space from one face
to displace the water saturating the medium, and the pressure at the
opposite face was held constant. Periodic boundary conditions were used in
the other four faces.
In terms of their accuracy, all the results for the Regenerated sandstone
are completely similar to those for the Berea sandstone. Therefore, we
present only two main results. Figure 5.12 presents the saturation of oil
versus the dimensionless time 𝑡 𝐷 , defined earlier. The agreement between
the three profiles is excellent.
Figure 5.13 presents the calculated relative permeabilities to oil and water.
Once again, similar to the Berea sandstone, the agreement between the
relative permeabilities, computed for the original image, and those for the
two coarser images, is excellent.

101

Figure 5.12: Comparison of time-dependence of saturation of oil, where 𝑡 𝐷 =
𝑡 /𝑡 𝑠 , with 𝑡 𝑠 being the time at which the system reaches steady state, in the
original image of the Regenerated sandstone and its curvelet-transformed
images

Figure 5.13: Relative permeabilities of oil and water in the Regenerated
sandstone as functions of the brine saturation 𝑆 𝑏 in the sandstone's original
image, as well as its curvelet-thresholded images.
0
0.2
0.4
0.6
0.8
1
0 0.2 0.4 0.6 0.8 1
S
o
t
D
Original
= 0.5
= 0.9
ε
ε
0
0.2
0.4
0.6
0.8
1
0.0 0.2 0.4 0.6 0.8 1.0
Kr
Sw
Original,
Original,
= 0.5,
= 0.5,
= 0.9,
= 0.9,
ε
ε
ε
ε
Krw
Krw
Krw
Kro
Kro
Kro

102

Table 5.2: Comparison of the number of the lattice nodes M, the
computation times (in CPU seconds), the speed-up factor S of the
computations, and the brine residual saturation Sr, computed
with the original and thresholded images.

D. Efficiency of the Computation
It took 196 CPU seconds to compute the CT of the 3D image of the Berea
sandstone, while inverse transforming of the coarse images took 203 CPU
seconds. The computation time for thresholding the image is insignificant.
The corresponding times for the Regenerated sandstone, the carbonate
medium, and the fractured sample were, respectively, 641, 23, and 8 CPU
second. In Table 5.2 we compare the total computation times, including the
processing of the images, for all the cases, as well as the speed-up factor of
the computations.
Table 5.2 indicates that computations with the coarse image of the Berea
sandstone that was obtained with the very high threshold of 𝜀 = 0.9 still
produce very accurate results, with a speed-up factors of 35, while those with
the thresholds ε = 0.5 are accelerated by a factor of 25.
On the other hand, the speed-up factor for the Regenerated sandstone is
between 8 and 11, still very significant. There are two reasons for a smaller
Berea Sandstone
M Threshold ε Sr Ke Time S
Original image
1691443
0.410 1,273 140,422 -
Curvelet-transformed image 72978 0.5 0.433 1,356 5,546 25
Curvelet-transformed image 58421 0.9 0.441 1,383 3,953 35
Regenerated Sandstone

Original image 4140157 0.1487 125 197969 -
Curvelet-transformed image 609748 0.5 0.1517 131 25183 8
Curvelet-transformed image 396901 0.9 0.1534 134 18114 11

103

speed-up factor for this porous medium. One is that even though the
porosities of the two porous media are virtually the same, the Regenerated
sandstone is more tortuous, hence necessitating larger number of lattice cells
to accurately represent the pore space. The second reason is the size of the
lattice in the original image of the sandstone, which is roughly 2.5 times
larger than that of the Berea sample. Indeed, the speed-up factor for the
Regenerated sandstone is also reduced by a factor of roughly 3, close to ratio
of the sizes of the two lattices.

5.7 Summary
Lattice-Boltzmann simulation of multiphase flow in heterogeneous porous
media, represented by their high-resolution images, has become practical,
but the computations are still intensive and take a long time. Thus, one must
address the issue of high computational cost of such simulations. We
proposed a method by which one preprocesses the 3D image by taking it to
the curvelet transform space and eliminating those details that do not
contribute significant to multiphase fluid flow. After inverting the coarser
image, it is utilized in the simulation of multiphase flow by the CFLB-MRT
approach. We demonstrated that the method accelerates the computations
very significantly by a speed-up factor of up to 35.
In a recent paper (Aljasmi and Sahimi, 2020) in which we solved the
Stokes' and the diffusion equations in the images of porous media in order to
compute their effective permeability and diffusivity, as well as the flow and
concentration fields, we showed that if the two equations are curvelet
transformed and solved in the curvelet-transformed image of the porous
medium, one obtains a speed-up factor of four or better, without

104

preprocessing and thresholding the image. The reason is that a coarser
computational grid suffices in the curvelet space because many of the CCs of
the image are either very small and, thus, do not contribute significantly to
fluid flow, and due to the aforementioned correlations between the voxels the
CCs are close to those of their neighbors, hence giving rise to a relatively
smooth local environment in the curvelet space that makes it possible to use
larger lattice blocks.
The same approach may be taken for the LB simulation of fluid flow and
transport in the images of porous media. One first determines the CT of the
Boltzmann equation that, for example, for a single component in the
Bhatnagar, Gross and Krook approximation (Bhatnagar, et. al., 1954) is
given by (He and Luo, 1997),

𝜕𝑓
𝜕𝑡
+𝐮 ·𝛻𝑓 +
1
τ
𝑓 =
1
τ
𝑔 , (5.11)

and similarly for the multicomponent version of the same equation. Here, τ
is the relaxation time due to particle collision, and 𝑔 is the Maxwell-
Boltzmann distribution function that in the 3D space is given by,

𝑔 =
𝜌 (2𝜋𝑅𝑇 )
3/2
exp[−
(𝐮 −𝐯 )
𝟐 2𝑅𝑇
] , (5.12)

with R and T being, respectively, the gas constant and the temperature. The
discretized Boltzmann equations used in the simulations can also be curvelet
transformed. The curvelet-transformed equation(s) are then solved
numerically in the curvelet space, and the numerical results are then
curvelet-transformed back to real space.

105

Bibliography
Ahrenholz, B. , Tölke, J. , Lehmann, P. A. Peters, Kaestner, A., Krafczyk, M.
and Durner, W. (2008) Prediction of capillary hysteresis in a porous
material using lattice-Boltzmann methods and comparison to
experimental data and a morphological pore network model, Adv. Water
Resour. 31, 1151.

Aljasmi, A., Sahimi, M. (2020). Efficient image-based simulation of flow and
transport in heterogeneous porous Media: Application of curvelet
transforms. Geophys. Res. Lett. 47, e2019GL085671.

Al-Zubi, S., Islam, N., Abbod, M. (2011): Multiresolution analysis using
wavelet, ridgelet, and curvelet transforms for medical image
segmentation. Int. J. Biomed. Imaging 2011, 136034.

An, S., Yu, H. and Yao, J., GPU-accelerated volumetric lattice Boltzmann
method for porous media flow. J. Pet. Sci. Eng. 156, 546 (2017).

An, S., Zhan, Y., Yao, J., Yu, H.W., Niasar, V. (2020). A greyscale volumetric
lattice Boltzmann method for upscaling pore-scale two-phase flow. Adv.
Water Resour. 144, 10711.

Andrew, M., Menke, H., Blunt, M.J., and Bijeljic, B. (2015). The imaging of
dynamic multiphase fluid flow using synchrotron-based X-ray
microtomography at reservoir conditions, Transp. Porous Media 110,
124.

Arns, C.H., Knackstedt, M.A., Pinczewski, W.V., and Lindquist,W.B. (2001).
Accurate computation of transport properties from microtomographic
images, Geophys. Res. Lett. 28, 33613364.

106

Arns, C.H., Knackstedt, M.A., Pinczewski, W.V., and Garboczi, E. (2002).
Computation of linear elastic properties from microtomographic images:
Methodology and agreement between theory and experiment,
Geophysics 67, 13481672.

Aslannejad, H., Hassanizadeh, S.M. (2017). Study of hydraulic properties of
uncoated paper: Image analysis and pore-scale modeling. Transp. Porous
Media 120, 67.

Aslannejad, H., Hassanizadeh, S.M., Celia, M.A. (2019). Characterization of
the interface between coating and fibrous layers of paper. Transp. Porous
Media 127, 143.

Babaei, M., King, P.R. (2011). A comparison between wavelet and
renormalization upscaling methods and iterative upscaling-downscaling
scheme. SPE Reservoir Simulation Symposium 1, 469.

Bakhshian, S., Shi, Z., Sahimi, M., Tsotsis, T.T., and Jessen, K. (2018).
Image-based modeling of gas adsorption and swelling in high-pressure
porous formations, Sci. Rep. 8, 8249.

Bear, J. (1972). Dynamics of Fluids in Porous Media. Dover, Mineola.

Berg, S., Ott, H., Klapp, S.A., Schwing, A., Neiteler, R., Brussee, Makurat,
A., Leu, L., Enzmann, F., Schwarz, J.-O., Kersten, M., Irvine,
Stampanoni, M. (2013). Real-time 3D imaging of Haines jumps in porous
media flow. Proc. Natl. Acad. Sci. U.S.A. 110, 3755.

Bermejo-Moreno, I., Pullin, D. (2008). On the non-local geometry of
turbulence, J. Fluid Mech., 603, 101-135.

Beylkin, G., Coifman, R., and Rokhlin, V. (1991). Fast wavelet transforms
and numerical algorithms. Comm. on Pure and Appl. Math. 44, 141–183.

Bhatnagar, P.L. Gross, E.P. and Krook, M. (1954). A model for collision
processes in gases. I. Small amplitude processes in charged and neutral
one-component systems, Phys. Rev. 94, 511.

107

Blunt, M.J., King, P.R (1991). Relative permeabilities from two- and three-
dimensional pore-scale network modelling. Transp. Porous Media 6, 407.

Blunt, M.J., King, M.J., Scher, H. (1992). Simulation and theory of two-phase
flow in porous media. Phys. Rev. A 46, 7680.

Blunt, M.J., Scher, H (1995). Pore-level modeling of wetting. Phys. Rev. E
52, 6387.

Blunt, M.J. (1997). Effects of heterogeneity and wetting on relative
permeability using pore level modeling. SPE J. 2, 70.

Blunt, M.J., Bijeljic, B., Dong, H., Gharbi, O., Iglauer, S., Mostaghimi, P.,
Paluszny, A., and Pentland, C. (2013). Pore-scale imaging and modelling.
Advances in Water Resources, 51, 197-216.

Blunt, M.J. (2017). Multiphase Flow in Permeable Media. Cambridge
University Press, Cambridge.

Boubchir, L., and Fadili. J. (2005). “Multivariate statistical modelling of
images with the curvelet transform,” in Proceedings of the 8th
International Conference on Signal Processing, Pattern Recognition, and
Applications, pp. 747–750.

Cand´es, E., Demanent, L., Donoho, D.L., and Ying, L. (2005). Fast discrete
curvelet transforms, Multiscale Model. Simul. 5, 861.

Chandler, R., Koplik, J., Lerman, K., Willemsen, J. (1982). Capillary
displacement and percolation in porous media. J. Fluid Mech. 119, 249.

Chen, S. and Doolen, G.D. (1998). Lattice Boltzmann method for fluid flows,
Annu. Rev. Fluid. Mech. 30, 329.

Chen, H. Chen, S. and Matthaeus, W.H. (1992). Recovery of the Navier-
Stokes equations using a lattice-gas Boltzmann method, Phys. Rev. A 45,
R5339.

108

Chen, Y. Li, Y. Valocchi, A.J. and Christensen, K.T. (2018) Lattice
Boltzmann simulations of liquid CO2 displacing water in a 2D
heterogeneous micromodel at reservoir pressure conditions, J.
Contaminant Hydro. 212, 14.

Cooper, L.J., Daly, K.R., Hallett, P.D., Naveed, M., Koebernick, N.,
Bengough, A.G., George, T.S., and Roose, T. (2017). Fluid flow in porous
media using image-based modelling to parametrize Richards’ equation,
Proc. Roy. Soc. London A 473, 20170178.

Dashtian, H., and Sahimi, M. (2014). Coherence index and curvelet
transformation for denoising geophysical data, Phys. Rev. E 90, 042810.

Daubechies, I. (1988). Orthonormal basis of compactly supported wavelets.
Commun. Pure Appl. Math. 41, 901.

Daubechies, I. (1992). Ten Lecture on Wavelets. SIAM, Philadelphia.

Davis, M. (2011). Finite Difference Methods, Department of Mathematics,
Imperial College London.

d'Humieres, D. Ginzburg, I., Krafczyk, M., Lallemand, P., and Luo, L.-S.
(2002). Multiplerelaxation-time lattice Boltzmann models in three
dimensions, Phil. Trans. R. Soc. Lond. A 360, 437.

Donoho, D.L. (1999). Wedgelets: nearly minimax estimation of edges, Ann.
Statist. 27, 859.

Ebrahimi, F. (2010). Invasion percolation: A computational algorithm for
complex phenomena. Comput. Sci. Eng. 12(2), 84.

Ebrahimi, F., and Sahimi, M. (2002). Multiresolution wavelet coarsening
and analysis of transport in heterogeneous media, Physica A 316, 160.

Ebrahimi, F., Sahimi, M. (2004). Multiresolution wavelet scale up of
unstable miscible displacements in flow through porous media. Transp.
Porous Media 57, 75.

109

Esmaeili, M. Dehnavi, A.M. and Rabbani, H. (2017). 3D curvelet-based
segmentation and quantification of drusen in optical coherence
tomography images, Hindawi J. Elect. Comput. Eng. 2017, Article
4362603.

Etemad, S., Behrang, A., Mohammadmoradi, P., Hejazi, S.H., Kantzas, A.
(2017). Effects of surface roughness and mineral heterogeneity on pore-
scale steam condensation. J. Pet. Sci. Eng. 159, 624.

Fatt, I. (1956) The network model of porous media. I. Capillary pressure
Characteristics, Trans. AIME 207, 155.

Fatt, I. (1956) The network model of porous media. III. Dynamical properties
of networks with tube radius distributions, Trans. AIME 207, 164.

Francois, M.M., Cummins, S.J., Dendy, E.D., Kothe, D.B., Sicilian, J.M.,
Williams, M.M. (2006). A balanced-force algorithm for continuous and
sharp interfacial surface tension models within a volume tracking
framework. J. Comput. Phys. 213, 141.

Friedlingstein, P., Solomon, S. (2005). Contributions of past and present
human generations to committed warming caused by carbon dioxide.
Proc. Natl. Acad. Sci. U.S.A. 102, 10832.

Ghanbarian, B., Sahimi, M., Daigle, H. (2016). Modeling relative
permeability of water in soil: Application of effective-medium
approximation and percolation theory. Water Resour. Res. 52, 5025.

Ghassemzadeh, J., Hashemi, M., Sartor, L., Sahimi, M. (2001). Pore network
simulation of fluid imbibition into paper during coating processes: I.
Model development. AIChE J. 47, 519

Ghassemzadeh, J., Sahimi, M. (2004). Pore network simulation of fluid
imbibition into paper during coating III: Modeling of the two-phase flow.
Chem. Eng. Sci. 59, 2281.

110

Ginzburg, I. (2005) Equilibrium-type and link-type lattice Boltzmann models
for generic advection and anisotropic-dispersion equation, Adv. Water
Resour. 28, 1171.

Ginzburg, I.(2006). Variably saturated flow described with the anisotropic
lattice Boltzmann methods, Comput. Fluids 35, 831.

Glimm, J., Grove, J.W., Li, X.L., and Zhao, N. (1999).Simple front tracking,
Contemp. Math. 238(2) 133–149.

Grinsted, A., Moore, J.C., Jevrejeva, S. (2004). Application of the cross
wavelet transform and wavelet coherence to geophysical time series.
Nonlinear Processes Geophys. 11. 561.

Gueyffier, D., Li, J., Nadim, A., Scardovelli, R., Zaleski, S. (1999). Volume-
of-fluid interface tracking with smoothed surface stress methods for
three-dimensional flows. J. Comput. Phys. 152, 423.

Gunstensen, A.K. Rothman, D.H. Zaleski, S. and Zanetti, G. (1991). Lattice
Boltzman model of immiscible fluids, Phys. Rev. A 43, 4320.

He, X. and Luo, L.-S. (1997).A priori derivation of the lattice Boltzmann
equation, Phys. Rev. E 55, R6333.

Heiba, A.A., Sahimi, M., Scriven, L.E., Davis, H.T. (1992). Percolation theory
of two-phase relative permeability. SPE Reservoir Eng. 7, 123.

Helmig, R., Schulz, P. (1997). Multiphase Flow and Transport Processes in
the Subsurface (Springer, Berlin).

Herrmann, F.J., Wang, D., Hennenfent, G., and Moghaddam, P.P. (2007).
Curvelet-based seismic data processing: A multiscale and nonlinear
approach, Geophysics 73, A1.

Herman, G.T. (2009). Fundamentals of Computerized Tomography: Image
Reconstruction from Projections, 2nd ed., Springer, New York.

111

Hirt, C.W., Nichols, B.D. (1981). Volume of fluid (VOF) method for the
dynamics of free boundaries, J. Comput. Phys. 39(1) 201–225.

Huang, H., Sukop, M., and Lu, X. (2015). Multiphase Lattice Boltzmann
Methods: Theory and Application. Wiley, New York.

Hunt, A.G., and Sahimi, M. (2017). Flow, transport, and reaction in porous
media: percolation scaling, critical-path analysis, and effective-medium
approximation, Rev. Geophys. 55, 993.

Iglauer, S., Favretto, S., Spinelli, G., Schena, G., Blunt, M.J. (2010). X-ray
tomography measurements of power-law cluster size distributions for the
nonwetting phase in sandstones. Phys. Rev. E. 82, 056315.

Kak, A.C., and Slaney, M. (2001). Principles of Computerized Tomographic
Imaging (SIAM, Philadelphia.

Kantzas, A., Chatzis, I. (1988). Network simulation of relative permeability
curves using a bond correlated-site percolation model of pore structure.
Chem. Eng. Commun. 69, 191.

Kantzas, A., Marentette, D.F., Jha, K.N.N. (1992). Computer-assisted
tomography: from qualitative visualization to quantitative core analysis.
J. Canadian Pet. Technol. 31, PETSOC-92-09-06.

Kelly, K., Ward, R. and Treitel, S. (1976). Synthetic seismograms: a finite
difference approach, Geophysics, 1, 2-27.

Knackstedt, M.A., Sheppard, A.P., Sahimi, M. (2001). Pore network
modeling of two-phase flow in porous rock: The effect of correlated
heterogeneity. Adv. Water Resour. 24, 257.

Kohanpur, A.H., Rahromostaqim, M., Valocchi, A.J., Sahimi, M. (2020). Two-
phase flow of CO2-brine in a heterogeneous sandstone: characterization
of the rock and comparison of the lattice-Boltzmann, pore-network, and
direct numerical simulation methods. Adv. Water Resour. 135, 103439.

112

Krämer, R., and Ziegler, C. (2013). Membrane Transport Mechanism: 3D
Structure and Beyond.

Kr𝑢 ̈ ger, T., Kusumaatmaja, H., Kuzmin, A., Shardt, O., Silva, G., and
Viggen, E.M. (2017). The Lattice Boltzmann Method. Springer, New
York.

Kumar, P., Foufoula-Georgiou, E., (1997). Wavelet analysis for geophysical
applications. Rev. Geophys. 35, 385-412.

Lallemand, P. and Luo, L.-S. (2000). Theory of the lattice Boltzmann
method: dispersion, dissipation, isotropy, Galilean invariance and
stability, Phys. Rev. E 61, 6546
.
Larson, R.G., Scriven, L.E., Davis, H.T. (1977). Percolation theory of residual
phases in porous media. Nature 268, 409.

Latva-Kokko, M. and Rothman, D.H. (2005). Static contact anglein lattice
Boltzmann models of immiscible fluids, Phys. Rev. E 72, 046701.

Lemmens, H.J., Butcher, R., Botha, P.W.S.K. (2011). FIB/SEM and
SEM/EDX: a new dawn for the SEM in the core lab? Petrophysics 52,
452.

Liang, Y. (2016) Rock fracture skeleton tracing by image processing and
quantitative analysis by geometry features, J. Geophys. Eng. 13, 273.

Lysmer, J. and Draker. L. (1972). A finite element method for seismology in
methods in computational physics, Vol. 11. Academic Press, New York.

Ma, J., and Hussaini, M. (2007). Three-dimensional curvelets for coherent
vortex analysis of turbulence, Appl. Phys. Lett., 91, 184101.

Ma, J., Hussaini, M., Vasilyev, O., and Le Dimet, F.-X. (2009). Multiscale
geometric analysis of turbulence by curvelets, Physics of Fluids, 21,
075104.

113

Ma, J., and Plonka, G. (2009). Computing with curvelets: From image
processing to turbulent flows, Comput. Sci. Eng. 11 (2), 72.

Mallat, S.G. (1989a). A theory for multiresolution signal decomposition: the
wavelet representation. IEEE Trans. Pattern Recog. Machine Intell. 11,
674.

Mallat, S.G. (1989b). Multiresolution approximations and wavelet
orthonormal bases of L2(R). Trans.Am. Math. Soc. 315, 69.

Mehrabi, A.R., and Sahimi, M. (1997). Coarsening of heterogeneous media:
Application of wavelets, Phys. Rev. Lett. 79, 4385.

Metz, B., Davidson, O., De Coninck, H., (2005). Carbon Dioxide Capture and
Storage: Special Report of the Intergovernmental Panel on Climate
Change. Cambridge University Press, London.

Meyer, H.I. (1953). Pore distribution in porous media, J. Appl. Phys. 24, 510

Meyer, Y. (1992). Wavelet and Operators, Cambridge University Press,
Cambridge.

Mohamad, A.A. (2011). Lattice Boltzmann Method (Springer, New York,).

Neelamani, R., Baumstein, A.I., Gillard, D.J., Hadidi, M.T., and Soroka,
W.L. (2008). Coherent and random noise attenuation using the curvelet
transform, The Leading Edge 27, 129.

Nordbotten, J.M., Celia, M.A. (2011). Geological Storage of CO2: Modeling
Approaches for Large-Scale Simulation (Wiley, New York).

Oak, M., Baker, L., Thomas, D. (1990). Three-phase relative permeability of
Berea sandstone. J. Pet. Technol. 42, 1054.

Olhede, S., and Walden, A.T. (2004). The Hilbert spectrum via wavelet
projections, Proc. R. Soc. Lond. A 460, 955.

114

Pancaldi, V., Christensen, K., King, P.R. (2007). Permeability up-scaling
using Haar wavelets. Transp. Porous Media 67, 395.

Piri, M., Blunt, M.J. (2005a). Three-dimensional mixed-wet random pore-
scale network modeling of two- and three-phase flow in porous media. I.
Model description. Phys. Rev. E 71, 026301.

Piri, M., Blunt, M.J. (2005b). Three-dimensional mixed-wet random pore-
scale network modeling of two- and three-phase flow in porous media. II.
Results. Phys. Rev. E 71, 026302.

Piller, M., Schena, G., Nolich, M., Favretto, S., Radaelli,F., and Rossi, E.
(2009). Analysis of hydraulic permeability in porous media: From high
resolution X-ray tomography to direct numerical simulation, Transp.
Porous Media 80, 57.

Popinet, S., Zaleski, S. (1999) A front-tracking algorithm for accurate
representation of surface tension, Int. J. Numer. Methods Fluids 30(6)
775–793.

Porter, M.L., Wildenschild, D., Grant, G., Gerhard, J.I. (2010). Measurement
and prediction of the relationship between capillary pressure, saturation,
and interfacial area in a NAPL-water-glass bead system. Water Resour.
Res. 46, W08512.

Qian, Y., d'Humieres, and Lallemand, P. (1992). Lattice BGK models for
Navier-Stokes equations, Europhys. Lett. 17, 479.

Rabbani, H. Vafadust, M. and Gazor, S. (2007). Image denoising in curvelet
transform domain using gaussian mixture model with local parameters
for distribution of noise-free coefficients, in 2007 4th IEEE/EMBS
International Summer School and Symposium on Medical Devices and
Biosensors; https://ieeexplore.ieee.org/document/4338291

Raeini, A.Q., Blunt, M.J., and Bijeljic, B. (2012). Modelling two-phase flow
in porous media at the pore scale using the volume-of-fluid method, J.
Comput. Phys. 231(17):5653–5668.

115

Raeini, A.Q. (2013). Modelling Multiphase Flow Through Micro-CT Images
of the Pore Space, Ph.D. Thesis, Imperial College of London.

Raeini, A.Q., Blunt, M.J., and Bijeljic, B. (2014). Direct simulations of two-
phase flow on micro-CT images of porous media and upscaling of pore-
scale forces, Adv. Water Resource. 74, 116.

Raeini, A.Q., Bijeljic, B., and Blunt, M.J. (2015). Modelling capillary trapping
using finite volume simulation of two-phase flow directly on micro-CT
images, Adv. Water Resource. 83, 102.

Rahaman, M., Sikdar, M. H., Hossain, M. B. (2015). Numerical Solution of
Diffusion Equation by Finite Difference Method, Journal of
Mathematics, p-ISSN: 2319-765X. 11.

Ramstad, T., Idowu, N., Nardi, C., Øren, P.-E. (2012). Relative permeability
calculations from two-phase flow simulations directly on digital images
of porous rocks. Transp. Porous Media 94, 487-504.

Rasaei, M.R., and Sahimi, M. (2008). Efficient simulation of water flooding
in three-dimensional heterogeneous reservoirs using wavelet
transformations: application to the SPE-10 model, Transp. Porous Media
72, 311.

Rasaei, M.R., Sahimi, M. (2008). Upscaling and simulation of waterflooding
in heterogeneous reservoirs using wavelet transformations: Application
to the SPE-10 model. Transp. Porous Media 72, 311.

Rasaei, M.R., Sahimi, M. (2009). Upscaling of the permeability by multiscale
wavelet transformations and simulation of multiphase flows in
heterogeneous porous media. Comput. Geosci. 13, 187.

Rezapour, A., Ortega, A., Sahimi, M. (2019). Upscaling of geological models
of oil reservoirs with unstructured grids using lifting-based graph
wavelet transforms. Transp. Porous Media 127, 661.

Sahimi, M., Hashemi, M., Ghassemzadeh, J. (1998). Site-bond invasion
percolation with fluid trapping. Physica A 260, 231-243.

116

Sahimi, M. (2011). Flow and Transport in Porous Media and Fractured Rock,
2nd ed., Wiley-VCH, Weinheim.

Sahimi, M., Heiba, A.A., Davis, H.T., Scriven, L.E. (1986). Dispersion in flow
through porous media: II. Two-phase flow. Chem. Eng. Sci. 41, 2123.

Sankey, M.H., Holland, D.J., Sederman, A.J., Gladden, L.F. (2009). Magnetic
resonance velocity imaging of liquid and gas two-phase flow in packed
beds. J. Magn. Reson. 196, 142.

Scheidegger, A.E. (1974). The Physics of Flow Through Porous Media, 3rd
ed. (University of Toronto Press, Toronto,).

Schindelin, J., Arganda-Carreras, I., Frise, E., Kaynig, V., Longair, M.,
Pietzsch, T., Preibisch, S., Rueden, C., Saalfeld, S., Schmid, B. (2012).
Fiji: an open-source platform for biological image analysis. Nature
Methods, 9, 676.

Sheppard, S., Mantle, M.D., Sederman, A.J., Johns, M.L., Gladden, L.F.
(2003). Magnetic resonance imaging study of complex fluid flow in porous
media: flow patterns and quantitative saturation profiling of amphiphilic
fracturing fluid displacement in sandstone cores. Magn. Reson. Imaging.
21, 365.

Shinde, A.A., Rahulkar, A.D., Patil, C.V. (2017). Fast discrete curvelet
transform-based anisotropic feature extraction for biomedical image
indexing and retrieval. Int. J. Multi-media Information Retrieval 6, 281-
288.

Shokri, N., Sahimi, M., Or, D. (2012). Morphology, propagation dynamics
and scaling characteristics of drying fronts in porous media. Geophys.
Res. Lett. 39, L09401.

Siau, J.F. (1984). Transport Processes in Wood. Springer, Berlin.

Starck, J.-L., Cand´es, E.J. and Donoho, D.L. (2002). The curvelet transform
for image denoising, IEEE Trans. Image Process. 11 (6), 670.

117

Succi, S. (2001). The Lattice Boltzmann Equations for Fluid Dynamics and
Beyond (Oxford University Press, Oxford).

Sun, B., Ma, J., Chauris, H., and Yang, H. (2009). Solving the wave equation
in the curvelet domain: a multiscale and multidirectional approach, J.
Seismic Exploration, 18, 385-399.

Swerin, A. (2018). Dimensional scaling of aqueous ink imbibition and inkjet
printability on porous pigment coated paper { A revisit. Ind. Eng. Chem.
Res. 57, 49.

Tahmasebi, P., and Sahimi, M. (2012). Reconstruction of three-dimensional
porous media using a single thin section, Phys. Rev. E 85, 066709.

Tahmasebi, P., Javadpour, F., and Sahimi, M. (2015). Multiscale and
multiresolution modeling of shales and their flow and morphological
properties, Sci. Rep. 5, 16373.

Tahmasebi, P., Sahimi, M., Kohanpur, A.H., and Valocchi, A. (2017). Pore-
scale simulation of flow of CO2 and brine in reconstructed and actual 3D
rock cores, J. Pet. Sci. Eng. 155, 21.

Tölke, J., Krafczyk, M. Schulz, M. and Rank, E. (2002). Lattice Boltzmann
simulations of binary fluid flow through porous media, Phil. Trans. R.
Soc. Lond. A 360, 535.

Tölke, J., Freudiger, S. and Krafczyk, M. (2006). An adaptive scheme using
hierarchical grids for lattice Boltzmann multi-phase flow simulations,
Comput. Fluids 35, 820.

Tölke, J., Prisco, G. D. & Mu, Y. (2013). A lattice Boltzmann method for
immiscible two-phase stokes flow with a local collision operator. Comput.
Math. Appl. 65, 864–881.

Truskey, G. A., Yuan, f., and Katz, D.F. (2009). Transport Phenomena in
Biological Systems, 2nd ed. Pearson Prentice Hall, Upper Saddle River.

118

Ubink, O. (1997). Numerical Prediction of Two Fluid Systems with Sharp
Interfaces, Ph.D. Thesis, Imperial College of London.

Valvatne, P.H., Blunt, M.J. (2004). Predictive pore scale modeling of two
phase flow in mixed wet media. Water Resour. Res. 40, W07406.

Verma, A. (2021). Curvelet transform analysis and denoising of images using
MATLAB; https://www.matlabcoding.com/2020/04/curvelet-transform-
analysis and.html.
Wildenschild, D., Armstrong, R.T., Herring, A.L., Young, I., Young, I.M.,
Carey, J.W. (2014). Exploring capillary trapping efficiency as a function
of interfacial tension, viscosity, and flow rate. Energy Procedia 4, 4945.

Wilkinson, D., Willimsen, J.F. (1983). Invasion percolation: A new form of
percolation theory. J. Phys. A 16, 3365.

Woiselle, J.-L. Starck, and J. Fadili. (2010). 3D curvelet transforms and
astronomical data restoration, Appl. Comput. Harmonic Anal. 28, 171.

Ying, L., Demanet, L., and Cand´es, E. (2005). 3D discrete curvelet
transform, in Proc. SPIE 5914, Wavelets XI, 591413;
https://doi.org/10.1117/12.616205.

119

APPENDIX: PARAMETERS OF THE MULTIRELAXATION
TIME LATTICE –BOLTZMANN SCHEME

We partition the 19 × 19 matrix M into four submatrices:

𝑀 = [
𝑀 1
𝑀 2
𝑀 3
𝑀 4
] . (A.1)

Then, the four submatrices are given by

𝑀 1
=
[

1 1 1 1
−30 −11 −11 −11
12 −4 −4 −4
0 1 −1 0
0 −4 4 0
0 0 0 1
0 0 0 −4
0 0 0 0
0 0 0 0
0 2 2 −1

1 1 1 1 1
−11 −11 −11 8 8
−4 −4 −4 1 1
0 0 0 1 −1
0 0 0 1 −1
−1 0 0 1 1
4 0 0 1 1
0 1 −1 0 0
0 −4 4 0 0
−1 −1 −1 1 1
]

, (A.2)

𝑀 2
=
[

1 1 1 1
8 8 8 8
1 1 1 1
1 −1 1 −1
1 −1 1 −1
−1 −1 0 0
−1 −1 0 0
0 0 1 1
0 0 1 1
1 1 1 1

1 1 1 1 1 1
8 8 8 8 8 8
1 1 1 1 1 1
1 −1 0 0 0 0
1 −1 0 0 0 0
0 0 1 −1 1 −1
0 0 1 −1 1 −1
−1 −1 1 1 −1 −1
−1 −1 1 1 −1 −1
1 1 −2 −2 −2 −2
]

, (A.3)

120

𝑀 3
=
[

0 −4 −4 2
0 0 0 1
0 0 0 −2
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0

2 2 2 1 1
1 −1 −1 1 1
−2 2 2 1 1
0 0 0 1 −1
0 0 0 0 0
0 0 0 0 0
0 0 0 1 −1
0 0 0 −1 −1
0 0 0 0 0
]

, (A.4)

and
𝑀 4
=
[

1 1 1 1
1 1 −1 −1
1 1 −1 −1
−1 1 0 0
0 0 0 0
0 0 1 −1
1 −1 −1 1
1 1 0 0
0 0 1 1

1 1 −2 −2 −2 −2
−1 −1 0 0 0 0
−1 −1 0 0 0 0
0 0 0 0 0 0
0 0 1 −1 −1 1
−1 1 0 0 0 0
−1 1 0 0 0 0
0 0 1 −1 1 −1
−1 −1 −1 −1 1 1
]

, (A.5)

Deleting the superscript eq for convenience and brevity, the nineteen entries
of the equilibrium moments m
eq
are given by

m0 = ρ , m1 = −11ρ + 19ρ0|v|2 , m2 = 3ρ − 11ρ0|v|2/2 , (A.6)
m3 = ρ0vx , m4 = −2m3/2 , m5 = ρ0vy , (A.7)
m6 = −2m5/3 , m7 = ρ0vz , m8 = −2m7/3 , (A.8)
m9 = ρ0(2 vx
2
− vy
2
– vz
2
) , m10 = −3m9/2 , m11 = ρ0(vy
2
– vz
2
) , (A.9)
m12 = −m11/2 , m13 = ρ0vxvy , m14 = ρ0vyvz , (A.10)
m15 = ρ0vzvx , m16 = m17 = m18 = 0 , (A.11)

with ρ0 being the mean fluid density in the system, which is taken to be unity
in simulations.
Since the relaxation rates matrix S is diagonal, we set Sii = Si for
brevity. Then, the nineteen entries of S are given by, S1 = S2 = S9 = S10 =
S11 = S12 = S13 = S14 = S15 = −Sν,
S4 = S6 = S8 = S16 = S17 = S18 = −Sm, and S0 = S3 = S5 = S7 = 0, with

121

1
𝑆 𝑣 =
3𝑣 |v|
2
𝛥𝑡
+0.5 (A.12)

where ν is the kinematic viscosity. Sm is given by (Ginzburg, I., 2005;2006):

𝑆 𝑚 =8
2−𝑆 𝑣 8−𝑆 𝑣 . (A.13)

The color gradient C, mentioned in the main text, is defined by

𝐶 =
3
|v|
2
𝛥𝑡
∑ 𝜔 𝑖 𝜑 (𝑥 +v
𝑖 𝛥𝑡 ,𝑡 )
𝑖 , (A.14)

where 𝜔 𝑖 are the weight coefficients in the D3Q19 LB model [38], which are
given by

𝜔 𝑖 = (𝑥 )=
{

1
3
𝑖 =0
1
18
𝑖 =1,2,···,6
1
36
𝑖 =7,8,···,18
(A.15)

Defining a vector n by, n = (nx, ny, nz) = C/|C|, the nonzero entries of m
s

added to the stress-related entries of equilibrium moments m
eq
are given by
(Chen, Y., et. al., 2018)

𝑚 1
𝑠 =−19𝜎 |𝐂 |, 𝑚 9
𝑠 = 0.5𝜎 |𝑪 |(2𝑛 𝑧 2
−𝑛 𝑦 2
−𝑛 𝑧 2
),𝑚 11
𝑠 = 𝜎 |𝐂 |(𝑛 𝑦 2
−𝑛 𝑧 2
) (A.16)
𝑚 13
𝑠 = 𝜎 |𝑪 |(𝑛 𝑥 𝑛 𝑦 ), 𝑚 14
𝑠 = 𝜎 |𝑪 |(𝑛 𝑦 𝑛 𝑧 ),𝑚 15
𝑠 = 𝜎 |𝑪 |(𝑛 𝑧 𝑛 𝑧 ) . (A.17)

Abstract (if available)

Abstract Flow, transport, reaction, adsorption and deformation (FTRAD) are phenomena that occur in a wide variety of heterogeneous materials and media, ranging from biological membranes, to composite solids, and porous formations. They have been studied for decades based on models of heterogeneous media that involved approximations, simplifications, and even unjustified assumptions. With advances in instrumentation, obtaining highly resolved images of such media has become feasible, opening the way to direct image-based simulation of the FTRAD. Due to complexity of such images, however, the simulations still require long intensive computations, limiting the size of the image that can be used. ? We propose a novel methodology based on curvelet transformation (CT), a multiscale and multidirectional transformation that provides an optimal sparse representation of objects?images in our case. Our methodology is described as following, the image (i) is taken to the curvelet space; (ii) all the curvelet coefficients that are smaller than a threshold are set to zero, and (iii) the image is reconstructed. We then demonstrate that simulation of the FTRAD in the reconstructed image yields as accurate results as those in the original one, while speeding up the computation by two orders of magnitude or better, hence taking a very significant step toward achieving the goal of making direct image-based simulation of the FTRAD a standard practice. ? The research presented in this dissertation demonstrates the virtues of the method, by testing our proposed method on several different applications. First, to demonstrate the increase in the computational efficiency of simulation of the flow and transport processes, we solve for unsteady-state diffusion problem and a single-phase fluid flow in complex images of porous media using our method and compare the results with the original image solution. ? We then extend our proposed method for speeding-up simulation of multiphase fluid flow in images of porous media by curvelet transformations, which are specifically designed for processing of images that contain complex curved surfaces. Most porous media contain correlations in their morphology and, therefore, their images carry redundant information that, in the curvelet transform space, can be removed efficiently and accurately in order to obtain a coarser image with which the computations are far less intensive. We utilize the methodology to simulate two-phase flow of oil and water in two-dimensional digital images of sandstone and carbonate samples, and demonstrate that while the results with the curvelet-processed images are as accurate as those with the original ones, the computations are speeded up by a factor of 110-150. Thus, the methodology opens the way toward achieving the ultimate goal of simulation of multiphase flow in porous media. ? Next, we use the same approach that utilizes curvelet transformations (CTs) for speeding-up simulation of multiphase flow on 3D images of porous media. As described previously, the CTs are designed for denoising of images that contain complex curved surfaces, such as those of heterogeneous porous media. This is possible because the morphology of most porous media contain extended correlations, implying that their images carry redundant information that can be eliminated by a suitable CT. As a result, simulation of multiphase flow in the coarser images are sped up very significantly. The new approach is used to simulate two-phase flow of CO? and brine in a 3D image of Berea sandstone. We demonstrate that the results with the CT-processed image are as accurate as those with the original one, while computations are significantly faster by a speed-up factor between one order of magnitude and more than two orders of magnitude. ? We also use the same methodology on a prime method for carrying out multiphase fluid flow simulations, which is the color-fluid lattice Boltzmann method with multirelaxation time (CFLB-MRT) collision operator. The simulations are, however, time consuming and intensive. We propose a method to accelerate image-based computations with the CFLB-MRT method. In the proposed method, the 3D image is preprocessed by curvelet transforming it and eliminating those details that do not contribute significantly to multiphase flow, which is done by thresholding the image. After inverting the coarser image back to the real space, it is utilized in the simulation of multiphase flow by the CFLB-MRT approach. As the test of the method, we carry out simulation of a two-phase flow problem in which the porous media are initially saturated by brine or water, which is then displaced by CO? or oil, injected into the pore space. The simulations are carried out with two types of sandstone. We show that the method accelerates the computations significantly by a speed-up factor of up to 35. ? The results in all applications indicate that excellent agreement between the solutions in the original images, and the less noisy and coarser images by two thresholds, while the speed of the computations increases significantly.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Machine-learning approaches for modeling of complex materials and media

PDF

Chemical and mechanical deformation of porous media and materials during adsorption and fluid flow

PDF

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

PDF

Effective flow and transport properties of deforming porous media and materials: theoretical modeling and comparison with experimental data

PDF

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

PDF

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

PDF

Flow and thermal transport at porous interfaces

PDF

Efficient stochastic simulations of hydrogeological systems: from model complexity to data assimilation

PDF

Investigation of adsorption and phase transition phenomena in porous media

PDF

Inverse modeling and uncertainty quantification of nonlinear flow in porous media models

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

PDF

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

PDF

Asphaltene deposition during co₂ injection in enhanced oil recovery applications

PDF

Modeling and simulation of complex recovery processes

PDF

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

PDF

Efficient connectivity assessment of heterogeneous porous media using graph theory

PDF

Molecular dynamics studies of protein aggregation in unbounded and confined media

PDF

On the transport dynamics of miscible fluids in porous media under different sources of disorder

PDF

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

PDF

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

PDF

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

Asset Metadata

Creator Aljasmi, Abdullah (author)

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

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Chemical Engineering

Degree Conferral Date 2021-08

Publication Date 07/08/2021

Defense Date 06/17/2021

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

Tag curvelet transformation,fluid flow and transport,OAI-PMH Harvest,porous media

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Sahimi, Muhammad (committee chair), Jessen, Kristian (committee member), de Barros, Felipe (committee member)

Creator Email aljasmi@usc.edu,toaljasmi@gmail.com

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

Unique identifier UC15277059

Identifier etd-AljasmiAbd-9710.pdf (filename)

Legacy Identifier etd-AljasmiAbd-9710

Document Type Dissertation

Format application/pdf (imt)

Rights Aljasmi, Abdullah

Type texts

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

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

Repository Name University of Southern California Digital Library

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

Repository Email cisadmin@lib.usc.edu