issertationThe Pennsylvania State University
The Graduate School
College of Earth and Mineral Sciences

PORE-SCALE IMAGING AND LATTICE BOLTZMANN MODELING OF SINGLEAND MULTI-PHASE FLOW IN FRACTURED AND MIXED-WET PERMEABLE
MEDIA

A Dissertation in
Energy and Mineral Engineering by Christopher James Landry

© 2013 Christopher James Landry

Submitted in Partial Fulfillment of the Requirements for the Degree of

Doctor of Philosophy

May 2013

The dissertation of Christopher James Landry was reviewed and approved* by the following:

Zuleima T. Karpyn
Associate Professor of Petroleum and Natural Gas Engineering
Dissertation Adviser
Chair of Committee

Li Li
Assistant Professor of Energy and Mineral Engineering

Russell T. Johns
Professor of Petroleum and Natural Gas Engineering

Maria Lopez de Murphy
Associate Professor of Civil Engineering

Luis Ayala
Associate Professor of Petroleum and Natural Gas Engineering
Associate Department Head for Graduate Education

*Signatures are on file in the Graduate School

ii

ABSTRACT

Three investigations of pore-scale single-phase and multiphase flow in fractured porous media and mixed-wet porous media are presented here. With an emphasis on validating and utilizing lattice Boltzmann models in conjunction with x-ray computed microtomography.

The objective of the first study is to investigate fracture flow characteristics at the pore-scale, and evaluate the influence of the adjacent permeable matrix on the fracture’s permeability.
We use X-ray computed microtomography to produce three-dimensional images of a fracture in a permeable medium. These images are processed and directly translated into lattices for single-phase lattice Boltzmann simulations. Three flow simulations are presented for the imaged volume, a simulation of the pore space, the fracture alone and the matrix alone. We show that the fracture permeability increases by a factor of 15.1 due to bypassing of fracture choke points through the matrix pore space. In addition, pore-scale matrix velocities were found to follow a logarithmic function of the distance from the fracture. Finally, our results are compared against previously proposed methods of estimating fracture permeability from fracture roughness, tortuosity, aperture distribution and matrix permeability.

In the second study we present a pore-scale study of relative permeability dependence on the strength of wettability of homogenous-wet porous media, as well as the dependence of relative permeability on the distribution and severity of wettability alteration of porous media to a mixedwet state. A Shan-Chen type multicomponent multiphase lattice Boltzmann model is employed to determine pore-scale fluid distributions and relative permeability.

iii

Mixed-wet states are

created – after pre-simulation of homogeneous-wet porous medium – by altering the wettability of solid surfaces in contact with the non-wetting phase. To ensure accurate representation of fluid-solid interfacial areas we compare LB simulation results to experimental measurements of interfacial fluid-fluid and fluid-solid areas determined by x-ray computed microtomography imaging of water and oil distributions in bead packs (Landry et al. 2011). The LB simulations are found to match experimental trends observed for fluid-fluid and fluid-solid interfacial areasaturation relationships. The relative permeability of both fluids in the homogenous-wet porous media is found to decrease with a decreasing contact angle. This is attributed to the increasing disconnection of the non-wetting phase and increased fluid-solid interfacial area of the wetting phase. The relative permeability of both fluids in the altered mixed-wet porous media is found to decrease. However the significance of the decrease is dependent on the connectivity of the unaltered solid surfaces, with less dependence on the severity of alteration.

In the third study we present sequential x-ray computed microtomography (CMT) images of matrix drainage in a fractured sintered glass granule pack. Sequential imaging captured the capillary-dominated migration of the non-wetting phase front from the fracture to the matrix in a brine-surfactant-Decane system. The sintered glass granule pack was designed to have minimal pore space beyond the resolution of CMT imaging, so that the pore space of the matrix connected to the fracture could be captured in its entirety. The segmented image of the pore space was then directly translated to a lattice to simulate the transfer of fluids between the fracture and the matrix using lattice Boltzmann (LB) modeling. This provided us an opportunity to validate the modeling technique against experimental images at the pore-scale. Although the surfactant was found to alter the wettability of the originally weakly oil-wet glass to water-wet,

iv

the fracture-matrix fluid transfer is found to be a drainage process, showing little to no countercurrent migration of the oil-phase.

The LB simulations were found to closely match

experimental rates of fracture-matrix fluid transfer, equilibrium saturation, irreducible wetting phase saturation and fluid distributions.

v

TABLE OF CONTENTS

List of Tables …………………………………………………………………..

ix

List of Figures ………………………………………………………………….

ix

Acknowledgements ……………………………………………………………

xv

CHAPTER 1: INTORDUCTION TO PORE-SCALE STUDIES …………….

1

CHAPTER 2: SINGLE-PHASE LATTICE BOLTZMANN …………………
SIMULATIONS OF PORE-SCALE FLOW IN FRACTURED
PERMEABLE MEDIA

9

2.1 Summary .……………………………………………………………….

9

2.2 Introduction …………………………………………………………….

9

2.3 Materials and Methods …………………………………………………

17

2.3.1 X-ray Computed Microtomography Scanner …..………………..

17

2.3.2 Fractured Permeable Media ….…...………………………………

18

2.3.3 Image Processing ….……………………………………………..

20

2.3.4 Fracture Pore Space Identification ...……………………………..

22

2.3.5 Surface Area and Volume ….……………………………………..

24

2.3.6 Aperture and Fracture Roughness .……………………………….

24

2.3.7 Lattice Boltzmann Simulations .………………………………….

28

2.4 Results ..………………………………………………………………..

31

2.4.1 Fracture and Matrix Pore Space Geometry .………………………

31

2.4.2 LB Simulation Results ……..……………………………………..

33

2.4.3 Fracture Permeability Estimation ………………………………..

36

2.5 Conclusions ….………………………………………………………..

43

vi

2.6 References …………..…………………………………………………

46

CHAPTER 3: RELATIVE PERMEABILITY OF HOMOGENOUS-WET ….
AND MIXED-WET POROUS MEDIA AS DETERMINED BY
PORE-SCALE LATTICE BOLTZMANN MODELING
3.1 Summary ..………………………………………………………………

55

3.2 Background ….…………………………………………………………..

56

3.3 Materials and Methods ….……………………………………………….

64

3.3.1 Experimental Measurements ……………………………………....

64

3.3.2 Single-Phase BGK Lattice Boltzmann Model ……………………..

65

3.3.3 Shan-Chen Multicomponent Lattice Boltzmann Model ………….

67

3.3.4 LB Model Implementation …………………………………………

69

55

3.4 Results …………………………………………………………………… 72
3.4.1 Fluid-Fluid and Fluid-Solid Interfacial Areas ……………………..

72

3.4.2 Pore Aperture and Fluid Distribution ………………………………

77

3.4.3 Two-Phase Flow in a Circular Pore ………………………………..

79

3.4.4 Relative Permeability of Homogenous-Wet States ……………….

81

3.4.5 Relative Permeability of Mixed-Wet States ……………………….

84

3.5 Conclusions ……………………………………………………………..

92

3.6 References ……..…………………………………………………………

95

CHAPTER 4: LATTICE BOLTZMANN MODELING AND ………………....
4D X-RAY COMPUTED MICROTOMOGRAPHIC IMAGING OF
FRACTURE-MATRIX FLUID TRANSFER

106

4.1 Summary ……………………………………………………………….

106

4.2 Introduction ……………………………………………………………...

107

4.3 Materials and Methods …………………………………………………..

110

vii

4.3.1. Fractured Porous Medium …………………………………………

110

4.3.2 Procedure …………………………………………………………..

111

4.3.3 CMT Imaging ……………………………………………………....

113

4.3.4 Image Processing …………………………………………………...

113

4.4 Results ………………………………………………………………….... 115
4.4.1 Pore Space Analysis ………………………………………………... 115
4.4.2 Immiscible Displacement …………………………………….…….

121

4.4.3 LB Simulation of Fracture-Matrix Fluid Transfer ………….……...

127

4.4.4 Pore-Scale Experimental and Model Results ………………..……..

130

4.5 Conclusions ……….……………………………………………………... 137
4.6 References …….…………………………………………………………. 140
CHAPTER 5: SUMMARY OF FINDINGS AND RECOMMENDATIONS ….
FOR FUTURE WORK

viii

148

LIST OF TABLES
Table 2.1.

Pore space geometry and fracture characteristics.

…………………33

Table 2.2.

Summary of LB permeability measurements and fracture permeability estimates.

…………………41

Table 4.1.

Temperature cycles for sintering of glass granules.

…………………110

LIST OF FIGURES
Figure 1.1

An example segmented CMT image showing the pore space voxels in white (A). Lattice Boltzmann models simulate fluids as swarming particles (B) on a discrete lattice, where each pore voxel is represented by a LB fluid node (C). The complex solid boundaries are also retained by translating the CMT image 1-to-1 voxel-tonode with the solid boundaries being translated to bounce-back nodes, shown in yellow (D).

…………………4

Figure 2.1.

Illustration of the double-sided ‘log-splitting’ approach used to propagate a rough fracture in the porous polyethylene rod.

…………………19

Figure 2.2.

Demonstration of simple thresholding image segmentation. At the left is a raw CMT image slice, the raw data is cropped (shown as dashed lines) to remove portions of fracture that were in contact with the splitting tool. The middle plot shows the frequency of CT registration numbers, and Gaussian fits to the solid and pore phases in the image, the intercept of these Gaussian fits is used to segment the image. At the right is a segmented image, in which gray represents the pore space and white represents the solid or grain.

…………………21

ix

Figure 2.3.

Example slices from the fracture space identification procedure. The left image shows the pore space (light gray) after 8 cycles of erosion, at which point all of the matrix pore space is removed, this locates the fracture
(dark gray). The middle image shows the pore space in contact with the fracture space after 2 cycles of reversedout erosions, the erosion cycle one short of the fracture pore space becoming connected with the matrix pore space; this identifies the starting point for expanding the fracture pore space out. The right image shows the fracture pore space after 6 cycles of expansion, the expansion cycle one short of the fracture pore space making contact with the grains or solid, this defines the fracture. …………………23

Figure 2.4.

Illustration of the measure of the vertical aperture of the fracture ( ), the perpendicular aperture ( ), the fracture profile distances from edge of image (
) and the bandwidth window ( ). The light gray is the matrix pore space, and the dark gray is the fracture.

…………………24

Figure 2.5.

Fracture aperture distribution. The small initial peak is a result of the matrix pore aperture distribution impacting the aperture distribution of the fracture at fracture boundaries. The distribution is overall log-normal, consistent with fractures in natural permeable media.

…………………26

Figure 2.6.

Poiseuille flow through a channel with a radius velocity profile, comparison of LB results and analytical solution. These results show the LB simulation is closely approximating flow in the small channel. …………………31

x

Figure 2.7.

Determination of Hurst coefficient using the variable bandwidth method. Data points are shown for four slices, the linear function shown (
) was determined from 16 equidistant slices spanning the range of the simulated pore volume. The Hurst exponent is the slope of the fitted linear function shown.

…………………32

Figure 2.8.

Normalized x-velocity glyph visualization of the fracture-alone (above) and the all-pore (below) simulations; the glyphs are larger and farther to the red side of the spectrum for greater velocities.
Th28maximum normalized velocity shown in the gly29s is .

…………………35

Figure 2.9.

̅̅
Mean normalized velocities (̅̅̅ ̅ ̅ ̅̅̅) as a function of distance from the fracture ( ), and logarithmic function fits. …………………27

Figure 3.1.

Segmented CMT image of the bead pack (A) with solid voxels colored white, and the corresponding lattice of the
LB model (B) with the bounce-back nodes in blue.

…………………71

Figure 3.2.

Specific fluid-fluid interfacial areas as a function of wetting phase saturation for the LB simulations and experimental CMT images.

…………………76

Figure 3.3.

Fractional fluid-solid interfacial areas as a function of wetting phase saturation for the LB simulations and experimental CMT images, the wetting fluid-solid interfacial areas are shown in black and the non-wetting fluid-solid interfacial areas are shown in gray.

…………………76

Figure 3.4.

Pore aperture and fluid distribution of initial

…………………78

xi

homogenous-wet states.

Figure 3.5.

Relative permeability in a circular pore with radius = 11
; LB simulation results and semi-analytical solution curves. …………………81

Figure 3.6.

Relative permeability determined by LB simulations for the initial homogenous-wettability states. The relative permeability of the wetting and non-wetting phase are shown in black and gray, respectively.

…………………82

Figure 3.7.

Images of the LB lattice with initial homogenouswettability (A), the fluid distribution of the non-wetting phase at the end of the simulation of the initial homogenous-wettability with and (B), and the alteration of the bounce-back densities of the nodes in contact with the non-wetting fluid (C). This image also summarizes the three parameters determining the mixed-wet state, initial wettability (A), saturation of alteration (B), and severity of alteration (C).

…………………86

Figure 3.8.

Relative permeability determined by LB simulations for the mixed-wet states of the originally homogenous-wet states, (A)
, (B)
, (C)
, and (D)
, with increasing severity of alteration. The relative permeability of the wetting and non-wetting phase are shown in black and gray, respectively.

…………………88

Figure 4.1.

Image of the fracture with the ROI highlighted (A), and the corresponding mean fracture aperture and contact area profiles with the ROI indicated by arrows.

…………………116

xii

Figure 4.2.

Image of the matrix pore space showing the fracture saturated with water phase and the matrix saturated with the oil phase, the ROI volume is highlighted (A). Also the corresponding porosity (B) and mean aperture (C) profiles with the ROI indicated with arrows.

…………………118

Figure 4.3.

Example 0.1183 cm3 image of sphere-packing method used to measure pore aperture distribution. The individual pore apertures are measured as the radius of the spheres seen in this image. The red, purple, teal, green, yellow, and light green represent pore aperture measurements between
,
,
,
,
, and mm, respectively. The background shows the solid in black and the pore space in white.

…………………120

Figure 4.4.

The fracture aperture distribution (left) and the pore aperture distribution (right) for the portion of the core to be imaged sequentially during water invasion of the matrix. …………………121

Figure 4.5.

Image of oil phase (white), water phase (black), glass
(gray) contact in the fracture after initial water phase injection. The contact angle appears in this image to be large, near neutral wetting conditions.

…………………122

Figure 4.6.

CMT images from the sequential scanning of the surfactant water phase injection showing a volume rendering (left) and a slice from a height of 1.0704 cm
(right). In the volume renderings the water phase is dark orange, the solid is light orange and the oil phase is white. In the slices the water phase is light gray, the glass is dark gray, and the oil phase is white.

…………………124

xiii

Figure 4.7.

Non-wetting phase saturation of the right and left matrix as determined from sequential CMT imaging. Early time refers to displacement before the interfacial tension reduction has stabilized. The linear fits shown represent a displacement following the rate function,
⁄

…………………126

.

Figure 4.8.

Volume images of the 270x150x150 voxel3 (0.4256 x
0.2365 x 0.2365 cm3) sample volumes of the pore space for the left and right matrix to be translated 1-to-1, voxel to lattice node, for lattice Boltzmann simulation (A). LB simulation setup showing the 5 lattice layer thick pressure boundaries, with the fracture represented as the non-wetting phase occupied pressure boundary (B).

…………………128

Figure 4.9.

Non-wetting phase saturation of the right and left matrix as determined from sequential CMT imaging and LB drainage simulations. The only variable between the LB simulations for the right and left matrix is the pore space geometry. …………………131

Figure 4.10.

Comparison of experimental (left) and LB simulation
(right) vertical saturation profiles for the right (A, B) and left (C, D) matrix, and horizontal saturation profiles for the right (E, F) and left (G, H) matrix.

…………………133

Figure 4.11.

Sample images of the LB simulations of fracture-matrix transfer. The non-wetting phase saturation front and simulation pore space are shown, the wetting phase is occupies the pore space not occupied by the non-wetting phase. …………………135

xiv

ACKNOWLEDGEMENTS
This material is based upon work supported by the National Science Foundation under
Grant No. 0747585 and the Marathon Alumni Centennial Graduate Fellowship. I would like to thank Dr. Orlando Ayala of the University of Delaware for providing computational resources, and for offering his expertise and time, many of the simulations shown here would not have been possible without his help. I would also like to thank Soheil Saraji of Piri Research Group at the University of Wyoming for providing interfacial tension measurements. I would also like to thank my doctoral committee for their thoughtful comments, and insightful discussion. Last I would like to thank my advisor, Dr. Zuleima
Karpyn, who had the unenviable task of keeping me on task.

xv

CHAPTER 1: INTRODUCTION TO PORE-SCALE STUDIES

The goal of studying multiphase fluid flow at the pore-scale in porous media and fractured porous media is to improve fluid transport prediction at the field-scale. Generally prediction of fluid transport at the field scale is determined by models that discretize the formation of interest
– attributing macroscale flow properties to individual gridblocks – and simulating flow under varying scenarios. Macroscale flow properties are conventionally determined by measuring fluid flow in cores take from the formation of interest. Pore-scale studies attempt to elucidate the connection between pore space characteristics and fluid mobility, with the intention of determining macroscale flow properties to be utilized in field-scale fluid transport prediction.

Accurate determination of macroscale flow properties such as permeability, relative permeability
– saturation relationships and capillary pressure – saturation relationships are fundamental to accurate field-scale fluid transport prediction. There are numerous techniques used to measure macroscale flow properties of porous media and fractured porous media, all of which determine macroscale flow properties by measuring fluid flux through cores taken from the formation of interest under varying applied boundary conditions. Although direct measurement of macroscale flow properties is straight forward there are downsides to these methods. For one they require cores that can be expensive to procure, and often these cores only offer a small sampling of a vast heterogeneous formation. The measurements themselves are also expensive, and in the case of fractured porous media, measuring flow in the fracture and porous matrix separately can be very difficult, if not impossible. Most importantly, these methods provide measurements, but no

1

insight into the cause for resulting macroscale flow measurements. The porous media is treated as essentially a “black-box”.

Pore-scale studies peer into this black-box to elucidate the relationship between pore space characteristics and resulting macroscale fluid flow. Although the physics of multiphase fluid flow is fairly well understood, the complexity of the solid boundary conditions encountered in pore spaces render any possibility of analytical solutions of fluid flow to be beyond practicality.
Complex non-invasive imaging techniques, such as x-ray computed microtomography (CMT) allow us to not only image pore spaces at resolutions great enough to capture the shape of individual pores, but also the distribution of fluids within these pores. Being able to image pore spaces and fluid distributions at the pore-scale provides insight into multiphase fluid flow, as in we are able to peer into the black-box, but this does not alleviate some of the issues we encountered with conventional measurements of multiphase fluid flow. It is still expensive, and requires cores.

Pore-scale models offer us a powerful tool to not only gain insight into the connection between pore space characteristics and resulting macroscale fluid flow. Pore-scale models use either imaged pore spaces or synthesized pore spaces to simulate fluid flow at the pore-scale. There are many different types of pore-scale models in existence, including classic computational fluid dynamics (discrete solutions of Navier-Stokes equation), pore-network models, smoothed particle hydrodynamics, the level-set method and lattice Boltzmann (LB) models. The best choice of model to use when investigating fluid flow at the pore-scale is dependent on the objective of the investigation, here we employ lattice Boltzmann models. The complexity of the

2

solid boundaries simulated in LB models is only limited by computational power, or in the case of imaged pore spaces, the resolution of the image.

This is an improvement over more

commonly used pore-network models which represent the pore space as a network of discrete geometric shapes (i.e. cylinders, triangular tubes, square tubes, spheres etc.). Briefly, LB models can be described as simulating fluids as swarming particles on a discrete lattice. A detailed description of the model used in the investigations presented here can be found in Chapter 3.
The complex solid boundaries of pore spaces imaged by CMT can be translated 1-to-1 voxel-tonode onto LB lattices To better visualize what is meant by maintaining the complexity of the solid boundaries of pore space in LB models please refer to figure 1.1. In this figure we show a
CMT image of a pore space and the corollary LB lattice. LB models simulate fluid flow in a purely local manner, no long range information is required during the evolution of the model, each node only requires knowledge of adjacent nodes. This means LB models have an inherent parallelism, and simulations can carried out across many processors using parallel computing techniques. One of the greatest challenges encountered in the implementation of LB models is the computational expense. Although LB models are inherently parallel and simulations can be spread across hundreds, if not thousands of processor cores, the computational demand is so great that multiphase simulations on 1003 node3 lattices can require upwards of a week on 200 cores to reach convergence. However, it should be noted the only expense is computational. It is possible to run numerous scenarios at only the cost of computer time.

3

A)

B)

Swarming
Particles

C)

Discrete LBE

D)

Figure 1.1.

An example segmented CMT image showing the pore space voxels in white (A). Lattice
Boltzmann models simulate fluids as swarming particles (B) on a discrete lattice, where each pore voxel is represented by a LB fluid node (C). The complex solid boundaries are also retained by translating the CMT image 1-to-1 voxel-to-node with the solid boundaries being translated to bounce-back nodes, shown in yellow (D).

There are two main objectives of the three investigations presented here in chapters 2, 3, and 4.
The first is to validate multi-phase pore-scale fluid distributions predicted by LB models against those imaged experimentally. One of the main goals of pore-scale studies is to better understand the connection between pore space characteristics and resulting macroscale flow, considering

4

this, it is imperative that we validate LB prediction of pore-scale fluid distribution. LB models have in many previous investigations been used to determine macroscale flow properties in the same way as is done conventionally, by imposing pressure boundaries and monitoring fluid flux.
In these studies LB models have been shown to accurately predict macroscale flow properties.
For further information regarding previous investigations into prediction of macroscale flow properties by LB please refer to subsection 3.2 Background. There are however only a handful of studies (subsection 3.2 Background) that have attempted to compare LB model pore-scale fluid distributions against experimentally imaged pore-scale fluid distributions. The second main objective is to take advantage of the intimate knowledge of fluid velocity and distribution offered by the LB model to perform numerical investigations into single-phase and two-phase fluid flow that would be difficult if not impossible to perform in a physical laboratory setting.

With these objectives in mind we present three investigations. In chapter 2 we take advantage of the intimate knowledge of velocity (velocity is known at every LB fluid node) offered by the LB model to investigate single-phase flow in fractures bounded by a porous matrix.

We simulate

flow in a fracture bounded by porous walls and the same fracture bounded by impermeable walls, to investigate the effect porous walls have on fracture permeability. Because we have intimate knowledge of velocity we simulate the flow in our fracture-matrix coupled system, and measure the flow in the fracture exclusively. This would be nearly impossible in a conventional, physical experimental system, because it is extremely difficult to separate the measurement of fluid flow in the fracture from that of the matrix. In this investigation we also compare the our results to existing fracture permeability estimations based on fracture characteristics, such as the mean aperture, standard deviation of the aperture distribution, wall roughness, tortuosity and

5

contact area. Since many formations of interest are fractured, and fluid flow is often dominated by flow in fractures in these formations, accurate estimations of fracture permeability are imperative to the success of field-scale fluid transport prediction.

In the third chapter we engage in validation of LB model prediction of fluid distribution by comparison to experimental results, and investigate fluid flow in mixed-wet states using the LB model. Experimental pore-scale fluid distribution data sets are sparse, which is one of the reasons so few investigations have attempted to compare model and experimental pore-scale fluid distributions. We use an experimental data set we previously published (Landry et al.
2011) to compare to our LB model results. And after presenting an argument for validation of our methods through comparison to experimental measurements, we engage in an investigation of wettability alteration to mixed-wet states and the effect this has on the relative permeability of the porous medium. We take advantage of the control of solid surface wettability (the adhesion force of every LB bounce-back node is tunable) offered by the LB model to investigate twophase flow in mixed-wet porous media.

Although numerical studies into the relative

permeability of mixed-wet porous media using pore-network models has existed for over twenty years (subsection 3.2 Background), this is the first thorough investigation of two-phase flow in mixed-wet porous media using LB to the author’s knowledge. Also much like in the first investigation, this investigation would be very difficult to perform in a conventional, physical experimental system, due to the difficulty in engineering the wettability of physical porous media. 6

In the fourth chapter we again engage in validation of LB model prediction of pore-scale fluid distribution by comparison to experimental results. In this investigation we image and simulate fracture-matrix fluid transfer at the pore-scale in a surfactant-brine-oil fluid system. Many formations of interest are fractured and the transfer of fluids from the high-permeability, lowstorage fractures to the low-permeability, high-storage matrix is of great interest. Predicting this migration requires a rigorous investigation into the mechanics that determine the transfer of fluids to and from the matrix. Understanding the physics of this process can be greatly enhanced by investigating fluid transfer at the interface of the fracture and matrix at the pore-scale. Many previous investigations at the pore-scale have used two-dimensional etched-glass micromodels to investigate fracture-matrix fluid transfer (see subsection 4.2 Introduction). We present a first of its kind three-dimensional study of fracture-matrix fluid transfer at the pore-scale.

We

sequentially image fracture-matrix fluid transfer at the pore-scale in a three-dimensional fractured porous medium using CMT. This is a novel use of three-dimensional CMT imaging, and the first time to the author’s knowledge that fracture-matrix fluid transfer has been captured in four dimensions at the pore-scale.

We compare pore-scale fluid distributions from this

experiment to LB model results to both validate the method of simulation and gain insight into the physics of fracture-matrix fluid transfer in a surfactant-brine-oil system.

The porous media used in these investigations is all synthetic, there are two reasons we do not use natural porous media in our investigations. First, the CMT imaging system has limited resolution. To compare LB and experimental results it is imperative that the model system is representative of the porous medium imaged, if there are pore spaces beyond the resolution of the CMT imaging system, they may affect the fluid flow in the experimental system in a way that

7

cannot be represented in the simulation.

Second, natural porous media can contain

heterogeneous surfaces that we cannot differentiate using the imaging techniques at our disposal.
Small differences in the wettability of natural porous media due to heterogeneous composition or surface precipitation may affect fluid flow in the experimental system in a way we cannot translate to the simulation. Investigations like those presented here, help to build a foundation of validation of pore-scale modeling techniques, and highlight the strengths and limits of the porescale models used. The ultimate goal is to apply these techniques to natural porous media, and free the determination of macroscale flow properties from the “black-box” of conventional measurement techniques.

8

CHAPTER 2: SINGLE-PHASE LATTICE BOLTZMANN SIMULATIONS OF PORESCALE FLOW IN FRACTURED PERMEABLE MEDIA

2.1 SUMMARY

The objective of this work is to investigate fracture flow characteristics at the pore-scale, and evaluate the influence of the adjacent permeable matrix on the fracture’s permeability. We use
X-ray computed microtomography to produce three-dimensional images of a fracture in a permeable medium. These images are processed and directly translated into lattices for single-phase lattice Boltzmann simulations. Three flow simulations are presented for the imaged volume, a simulation of the pore space, the fracture alone and the matrix alone. We show that the fracture permeability increases by a factor of 15.1 due to bypassing of fracture choke points through the matrix pore space. In addition, pore-scale matrix velocities were found to follow a logarithmic function of the distance from the fracture. Finally, our results are compared against previously proposed methods of estimating fracture permeability from fracture roughness, tortuosity, aperture distribution and matrix permeability.

2.2 INTRODUCTION

Current understanding of the physical phenomena concerning fluid flow through fractures is limited. Fluid flow in fractured formations depends on both the fracture system and the adjacent permeable rock, also called the rock matrix. Fractures control the overall conductivity

9

of the rock while the permeable matrix provides fluid storage capacity. Studying fracture structures, flow through fractures and fracture-matrix interactions presents significant challenges with important applications in areas such as geo-environmental remediation, geohazard mitigation, geothermal exploitation, and hydrocarbon production. In particular, development of unconventional hydrocarbon resources, such as coalbed methane, gas shales and tight sands, depends on our ability to represent fracture flow and fracture-matrix interactions for the design of recovery strategies and performance predictions. The deliverability of these unconventional hydrocarbon resources is highly controlled by the connectivity of natural fracture network, the accessibility of fractures to the matrix storing hydrocarbons and the fracture permeability. Of fundamental interest to the study of flow in fractured permeable media is the determination of fracture permeability from fracture and matrix pore space geometry.

The study of flow in fractures is, for the most part, limited to fractures in rock with impermeable fracture walls. The simplest approximation of single phase fluid flow in fractures with impermeable walls is given by the cubic law, the analytical solution to the NavierStokes equation for viscous flow through wide smooth parallel plates (no-slip walls). Given the wall is wide enough (Width >> aperture) and the length across which a pressure head is applied is long enough (Length > Width), the cubic law is defined as:

(

where aperture, )

is the flow rate through the fracture, w is the fracture width, is the dynamic viscosity,

is the fracture

is the pressure gradient, and in the context of Darcy’s

10

equation fracture permeability ( ), is defined as,

⁄

. The simplest modification

to the cubic law is to replace the aperture, with the mean aperture ( ̅ ), this is referred to as the alternate cubic law (ACL). This simple definition of flow through fractures does not include important factors influencing fluid flow in natural fracture geometries, such as, fracture roughness, asperities, anisotropies and aperture distribution. These factors have been shown to greatly influence flow in fractures (Auradou et al. 2005, Keller et al. 1995, Konzuk and
Kueper 2004, Nolte et al. 1989, Tsang and Witherspoon 1981, Tsang and Witherspoon
1983, Witherspoon et al. 1980). As the effect of these factors became apparent, empirical and semi-analytical modifications of the cubic law were made in the hopes of developing an estimation of fracture permeability from fracture geometry statistics. Witherspoon et al. 1980 studied flow across smooth and rough non-contacting surfaces and introduced a modification to the fracture permeability that included the fracture roughness,

where is

a surface roughness factor ranging from 1.04 to 1.65. After the introduction of fractal

geometry, it became apparent fractures possess self-affine fractal properties regardless of the host rock (Brown and Scholz 1985, Crandall et al. 2010, Poon et al. 1992, Schmittbuhl et al.
1995). The fractal property of fractures is often described by its Hurst exponent ( ), related to the fractal dimension (

) by

. With this knowledge, investigators began to study

flow in synthetically developed fractures that followed fractal geometries. The study of flow through rough fractures can be difficult to quantify and observe experimentally, thus investigators turned to recent developments in pore-scale fluid flow models, and computational

11

fluid dynamics (CFD) to numerically simulate flows in complex pore spaces. A variety of pore scale models have been developed, here we focus on lattice Boltzmann (LB) methods
(McNamara and Zanetti 1988), and to a lesser extent its forbearer, lattice gas cellular automata
(LGCA) (Frisch et al. 1986) methods.

Lattice techniques hold great promise in their application to the study of flow in permeable media; they are capable of simulating flow in any complex media – given a lattice of adequate resolution is used – as well as being easily used in parallel processing. LB and LGCA methods have been used to study single-phase flow in both permeable media (e.g., Boek and Venturoli
2010, Ferreol and Rothman 1995, Genabeek and Rothman 1996, Koponen et al. 1997, Maier et al. 1998, Manz et al. 1999, Pan et al. 2004) and fractures (e.g., Auradou et al. 2005, Boutt et al. 2006, Drazer and Koplik 2002, Eker and Akin 2006, Kim et al. 2003, Madadi and Sahimi
2003, Nazridoust et al. 2006, Stockman et al. 1997, Zhang and Kang 2004, Zhang et al. 1996).
Zhang et al. 1996 studied single-phase flow through synthesised self-affine fractal fracture using LCGA; they proposed an empirical estimation of fracture permeability as a function of mean aperture:

̅

with a value of
A

of 2.67, 3.15 and >4 for fractures with a Hurst exponent of 0.8, 0.3 and -0.5.

value of 2 would be expected for flow between parallel plates, the value of

is expected to

increase monotonically with roughness. Eker and Akin 2006 also investigated single phase flow in synthesized self-affine fractal fractures with Hurst exponents ranging from -0.21 to -0.50

12

using LB methods, and following the work of Zhang et al. 1996 proposed

to be a linear

function of fractal dimension,

These empirical predictions only take the fracture roughness into consideration, other factors governing fracture geometry such as asperities and aperture distribution must also be included in a complete prediction for the estimation of fracture permeability. Zimmerman and Bodvarsson
1996 attempted to relate the effective hydraulic fracture aperture ( ̅ ) to the statistics of aperture distribution. The effective hydraulic fracture aperture is a measured property of flow through a fracture, it is the equivalent aperture for flow between smooth plates for a measured flow rate
(Q),
(

̅

)

The relation between ( ̅ ) and the statistics of aperture distribution was derived from the application of Reynold’s lubrication equation. As they noted, asperities can be viewed as a part of the aperture distribution, however, it is convenient to treat the contact areas or asperities separately from the aperture distribution. The ̅ is proportional to the fracture permeability
(

̅ ⁄

), thus their prediction for the relation of the ̅ to fracture statistics is also a

definition of fracture permeability:

̅

̅ (

̅ )

13

is the standard deviation of the fracture aperture distribution and C is the fractional contact between fracture walls. The applicability of numerical results from synthesised fractures has merit, but validation of these methods requires studying real fracture geometries.
Imaging of fractures has become simplified by the employment of sophisticated three dimensional imaging techniques, such as X-ray computer microtomography (CMT).

The use of X-ray CMT in the study of flow in permeable media and fractures has become widespread, offering a non-invasive method to image complex pore geometries and fractures in three dimensions. High resolution images of pore space are often used in conjunction with porescale models to investigate pore scale flow mechanisms that are difficult to measure and observe. A large variety of permeable media has been imaged using CMT, including macro permeable clays (Heijs et al. 1995) soils (Perret et al. 1999), bead packs (Culligan et al.
2006), unconsolidated sands (Brusseau et al. 2006, Wildenschild et al. 2005), sandstones
(Auzerais et al. 1996, Nakashima et al. 2004) and limestones (Okabe and Blunt 2004). Many fractures have also been studied using CMT, including single brittle fractures in Berea sandstone (Karpyn et al. 2007), single planar artificial fractures of increasing aperture in limestone (Ketcham et al. 2010), lab-scale hydraulic fractures in limestone (Renard et al.
2009), natural fractures in limestone (Wennberg et al. 2009), and single brittle fractures in granite (Johns et al. 1993, Keller et al. 1999), to name a few.

Nazridoust et al. 2006 studied flow in a single tensile fracture in Berea sandstone imaged by
Karpyn et al. 2007. They proposed studying flow in fractures in a manner analogous to the

14

study of pressure losses due to friction in pipes and ducts. The resistance to flow is measured as a friction factor, analogous to fracture permeability. For flow between smooth parallel plates the friction factor (f) is defined as,

They introduced a new friction factor dependant on fracture tortuosity and aperture distribution, and fitted to CFD results for flows of varying velocity,

̅

where

̅

⁄

(

̅

is the Reynolds fracture number and

)

is the kinematic viscosity; this can be

related to the fracture permeability by the pressure drop,

(

̅

̅

where ̅ is the corrected mean fracture aperture, defined as, ̅ viscosity, and

is the tortuosity,

15

)

̅

,

is the dynamic

where

is the length of the tortuous actual path of flow and

is the length of the fracture,

arriving at a semi-analytical estimation of fracture permeability,

̅
(

)

̅

Crandall et al. 2010 also studied flow in a single fracture in Berea sandstone imaged by CMT using CFD. They used the same image of a fracture as Nazridoust et al. 2006, and attempted to include matrix effects on flow by making the CFD grid adjacent to the fracture permeable.
Following the work of Nazridoust et al. 2006 they redefined the friction factor for the fracture to include the flow through the matrix,

(

)

̅

̅

⁄̅

̅

and equivalently the fracture permeability,

̅ (

(
(

where

is the matrix permeability and

⁄̅ )

)

̅
̅

)

is the width of the matrix. Overall they found the

friction factor decreased with increasing matrix permeability, and the fracture permeability is directly proportional to the matrix permeability.

16

As of yet, there is no agreed upon approach for estimation of fracture permeability from fracture geometry and the flow behavior of the rock matrix. This is a fundamental issue, and a key building block for adequate description and simulation of fluid flow through fractured reservoirs. To the author’s knowledge, this is the first pore-scale study of flow in fractured permeable media in which a three-dimensional rough fracture is explicitly interconnected with the pore network of the adjacent permeable rock. We use CMT to image the fractured permeable media at resolutions high enough to capture the solid surface curvature of the grains; hence we are able image the entirety of the pore space including that of the matrix.
This image is translated directly into a lattice to be used by LB simulations of single phase flow. To study the effect of the matrix on flow in the fracture, and vice versa, three simulations are presented. One in which the entirety of the pore space is simulated (fracture in permeable rock), one in which the matrix pore space is removed and the flow is simulated in the fracture alone (analogous to a fracture in impermeable rock), and lastly five simulations are run on portions of the pore space that do not contain any part of the fracture and averaged to determine the permeability of the matrix alone. The results of these simulations are compared to the empirical and semi-analytical estimations of fracture permeability discussed above.

2.3 MATERIALS AND METHODS
2.3.1 X-ray Microtomography Scanner

The pore space of the fractured permeable medium was imaged using X-ray (CMT).
Scanning took place at Penn State’s Center for Quantitative X-ray Imaging (CQI). X-ray CMT imaging is a non-invasive imaging technique that produces a three-dimensional grid of voxels

17

– much like pixels for a two-dimensional image – each of which contain an integer value CMT registration number. The scan of the core in the scan presented here is composed of 800 slices, each of these slices being 1024 x 1024 voxels. Each slice of the three-dimensional grid of CMT registration numbers is produced by a mathematical reconstruction of 1440 two-dimensional images taken perpendicular to the axis of the core over the course of a full 360° core rotation.
CMT registration numbers produced by the mathematical reconstruction are dependent on the relative X-ray attenuation, or X-ray opaqueness, of the material occupying each voxel. X-ray attenuation by a material is a function of the density and apparent atomic number of the material. These images are segmented to identify the material occupying each voxel, to be further elaborated on in the image processing section. The half inch diameter cores in this investigation result in CMT scans with a voxel resolution of 13.606 x 13.606 x 12.405 μm3, the latter represents the thickness of the slice. However, to directly translate the CMT image of the pore space to a LB simulation requires a cubic lattice, thus the images are treated as having cubic voxels with side lengths of 13.606 μm.

2.3.2 Fractured Permeable Media

Porous polyethylene rods manufactured by Small Parts Inc., were used as the permeable medium. These rods are 0.127 cm in diameter and formed by sintering medium/coarse, angular/subangular polyethylene grains into rod form. This creates a permeable structure similar to a medium/coarse, well-sorted unconsolidated sand or sandstone. A synthetic permeable medium is used in lieu of a natural permeable medium to ensure the pore space of the matrix is resolvable. Polyethylene is a malleable material, thus a double- sided ‘log-splitting’

18

approach was used to create a fracture in the porous polyethylene rods (Figure 2.1), instead of more conventional techniques (i.e., modified Brazilian test) used to induce fractures on natural rock. Two splitting tools are driven along the desired fracture plane, not down the center of the rod, but along the sides. This results in a rough fracture propagating between the splitting tools.
The outer portions of the fracture that were in contact with the spitting tool are cropped from the image and not used in the LB simulation. This crop results in a total pore space lattice volume of 520 x 520 x 800 voxel3.

Figure 2.1.

Illustration of the double-sided ‘log-splitting’ approach used to propagate a rough fracture in the porous polyethylene rod.

19

2.3.3 Image processing

The raw, cropped CMT images are processed and segmented for analysis. The raw images are first filtered using a 3 x 3 x 3 kernel median filter to remove salt-and-pepper type noise. To quantify measurements of the pore space, the filtered images are segmented into the phases present, pore and solid. The contrast between the solid phase and pore space is sufficient for application of simple thresholding. Simple thresholding segments the image by declaring all voxels below a threshold value as belonging to one population and all those above to the other.
The threshold value is determined by fitting Gaussian curves to a frequency plot of the CT registration numbers (Figure 2.2) (Mees et al. 2003). The frequency plot contains three peaks, the first is associated with the pore space, and the second two are associated with the solid, with respect to increasing CT registration number. The intercept of the two Gaussian curves for the respective populations is used as the threshold value, here at a CT registration number of 17323. The two peaks of the solid represent the two different polyethylene particulates used in the construction of the permeable rod.

20

Figure 2.2.

Demonstration of simple thresholding image segmentation, at the left is a raw CMT image slice, the raw data is cropped (shown as dashed lines) to remove portions of fracture that were in contact with the splitting tool. The middle plot shows the frequency of CT registration numbers, and Gaussian fits to the solid and pore phases in the image, the intercept of these Gaussian fits is used to segment the image. At the right is a segmented image, in which grey represents the pore space and white represents the solid or grain.

The segmented images contain some sporadic solid voxels misidentified by the segmentation procedure in the pore space. To remove these sporadic voxels, the segmented 3D image of the pore space undergoes a percolation. A percolation finds all like voxels in contact. First a percolation is conducted on the solid voxels. The contact points between individual grains are beyond the resolution of the CMT images, thus this percolation finds all the connected solid voxels or ‘true’ solid voxels. The ‘negative’ of the connected solid voxels is then considered the pore space. To remove any isolated pores from the image of the pore space a second percolation is carried out on the negative of the connected solid voxels. The second percolation finds all pore space voxels in contact, ensuring the pore space to be simulated does not contain any unnecessary, non-interacting (isolated pores) pore space. This prepares the pore space for direct interpretation into a lattice for use by LB simulations.

21

2.3.4 Fracture Pore Space Identification

The identification of the fracture is non-trivial; unlike previous investigations of fractures, the distinction between the fracture and matrix is not readily apparent. To identify the pore space associated with the fracture pore space – also simply referred to as the fracture – a sequence of erosions, expansions and percolations were performed on the segmented images of the pore space. A cycle of erosion substitutes all the pore space edge voxels with solid voxels, thus removing a single layer of voxels from the pore space volume. Edge voxels are defined as those with voxel face contact between the two phases present, namely solid and pore. Expansion is the opposite of erosion.

The first step in the identification of the fracture pore space is the isolation of the locality of the fracture. If the fracture aperture is generally greater than the greatest distance between grains in the pores of the matrix, the fracture can be located by eroding the pore space until a single pore is remaining. This point was reached after eight cycles of erosion for the pore space presented here. This single pore only captures the core of the fracture pore space, to identify the entirety of the fracture this pore must be dilated out. Each cycle of erosion results in a loss of detail, to minimize this effect, the core of the fracture pore space is reversed-out by two cycles of erosion. It is important to note that this is not two cycles of expansion. We have the images preceding each erosion; we can retain loss of detail to these two cycles of erosions. The dilation cannot continue in this reversing-out manner. The next reversed-out erosion opens up the pore throats connecting to the matrix pore space, to continue in this manner would erroneously begin to identify the matrix pore space as fracture pore space. Instead, at this point

22

we perform cycles of expansion until the fracture pore space contacts the solid surface. After seven cycles of expansion the dilated fracture pore space overlaps solid voxels, thus six cycles of expansion defines the fracture. This sequence is summarized in Figure 2.3.

Figure 2.3.

Example slices from the fracture space identification procedure. The left image shows the pore space (light grey) after eight cycles of erosion, at which point all of the matrix pore space is removed, this locates the fracture (dark grey). The middle image shows the pore space in contact with the fracture space after two cycles of reversed-out erosions, the erosion cycle one short of the fracture pore space becoming connected with the matrix pore space; this identifies the starting point for expanding the fracture pore space out. The right image shows the fracture pore space after six cycles of expansion, the expansion cycle one short of the fracture pore space making contact with the grains or solid, this defines the fracture.

The flow in the matrix adjacent to the fracture is influenced by the high rate of flow in the fracture, and differs significantly from the flow in the rest of the matrix. To investigate the extent of this effect, the flow in the pore space adjacent to the fracture pore space is analyzed.
To measure this flow, the voxels adjacent to the fracture must be identified as a function of distance from the fracture. To do so, the fracture pore space is expanded to identify voxels at equal spatial distance from the fracture throughout the entirety of the pore space. These expansions identify voxels as a function of distance from the fracture, and provide a measure of flow velocity as a function of distance from the fracture.

23

2.3.5 Surface Area and Volume

The segmented images are used to measure surface areas and volumes of grains and pore space.
Volumes are calculated as the sum of all solid and pore voxels ( matrix (

⁄(

) is defined as,

⁄ , where

). The porosity of the

) and the fracture porosity (

) is defined as,

is the facture volume. The surface areas are determined using the image

processing software Avizo Fire 6.3. This software uses a modified marching cubes algorithm to
“wrap” a tetrahedral surface around the segmented image data for the solid phase. This surface quantifies the solid surface area (

). This is normalized to the bulk permeable volume of the

matrix to find the specific solid surface area (

) defined as,

bulk volume of the matrix defined as,

⁄

, where

is the

.

2.3.6 Aperture and Fracture Roughness

It is assumed that the slices of the CT scan run orthogonal to the vertical plane of the fracture; the measurement of apertures and fracture roughness are taken per slice of the scan.
There are two manners in which aperture can be measured, the perpendicular aperture (hp) and the vertical aperture (h) (Figure 2.4). The apertures measured here are vertical apertures, albeit the perpendicular aperture is not subject to fracture orientation, as the vertical aperture is, considering the velocities reported here follow a Cartesian coordinate system defined by the images, the vertical aperture was considered a more appropriate method of measurement.
The aperture distribution is shown in figure 2.5; the small initial peak is a result of the matrix pore aperture distribution impacting

the aperture distribution at fracture boundaries. The

24

fracture boundaries are defined by the matrix pore space, this will define the lower end of the fracture aperture distribution. This effect on the aperture distribution becomes negligible as the fracture aperture to matrix pore aperture ratio increases, but for our fracture the effect is readily apparent. The distribution is overall log-normal, consistent with fractures in natural permeable media (Karpyn et al. 2007, Keller et al. 1999, Ketcham et al. 2010). In this study, fracture apertures below the maximum pore aperture found in the matrix (~0.33 mm) are considered fracture contact areas. Pore apertures are measured in the same manner as fracture apertures. The definition of contact area used here is,

where

is the number of fracture apertures measured that are less than the maximum vertical

aperture of the matrix pores and

is the total number of measured fracture apertures.

25

Figure 2.4.

Illustration of the measure of the vertical aperture of the fracture (h), the perpendicular aperture (hp), the fracture profile distances from edge of image (Bj, Bj+s) and the bandwidth window (s). Note, The light gray is the matrix pore space, and the dark gray is the fracture.

0.03

Frequency

0.025
0.02
0.015
0.01
0.005
0
0

0.5

1

1.5

2

2.5

Fracture Aperture (h) (mm)

Figure 2.5.

The fracture aperture distribution, the small initial peak is a result of the matrix pore aperture distribution impacting the aperture distribution of the fracture at fracture boundaries. The distribution is overall log-normal, consistent with fractures in natural permeable media.

26

The fracture profile is defined by the fracture pore space in individual two-dimensional slices of the image. The roughness of the fracture profiles of individual slices is characterized by the
Hurst exponent ( ) also referred to as the roughness exponent or self-affine exponent. We use the variable bandwidth method to determine the Hurst exponent ( ) of the fracture profile for individual slices. We assume our fracture is well described by fractional Brownian motion, such that a scaling relationship exists between the size of the sampling window ( ) and the standard deviation of the fracture profile (

):

The Hurst exponent can be determined from the slope of a log-log plot of standard deviation
(

) vs. sampling window size ( ).

The standard deviation is calculated according to a

modified form of (Dougan et al. 2000):

√

where

∑ (

)

is the number of datapoints contained in the dataset, j is the index of the fracture profile

height ( ) (Figure 2.4), and wherein a measurement of

is a count of the number of samples that contain an asperity, is not possible. The assumption that a fracture is well described by

fractional Brownian motion and the variable bandwidth method is an appropriate method for the measurement of the Hurst exponent has been found to be generally true for fractures in natural rocks by previous investigators (Crandall et al. 2010, Renard et al. 2009).

The

tortuosity ( ), previously defined, requires a measurement of the mean length of a streamline flowing through the fracture ( ), this is found using the programme Paraview. The stream tracer function seeds a particle at the inlet and progress the particle in the velocity field. The

27

mean length of a streamline is determined from 50 streamlines seeded at regular intervals along the fracture inlet.

2.3.7 Lattice Boltzmann Simulations

Most petroleum engineers are accustomed to the ‘top-down’ methods of simulating fluid flow, in which partial differential equations defining fluid flow are discretized using finitedifference, finite-element and finite-volume techniques. The LB method has its roots in
LCGA methods, which simulate flow using a ‘bottom-up’ approach, in which, fluids are simulated by a swarm of particles moving along discrete directions on a lattice. The streaming and collision of these particles are governed by a set of rules such that the time-averaged motion of the particles is consistent with the Navier-Stokes equations. The LCGA method is subject to numerical instability, to avoid instability the LB method replaces these particles with averages referred to as particle distribution functions. For a review of LB methods, see Succi 2001 and
Sukop and Thorne 2006. A detailed description of the LB model is given in chapter 3.

The lattice of LB methods is well suited for direct pore space lattice construction from CMT images. Each voxel is directly translated to a single node in the Boltzmann lattice. The pore voxels are represented by nodes on which ‘particles’, represented by particle distribution functions, can freely move upon, the solid voxels are represented by non-interacting nodes, and the solid boundary voxels are represented by bounce-back nodes. Bounce-back nodes are essentially algorithm devices that return any exiting ‘particles’ back into the pore lattice domain with the opposite momentum. Here, we use the half-way bounce-back condition on the nodes at the solid surface, this ensures a microscopic, and hence, a macroscopic no-slip condition

28

halfway between the pore nodes and the bounce-back nodes (He et al., 1997). A density is associated with each node, this gives the system mass, and as a result pressure. The LB simulation presented here simulates a quasi-incompressible fluid, wherein the lattice node density ( ) is directly proportional to the lattice pressure (

), following from the ideal

gas equation of state. The ratio of the determining forces involved in flow, namely the inertial and viscous forces, are conventionally defined by Reynolds numbers. The physical and lattice Reynolds numbers (

) are defined as,

and,

where simulated, is the reference velocity, taken as the average velocity of the flow to be is the reference length taken as the length of the simulated volume,

physical dynamic viscosity,

is the maximum lattice velocity,

is the lattice kinematic viscosity and

is the

is the lattice resolution,

is the lattice density. Here, a water flow with a

is simulated, a stokes flow. The lattice kinematic viscosity was set to 1/6
(generally accepted to be the most numerically stable lattice kinematic viscosity), and the resolution was set equal to 520, corresponding to CMT image resolution. To impose a flow on the lattice the pressure at the inlet and outlet are set to constant values. To set the pressure at the inlet and outlet, the densities of the nodes at these boundaries are maintained throughout the simulation. Since the density is related to the pressure by the equation of state, this ensures the pressures will remain constant at these boundaries (Zou and He 1997). All other boundaries are treated as periodic, as in, when a particle exits a boundary it directly enters the boundary exactly

29

opposite of where it exited. The ‘particles’, driven by the pressure (lattice density) gradient between the inlet and outlet, undergo a single stream and collision step during each iteration of the LB simulation until steady state flow is reached. Each lattice site has an associated macroscopic velocity vector with velocity components reported in the three Cartesian directions (

, where

). Steady-state flow is reached when the mean

lattice velocity magnitude converges.

The open source LB code Palabos v0.7r3 was used to perform the LB simulations. Simulations were carried out on Penn State’s HPC lion-X cluster. The LB simulations used a D3Q19 lattice with Bhatnagar-Gross-Krook (BGK) dynamics. Pan et al. (2006) investigated the accuracy of LB simulations of flow in sphere packs using multiple-relaxation-time (MRT), single-relaxation-time (SRT) and BGK dynamics. They found for a lattice viscosity of 1/6
(used here) the BGK model was as accurate as the more computationally expensive MRT model. To validate the simulation, a Poiseuille flow was simulated through a narrow cylindrical channel with a radius ( ) of 0.33 mm. The analytical solution to the velocity profile for Poiseuille flow is,

The velocity profiles for the analytical solution and LB simulation results for Poiseuille flow are shown in figure 2.6.

30

0.8
0.7
0.6 u (m/s)

0.5
0.4
0.3

LB sim
0.2

Analytical

0.1

0.0
-0.33 -0.27 -0.21 -0.15 -0.09 -0.03 0.03 0.09 0.15 0.21 0.27 0.33 x (m)

Figure 2.6.

Poiseuille flow through a channel with a radius a = 0.33(mm) velocity profile, comparison of LB results and analytical solution.

We observe that the LB simulation is closely

approximating flow in the small channel.

2.4 RESULTS
2.4.1 Fracture and Matrix Pore-Space Geometry

A summary of all measured properties of the pore space are summarised in Table 2.1. The matrix porosity (φm) was found to be ~0.32, significantly lower than the values of 0.40–0.50 reported by Prodanovic et al. 2006, which suggests inconsistency in the manufacturing process of these rods. The fracture porosity (φf) of the pore space is 0.243; the volume of the fracture is significant, and the fracture is expected to dominate flow. The roughness of the fracture was measured for 16 slices evenly distributed along the vertical extent of the fracture
(every 0.652 mm). The Hurst exponents were determined for each of these slices; figure 2.7 shows a log-log plot of deviation vs. window size for some of these slices. The Hurst exponent ranges from 0.51 to 0.69, with a mean Hurst exponent (ε) of 0.612; this is a nominal value of

31

roughness for permeable media. The mean Hurst exponent of 0.612 is similar to Hurst exponents measured for tensile fractures in Berea sandstone (0.38 to 0.45) (Crandall et al. 2010), and single hydraulically induced fractures in annular cores of Fountainbleu sandstone (0.4 to 0.5) (Renard et al. 2009). The mean tortuosity was determined from the tortuosity measurement of 50 streamlines, and found to be 0.52. This is a relatively high tortuosity, the tortuosities measured for a single tensile fractures in Berea sandstone ranged from 0.29 to 0.30 (Nazridoust et al.
2006) and 0.034 to 0.064 (Crandall et al. 2010). The differences in tortuosity measured by
Nazridoust et al. 2006 and Crandall et al. 2010 are a result of different methods for determining actual length of streamlines. Crandall et al. 2010 determined the actual length from seeded streamlines in the same manner used here, Nazridoust et al. 2006 assumed the fracture profile was representative of the tortuosity.

2.5
Slice 200

log (σfp)

2
1.5
1
0.5
0
0

0.5

1

1.5

2

2.5

log (s)
Figure 2.7.

Determination of Hurst coefficient using the variable bandwidth method, the data points are shown for four slices, the linear function shown (ε = 0.612) was determined from 16 equidistant slices spanning the range of the simulated pore volume. The Hurst exponent is the slope of the fitted linear function shown.

32

Table 2.1

Pore Space Geometry and Fracture Characteristics

̅

2.4.2 LB Simulation Results

The objective of running the three LB simulations is to study flow in a fractured permeable medium, a fractured impermeable medium and an unfractured permeable medium, three LB simulation results are presented respectively, one of the entire pore space (fractured permeable media), one of the fracture alone (fractured impermeable media), and an average of five simulations of the matrix alone (permeable media). The velocities are normalized to the mean xvelocity (̅̅̅), also the Darcy velocity, for each voxel with index ,

|
̅̅̅

|
̅̅̅

33

|

|
̅̅̅

A three-dimensional velocity glyph visualization is presented in figure 2.8, this visualization illustrates the channeling of flow in the fractures for both simulations. In general, the flow through the fracture is significantly increased by the permeable matrix. Also, the flow in the fracture-alone simulation shows a single readily apparent high velocity choke point, while the flow in the fracture – aided by the permeable matrix – is more widely dispersed over the fracture.
The matrix only LB simulations were run on 200 x 200 x 400 voxel portions of the pore space that did not contain any part of the fracture; the average of five of these simulations was used to determine matrix flow properties. The permeability was determined from the mean x-, or Darcy velocities, ̅̅̅

The flow in the fracture for the all-pore space simulation dominates the flow, with 88% of the flow occurring in the fracture with a velocity 4.24 times the mean x-velocity. The permeability of the fracture with a permeable matrix is 15.1 times greater than the fracture alone. The permeability of the matrix in the all-pore space simulation is also increased by the proximity of the fracture by 2.14 times. However, if the simulation included a greater extent of the matrix – as in a greater distance between fractures – this effect would dissipate. As in, if a greater extent of the matrix was given this apparent increase in permeability would drop to unity. Keep in mind the periodic boundaries imply that the flow simulated is the equivalent of imposing a distance between ‘fractures’ equal to the dimension of the simulated volume orthogonal to flow.

34

Figure 2.8.

Normalized x-velocity glyph visualization of the fracture-alone (above) and the all-pore
(below) simulations; the glyphs are larger and farther to the red side of the spectrum for greater velocities. The maximum normalized velocity shown in the glyphs is

.

To explore the flow adjacent to the fracture, the velocities are plotted as a function of distance

35

from the fracture in figure 2.9. The velocities appear to decrease logarithmically as a function of distance from the fracture following the form,

( )

where

( )

is the velocity in the x-, y- or z-direction and

is the distance from the fracture. This

is a purely empirical fitting of data. The high velocities in the matrix immediately adjacent to the fracture persist for only a short distance, a distance on the order of a single grain diameter.
The previously mentioned apparent increase in the matrix permeability in the all-pore simulation is a result of these high velocities.

36

Figure 2.9.

Mean normalized velocities (

) as a function of distance from the fracture ( ),

and logarithmic function fits.

2.4.3 Fracture Permeability Estimation

A summary of the LB simulation results and fracture permeability estimates from previously mentioned investigations are presented in Table 2.2. The fracture permeability determined from
LB simulation results for the all-pore simulation and the fracture-alone simulation are 31.9e-9 and 2.11e-10 m 2 , respectively. The ACL determines an estimate of fracture permeability from the mean aperture, ignoring all other aspects of the fracture geometry and matrix permeability.

37

The ACL overestimates the fracture permeability in the all-pore simulation and the fracturealone simulation by 11.1 and 168 times, respectively. The ACL is not expected to give a good estimate of permeability for our fracture. The ACL estimate decreases in accuracy with the greater influence of unaccounted geometry determining flow in the fracture, such as the aperture distribution, tortuosity and fracture roughness. If the fracture wall separation was significantly increased, the tortuosity would decrease and the roughness would interfere with a smaller percentage of the flow. Witherspoon et al. 1980 fit an empirical prediction of fracture permeability to measurements from flow experiments. They determined the fracture permeability from the mean aperture and a surface roughness factor, here set to the highest value
(

) of the range of values suggested. Witherspoon et al. 1980 overestimates the

fracture permeability in the all-pore simulation and the fracture-alone simulation by 6.71 and
102 times, respectively. Witherspoon et al. 1980 was determined empirically from flow experiments, thus it is subject to the nature of the fractures used to fit it. If these fractures were not as rough and tortuous as ours we would expect these characteristics to have less of an effect on their flow experiments, and more closely approximated by the ACL. Zhang et al.
1996 determines the fracture permeability from the mean aperture and a power factor, here set to a value (

) linearly interpolated from the

Hurst exponent of 0.80 (

) and 0.30 (

values determined for fractures with a
). This estimate, like the estimate of

Witherspoon et al. 1980, was determined empirically, and unlike Witherspoon et al. 1980, from LCGA simulations on synthetic fractures, instead of flow measurements. Zhang et al. 1996 underestimates the fracture permeability in the all-pore simulation by 3.85 times, and overestimates fracture permeability in the fracture-alone simulation by 3.94 times. The simulated flow in Zhang et al. 1996 took place in a fracture with impermeable walls, thus their

38

estimate of fracture permeability is most appropriately tested by the fracture-alone LB simulation, and the overestimate is a result of the differences in the fractures simulated by them and the fracture imaged here. The fractures simulated by Zhang et al. 1996 had no induced tortuosity, nor varying aperture. Tortuosity and varying aperture within a fracture will increase the resistance to flow, without consideration of these fracture characteristics we would expect an overestimation of permeability. However, considering the simplicity of their proposed estimation of fracture permeability it is a vast improvement over ACL and Witherspoon et al. 1980. Eker and Akin 2006 proposed the β value of the Zhang et al. 1996 is a linear function of the fractal dimension, as we see in table 2.2 their estimate is by far the most inaccurate. This is not surprising, considering this linear function was fit to synthesised (

fractures

in

impermeable

rocks

). Their linear function for

nor ours, suggesting a linear function of

numerical

flow

simulations

on

with abnormally high fractal dimensions

does not fit the results of Zhang et al. 1996,

determined from the fractal dimension alone will

not accurately estimate fracture permeability outside from those for which it was fit.

The fracture permeability estimate of Zimmerman and Bodvarsson 1996, unlike the above estimations of fracture permeability, is semi-analytically derived. Their estimate is only applicable to fractures in impermeable rock, and underestimates the fracture permeability for the fracture-alone LB simulation by approximately 48%. The underestimation indicates that characterizing the fracture by the mean and standard deviation of the fracture aperture alone does not provide enough information to accurately predict the fracture permeability.

Nazridoust et al. 2006 derived a semi-analytical prediction for fracture permeability that assumes

39

the flow in the fracture can be treated in much the same manner as flow in pipes.
Conventionally in pipe and duct flow, pressure losses are a function of friction factors, which are dependent on flow rate. The cubic law, which was the basis for previously evaluated fracture permeability estimations, assumes Stoke’s or laminar flow. The smooth, parallel plate solution used by Nazridoust et al. 2006 includes inertial effects, thus the resistance to flow is dependent on flow rate. The fracture permeability estimation of Nazridoust et al. 2006 is dependent on the flow rate, fracture aperture mean, standard deviation and tortuosity. They fitted their function for the friction factor to CFD simulations performed on two-dimensional
CMT images of fractures in Berea sandstone. These two-dimensional fracture profiles contained no asperities. Their derivation does not account for permeable fracture walls, thus only the fracture-alone LB simulation is comparable. Nazridoust et al. 2006 overestimate the permeability in the fracture-alone LB simulation by 1.33 times, but it is the best estimate of those tested. This overestimation can be attributed to the significant amount of asperities in our fracture, which are not accounted for by Nazridoust et al. 2006, and in particular the effect these asperities have on a three-dimensional fracture.

40

Table 2.2.

Notes:

Summary of LB Permeability Measurements and Fracture Permeability Estimates

Summary of LB results ( estimates, where

- Alternate Cubic Law,

(X. Zhang et al., 1996),
1996),

) and fracture permeability
- (Witherspoon et al., 1980),

-(Eker and Akin, 2006),

- (Nazridoust et al., 2006) and

-

- (Zimmerman and Bodvarsson,

- (Crandall et al., 2010);

is the only

estimate that considers matrix permeability, and the only estimate appropriate for the allpore LB simulation.

Crandall et al. 2010 built upon the work of Nazridoust et al. 2006 by adding the matrix permeability to the prediction of fracture permeability. They used the same CMT images of a fracture in Berea sandstone as Nazridoust et al. 2006 and performed two-dimensional CFD simulations with a homogenous, isotropic fracture wall with permeability ranging from 2E-16 to
2E-12 mm 2 (0.2 to 2000 mD). Since their derivation accounts for a permeable fracture wall,

41

only the all-pore LB simulation is comparable. Crandall et al. 2010 underestimates the fracture permeability in the all-pore simulation by 6.72 times. Although the prediction for the LB simulation of the fracture alone (Nazridoust et al. 2006) is off by a relatively small percentage
(33%), the prediction for the all-pore space fracture permeability (Crandall et al. 2010) is off by
672%. Here are some explanations for this large difference. First, Crandall et al. 2010 validated their derivation of fracture permeability from CFD results for matrix permeability well below our matrix permeability of 4.92E-11mm 2 (49.7 D). Also, there is a fundamental difference in the nature of the matrix for our simulations and theirs. Our simulation directly simulates an actual matrix pore space, this allows for fluid to channel in and out of the fracture through the matrix, helping to bypass fracture choke points. This will greatly increase the flow in the fracture, as we observe from our results. The matrix in Crandall et al. 2010 does not allow for this mechanism of bypass. In CFD the rock matrix is divided up into a finegrid of finite-volumes, and each of these volumes are homogenous. This form of simulation assumes the finite-volumes are on the order of the representative elementary volumes (REV) of flow for the permeable medium being simulated. The fracture pore space does not have a REV, but the matrix does, and it’s REV is dependent on the grain size of the permeable media. The
REVs of sandstone, and the matrix simulated in the all-pore LB simulation presented here, are many orders of magnitude larger than the finite volumes used by Crandall et al. 2010.
Crandall et al. 2010 use of CFD is arguably a fair representation of flow in fractured permeable media with much smaller REVs.

42

2.5 CONCLUSIONS

A single rough fracture was propagated along the axis of a porous polyethylene rod; this fractured permeable medium was imaged using CMT. The fracture pore space identification process, consisting of a succession of erosions, percolations and expansions on the segmented processed CMT images of the pore space, presented here appears to capture the fracture well.
The fracture was found to have a log-normal aperture distribution and a Hurst exponent of
0.612, similar to tensile fractures in Berea and Fountainbleu sandstone. Three LB simulations were performed on the imaged pore space, one of all-pore space imaged (fracture in permeable media), one of the fracture alone (fracture with impermeable walls), and an average of five simulations on portions of the pore space that did not contain any part of the fracture to be representative of flow in the matrix only. The fracture permeability was increased by a factor of
15.1 by the presence of the highly permeable matrix, and the matrix permeability was increased by a factor of 2.14. However, this increase in the matrix permeability is a result of the presence of the fracture, and is dependent on fracture spacing. The x-, y-, and z-velocities were found to be a logarithmic function of the distance from the fracture. Although developing a generalised form of this equation was beyond the scope of this paper.

Fracture permeability estimates from a variety of previous investigations were compared to our
LB simulation results; all of these estimates assume an impermeable fracture wall, except for
Crandall et al. 2010. The ACL does not account for fracture roughness, tortuosity, asperities and aperture distribution, and has been widely accepted to grossly overestimate the fracture permeability. The ACL overestimated the fracture permeability by two orders of magnitude.

43

Zhang et al. 1996 suggest permeability is a power function of mean aperture and a fitted power factor β. Although this will provide a good estimate of permeability if a value of β is found that fits the fracture flow measurements, it has no power of prediction without a function for β dependent on fracture characteristics. Eker and Akin 2006 proposed a linear function for β dependent on the fractal dimension of the fracture wall roughness; however their prediction for β did not fit the results of Zhang et al. 1996, nor ours. A linear function for β dependent only on the fracture wall fractal dimension is inadequate. Zimmerman and Bodvarsson 1996 derived a semi-analytical prediction of fracture permeability dependent on the mean fracture aperture, fracture aperture deviation and

fractional contact

area from the Reynold’s lubrication

equation. The permeability predicted by this function underestimated fracture permeability approximately 48%, indicating accurate prediction of fracture permeability requires more information than the mean and standard deviation of the fracture aperture provide. Nazridoust et al. 2006 derived a semi-analytical function of fracture permeability dependent on mean fracture aperture, fracture aperture standard deviation and flow rate from the parallel-plate flow solution for flow with significant inertial effects. Their prediction of fracture permeability provided the closest match to LB simulations, overestimating the fracture permeability by a mere
33%. The slight overestimation can be attributed to neglecting the effect of asperities.
Crandall et al. 2010 built upon the work of Nazridoust et al. 2006 to include the effect of matrix permeability on flow in fractured permeable media. Their fracture permeability prediction underestimated the fracture permeability in the LB simulation of the fractured permeable medium by a factor of 6.72 times. As discussed, the simulations used in Crandall et al. 2010 assumed the flow in the matrix can be represented by homogenous and isotropic finite volume elements. However, they used REVs that are much smaller than the REVs of

44

sandstones and the permeable media used in our simulations. Our permeable medium is simulated at the pore-scale without any assumptions made regarding the matrix pore space.
This allows for bypassing of fracture choke points through the pore spaces of the matrix adjacent to the fracture. This dynamic is not possible in CFD simulations. It is for these reasons flow in the fracture was so much greater in our simulations.

LB methods, and other pore-scale models are tools that give us insight into fluid flow through permeable media at scales that are very difficult to observe; particularly in the study of fractured permeable media, for which measuring the flow in the fracture alone would prove very difficult. It is for this reason validating our results would require more indirect methods, such as measuring the mean permeability of the entire matrix-fracture system and comparing this permeability to simulation results. However, this investigation was not meant to be an evaluation of the accuracy of LB methods, the evaluation of the accuracy of single-phase LB models has been addressed by many of the previous researchers mentioned. A comprehensive prediction of fracture permeability from fracture and matrix characteristics for fractures in permeable rock that includes all determining factors, including the flow behavior of the matrix, remains elusive. However, pore-scale studies such as this one provide a tool to help minimize unrealistic assumptions made in more simplistic representations of flow in fractured permeable media.

45

2.6 REFERENCES

Auradou, H., Drazer, G., Hulin, J.P. and Koplik, J. (2005) ‘Permeability anisotropy induced by the shear displacement of rough fracture walls’, Water Resources Research, Vol. 41, No. 9,
W09423, DOI: 10.1029/2005WR003938.

Auzerais, F.M., Dunsmuir, J., Ferreol, B.B., Martys, N., Olson, J., Ramakrishnan, T.S.,
Rothman, D.H. and Schwartz, L.M. (1996) ‘Transport in sandstone: a study based on three dimensional microtomography’, Geophysical Research Letters, Vol. 23, No. 7, pp.705–708.

Boek, E. S. and Venturoli, M. (2010) ‘Lattice-Boltzmann studies of fluid flow in porous media with realistic rock geometries’, Computers and Mathematics with Applications, Vol. 59, No.
7, pp.2305–2314.

Boutt, D.F., Grasselli, G., Fredrich, J.T., Cook, B.K. and Williams, J.R. (2006) ‘Trapping zones: the effect of fracture roughness on the directional anisotropy of fluid flow and colloid transport in a single fracture’, Geophysical Research Letters, Vol. 33, L21402, DOI: 10.1029/
2006GL027659.

Brown, S.R. and Scholz, C.H. (1985) ‘Broad bandwidth study of the topography of natural rock surfaces’, Journal Geophysical Research, Vol. 90, No. B14, pp.12575–12582.

46

Brusseau, M.L., Peng, S., Schnaar, G. and Costanza-Robinson, M.S. (2006) ‘Relationships among air-water interfacial area, capillary pressure, and water saturation for a sandy porous medium’, Water Resources Research, Vol. 42, W03501, DOI: 10.1029/2005WR004058.

Crandall, D., Bromhal, G. and Karpyn, Z.T. (2010) ‘Numerical simulations examining the relationship between wall-roughness and fluid flow in rock fractures’, International Journal of
Rock Mechanics and Mining Sciences, Vol. 47, No. 5, pp.784–796.

Culligan, K.A., Wildenschild, D., Christensen, B.S.B., Gray, W.G. and Rivers, M.L. (2006)
‘Pore-scale characteristics of multiphase flow in porous media: a comparison of air-water and oil-water experiments’, Advances in Water Resources, Vol. 29, No. 2, pp.227–238.

Dougan, L.T., Addison, P.S. and McKenzie, W.M.C. (2000) ‘Fractal analysis of fracture: a comparison of dimension estimates’, Mechanics Research Communications, Vol. 27, No. 4, pp.383–392. Drazer, G. and Koplik, J. (2002) ‘Transport in rough self-affine fractures’, Physical Review
E, Vol. 66, 026303, DOI: 10.1103/PhysRevE.66.026303.

Eker, E. and Akin, S. (2006) ‘Lattice Boltzmann simulation of fluid flow in synthetic fractures’,Transport in Porous Media, Vol. 65, No. 3, pp.363–384.

47

Ferreol, B. and Rothman, D.H. (1995) ‘Lattice-Boltzmann simulations of flow-through
Fontainebleau sandstone’, Transport in Porous Media, Vol. 20, No. 1, pp.3–20.

Frisch, U., Hasslacher, B. and Pomeau, Y. (1986) ‘Lattice-gas automata for the Navierstokes equation’, Physical Review Letters, Vol. 56, No. 14, pp.1505–1508.

Genabeek, O.V. and Rothman, D.H. (1996) ‘Macroscopic manifestations of microscopic flows: phenomenology from simulation’, Annual Review of Earth and Planetary Sciences, Vol.
24, No. 1, pp.63–87.

He, X., Zou, Q., Luo, L. and Dembo, M. (1997) ‘Analytic solutions of simple flows and analysis of nonslip boundary conditions for the lattice Boltzmann BGK model’, Journal of Statistical
Physics, Vol. 87, Nos. 1/2, pp.115–136.

Heijs, A.W.J., Delange, J., Schoute, J.F.T. and Bouma, J. (1995) ‘Computed-tomography as a tool for nondestructive analysis of flow patterns in macroporous clay soils’, Geoderma, Vol.
64, Nos. 3/4, pp.183–196.

Johns, R.A., Steude, J.S., Castanier, L.M. and Roberts, P.V. (1993) ‘Nondestructive measurements of fracture aperture in crystalline rock cores using X-ray computedtomography’, Journal of Geophysical Research-Solid Earth, Vol. 98, No. B2, pp.1889–1900.

48

Karpyn, Z.T., Grader, A.S. and Halleck, P.M. (2007) ‘Visualization of fluid occupancy in a rough fracture using micro-tomography’, Journal of Colloid and Interface Science, Vol. 307,
No. 1, pp.181–187.

Keller, A.A., Roberts, P.V. and Blunt, M.J. (1999) ‘Effect of fracture aperture variations on the dispersion of contaminants’, Water Resources Research, Vol. 35, No. 1, pp.55–63.

Keller, A.A., Roberts, P.V. and Kitanidis, P.K. (1995) ‘Prediction of single-phase transport parameters in a variable aperture fracture’, Geophysical Research Letters, Vol. 22, No. 11, pp.1425–1428. Ketcham, R.A., Slottke, D.T. and Sharp, J.M. (2010) ‘Three-dimensional measurement of fractures in heterogeneous materials using high-resolution X-ray computed tomography’,
Geosphere, Vol. 6, No. 5, pp.499–514.

Kim, I., Lindquist, W.B. and Durham, W.B. (2003) ‘Fracture flow simulation using a finitedifference lattice Boltzmann method’, Physical Review E, Vol. 67, No. 4, 046708, DOI:
10.1103/PhysRevE.67.046708.

Konzuk, J.S. and Kueper, B.H. (2004) ‘Evaluation of cubic law based models describing single-phase flow through a rough-walled fracture’, Water Resources Research, Vol. 40, No.
2, W02402, DOI: 10.1029/2003WR002356.

49

Koponen, A., Kataja, M. and Timonen, J. (1997) ‘Permeability and effective porosity of porous media’, Physical Review E, Vol. 56, No. 3, pp.3319–3325.

Madadi, M. and Sahimi, M. (2003) ‘Lattice Boltzmann simulation of fluid flow in fracture networks with rough, self-affine surfaces’, Physical Review E, Vol. 67, No. 2, 026309, DOI:
10.1103/PhysRevE.67.026309.

Maier, R.S., Kroll, D.M., Kutsovsky, Y.E., Davis, H.T. and Bernard, R.S. (1998) ‘Simulation of flow through bead packs using the lattice Boltzmann method’, Physics of Fluids, Vol. 10,
No.1, pp.60–74.

Manz, B., Gladden, L.F. and Warren, P.B. (1999) ‘Flow and dispersion in porous media: lattice-Boltzmann and NMR studies’, AICHE Journal, Vol. 45, No. 9, pp.1845–1854.

McNamara, G.R. and Zanetti, G. (1988) ‘Use of the Boltzmann equation to simulate latticegas automata’, Physical Review Letters, Vol. 61, No. 20, pp.2332–2335.

Mees, F., Swennen, R., Van Geet, M. and Jacobs, P. (2003) Applications of X-ray
Computed Tomography in Geosciences, p.57, Geological Society Special Publications, London.

Nakashima, Y., Nakano, T., Nakamura, K., Uesugi, K., Tsuchiyama, A. and Ikeda, S.
(2004) ‘Three-dimensional diffusion of non-sorbing species in porous sandstone: computer

50

simulation based on X-ray microtomography using synchrotron’, Journal of Contaminant
Hydrology, Vol. 74, Nos. 1/4, pp.253–264.

Nazridoust, K., Ahmadi, G. and Smith, D.H. (2006) ‘A new friction factor correlation for laminar, single-phase flows through rock fractures’, Journal of Hydrology, Vol. 329, Nos. 1/2, pp.315–328. Nolte, D.D., Pyraknolte, L.J. and Cook, N.G.W. (1989) ‘The fractal geometry of flow paths in natural fractures in rock and the approach to percolation’, Pure and Applied Geophysics,
Vol. 131, Nos. 1/2, pp.111–138.

Okabe, H. and Blunt, M.J. (2004) ‘Prediction of permeability for porous media reconstructed using multiple-point statistics’, Physical Review E, Vol. 70, No. 6, 066135,
DOI: 10.1103/ PhysRevE.70.066135.

Pan, C., Luo, L-S. and Miller, C.T. (2006) ‘An analysis of lattice Boltzmann schemes for porous medium flow simulation’, Computers and Fluids, Vol. 35, Nos. 8–9, pp.898–909.

Pan, C.X., Prins, J.F. and Miller, C.T. (2004) ‘A high-performance lattice Boltzmann implementation to model flow in porous media’, Computer Physics Communications, Vol.
158, No. 2, pp.89–105.

51

Perret, J., Prasher, S.O., Kantzas, A. and Langford, C. (1999) ‘Three-dimensional quantification of macropore networks in undisturbed soil cores’, Soil Science Society of
America Journal, Vol. 63, No. 6, pp.1530–1543.

Poon, C.Y., Sayles, R.S. and Jones, T.A. (1992) ‘Surface measurement and fractal characterization of naturally fractured rocks’, Journal of Physics D, Applied Physics, Vol. 25,
No. 8, pp.1269–1275.

Prodanovic, M., Lindquist, W.B. and Seright, R.S. (2006) ‘Porous structure and fluid partitioning in polyethylene cores from 3D X-ray microtomographic imaging’, Journal of Colloid and
Interface Science, Vol. 298, No. 1, pp.282–297.

Renard, F., Bernard, D., Desrues, J. and Ougier-Simonin, A. (2009) ‘3D imaging of fracture propagation using synchrotron X-ray microtomography’, Earth and Planetary Science
Letters, Vol. 286, Nos. 1/2, pp.285–291.

Schmittbuhl, J., Schmitt, F. and Scholz, C. (1995) ‘Scaling invariance of crack surfaces’, Journal of Geophysical Research, Vol. 100, No. B4, pp.5953–5973.

Stockman, H.W., Li, C. and Wilson, J.L. (1997) ‘A lattice-gas and lattice Boltzmann study of mixing at continuous fracture junctions: importance of boundary conditions’, Geophysical
Research Letters, Vol. 24, No. 12, pp.1515–1518.

52

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

Sukop, M. and Thorne, D.T. (2006) Lattice Boltzmann Modelingan Introduction for
Geoscientists and Engineers, Springer-Verlag, Berlin, Heidelberg.

Tsang, Y.W. and Witherspoon, P.A. (1981) ‘Hydromechanical behavior of a deformable rock fracture subject to normal stress’, Journal of Geophysical Research, Vol. 86, No. 10, pp.9287–9298. Tsang, Y.W. and Witherspoon, P.A. (1983) ‘The dependence of fracture mechanical and fluidflow properties on fracture roughness and sample-size’, Journal of Geophysical Research, Vol.
88, No. 3, pp.2359–2366.

Wennberg, O.P., Rennan, L. and Basquet, R. (2009) ‘Computed tomography scan imaging of natural open fractures in a porous rock; geometry and fluid flow’, Geophysical Prospecting,
Vol. 57, No. 2, pp.239–249.

Wildenschild, D., Hopmans, J.W., Rivers, M.L. and Kent, A.J.R. (2005) ‘Quantitative analysis of flow processes in a sand using synchrotron-based X-ray microtomography’, Vadose Zone
Journal, Vol. 4, No. 1, pp.112–126.

53

Witherspoon, P.A., Wang, J.S.Y., Iwai, K. and Gale, J.E. (1980) ‘Validity of cubic law for fluidflow in a deformable rock fracture’, Water Resources Research, Vol. 16, No. 6, pp.1016–1024.

Zhang, D.X. and Kang, Q.J. (2004) ‘Pore scale simulation of solute transport in fractured porous media’, Geophysical Research Letters, Vol. 31, No. 12, L12504, DOI: 10.1029/
2004GL019886.

Zhang, X., Knackstedt, M.A. and Sahimi, M. (1996) ‘Fluid flow across mass fractals and selfaffine surfaces’, Physica A: Statistical Mechanics and its Applications, Vol. 233, Nos. 3/4, pp.835–847. Zimmerman, R.W. and

Bodvarsson, G.S. (1996) ‘Hydraulic

conductivity of rock

fractures’, Transport in Porous Media, Vol. 23, No. 1, pp.1–30.

Zou, Q. and He, X. (1997) ‘On pressure and velocity boundary conditions for the lattice
Boltzmann BGK model’, Physics of Fluids, Vol. 9, No. 6, pp.1591–1598.

54

CHAPTER 3: RELATIVE PERMEABILITY OF HOMOGENOUS-WET AND MIXEDWET

POROUS

MEDIA

AS

DETERMINED

BY

PORE-SCALE

LATTICE

BOLTZMANN MODELING

3.1 SUMMARY

We present a pore-scale study of the dependence of relative permeability dependence on the strength of wettability of homogenous-wet porous media, as well as the dependence of relative permeability on the distribution and severity of wettability alteration of porous media altered to a mixed-wet state. A Shan-Chen type multicomponent multiphase lattice Boltzmann model is employed to determine pore-scale fluid distributions and relative permeability. Mixed-wet states are created – after pre-simulation of homogeneous-wet porous medium – by altering the wettability of solid surfaces in contact with the non-wetting phase.

To ensure accurate

representation of fluid-solid interfacial areas we compare LB simulation results to experimental measurements of interfacial fluid-fluid and fluid-solid areas determined by x-ray computed microtomography imaging of water and oil distributions in bead packs (Landry et al. 2011). The
LB simulations are found to match experimental trends observed for fluid-fluid and fluid-solid interfacial area-saturation relationships.

The relative permeability of both fluids in the

homogenous-wet porous media is found to decrease with a decreasing contact angle. This is attributed to the increasing disconnection of the non-wetting phase and increased fluid-solid interfacial area of the wetting phase. The relative permeability of both fluids in the altered mixed-wet porous media is found to decrease. However the significance of the decrease is

55

dependent on the connectivity of the unaltered solid surfaces, with less dependence on the severity of alteration.

3.2 BACKGROUND

Macroscale multiphase flow properties in porous media (i.e. relative permeability and capillary pressure) are strongly dependent on wettability. The wettability of a porous medium is initially determined by its mineral composition. However, brines and oils can alter the wettability of porous media, such that the multiphase flow properties of a rock can be significantly dependent on saturation history. Subsequently, the majority of consolidated and unconsolidated porous media have a mixed or fractional wettability. The process occurs in three steps:
1) Rocks initially develop wettability as a result of mineral-water interactions.
2) Migration of non-aqueous phases into the rock alters the wettability of the rock.
3) Subsequent cycles of displacement further alter the wettability of the rock.
Wettability is conventionally measured using the Amott-Harvey method or the USBM method.
Both measure an index of wettability from displacement experiments. Indices range from -1-to-1 for perfectly oil-wet to perfectly water-wet conditions. At the pore-scale, wettability is described by contact angles. The contact angle is measured through the aqueous phase and ranges from 0° to 180°, for perfectly water-wet to perfectly oil-wet, respectively. Common reservoir porous media includes siliclastic rocks (i.e. sandstone) primarily composed of quartz, but also containing significant amounts of feldspars, calcite, dolomite and clays; and carbonate rocks (i.e. limestone) composed primarily of calcite, dolomite, clays and organics.

56

Conventional measurement of oil recovery and macroscale flow properties from mixed-wet reservoir rock cores demonstrate some of the important trends observed in field-scale variations of saturation history and wettability. Janhunandan and Morrow 1995 investigated changes in wettability and oil recovery from Berea cores aged with two different crude oils and varying initial water saturations (Swi), aging temperatures and brine compositions. Regardless of oil type and brine composition the Amott index of the rock increased monotonically with S wi. Jerauld and Rathmell 1997 observed the same trend in oil-wetness and connate water saturation for
Prudhoe Bay sandstone cores composed primarily of quartz. Core wettability trended more oilwet with the decrease in connate water saturation above the oil-water contact line of the reservoir. With more oil present in the rock, a greater portion of the pore space will be in contact with the oil, resulting in larger surface area of the pore space subject to wettability alteration.
Both studies also demonstrated the dependence of wettability alteration on the brine composition, oil composition and aging temperature.

Schembre et al. 2006 measured the

wettability and oil recovery of cores from a diamtomite (carbonate) reservoir. The Amott indices were found to increase significantly from 0.2-0.55 to 0.55-0.85 with temperature increases from
45 to 230 C, and correlate with the percent content of water-wet clays. The increase in the
Amott indices were also found to increase oil recovery. The temperature dependence indicates variation in the temperature of a reservoir will affect the wettability of the rock.

Once a fluid is in contact with a mineral surface in the porous media the alteration in wettability is controlled by fluid-mineral surface chemistry. As summarized by Buckley and Liu 1998 the mechanisms involved in crude oil/brine/solid interaction, include polar interactions between the polar functional groups in the oil and polar surface sites on the mineral, surface precipitation,

57

acid/base reactions and ion binding. One of the most common minerals found in siliclastic rocks is quartz. The surface interactions between quartz and brine depend on both the pH and ionic strength of the brine.

The affinity of brine to quartz in the presence of an oil phase is

proportional to the amount of ion adsorption. Ion adsorption increases with increasing ionic strength and decreases for pHs in the proximity of the pH of zero-point surface charge (pHzpc =
1.5-4.0) (Barranco et al. 1997). The same behavior is observed in coal, with the strongest water wettability found at alkaline and acidic conditions. Chaturvedi et al. 2009 demonstrated the effect pH has on the wettability, and subsequently, the relative permeability of coal. At brine pHs of 10, and 7, the intersection point of relative permeability curves were S w = 0.65 (waterwet) and Sw = 0.43 (weakly oil-wet), respectively. Carbonate wettability is also sensitive to brine pH and composition, as demonstrated in Yu et al. 2009 chalk cores experienced very little spontaneous imbibition of reservoir brine, while seawater altered the wettability to water-wet and spontaneously imbibed the core. Evje and Hiorth 2011 modeled these results with the inclusion water-rock chemistry, specifically the dissolution and precipitation of calcite, magnesite and anhydrite in the presence of Mg2+ and SO42-. They attributed the increase in wettability to the presence of SO42- in the seawater brine, which increased the dissolution of calcite, exposing new water-wet surfaces. One of the most important mechanisms of wettability alteration by crude oils is the adsorption and deposition of asphaltenes on mineral surfaces. Asphaltenes can render initially water-wet mineral surfaces oil-wet.

Saraji et al. 2010 measured the adsorption of

asphaltene to calcite, dolomite and quartz by asphaltene containing oils under flow conditions.
The asphaltene was found to adsorb the most to calcite, and to lesser extent quartz and dolomite.
Wettability was not measured, but it was concluded that the amount of asphaltene adsorbed was sufficient to form a monolayer on all minerals. The alteration of mineral wettability is controlled

58

by mineral-fluid interactions. Knowledge of individual mineral-fluid surface chemistry provides insight into which minerals in a porous medium are being altered.

Tabriziky et al. 2011

measured changes to contact angles and wettability indices of quartz, kaolinite and calcite after aging in presence of oil with and without asphaltene content, and brines containing either MgCl 2 or Na2SO4. Wettability was measured using BET H2O surface adsorption, and contact angles measured on flat surfaces of calcite and quartz to investigate the dependence of wettability alteration on ionic species and the presence of asphaltene. Quartz was neutral wet for oil without asphaltenes and distilled water, and weakly oil wet for oil with asphaltenes and distilled water.
SO42- had a negligible effect; however, Mg- significantly altered quartz to a water wet state in the presence of asphaltenes, and an oil-wet state without the presence of asphaltenes. Kaolinite was neutral wet for oil with and without asphaltenes and distilled water. SO 42-, as with quartz, had a negligible effect. Mg - significantly altered kaolinite to a water wet state in the presence of oil with and without asphaltene. Calcite was water wet in presence of oil with and without asphaltene and distilled water. SO42- altered the wettability to a weakly water-wet state, and Mg altered the wettability further to a weakly oil-wet state. Studies such as these demonstrate the complexity involved in fluid/mineral interactions and the resulting alteration of wettability.

There are two factors to consider in the study of wettability alteration, the chemistry of alteration and the distribution of fluids in the pore geometry with respect to mineral constituents of the porous media. There are few studies that have attempted to measure contact angles at the pore scale within natural porous media, where mineral heterogeneity exists in individual pores. Robin et al. 1995 imaged the distribution of oil and water phases in a reservoir sandstone and carbonate before and after wettability alteration using cryo-SEM.

59

The images provided qualitative

evidence of individual mineral wettability alteration within individual pores. Many of the trends in mineral wettability alteration discussed above are found in these images. Although these types of images confirm the assumptions of bulk studies, they do not provide much insight into the 3D pore geometry that governs the movement of fluids. Non-destructive 3D imaging of pore geometry and fluid distributions using x-ray computed microtomography (CMT) provides us with a powerful tool to image and analyze pore geometries and fluid distributions in these geometries. Many investigators have used CMT to image a variety of fluid distributions in unconsolidated materials (Culligan et al. 2006, Brusseau et al. 2006, Brusseau et al. 2008,
Constanza-Robinson et al. 2008, Al-Raoush et al. 2009, Lebedeva and Fogden 2011) and rock
(Coles et al. 1998, Turner et al. 2004, Prodanovic et al. 2007, Iglauer et al. 2011,, Silin et al.
2011, Kumar et al. 2012). Armed with knowledge of the fluid/mineral wettability alteration and pore geometry, pore scale models can be employed to investigate the macroscale multiphase flow properties of mixed-wet porous media.

The focus of this paper is to further elucidate the relationship between mixed-wet states and relative permeability. In the past two decades pore network models have been successfully used to model wettability alteration, multiphase fluid flow of mixed-wet states, and saturation history dependent relative permeability and capillary pressure hysteresis (Blunt 1997a, Blunt 1998a,
Blunt et al. 2002, Dixit 1996, Dixit 1998, Al-Futaisi and Patzek 2004, Hui and Blunt 2000,
Jackson et al. 2003). Kovscek et al. 1993 introduced a physically based model for wettability alteration within individual pores; each pore is represented by a capillary tube with a star-shape cross-section. The tube is initially occupied by the wetting phase, during drainage the wetting phase is displaced by a piston-like mechanism by the non-wetting phase. The non-wetting phase

60

will now occupy the center of the pore, leaving the corners occupied by the wetting phase. For the surface of the pore to make contact with the non-wetting phase the capillary pressure must overcome a critical capillary pressure at which the wetting film destabilizes to molecular thickness. Once contact is made the surface of the pore in contact with the oil is altered in its strength of wettability.

This concept was further extended to squares (Blunt 1997a, Blunt

1997b), triangles and lenses (Oren et al. 1998). As is noted in Blunt 1997, the disjoining capillary pressure is not known for each individual pore, pore network models cope with this by assigning a wettability alteration to a fraction of the pores occupied by oil. These fractions can then be varied to observe resulting multiphase flow behavior. The extraction of representative pore networks from 3D images was introduced by Oren et al. 1998, and has led to numerous studies of multiphase flow in a variety of porous media imaged by CMT and modeled by pore network models.

In mixed-wet media there are three general possible modes of assigning

fractional wettability to pores - the largest pores, the smallest pores, or at random (Oren and
Bakke 2003, Hoiland et al. 2007). The choice of contact angles representing the wettability can also be uniform or distributed at random over a range. With these free parameters almost any wetting state can be modeled (Fenwick and Blunt 1998, Hui and Blunt 2000, Jackson et al. 2003,
Al-Futaisi and Patzek 2004, Speiteri et al. 2008, Svirsky et al. 2007, Zhao et al. 2010). There are concerns regarding the use of pore network models for the study of wettability alteration and its effect on multiphase flows. The simplification of the pore space into capillary tubes of varying shape is not predisposed to the inclusion of knowledge of the mineral composition of the pores. The surface chemistry is highly dependent on accurate representation of fluid/mineral interfacial areas, which may be difficult to ascertain without proper inclusion of the mineral composition of pores.

61

In this paper we use a more rigorous pore scale modeling approach. Lattice Boltzmann (LB) methods can directly use 3D images of pore space by a one-to-one correlation of image voxels to lattice sites. The accurate representation of any pore geometry is then only a question of lattice resolution. This makes the model very appealing over pore network models. However, it is far more computationally expensive, requiring on the order of 10 4-5 more cpu time than an equivalent pore network model (Vogel et al. 2005).

The computational needs of LB are

compensated by its inherent parallelism – in most models each node is only aware of its neighboring nodes.

Multiphase LB models have been validated against experimental

measurements of capillary pressure in bead packs (Pan et al. 2004, Schaap et al. 2007, Porter et al. 2009) and sandstones (Ramstad et al. 2010, Ramstad et al. 2012), as well as experimental measurements of relative permeability in sphere packs (Ghassemi and Pak 2009, Hao and Cheng
2010) and sandstones (Boek and Venturoli 2010, Ramstad et al. 2010). Investigations of mixed wettability using LB are limited to the work of Hazlett et al. 1998. In this investigation a CMT image of a reservoir sandstone was translated directly into a lattice for LB multiphase displacement simulations for a strongly water-wet sandstone and a mixed-wet sandstone. The mixed wet state was determined by simulating a drainage to irreducible water saturation. At the end of the drainage, the pore walls in contact with the non-wetting phase are altered in wettability. The measurements of capillary pressure and relative permeability matched

experimentally observed trends. Although there has been a fair amount of validation of LB methods via matching experimental measurements of macroscale flow properties, few have compared fluid distributions imaged using CMT with LB model results. Sukop et al. 2008 compared two phase fluid distributions in a sand pack imaged using CMT to distributions

62

resulting from LB simulations. The saturation per slice was found to match experimental results well, but there was no analysis presented of interfacial area measurements. Porter et al. 2009 compared interfacial fluid-fluid areas measured by CMT to LB displacement simulations and found a good match for drainage simulations.

Our objective is to investigate the evolution of wettability alteration as a result of saturation history and the initial state of wettability, as well as the sensitivity of relative permeability to this evolution. Proper pore-scale simulations of wettability alteration require accurate predictions of fluid-solid interfacial areas. To test the ability of the LB model to predict these interfacial areas we will first compare CMT images of oil/brine fluid distributions in water-wet and weakly oilwet bead packs from Landry et al. 2011 to the results of LB simulations. The bead packs can be viewed as simple analogs of water-wet siliclastic porous media and weakly oil-wet carbonate porous media. There are three variables of interest in our investigation of relative permeability
– saturation relationships and their dependence on wettability, the initial wetting state of the porous media, the saturation at which wettability alteration takes place, and the degree of alteration. As was done in the work of Hazlett et al. 1998 solid surfaces in contact with the nonwetting fluid, after fluid distribution is established by homogenous-wet LB simulations, will be altered to create a mixed-wet porous medium.

Relative permeability measurements on the

initial and mixed-wet states will then be carried out using the LB model. We will be able to evaluate the extent to which the LB model can recreate pore scale distributions of immiscible fluids, and compare relative permeability measurements of mixed-wet states. In this simplified scenario, the initial wettability, saturation of alteration and severity of alteration control the mixed-wet state of the porous media. The only variable in our simulations is wettability, as a

63

result of these three parameters.

Natural porous media is subject to a far more complex

wettability alteration process, dependent on numerous variables related to the mineral-fluid interactions. However, the resulting evolution of relative permeability – saturation relationships can be investigated without such considerations.

3.3 MATERIALS AND METHODS
3.3.1 Experimental Measurements

The validation of pore-scale models at the pore-scale is limited in the literature. Here we compare LB simulation results to CMT images of brine/oil distributions in strongly water-wet
(glass) and weakly oil-wet (polyethylene) bead packs from Landry et al. 2011. Both bead types were spherical with a size range between 0.425 and 0.600 mm, and packed in a column 25.4 mm in diameter and ~80-90 mm in length. The fluids used were a 4-8% NaI brine and kerosene, with dynamic viscosities of 1.00 and 2.43 cP, respectively. The bead packs were initially fully saturated with brine, followed by capillary-dominated displacement by kerosene, then brine. At the end of each displacement CMT images were taken with a voxel resolution of 0.0260 x 0.0260 x 0.0292 mm3 and 0.0259 x 0.0259 x 0.0274 mm3, for the glass and polyethylene bead packs, respectively. Due to the relatively large pore aperture of these porous media gravity plays a significant role in the distribution of fluids in these bead packs.

This results in a classic

transition-zone saturation profile in both bead packs at the end of drainage, from which fluid distributions can be observed at varying fluid saturations. From these images specific interfacial fluid areas were measured, as well as blob size and surface area distributions.

64

These

measurements of volume and area will be compared in this paper to the results of LB simulations. 3.3.2 Single-Phase BGK Lattice Boltzmann Model

The LB method is a mesoscopic method based on microscopic particle dynamics that provides numerical solutions to macroscopic hydrodynamics. In short, fluids are modeled as swarms of streaming and colliding particles. function, the lattice, and

Each node on the lattice contains a particle distribution

, where is the index of each discrete velocity,

is the location of the node on

is time; here we use a D3Q19 lattice (Qian et al 1992). The D3Q19 particle

distribution function has 19 discrete velocities,

, including a zero velocity and 18 velocities

pointing to neighboring nodes. The particle distribution function is also referred to as the density function – the macroscopic density at each node is represented by a distribution of densities of particles moving along each discrete velocity. The particle distribution function evolves in time according to the lattice Boltzmann equation,
[
where

is the relaxation parameter,

lattice units 1 t.u.) and

is a discrete time step (

]
, for one time iteration, in

is the equilibrium particle distribution function. The left hand

side of the equation represents the streaming of particles by passing the particle distribution at each discrete velocity to the respective neighboring node. The right-hand side of the equation represents the collision of particles as a partial relaxation to the equilibrium particle distribution.
The relaxation parameter

(also known as the collision interval) is representative of the rate of

particle collisions, and is related to the kinematic viscosity, , of the lattice by,
65

,

where

is the speed of sound of the lattice, or propagation speed, equal to

√ . The discrete

is the ratio of the lattice spacing – equal to 1 lattice length unit

lattice speed unit

–

and the time step. Here we use Bhatnar-Gross-Krook (BGK) dynamics, meaning the relaxation parameter is defined by a single value. Although it has been shown that BGK dynamics results in viscosity-dependent permeability measurements (Pan et al. 2004), as opposed to more computationally demanding dynamics (e.g. Multiple-Relaxation Time dynamics), for our purposes BGK dynamics will suffice. The D3Q19 equilibrium particle distribution function can be calculated as
[
where

]

is the weight of each discrete velocity,

macroscopic velocity moment of

is the density of

. The discrete velocities

, and

are defined as follows,

[

]
[

with weights density is the

]
⁄

,

⁄

⁄

,

is obtained by summing the particle densities,

∑

, and the macroscopic velocity

is obtained by summing the particle momentum and dividing by density , impose external forces,

. The

∑

⁄ . To

, on the lattice (i.e. gravity), momentum is added to the macroscopic

velocity in the following way,

where

is the macroscopic velocity calculated from the particle distribution function prior to

collision. Each time iteration of the LB equation proceeds in three steps. First, the particle
66

distributions are streamed to respective nodes. Second, the macroscopic density and velocity are determined from these new particle distributions, and when present, with the inclusion of an exterior force as stated above. Third, the equilibrium particle distribution is obtained and the particle distribution undergoes collision. There are three types of lattice nodes, fluid nodes, bounce-back nodes and solid nodes. The nodes on which fluids move are fluid nodes, the nodes on the surface of solids are bounce-back nodes, and nodes on which no fluids exist are solid nodes. The bounce-back nodes act as algorithmic devices that return incoming streamed particle densities with the opposite momentum, producing a no-slip condition halfway between the fluid node and the bounce-back node.

3.3.3 Shan-Chen Multicomponent Lattice Boltzmann Model

We employ the Shan-Chen Multicomponent (MC) LB model to simulate two-phase immiscible fluid flow (Shan and Doolen 1995, Martys and Chen 1996). The Shan-Chen MC LB model is one of the least computationally demanding multiphase LB models, and for this reason one of the most commonly used multiphase LB models. Two immiscible fluids are simulated on the lattice by representing each fluid phase with its own particle distribution function,

, where

is the index for each component particle distribution function. The component particle distributions interact via a pseudo-potential interparticle force (
∑

where

) defined as,

̅

is the parameter controlling the strength of the interparticle force, when

is positive

the force is repulsive. Only the nearest-neighbor nodes are considered in the calculation of the

67

interparticle force. The adhesion force between components and solid surfaces is created by imposing fictitious component densities,

, on bounce-back nodes,
∑

where

is an indicator function that denotes the existence of a bounce-back node, as in, when of

is occupied by a bounce-back node. Conventionally, negative values

are used for wetting fluids and positive values of

fluids, with

are used for the non-wetting

. The fictitious bounce-back densities are not densities in the

same manner as fluid nodes, they are simply values that control the strength of adhesion. These forces are included in the lattice Boltzmann equation in the same manner as the exterior force,

The total momentum of all particles at each lattice node most be conserved by the collision operator, thus the particle distribution functions for both components must be included in calculations, ∑
∑
The surface densities of bounce back nodes (

) in a Shan-Chen D3Q19 two

component system (wetting and non-wetting) can be estimated for a given contact angle by the following equation proposed by Huang et al. 2007,

68

where

is the density of the wetting fluid and

wetting fluid in the wetting fluid (for values of

is the dissolved density of the nonsufficient to segregate component fluids the

dissolved density is very low). The lattice pressure at each node is determined by the D3Q19
Shan-Chen MC LB model equation of state,
[

]

[

]

All LB simulations presented here were executed using the open source code Palabos.

3.3.4 LB Model Implementation

There are a collection of parameters that are chosen by the user to define a LB model simulation of fluid flow.

The resolution of the lattice is an important factor when designing a LB

simulation, optimally the resolution should be as fine as computing resources allow. However, due to the immense computational needs of LB simulations it is often the case that computing resource limits and time constrictions determine the resolutions of LB fluid flow simulations in porous media. Also, the necessary resolution to adequately represent flow in a porous media is dependent on the pore structure, thus determining an adequate resolution for simulations merits further investigation. Given our computing resources we chose to translate our CMT images with a 1-to-1 correlation, voxel-to-node. This ensures we are using a resolution great enough to capture the fluid-fluid and fluid-solid interfacial areas measured by CMT, and results in lattice sizes within the limits of our computing resources. The imaged voxel resolution of the glass bead pack is 0.0260 x 0.0260 x 0.0292 mm3, to impose a 1-to-1 correlation of voxel-to-node the raw image of the pore geometry is resampled to 0.0260 x 0.0260 x 0.0260 mm3 and segmented into pore and solid voxels. The pore voxels of the segmented image volume are translated as
69

lattice nodes, and the solid voxels that are in contact with pore voxels are translated as bounceback nodes, with all other solid voxels being ignored in the lattice (Figure 3.1). The resolution of the imaged voxels determines the physical size of a lattice length unit, length of a CMT voxel, and the physical mass of lattice mass unit,

, here being the side
, being the volume of a

CMT voxel times the density of the physical fluid being simulated.

Other parameters control the density and viscosity ratio of the fluids and the wettability of the solid surfaces. Throughout the LB simulations presented here the sum of the lattice density of the components,

, is set to 1.0

. The time relaxation parameter , is set to

1 for both fluids; the purpose of this being two-fold. The relaxation parameter determines the viscosity of each component, and thus the viscosity ratio of the fluids. Our experimental fluids
⁄

are kerosene and water with a viscosity ratio,

, of 2.3 and 0.43 for the glass bead

pack and polyethylene bead pack (Note: in the glass bead pack the kerosene is the non-wetting phase, in the polyethylene bead pack the water is the non-wetting phase) displacements, respectively. The capillary number,

, is a dimensionless measure of the ratio of viscous to

capillary forces, and is defined as,

where

is the dynamic velocity,

is the fluid velocity in the direction of flow, and

is the

interfacial tension. The experimental displacements were very low, on the order of
. Given the interfacial tension forces are so much greater than the viscous forces; we assume simulating fluids with a viscosity ratio of 1 will suffice for our needs. Also the use of
BGK dynamics has been shown to result in viscosity-dependent permeability outside timerelaxation parameters of 1 (Pan et al. 2006).
70

A time relaxation parameter of 1 is both

numerically stable and representative of the physical system to be modeled. The interfacial tension of the lattice is determined by the value of

, we use a value of 1.8, as suggested by

Huang et al. 2007, giving a lattice interfacial tension of 0.100. In the SC MC LB model, values of lower than 1.8 decrease the contrast of the component fluids, while higher values result in

an undesirable compression of fluid components (Schaap et al. 2007, Huang et al. 2007). As was previously mentioned, the desired contact angle can be determined by equation proposed by
Huang et al 2007 for given bounce-back node densities

and

. We do not know

exactly what the contact angles are in our experimental system, we are only aware of the fact that the glass bead pack is water-wet and the polyethylene bead pack is weakly oil-wet. To compare
LB results to experimental results we will simulate a range of contact angles.

Figure 3.1:

Segmented CMT image of the bead pack (A) with solid voxels colored white, and the corresponding lattice of the LB model (B) with the bounce-back nodes in blue.

71

3.4 RESULTS
3.4.1 Fluid-Fluid and Fluid-Solid Interfacial Areas

The objective of this investigation is to study flow in mixed-wet porous media, it is important that the model can adequately capture the interfacial areas between solids and fluids. To validate our method we compare the fluid surface and interfacial areas determined by LB simulations to experimental CMT measurements. As previously mentioned and illustrated in figure 3.1 a 1-to-1 correlation of CMT voxel to lattice node is used. The image chosen was 100 3 voxel3, or ~5 beads to a side. Larger lattices up to 2003 voxel3 were also used in initial simulations, but the results were similar to results using the smaller 100 3 voxel3 lattices at an 1/8 of the computational cost. In our steady-state simulations the initial density of each fluid component is uniformly set at a desired saturation with a sum of 1,
⁄

, and an exterior force of

is applied to both component fluids to initialize concurrent flow

resulting in capillary-dominated flow (

). This setup allows the fluids to segregate

without out further input within the pore space. The disadvantage of this setup is it allows fluids access to pores that would not be accessible in displacement processes, and also lead to a more scattered non-wetting fluid distribution than what would be observed experimentally. All other parameters of simulation are the same as those used in the simulation of the circular pore, with the exception of the bounce-back densities (

).

To compare the fluid distributions of the LB simulation to experimental images we measure the specific fluid surface areas and specific interfacial areas.

For details on the experimental

measurement of these areas refer to Landry et al. 2011. To measure these areas and volumes of

72

component fluids from the LB results the lattice must be first segmented to identify the fluid occupying each node, here the fluid occupying the node,

, is defined as the component

fluid with the greatest density,

The volume of the fluid component is,
∑
where the

. The volume of the solid is the product of the sum of all the

bounce-back and non-interacting nodes and solid surface areas,

. The wetting-phase, non-wetting phase, and

, are measured using a modified marching cubes algorithm,

this algorithm wraps a triangular-mesh around the segmented data (Dalla et al. 2002). The specific fluid and solid surface areas,

, are defined as,

⁄

, where

volume of the whole lattice. The specific solid surface area of the lattice,

is the
.

One of the drawbacks of using the LB model to simulate strongly-wetting porous media is the occurrence of a monolayer or 1-node thick wetting phase film surrounding all bounce-back nodes. This film is a consequence of using a mesoscopic model founded on particle dynamics; the pseudo-potential affinity of the wetting phase on nodes adjacent to the bounce-back nodes is too strong to be dislodged by the non-wetting phase.

This monolayer can exaggerate the

measurement of the wetting-phase surface areas, which will have an effect on interfacial areas.
Very little flow occurs in these monolayers, so although they may exaggerate wetting phase surface areas, they should not have a significant effect on flow.

73

We are interested in comparing three specific interfacial areas, the specific fluid-fluid interfacial area, the non-wetting fluid-solid interfacial area and the wetting fluid-solid interfacial area. The specific fluid-fluid interfacial area is the total area of the fluid-fluid interfaces (meniscus) normalized to bulk volume and is defined as (Dalla et al. 2002),

The specific fluid-fluid interfacial area is important to the study of multiphase flows; it is here that energy and mass are transferred between fluid phases. The specific wetting phase-solid interfacial area, as, , and the specific non-wetting phase-solid interfacial area,

are defined

. The solid surface area varies in the images, therefore we normalize the

specific fluid-solid interfacial areas to the specific solid surface area, giving us the fractional wetting fluid-solid interfacial area, area, , and the fractional non-wetting fluid-solid interfacial

, defined as
⁄

.

The wettability alteration that will be investigated here will occur on bounce-back nodes in contact with the non-wetting phase, and the fraction of the surface that will be altered is equal to the fractional non-wetting fluid-solid interfacial area.
,
respectively.

resulting

from

setting

Four contact angles are considered
,

A comparison of LB and experimental measurements of specific fluid-fluid

interfacial areas can be found in figure 3.2. Qualitatively, the LB simulations show a distinctive curve trending towards a maximum at wetting phase saturations between 0.3 and 0.5. This same curve is observed in the experimental measurements from the glass and polyethylene bead packs, and has also been reported by both previous experimental and numerical investigations (Reeves and Celia 1996, Brusseau et al. 2006, Culligan et al. 2006, Schaap et al. 2007, Joekar-Niaser et
74

al. 2009, Porter et al. 2009, Landry et al. 2011). We also observe an increase in specific fluidfluid interfacial area with decreasing contact angles, as would be expected, and is also reflected in the experimental measurements of the weakly oil-wet polyethylene bead pack and the moderately water-wet glass bead pack. A rough estimation by visual inspection from these curves would suggest that the contact angles of the glass and polyethylene bead packs are near and , respectively. However, we cannot make a quantitative comparison of

model and experimental results without knowledge of the experimental contact angles.

In

general, the LB simulations capture the trends of the relationship between saturation and interfacial fluid-fluid areas. The LB simulations also agree with the trends of the experimental measurements of the fractional fluid-solid interfacial areas, as shown in figure 3.3.

The

fractional non-wetting phase interfacial area decreases with decreasing contact angles, this is a result of the wetting phase dominating occupation of the smaller pore spaces and forcing the non-wetting phase into the largest pores where its fractional interface with the solid decreases.
This relationship affects the fraction of a solid surface that is in contact with a wettability altering fluid, and also the location of the fluid in the pore space (i.e. the small pores where clay minerals or organics may be found, or the large pores where multiple mineral constituents may be found.)

75

[
[mm-1] ]

2.4
2.2
2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

0.1
0.2
0.3
0.4
P
Poly beads
G
Glass beads

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Sw
Figure 3.2:

Specific fluid-fluid interfacial areas as a function of wetting phase saturation for the LB simulations and experimental CMT images.

1
0.9
0.8
0.7
0.1
0.2
0.3
0.4
Poly beads
P
Glass beads
G

Xfs[]

0.6
0.5

0.4
0.3

0.2
0.1

0
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Sw
Figure 3.3:

Fractional fluid-solid interfacial areas as a function of wetting phase saturation for the LB simulations and experimental CMT images, the wetting fluid-solid interfacial areas are shown in black and the non-wetting fluid-solid interfacial areas are shown in gray.

3.4.2 Pore Aperture and Fluid Distribution
76

To measure the pore aperture distribution of the pore space, first a skeleton or medial axis, of the pore space must be constructed. The skeleton is constructed by thinning algorithms (Lee et al.
1994) using the software Avizo Fire 6.3. The skeleton of the pore space is a connected network of voxels that represent the midpoint between solid walls in the pore space. Each voxel in the skeleton can be thought of as the center of a sphere that makes contact with the wall in at least two places. There is a significant amount of literature regarding the subject of pore space topology description (Thovert et al. 1993, Lindquist et al. 1996, Bakke and Oren 1997,
Prodanovic et al. 2006). The pore aperture is the maximum radius of a sphere that can occupy a pore. To determine the pore aperture distribution, a “sphere packing” of sequentially smaller spheres are strung along the skeleton. First, the voxels of skeleton are labeled with their distance to the pore wall. The labeled skeleton is then searched for voxels that fall within a user-defined range. These voxels then become the center of a sphere – with a radius equal to the distance to the closest pore wall – placed in the pore space, and all voxels of the pore space that fall within this sphere are labeled.

The pore aperture distribution can then be determined by simply

summing the labeled voxels. The fluid distribution can also be easily measured by masking the segmented fluid occupation (

) with the pore aperture-labeled pore space.

The pore

aperture distribution and fluid distribution are presented in figure 3.4. We can see from the fluid distribution that decreasing the contact angle pushes the non-wetting phase into the larger pores.
This is a result of the increasing affinity of the wetting phase for the solid surfaces, resulting in the wetting phase dominance of smaller pore spaces. This can have a substantial effect on the dependence of fluid mobility on wettability alteration. If only the larger pores that would be

77

occupied by the non-wetting phase in a homogenous-wet state are altered the effect on the relative permeability of the fluids will be minimal.

Volume-weighted fraction

0.3
0.25
0.2
0.15
0.1
0.05
0
0

0.026 0.052 0.078 0.104 0.13 0.156 0.182 0.208
Pore aperture [mm]

Non-wetting phase fluid fraction

1.2
0.1
0.2
0.3
0.4

1
0.8
0.6
0.4
0.2
0
0

0.026 0.052 0.078 0.104 0.13 0.156 0.182 0.208
Pore aperture [mm]

Figure 3.4:

Pore aperture and fluid distribution of initial homogenous-wet states.

78

3.4.3 Two-phase Flow in a Circular Pore

To test our model against a known semi-analytical solution for two-phase immiscible flow we simulate flow in a circular pore. Two-phase immiscible flow is commonly described using an empirical extension of Darcy’s law,

where the mean velocity of each fluid parallel to the direction of flow, dynamic viscosity,

, pressure gradient,

relative permeability defined as

, is a function of the

, absolute or single-phase permeability,
, where

, and the

is the effective permeability. In a

circular pore where the wetting fluid is distributed as an annulus along the walls of the pore and the viscosity ratio is one, the wetting phase relative permeability is a function of the wetting phase saturation,

, and the non-wetting phase relative permeability is a function of
⁄ (Ramstad et al. 2010). To simulate two-

non-wetting phase saturation,

phase flow in a circular pore at varying saturations the wetting phase fluid component is initially distributed as an annulus along the pore walls with the non-wetting phase fluid component occupying the nodes in the center in correlation with the desired saturation. The walls of the pore are strongly-wetting to maintain the annular distribution of the wetting phase. Otherwise, under weakly-wetting conditions, the wetting phase will separate and lose its annular shape, the annular shape being an assumption of the semi-analytical solution of immiscible two-phase flow in a circular pore. To determine the relative permeability of component fluids in the LB simulation the momentum in the direction of flow of each fluid is summed over all pore nodes at the given wetting phase saturations,

79

∑ and normalized to the fully saturated momentum,

To impose concurrent flow in a circular pore with a radius of 11 nodes an exterior force of
⁄

is applied to both component fluids, and periodic boundaries

are imposed at the inlet and outlet. The system is considered to have reached steady-state when the convergence criteria were met,

A choice of

was chosen, based on numerical experiments. In figure 3.5 are the results

of the LB simulation and the semi-analytical solutions for two-phase flow in a circular pore.
The difference between the semi-analytical solution and the LB simulations is due to the discrete nature of the rendering of the circular pore and the resolution. Previous investigations have shown multicomponent LB models are increasingly more accurate with finer resolutions
(Ramstad et al. 2010).

80

1
0.9
0.8

0.7

kr

0.6

0.5
0.4
0.3
0.2
0.1

0
0

0.2

0.4

0.6

0.8

1

S(w)
Figure 3.5:

Relative permeability in a circular pore with radius = 11

; LB simulation results and

semi-analytical solution curves.

3.4.4 Relative Permeability

The relative permeability of the porous media is determined by LB simulations in the same manner as described in the section Two-phase flow in a circular pore. With the simulation volume being bounded by walls perpendicular to flow, and as was used in the circular pore simulations, periodic boundaries are applied in the direction flow (i.e. the inlet and outlet).
However, as was described in Martys and Chen 1996, the total momentum of the component lattices did not display a linear relationship with exterior forcing below

.

To

determine the relative permeability exterior forces between
⁄

were used. The relative permeability curves for the initial homogenous-wettability
81

porous media are shown in figure 3.6. At low wetting phase saturations the non-wetting phase relative permeability exceeds one and the wetting phase relative permeability is slightly negative.
This has also been reported in the preliminary results of Boek and Ventourilli 2010. These saturations are at or below experimentally observed irreducible wetting phase saturations of and for the glass and polyethylene beads respectively, and are unlikely to

occur in displacement processes. We are able to measure the relative permeability at these saturations using the LB simulations due to the employment of a steady-state setup which allows

kr

for any initial saturation. However, some of these simulations may not have a physical corollary.

1.2
1.1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1

0.1
0.2
0.3
0.4

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Sw
Figure 3.6:

Relative permeability determined by LB simulations for the initial homogenous-wettability states. The relative permeability of the wetting and non-wetting phase are shown in black and gray, respectively.

82

Generally, figure 3.6 shows decreasing the contact angle also decreases the relative permeability of both phases at all saturations. The variation is small due to the high porosity and connectivity of the pore space. Also, as would be expected, the crossover point for the relative permeability curves occurs at

. The relative permeability of the fluids is governed by two competing

mechanisms that act to inhibit and promote fluid flow, namely, fluid-solid interfacial area and connectivity, respectively. Greater fluid-solid interfacial area leads to greater resistance to flow, while increased connectivity leads to greater mobility.

The dominant factor being the

connectivity of the phase, which results in the general trend of increasing phase relative permeability with increasing phase saturation. Unlike the fluid-solid interfacial area, we do not have a simple method of measuring the connectivity of a fluid phase, but given that we have measurements of the solid-fluid interfacial areas we can infer differences in relative permeability that cannot be attributed to trends in the solid-fluid interfacial areas are a result of connectivity.
Also the connectivity correlates with the fluid distribution, with the more connected wetting phase occupying smaller pores, as was shown in figure 3.4.

How the wetting phase and non-wetting phase respond to these competing mechanisms can be very different.

In figure 3.6 we can see the net effect on the non-wetting phase relative

permeability is an overall decrease as the fluids become more disconnected (occupy larger pores). Although there is less fluid-solid interfacial area this effect is dominated by the decrease in connectivity. Unlike the non-wetting phase, the wetting phase relative permeability shows very little dependence on the wetting strength of the porous medium for

. We know that

there will not be a decrease in the connectivity of the wetting phase as the contact angle is decreased, thus the small decrease in mobility between

83

and the other homogenous-wet

states can only be attributed to an increase in fluid-solid interfacial areas. It is not surprising that the non-wetting phase relative permeability is more greatly influenced by connectivity, being the less connected phase, and the wetting phase relative permeability is more greatly influenced by the fluid-solid interfacial area, being the generally more connected phase. However, regarding the non-wetting phase relative permeability, a similar previous investigation reported the opposite effect from what is reported here. Li et al. 2005 simulating two-phase flow in a homogenous sphere pack using a multiple-relaxation-time Shan-Chen type multicomponent LB model reported an increase in the non-wetting phase relative permeability for neutral wet states
(

) over strongly wet states (

). This could be attributed to differences in the pore

space geometry resulting in the prevalence of increased non-wetting phase mobility due to the decrease in solid-fluid interfacial area, over the decrease in mobility associated with decreased connectivity as was observed here. The homogenous sphere pack used by Li et al. 2005 may have a generally narrower distribution of pore apertures, meaning fewer large pore spaces for the non-wetting phase to occupy, increasing the dependence of the non-wetting phase relative permeability on the solid-fluid interfacial area. The interaction of these competing mechanisms is sensitive to small differences in pore space geometry. In the mixed-wet states individual pores can have mixed-wettability resulting in the competition of these mechanisms within each pore.

3.4.5 Relative Permeability of Mixed-Wet States

As was previously stated, the bounce-back density of the nodes in contact with the non-wetting phase after fluid distribution is established in the homogenous-wetting simulations, are altered to create a mixed-wet porous medium (figure 3.7). The mixed-wet state is the result of three

84

parameters, the initial homogenous-wettability, the saturation at which alteration takes place, and the severity of alteration. Here we will alter the wettability of each of four homogenous-wetting states at a wetting phase saturation near 0.5, with four increasingly severe alterations, resulting in sixteen mixed-wet states. The fictional bounce-back density of nodes in contact with the nonwetting phase are altered to preferentially wet the non-wetting phase with contact angles
(

). To avoid confusion, the initial

wetting phase is from here on referred to as the wetting phase of the mixed-wet states, and used as the wetting phase saturation in figures. according to the following indexing, is the

Also, the mixed-wet states will be referred to

, where

for the altered wettability (i.e.

is the

for the initial wettability, and

would refer to a state that was initially

wet, then altered as previously stated to

).

The initial homogenous-wet states will use the same indexing as the mixed-wet states absent a superscript. There are numerous mixed-wet states that could be considered using this

framework, we limit this investigation to studying the mixed-wet state resulting at

.

These mixed-wet states are the result of wettability alteration to the homogenous-wet states at
, and therefore the mixed-wet states are dependent on the fluid distribution of the phases at

. The wetting phase at this saturation is described as occupying the smaller

pore spaces (figure 3.4), and is the more well connected phase. Thus the surfaces altered by the non-wetting phase are disconnected, with an increasing fraction of the surface altered (increasing fractional solid-fluid interfacial area) with increasing contact angle.

85

Figure 3.7:

Images of the LB lattice with initial homogenous-wettability (A), the fluid distribution of the non-wetting phase at the end of the simulation of the initial homogenous-wettability with and (B), and the alteration of the bounce-back densities of the nodes in

contact with the non-wetting fluid (C). This image also summarizes the three parameters determining the mixed-wet state, initial wettability (A), saturation of alteration (B), and severity of alteration (C).

The relative permeability curves and the fractional solid-fluid interfacial areas for the mixed-wet states are shown in figure 3.8. Recapping the observations of the homogenous-wet states, the non-wetting phase was found to decrease with decreasing contact angle due to decreasing connectivity, while the wetting phase was found to decrease with decreasing contact angle due to increasing solid-fluid interfacial area. In general, the wettability alteration has little effect on the wetting phase relative permeability. The fractional wetting phase solid-fluid interfacial areas of the mixed-wet states are generally lower than those of the mixed-wet states. Unlike in the homogenous-wet states this decrease does not correlate with an increase of the wetting phase relative permeability.

We can conclude that any differences in the wetting phase relative

permeability of mixed-wet states cannot be attributed to differences in the solid-fluid interfacial area. Instead the insensitivity of the wetting phase relative permeability to wettability alteration can be attributed to the fluid distribution of the wetting phase at the saturation of alteration. At
86

the wetting phase has developed a well-connected pathway (occupies smaller pores, has a high fractional solid-fluid interfacial area). The solid in contact with the wetting phase remains unaltered, thus, this connected path remains in the mixed-wet states. For the mixed-wet states there is some decrease in the wetting phase relative permeability observed for the mixed-wet state. This decrease is limited to the

mixed-wet state as a result of its

poorer connectivity of the established unaltered wetting phase pathways in comparison to those of the

,

, and

mixed-wet states. But even this small decrease only exists for a

mixed-wet state suffering a severe alteration. Overall, the wetting phase relative permeability is unaffected by the wettability alteration, the established unaltered wetting phase pathways that encourage the wetting phase to occupy the smaller pores dominates the mobility of the wetting phase. 87

B)

1.1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1

1

0.9
0.8
0.7

0.10.0
0.1
0.2
0.3

0.6

0.10.0
0.1
0.2
0.3

kr

kr

A) 1.2

0.5
0.4
0.3
0.2
0.1
0

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

0

Sw
1

0.9
0.8
0.7

0.10.0
0.1
0.2
0.3

0.6

0.10.0
0.1
0.2
0.3

kr

kr

D)

1.1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1

0.5
0.4
0.3
0.2
0.1
0

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0

1

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

Sw

Sw

E) 1.2

F)

1.1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1

1

0.9
0.8
0.7

0.10.0
0.1
0.2
0.3

0.6

0.10.0
0.1
0.2
0.3

kr

kr

1

Sw

C) 1.2

0.5
0.4
0.3

0.2
0.1
0

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

0

Sw

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

Sw

G)1.2

H)

1.1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1

1

0.9
0.8
0.7

0.10.0
0.1
0.2
0.3

0.6

0.10.0
0.1
0.2
0.3

kr

kr

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.5
0.4
0.3
0.2
0.1
0

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0

1

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Sw

Sw

88

1

Figure 3.8:

Relative permeability and fractional fluid-solid interfacial areas determined by LB simulations for the mixed-wet states of the originally homogenous-wet states, (A, B)
, (C, D)

, (E, F)

, and (G, H)

,

with increasing severity of alteration. The relative permeability and fractional fluid-solid interfacial areas of the wetting and non-wetting phase are shown in black and gray, respectively. The non-wetting phase relative permeability is significantly altered in the mixed-wet states, and as in the case of the wetting phase, the most significant effect is seen for the states. For the

,

, and

mixed-wet

mixed-wet states the severity of the wettability alteration

has little effect on the mobility of the non-wetting phase, but the existence of a wettability alteration does. This can be understood as a result of solid-fluid interaction on altered solid surfaces. These surfaces are no longer dominated by the initial wetting phase, resulting in an increase of the non-wetting phase fluid-solid interfacial area and a decrease in the non-wetting phase mobility. The increase in the non-wetting phase relative permeability correlates with an increase in the non-wetting phase fluid-solid interfacial area. In general, the non-wetting phase becomes pinned by its affinity to these altered solid surfaces, the lubricating effect the wettingphase has on the non-wetting phase mobility is lost. The non-wetting phase responds to the strength of alteration much as the wetting phase responded in the homogenous-wet states to increasing strength of wettability. Unlike the

,

, and

mixed-wet states, for the

mixed-wet states the severity of alteration appears to be proportional to the decrease in the non-wetting phase relative permeability. This also can be considered in the context of the fluid distribution of the homogenous-wet

, being close to neutrally-wet the non-wetting phase

occupies smaller pore spaces and has a greater fluid-solid interfacial area than the more strongly
89

wetting scenarios. Thus, the alteration of the surface occurs not only on a larger fraction of the solid surface, but also in the smaller pore spaces – much like the wetting phase in the homogenous-wet states, the non-wetting phase relative permeability decreases due to the greater fluid-solid interfacial areas. This increase in fluid-solid interfacial areas is reflected in the measurements of non-wetting fluid-solid interfacial areas (figure 3.8). Incremental increases in the non-wetting fluid-solid interfacial area for the

mixed-wet states resulted in

incrementally decreasing non-wetting phase permeability. The non-wetting fluid-solid interfacial area of the

,

, and

mixed-wet states increases significantly with the existence of

an alteration, but the increase does not respond to increasing severity of alteration, which is reflected in the relative permeability curves. Overall the non-wetting phase responded to the wettability alteration with a decrease in relative permeability as a result of increasing nonwetting fluid-solid interfacial areas.

Had the wettability alteration occurred at high wetting phase saturations we would expect the same result to a lesser extent, simply due to less of the solid surface being altered. Wettability alteration at lower wetting phase saturations would have a somewhat increasing effect, however, not as much as one may at first think. The extreme example of altering the entire surface when fully saturated with the non-wetting phase would result in a simple role reversal of the wetting and non-wetting phases.

This would suggest that at lower wetting phase saturations of

wettability alteration we could expect the wetting phase relative permeability to “rebound” and begin to respond by increasing with decreasing fluid-solid interfacial area. However, saturations below irreducible wetting phase saturation are unlikely to occur.

At saturations near the

irreducible wetting phase saturation the wetting phase will dominate the smaller pore spaces.

90

And as occurred in our simulations, wettability alteration will not occur in these smaller pore spaces. The wetting phase will continue to be pinned to these smaller pore spaces, and generally will not experience any increase in mobility, as one might intuitively expect when considering the extreme example of the entire surface being altered. At these saturations the non-wetting phase will show a slightly greater effect than what we see here at

, however the

differences will not be great. In our simulations of mixed-wet states the non-wetting phase relative permeability has been significantly decreased by the wettability alteration.

There is

very little mobility left to be lost.

Comparable studies to this one are limited for the most part to numerical experiments using pore network modeling (Blunt et al. 2002, Valvatne and Blunt et al. 2004, Zhao et al. 2010, Gharbi and Blunt 2012). There are a few notable differences between PNMs and LB models. As was described in the background section, pore network models use networks of capillary tube elements to represent the pore space topology, while LB models use direct translations of CMT images to construct the lattice.

This also means that wetting films are simulated by very

different means. In PNMs the wetting phase is simulated as thin films, allowing for very low wetting phase saturations to maintain connectivity.

Consequently in mixed-wet scenarios it is

possible for a pore to contain a thin film of one phase sandwiched between the other phase occupying the center and corners of the pore (Blunt et al. 2002). The wetting phase films in a lattice Boltzmann simulation are limited in their thickness to the resolution of the lattice, and may or may not appear depending on the resolution of the lattice and the strength of the bounceback density (

) used.

Also PNMs simulate displacements, unlike here, where our

simulations are steady state. The PNM study of mixed-wettability in six different types of

91

limestone of Gharbi and Blunt 2012 offers a fair comparison to this work.

This work

investigated the mixed-wet states of PNMs constructed from CMT images of limestone cores.
They also found a significant reduction in the non-wetting phase (oil) relative permeability for wettability alteration fractions of 0.25 and 0.5 – similar to solid surface fractions altered at here. However, they also report a significant reduction in the wetting phase relative permeability, not seen here. This can be attributed to not only differences in the pore space geometry, but also the method in which fluid films are simulated. At remains well connected for homogenous-wet states

,

our wetting phase
, and

, and as was

previously stated there remains well-connected unaltered solid-surface for the wetting-phase to be established, ensuring an insignificant effect on the wetting phase relative permeability for the mixed-wet states. In the mixed-wet states of Gharbi and Blunt 2012 much of the wetting phase existed as thin films, resulting in very low wetting phase permeability. The inability to simulate fluid films beyond the resolution of the lattice is one of the weaknesses of the LB method.
However the ability of the LB method to simulate fluids in realistic pore geometries makes it a very attractive technique.

3.5 CONCLUSIONS

We have presented a study of the dependence of relative permeability on wettability for both homogenous-wet and mixed-wet states using LB modeling. Comparison between LB simulation and experimental results show the LB simulations capture the trends of fluid-fluid and fluid-solid interfacial areas well. Since our mixed-wet states are created by altering the wettability of the

92

solid in contact with the non-wetting phase it is important that these interfacial areas are accurately captured to create realistic mixed-wet states.

The homogenous-wet states display a decrease in relative permeability for both phases with decreasing contact angles. This is attributed here to the competition of two mechanisms that increase and decrease fluid flow. With decreasing contact angle the non-wetting phase becomes more disconnected and occupies larger pore spaces. This occupation results in less non-wetting fluid-solid interfacial areas which acts to increase the mobility of the non-wetting phase, however, the decreasing connectivity dominates this effect resulting in lower non-wetting phase relative permeability. The decrease in the wetting phase relative permeability with decreasing contact angle is a result of greater fluid-solid interfacial area which acts to decrease the mobility of the wetting phase; any increases in connectivity of the wetting phase cannot overcome this effect, due to the pinning of the wetting phase to the smaller pore spaces.

The pores of the mixed-wet state are themselves mixed, the interaction between the two competing mechanisms controlling flow will also compete within individual pores. This can lead to a rich assortment of fluid mobility behavior as a result of varying the wettability alone.
In this investigation we find that the connectivity of the wetting phase at the saturation of wettability alteration is critical to the significance of the effect the alteration to a mixed-wet state will have on the wetting phase mobility. In general, the severity of alteration is inconsequential if the unaltered solid surface is well connected. As is the case for initial homogenous-wet states with a

at

The non-wetting phase relative permeability demonstrates a

significant decrease in mobility proportional the severity of alteration when the solid surfaces

93

being altered includes a significant portion of the smaller pore spaces. As is the case for initial homogenous-wet state with a at at

For initial homogenous-wet states with a

the resulting mixed-wet states all show a significant decrease in the non-

wetting phase mobility, but no dependence on the severity of the alteration, even when the solid surface is altered to neutrally wet. Overall the most important factor governing non-wetting phase mobility in mixed-wet states is the fluid-solid interfacial areas, greater interfacial fluidsolid surface areas lead to decreases in mobility. The responses to a wettability alteration at a wetting phase saturation of 0.5 shown here, are comparable to what we would expect for wettability alteration at lower wetting phase saturations due to the pinning of the wetting phase to the smaller unaltered pore spaces.

The framework presented here to study the evolution of wettability alteration and its effect on the relative permeability of two-phase flow in a simple porous medium could easily be applied to natural consolidated and unconsolidated porous media.

Our initial wettability states are

homogenous, for natural porous media the initial wettability will be determined by the mineral constituents and their locations in the pore space geometry. We present simulations of varying alteration severity, for natural porous media the level of wettability alteration is a result of complex fluid-mineral interactions and would likely require bulk studies of wettability alteration for individual mineral constituents. Given this knowledge can be determined using imaging and laboratory techniques, LB simulations can be easily included to estimate fluid flow properties, such as relative permeability, for a wide assortment of scenarios.

94

3.6 REFERENCES

Al-Futaisi, and Patzek T.W. Secondary imbibition in NAPL-invaded mixed-wet sediments.
Journal of Contaminant Hydrology 2004; 74:61-81. doi:10.1016/j.jconhyd.2004.02.005
Hui M., and Blunt M.J. Wettability on three-phase flow in porous media. Journal of Physical
Chemistry B 2000; 104:3833-3845.

Alipour Tabrizky V., Hamouda A. A., and Denoyel R. Influence of magnesium and sulfate ions on wettability alteration of calcite, quartz, and kaolinite: Surface energy analysis. Energy and
Fuels 2011; 25:1667-1680. doi:10.1021/ef200039m.

Al-Raoush R.I. Impact of wettability on pore-scale characteristics of residual nonaqueous phase liquids. Environmental Science and Technology 2009; 43:4796-4801. doi:10.1021/es802566s

Bakke S., and Oren P.E. 3-D pore-scale modeling of sandstones and flow simulations in the pore networks. SPE Journal 1997; 2(2):136-149. doi:10.2118/35479-PA

Barranco Jr. F.T., Dawson H.E., Christener J.M., and Honeyman B.D. Influence of aqueous pH and ionic strength on the wettability of quartz in the presence of dense non-aqueous-phase liquids. Environmental Science and Technology 1997; 37:676-681.

95

Boek E.S., and Venturoli M. Lattice-Boltzmann studies of fluid flow in porous media with realistic rock geometries. Computers and Mathematics with Applications 2010; 59:2305-2314. doi:10.1016/j.camwa.2009.08.063 Blunt M.J.

Effects of heterogeneity and wetting on relative permeability using pore level

modeling. SPE Journal 1997; 2(1):70-87. doi:10.2118/36762-PA

Blunt M.J. Pore level modeling of the effects of wettability. SPE Journal 1997; 2(4):494-510. doi:10.2118/38435-PA Blunt M.J., Jackson M.D., Piri M., and Valvatne P.H. Detailed physics, predictive capabilities and macroscopic consequences for pore network models of multiphase flow. Advances in Water
Resources 2002; 25:1069-1089.

Brusseau M.L., Peng S., Schnaar G., and Constanza-Robinson M.S. Relationships among airwater interfacial area, capillary pressure, and water saturation for a sandy porous medium.
Water Resources Research 2006; 42:W03501. doi:10.1029/2005WR004058

Brusseau M.L., Janousek H., Murao A., and Schnaar G. Synchrotron X-ray microtomography and interfacial partitioning tracer test measurements of NAPL-water interfacial areas. Water
Resources Research 2008; 44:W01411. doi:10.1029/2006WR005517

96

Buckley J.S. and Liu Y. Some mechanisms of crude oil/brine/solid interactions. Journal of
Petroleum Science and Engineering 1998; 20:155-180. doi:10.2118/37230-PA

Coles M.E., Hazlett R.D., Spanne P., Soll W.E., Muegge E.L., and Jones K.W. Pore level imaging of fluid transport using synchrotron X-ray microtomography. Journal of Petroleum
Science and Engineering 1998; 19:55-63.

Constanza-Robinson M.S., Harrold K.H., and Lieb-Lappen R.M.

X-ray microtomography

determination of air-water interfacial area-water saturation relationships in sandy porous media.
Environmental Science and Technology 2008; 42: 2949-2956. doi:10.1021/es072080d

Culligan K.A., Wildenschild D., Christensen B.S.B., Gray W.G. and Rivers M.L. Pore-scale characteristics of multiphase flow in porous media: A comparison of air-water and oil-water experiments. Advances

in

Water

Resources

2006;

29:227-238.

doi:10.1016/j.advwatres.2005.03.021

Dalla E., Hilpert M., and Miller C.T. Computation of the interfacial area for two-fluid porous medium systems. Journal of Contaminant Hydrology 2002; 56:25-48.

Dixit A.B., McDougall S.R., and Sorbie K.S. A pore-level investigation of relative permeability hysteresis in water-wet systems. SPE Journal 1998; 3(2):115:123. doi:10.2118/37233-PA

97

Evje S., and Hiorth A. A model for interpretation of brine-dependent spontaneous imbibition experiments. Advances

in

Water

Resources

2011;

34:1627-1642.

doi:10.1016/j.advwatres.2011.09.003

Fenwick D.H., and Blunt M.J. Network modeling of three-phase flow in porous media. SPE
Journal 1998; 3(1):84-96. doi:10.2118/38881-PA

Gharbi O., and Blunt M.J. The impact of wettability and connectivity on relative permeability in carbonates: A pore network modeling analysis. Water Resources Research 2012; 48:W12513. doi:10.1029/2012WR011877 Ghassemi A., and Pak A. Numerical study of factors influencing relative permeabilities of two immiscible fluids flowing through porous media using lattice Boltzmann method. Journal of
Petroleum Science and Engineering 2011; 77:135-145. doi:10.1016/j.petrol.2011.02.007

Hao L., and Cheng P. Pore-scale simulations on relative permeabilities of porous media by lattice Boltzmann method. International Journal of Heat and Mass Transfer 2010; 53:1908-1913. doi:10.1016/j.ijheatmasstransfer.2009.12.066 Hazlett R.D., Chen S.Y. and Soll W.E. Wettability and rate effects on immiscible displacement:
Lattice Boltzmann simulation in microtomographic images of reservoir rocks. Journal of
Petroleum and Science Engineering 1998; 20:167-175

98

Hoiland L.K., Spildo K., and Skauge A.

Fluid flow properties for different classes of

intermediate wettability as studied by network modeling. Transport in Porous Media 2007;
70:127-146. doi:10.1007/s11242-006-9088-x

Huang H., Thorne D.T., Schaap M.G., and Sukop M.C. Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models. Physical
Review E 2007; 76:066701. doi: 10.1103/PhysRevE.76.066701

Iglauer S., Paluszny A., Pentland C.H., and Blunt M.J.

Residual CO 2 imaged with X-ray

microtomography.

Letters

Geophysical

Research

2011;

38:L21403.

doi:10.1029/2011GL049680

Jackson M.D., Valvatne P.H., and Blunt M.J. Prediction of wettability variation and its impact on flow using pore- to reservoir-scale simulations.

Journal of Petroleum Science and

Engineering 2003; 39:231-246. doi:10.1016/S0920-4105(03)00065-2

Jadhunandan P.P., and Morrow N.R. Effect of wettability on waterflood recovery for crudeoil/brine/rock systems. SPE Reservoir Engineering 1995; 10(1):40-46. doi:10.2118/22597-PA

Jerauld G.R., and Rathmell J.J. Wettability and relative permeability of Prudhoe Bay: A case study in

mixed-wet

reservoirs.

SPE

Reservoir

doi:10.2118/28576-PA

99

Engineering

1997;

12(1):58-65.

Joekar-Niaser V., Hassanizadeh, S.M., and Leijnse A. Insights into the relationships among capillary pressure, saturation, interfacial area and relative permeability using pore network modeling. Transport in Porous Media 2008; 74:201-219. doi:10.1007/s11242-007-9191-7

Kovscek A.R. Wong H., and Radke C.J. A pore-level scenario for the development of mixed wettability in oil reservoirs. AIChE Journal 1993; 39(6):1072-1085.

Kumar M., Fogden A., Senden T., and Knackstedt M.

Investigation of pore-scale missed

wettability. SPE Journal 2012; 17(1):20-30. doi: 10.2118/129974-PA

Landry C.J., Karpyn Z.T., and Piri M. Pore-scale analysis of trapped immiscible fluid structures and fluid interfacial areas in oil-wet and water-wet bead packs. Geofluids 2011; 11:209-227. doi: 10.1111/j.1468-8123.2011.00333.x

Lebedeva E.V., and Fogden A.

Micro-CT and wettability analysis of oil recovery from

sandpacks and the effect of waterflood salinity and kaolinite. Energy and Fuels 2011; 25:56835694. doi:10.1021/ef201242s

Lee T.C., Kashyap R.L. and Chu C.N. Building skeleton models via 3D medial surface/axis thinning algorithms. CVGIP, Graphical Models and Image Processing 1994; 54:462-478.

Li H., Pan C., and Miller C.T. Pore-scale investigation of viscous coupling effects for two-phase flow in porous media. Physical Review E 2005; 72:026705. doi:10.1103/PhysRevE.72.026705

100

Lindquist B., Lee S.M., Jones K.W., and Spanne P. Medial axis analysis of void structure in three-dimensional tomographic images of porous media. Journal of Geophysical Research 1996;
101(B4):8297-8310.

Martys N.S., and Chen H. Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method. Physical Review E 1996; 53(1):743-750.

Oren P.E., Bakke S., and Arntzen O.J. Extending predictive capabilities to network models.
SPE Journal 1998; 3(4):324-336. doi:10.2118/52052-PA

Oren P.E., and Bakke S. wettability effects.

Reconstruction of Berea sandstone and pore-scale modeling of

Journal of Petroleum Science and Engineering 2003; 39:177-199.

doi:10.1016/S0920-4105(03)00062-7

Pan C., Hilpert M., and Miller C.T. Lattice-Boltzmann simulation of two-phase flow in porous media. Water Resources Research 2004; 40:W01501. doi:10.1029/2003WR002120
Pan C., Luo L.S., and Miller C.T. An evaluation of lattice Boltzmann schemes for porous medium flow simulation. Computers & Fluids 2006; 35:898-909. doi:10.1016/j.compfluid.2005.03.008 101

Porter M.L., Schaap M.G., and Wildenschild D. Lattice-boltzmann simulations of the capillary pressure-saturation-interfacial area relationship for porous media. Advances in Water
Resources; 32:1632-1640. doi:10.1016/j.advwatres.2009.08.009

Prodanovic M., Lindquist W.B., and Seright R.S. 3D image-based characterization of fluid displacement in

a

Berea

core.

Advances

in

Water

Resources;

30:214-226.

doi:10.1016/j.advwatres.2005.05.015

Prodanovic M., Lindquist W.B., and Seright R.S. Porous structure and fluid partitioning in polyethylene cores from 3D x-ray microtomographic imaging. Journal of Colloid and Interface
Science 2006; 298:282-297. doi:10.1016/j.jcis.2005.11.053

Qian Y.H., d’Humieres D., and Lallemand P. Lattice BGK Models for Navier-Stokes Equations.
Europhysics Letters 1992; 17(6):479-484. doi:10.1209/0295-5075/17/6/001

Ramstad T., Oren P.E., and Bakke S. Simulation of two-phase flow in reservoir rocks using a lattice Boltzmann method. SPE Journal 2010; 15(4):917-927. doi: 10.2118/124617-PA

Ramstad T., Idowu N., Nardi C., and Oren P.E. Relative permeability calculations from twophase flow simulations directly on digital images of porous rock. Transport in Porous Media
2012; 94:487-504. doi:10.1007/s11242-011-9877-8

102

Reeves P.C., and Celia M.A. A functional relationship between capillary pressure, saturation, and interfacial area as revealed by a pore-scale network model. Water resources research 1996;
32(8):2345-2358.

Robin M., Rosenberg E., and Fassi-Fihri O. Wettability studies at the pore level: A new approach by use

of cryo-SEM.

SPE

Formation Evaluation 1995;

10(1):11-19.

doi:10.2118/22596-PA

Saraji S., Goual L., and Piri M.

Adsorption of asphaltenes in porous media under flow

conditions. Energy and Fuels 2010; 26:6009-6017. doi:10.1021/ef100881k.

Schaap M.G., Porter M.L., Christensen B.S.B, and Wildenschild D. Comparison of pressuresaturation characteristics derived from computed tomography and lattice Boltzmann simulations.
Water Resources Research; 43:W12S06. doi:10.1029/2003WR002120

Schembre J.M., Tang G.Q., and Kovscek A.R. Wettability alteration and oil recovery by water imbibition at elevated temperatures.

Journal of Petroleum Science and Engineering 2006;

52:131-148. doi:10.1016/j.petrol.2006.03.017

Shan X., and Doolen G. Multicomponent lattice-Boltzmann model with interparticle interaction.
Journal of Statistical Physics 1995; 81:379-393.

103

Silin D., Tomatsu L., Benson S.M., and Patzek T.W. Microtomography and pore-scale modeling of two-phase fluid distribution. Transport in Porous Media; 86:495-515. doi:10.1007/s11242010-9636-2

Spiteri E.J., Juanes R., Blunt M.J., and Orr Jr. F.M. A new model of trapping and relative permeability hysteresis for all wettability characteristics.

SPE Journal; 13(3):277-288.

doi:10.2118/96448-PA

Sukop M.C., Huang H., Chen L.L., Deo M.D., Oh K., and Miller J.D. Distribution of multiphase fluids in porous media: Comparison between lattice Boltzmann modeling and micro-x-ray tomography. Physical Review E 2008; 77:026710. doi: 10.1103/PhysRevE.77.026710

Svirsky D.S., van Dike M.I.J., and Sorbie K.S. Prediction of three-phase relative permeabilities using a pore-scale network model anchored to two-phase data. SPE Reservoir Evaluation and
Engineering; 10(5):527-538. doi:10.2118/89992-PA

Thovert J.F., Salles J., and Adler P.M. Computerized characterization of the geometry of real porous media: their discretization, analysis and interpretation. Journal of Microscopy 1993;
170(1):65-79.

Turner M.L., Knufing L., Arns C.H., Sakellariou A., Senden T.J., Sheppard A.P., Sok R.M.,
Limaye A., Pinczewski W.V., and Knackstedt M.A. Three-dimensional imaging of multiphase flow in porous media. Physica A 2004; 339:166-172. doi:10.1016/j.physa.2004.03.05

104

Vogel H.J., Tolke J., Schulz V.P., Krafczyk M., and Roth K. Comparison of a LatticeBoltzmann model, a full-morphology model, and a pore network model for determining capillary pressure-saturation relationships. Vadose Zone Journal 2005; 4:380-388. doi:10.2136/vzj2004.0114 Yu L., Evje S., Kleppe H., Karstad T., Fjelde I., and Skjaeveland S.M. Spontaneous imbibition of seawater into prefentially oil-wet chalk cores – Experiments and simulations. Journal of
Petroleum Science and Engineering 2009; 66:171-179. doi:10.1016/j.petrol.2009.02.008

Zhao X., Blunt M.J., and Yao J. Pore-scale modeling:effects of wettability on waterflood oil recovery. Journal

of

Petroleum

Science

doi:10.1016/j.petrol.2010.01.011

105

and

Engineering;

71:169-178.

CHAPTER 4: PORE-SCALE LATTICE BOLTZMANN MODELING AND 4D X-RAY
COMPUTED MICROTOMOGRAPHY IMAGING OF FRACTURE-MATRIX FLUID
TRANSFER

4.1 SUMMARY

We present sequential x-ray computed microtomography (CMT) images of matrix drainage in a fractured sintered glass granule pack. Sequential imaging captured the capillary-dominated migration of the non-wetting phase front from the fracture to the matrix in a brine-surfactantDecane system. The sintered glass granule pack was designed to have minimal pore space beyond the resolution of CMT imaging, so that the pore space of the matrix connected to the fracture could be captured in its entirety. The segmented image of the pore space was then directly translated to a lattice to simulate the transfer of fluids between the fracture and the matrix using lattice Boltzmann (LB) modeling. This provided us an opportunity to validate the modeling technique against experimental images at the pore-scale. Although the surfactant was found to alter the wettability of the originally weakly oil-wet glass to water-wet, the fracturematrix fluid transfer is found to be a drainage process, showing little to no counter-current migration of the oil-phase. The LB simulations were found to closely match experimental rates of fracture-matrix fluid transfer, equilibrium saturation, irreducible wetting phase saturation and fluid distributions.

106

4.2 INTRODUCTION

The displacement of fluids in fractured rock by transfer of fluids between the high-permeability low-storage fractures and low-permeability high-storage surrounding rock is a common process that occurs during hydrocarbon recovery, carbon dioxide injection, and groundwater contamination. During an injection into, or flooding of, a fractured rock, fluid migration is determined by the complex interaction of capillary and viscous forces that control flow in the fractures, as well as by the transfer of fluids to and from the matrix. Predicting this migration requires a rigorous investigation into the mechanics that determine the transfer of fluids to and from the matrix.

Understanding the physics of this process can be greatly enhanced by

investigating fluid transfer at the interface of the fracture and matrix at the pore-scale.

Experimental studies of fracture-matrix fluid transfer at the pore-scale have been limited, for the most part, to etched-glass micromodels and two-dimensional glass bead models. Micromodels are transparent, two-dimensional representations of pore networks (Wan et al. 1996,
Karadimitriou and Hassanizadeh 2012). Using these models, a wide variety of multiphase flow processes involving fracture-matrix transfer have been investigated in recent years including spontaneous imbibition into water-wet media (Rangel-German and Kovscek 2006, Hatiboglu and
Babdagli 2008, Hatiboglu and Babdagli 2011), CO2 flooding of naturally fractured oil-wet reservoirs at immiscible, near-miscible, and miscible pressure (Er et al. 2010), continuous gas, water and alternating gas-water injection in oil-wet naturally fractured reservoirs (Dehghan et al.
2012), and multiphase flow in multiple intersecting fractures (Buchgraber et al. 2012). These

107

studies have provided invaluable insight into the physics involved in multiphase fracture-matrix fluid transfer; however, their limitation to two-dimensions assumes flow in the missing third dimension will have minimal impact on fluid flow. In the investigation presented here, the third dimension is retained by utilizing X-ray computed microtomography (CMT), which is a noninvasive three-dimensional imaging technique. This allows us to image fracture-matrix fluid transfer of a three-dimensional porous media at the pore-scale. Also, images obtained of the pore space can be directly inputted into pore-scale lattice Boltzmann models for further investigation of the mechanics controlling flow.

X-ray computed microtomography (CMT) has become a commonly used method of imaging pore spaces and fluid distributions of unconsolidated sands and beadpacks (Culligan et al. 2006,
Brusseau et al. 2006, Brusseau et al. 2008, Constanza-Robinson et al. 2008, Al-Raoush et al.
2009, Lebedeva and Fogden 2011), rock (Coles et al. 1998, Turner et al. 2004, Prodanovic et al.
2007, Iglauer et al. 2011, Silin et al. 2011, Kumar et al. 2012) and fractures (Karpyn et al. 2007) at the pore-scale in three dimensions. Additionally, numerous investigations using CMT have imaged artificially induced (Ketcham et al. 2010) and naturally occurring fractures in rocks
(Keller et al. 1999, Renard et al. 2009, Wennberg et al. 2009). Due to advances in CMT imaging technology within the past two decades, it is now possible to image pore spaces and fluid distributions at the pore scale in three dimensions. However, carrying out experiments using
CMT is time-consuming and costly. Robust, pore-scale models offer the possibility of simulating a wide variety of scenarios in order to improve our predictive capabilities regarding fracturematrix fluid transfer, at the cost of only computational time. Here, we focus on the use of lattice
Boltzmann (LB) models, although numerous pore-scale modeling techniques currently exist

108

(Meakin and Tartakovsky 2009). Multiphase LB models have been validated against experimental measurements of capillary pressure in bead packs (Pan et al. 2004, Schaap et al.
2007, Porter et al. 2009) and sandstones (Ramstad et al. 2010, Ramstad et al. 2012), as well as experimental measurements of relative permeability in sphere packs (Ghassemi and Pak 2009,
Hao and Cheng 2010) and sandstones (Boek and Venturoli 2010, Ramstad et al. 2010).
Although there has been a fair amount of validation of LB methods via matching of experimental measurements of macroscale flow properties, few have compared fluid distributions imaged using CMT with LB model results. Sukop et al. 2008 compared two-phase fluid distributions in a sand pack imaged using CMT to distributions resulting from LB simulations. The saturation per slice was found to match experimental results well. Porter et al. 2009 compared interfacial fluid-fluid areas measured by CMT to LB displacement simulations, and found a good match for drainage simulations.

We are interested in investigating fracture-matrix fluid transfer at the pore-scale using experimental, sequential, three-dimensional CMT imaging of fluid transfer and lattice Boltzmann modeling. A synthetic core was manufactured in order to contain a minimal amount of pore space unresolvable by CMT imaging, and was subsequently fractured to serve as our fractured porous medium.

This fractured core is more closely analogous to a real rock than two-

dimensional micromodels, and it also has the previously stated advantage of allowing fluids to move in three directions. One of the advantages of using LB models in conjunction with CMT is the ease of translating images of the pore space into a simulation lattice. Here, we simply translate the CMT image of the pore space 1-to-1, image voxel to lattice node, thus retaining the complex solid boundary conditions of the pore space in the model. Also, by manufacturing the

109

core to have minimal pore space unresolvable by CMT, we can be confident that the modeled pore space is a close representation of the physical pore space observed during the experiment.

4.3 MATERIALS AND METHODS
4.3.1 Fractured Porous Medium

The fractured porous medium is a sintered glass granule pack, fractured using a modified
Brazilian test. Crushed glass was sieved for granules ranging from 115 to 210 μm in diameter.
These granules were placed in a 7.62 cm diameter, 10.2 cm long mullite tube and sintered following the temperature cycle summarized in table 4.1. Cores with a diameter of 1.27 cm were drilled from these sintered glass granule packs. A single, tensile fracture was induced along the axis of the core, and the two resulting halves of the core were glued along the edges of the fracture to ensure the two halves are immobile. The resulting fractured core is 3.12 cm in length and ~1.5 cm in diameter.

Table 4.1: Temperature cycles for sintering of glass granules.

Step
Increase
Hold
Increase
Hold
Decrease
Hold
Decrease

Time [h]
8
1.5
1
1
1
1.5
16

Cumulative Time
[h]
8
9.5
10.5
11.5
12.5
14
30

110

Temperature [C]
25 → 600
600
600 → 645
645
645 → 600
600
600 → 25

4.3.2 Procedure

The core was placed between the inlet and outlet of the core holder and encased in fuel-resistant
Teflon heat-shrink tubing to create an airtight seal at both the inlet and outlet of the core holder.
The core abuts the face of the inlet and outlet ports. This face was designed with crisscrossing channels in order to uniformly spread the injected/received fluid across the face of the inlet/outlet ports for even fluid delivery/recovery to the core. The core and core holder assembly was placed vertically in the CMT scanner under vacuum. The core was then scanned in its entirety to acquire a dry image of the pore space. Following the dry scan, the core was presaturated with the oil phase, 40% wt iodododecane in decane. The decane was doped with iodododecane in order to create contrast between phases in the CMT images. The fully oil-saturated core was then scanned in its entirety to acquire an image of the connected pore space. After scanning, the oil phase was displaced by injection of the water phase, 5%wt NaCl DI water. The interfacial tension between the water and oil phase was measured using the sessile drop method, and found to be 42.86 mN/m. The flow rates throughout the experiment were maintained at a rate of 0.1 to
1.0 ml/min.

The balance of viscous and capillary forces are commonly measured by the

capillary number,

⁄ , a dimensionless measure of viscous and capillary forces. The

displacements here have capillary numbers between 10 -5 and 10-7. All displacements in the experiment were considered to be capillary-dominated.

The CMT scanner is equipped to provide real-time radiographs (2D side-view images) of a 0.532 cm length of the core during water injection. During the water injection no water was observed entering the matrix. At the end of the initial water phase injection, the core was scanned in its

111

entirety, confirming that the water phase had only succeeded in displacing the oil phase from the fracture. Although glass is often considered a water-wet material, the presaturation with the oil phase and exposure to x-rays rendered the core weakly oil-wet. The water phase, being the nonwetting phase, occupies only the largest portions of the pore space necessary to connect the inlet with the outlet, i.e. the fracture. Due to the vertical orientation of the core and the difference in densities of the fluids, there is a difference in the pressure (

) of the heavier water-occupied

fracture and lighter oil-occupied matrix, equal to densities of water and oil,

, where

is the gravitational constant, and

and

are the

is the distance from the top or

inlet to the core. Because there are no viscous forces driving the water phase into the matrix, only capillary and gravitational forces are acting in the transfer of fluid to and from the matrix.
The capillary pressure ( ) is the difference in pressure between the non-wetting (

) and

wetting phase ( ). In our system the non-wetting phase occupies the fracture, while the wetting phase occupies the matrix, thus

. For the non-wetting phase to penetrate the matrix, the

difference in pressure due to gravity forces must overcome the capillary entry pressure (

).

The capillary entry pressure is the capillary pressure necessary for the non-wetting phase to invade a pore occupied by the wetting phase.

The capillary entry pressure is directly

proportional to the interfacial tension, and can be lowered by decreasing the interfacial tension.

In order to lower the interfacial tension of our fluid pair, 2 ml/l of dish soap was added to the water phase, lowering the interfacial tension to 0.30 mN/m. injected at the same rates as the water-phase injection.

The surfactant-water phase is

Once

the surfactant-water phase

injection is commenced the surfactant-water phase mixes with the water phase already occupying the fracture. Once the surfactant becomes established in the fractured porous medium, the

112

interfacial tension of the system will drop and the surfactant water-phase will begin to displace the oil phase (wetting phase) from the matrix. During the injection of the surfactant water phase, a 0.506 cm section of the core was sequentially scanned, producing 14 3D images of the water phase entering the matrix. This volume of the core is referred to as the region of interest (ROI) from this point on. Each scan takes approximately 10 minutes, with time in between scans ranging from 1 to 20 minutes. The injection is not stopped during the scans. The movement of the water phase into the matrix was slow enough that the moving interface did not appreciably blur our images. The sequential scanning of the water phase invasion was ended after 220 minutes, at which time the core was scanned in its entirety.

4.3.3 CMT Imaging

X-ray computed microtomography has become a popular technique to non-invasively image the pore spaces of porous media in 3D.

A CMT “scan” involves taking a sequence of 2-D

radiographs of a rotating sample and mathematically reconstructing them into a 3-D image composed of voxels. Each voxel is much like a pixel in a 2D image, and contains a CMT registration number, which is a relative measure of the x-ray attenuation of the material occupying a particular voxel.

Images reported here were taken at Pennsylvania State

University’s Center for Quantitative x-ray Imaging (CQI). Scans of the core in its entirety resulted in images 1024x1024x1980 voxel3, while sequential scans during the water phase invasion of the matrix were 1024x1024x321 voxel3.

Only a section of the core could be

monitored during transient imaging because imaging of a larger section would have come at a

113

loss to the overall number of images that could be acquired. All scans have a voxel resolution of
0.015762 x 0.015762 x 0.015767 mm3, with the latter representing height, i.e. slice thickness.

4.3.4 Image Processing

The CMT images must be further processed in order to determine meaningful information. First, all images are filtered using a 3x3x3 kernel median filter to remove salt and pepper type noise.
The CMT registration numbers of the three phases present, i.e. glass, oil phase, and water phase, are significantly different, allowing for identification of the phase occupying each voxel.
Segmentation here was accomplished through simple thresholding. Other, more complex thresholding techniques, such as bi-level thresholding with indicator kriging, were attempted, but did not appear to result in any visually apparent improvement in segmentation. In simple thresholding, two phases are segmented according to a threshold CMT value. Voxels with CMT registration numbers below the threshold are labeled as one phase and those with CMT registration numbers above the threshold are labeled as the other phase. To segment our three phase images, a two-step segmentation process was used. First, the glass and pore spaces were segmented in the presaturation image using a threshold determined by fitting a Gaussian curve to the peak in the CMT registration number frequency distribution corresponding to the glass. The threshold was taken as the CMT registration number corresponding to where the Gaussian curve has a frequency value below 1% of the peak value; all voxels with CMT registration numbers above this value are considered to belong to the pore space. This segmentation was then used to mask the three phase images. By subtracting out the glass, only two phases, the oil phase and water phase, are left.

These two phases were then segmented in the same manner as the

114

segmentation of the glass from the pore space. With the images subsequently segmented into glass, oil phase, and water phase, measurements of the pore space and fluid distribution could be made. However, before measurements associated with the matrix and fracture could be taken, the pore space had to be labeled as belonging to either the fracture or the matrix. The fracture and matrix pore spaces are connected, and the distinction between them is not obvious. To identify the pore space associated with the fracture and the matrix, a process of image erosion, isolation, and dilation, detailed in a previous study of ours (Landry and Karpyn 2012), was used. The matrix measurements are taken from the portion of the matrix parallel to the fracture and within 0.5 cm of the fracture face. After segmentation of the images and labeling of the pore space, volume measurements could be made from the image.

4.4 RESULTS
4.4.1 Pore Space Analysis

During the surfactant fluid penetration of the matrix, it became clear that the matrix on either side of the fracture was appreciably different, and that these differences resulted in significantly different rates of water phase invasion into the matrix. To measure the characteristics of the matrix separately, for what will be referred to from here on as the right and left sides, a line of division was determined from the center point of the fracture aperture measurements. In our images, the fracture is oriented perpendicular to image voxels, therefore we found it appropriate to measure the so-called “vertical aperture” ( ), i.e. the distance between the fracture walls. An

115

image of the fracture, and a corresponding mean fracture aperture profile, are shown in figure
4.1. The mean fracture aperture was measured for each slice. The ROI is also highlighted in figure 4.1.

B)

A)

Figure 4.1:

Image of the fracture with the ROI highlighted (A), and the corresponding mean fracture aperture and contact area profiles with the ROI indicated by arrows.

Overall, the fracture profile shows two wavelengths of decreasing and increasing mean aperture, with the two minimums corresponding to significant contact areas between the fracture walls.
The portion of the core to be imaged sequentially during the surfactant water phase injection of the matrix contains, at its center, a large “asperity”, meaning that the aperture is smaller than the average pore aperture and is thus indistinguishable from the matrix using the labeling method
116

described here. This portion of the core was chosen because it is located near the center and contains a portion of the fracture that has a significantly varying aperture and the occurrence of an asperity. One of the objectives of this study is to investigate the dependence of fracturematrix fluid transfer on the fracture aperture and the existence of asperities.

The porosity of the right and left matrix was calculated by summing the voxels of the segmented, labeled pore space images from the presaturation scan. Calculating the porosity for each slice, a profile of the porosity is presented in figure 4.2. The porosity of the core varies from ~0.1 to ~0.35, with the lowest porosities being found at the bottom and top of the core.
Also, the porosity of the left matrix is generally slightly greater than that of the right matrix as a result of the sintering process. The top, bottom, and sides, with the left matrix being closer to the center, of the sintered sample reached greater temperatures, resulting in greater compaction. The porosity of the left and right matrix of the portion of the core to be sequentially imaged is 0.3030 and 0.2698, respectively. Although the porosity gives a bulk measure of the volume of pore space present, and provides information in regards to the level of sintering compaction, the displacement of the wetting phase in the matrix by the non-wetting phase in the fracture is not dependent on the porosity. Instead, this displacement is strongly dependent on the pore aperture of the matrix. The porosity and the pore aperture are not necessarily related, it is both possible to have a high porosity, small pore aperture porous medium and a low porosity, large pore aperture porous medium. Another objective of this study was to design a synthetic 3D porous medium that could be imaged with minimum unresolvable pore space, thereby allowing very accurate and complete measurements of pore space geometry, such as the matrix pore aperture distribution.

117

A)

Figure 4.2:

B)

C)

Image of the matrix pore space showing the fracture saturated with water phase and the matrix saturated with the oil phase, the ROI volume is highlighted (A). Also the corresponding porosity (B) and mean aperture (C) profiles with the ROI indicated with arrows. A significant amount of literature regarding the description of pore space topology exists
(Thovert et al. 1993, Lindquist et al. 1996, Oren and Bakke 1997, Prodanovic et al. 2006).
Although there are numerous measurements that can be made regarding the pore space topology, we focused on the pore aperture distribution. In order to measure the pore aperture distribution of the pore space, first a skeleton, or medial axis, of the pore space must be constructed by thinning algorithms (Lee et al. 1994) using the software Avizo Fire 6.3. The skeleton of the pore space is a connected network of voxels that represents the midpoint between solid walls in the pore space. Each voxel in the skeleton can be thought of as the center of a sphere that makes contact with the wall in at least two places.

The pore aperture is the maximum radius of a
118

sphere that can occupy a pore. To determine the pore aperture distribution, a “sphere packing” of sequentially smaller spheres are strung along the skeleton. First, the voxels of skeleton are labeled with their distance to the pore wall. Next, the labeled skeleton is searched for voxels that fall within a user-defined range. These voxels then become the center of a sphere – with a radius equal to the distance to the closest pore wall – placed in the pore space. The algorithm treats all voxels that fall within this sphere as solids in order to avoid overlapping measurements of pore aperture between spheres from different aperture ranges.

An example image of the pore

apertures is shown in figure 4.3. The mean pore aperture is defined here as the volume-weighted mean pore aperture,
∑
where

represents all pore apertures measured,

is the pore aperture, and

is the volume

weight of each pore aperture measured,
∑
The weight of each pore aperture is equal to the volume ( ) of the associated sphere divided by the sum of all pore aperture sphere volumes. The mean pore aperture profile generally follows the trends seen in the porosity profile, with the largest pore apertures occurring in the middle of the core for both of the matrixes, and the left matrix pore apertures being generally greater than the right (figure 4.2). The fracture and pore aperture distributions for the ROI are given in figure
4.4. The fracture and pore aperture distributions are lognormal, similar to what is observed in
CMT images of Berea sandstone (Prodanovic et al. 2007) and fractures in Berea sandstone
(Karpyn et al. 2007). The mean pore aperture of the left and right matrix is 0.0770 and 0.0721
(mm), respectively.

119

Figure 4.3.

Example 0.1183 cm3 image of sphere-packing method used to measure pore aperture distribution. The individual pore apertures are measured as the radius of the spheres seen in this image. The red, purple, teal, green, yellow, and light green represent pore aperture measurements between
,

,

,

, and

, mm, respectively. The background shows the solid in black and the pore space in white.

120

Figure 4.4.

The fracture aperture distribution and the pore aperture distribution for the portion of the core to be imaged sequentially during water invasion of the matrix.

4.4.2 Immiscible Displacement

After the core is presaturated with the oil phase, the water phase without the addition of a surfactant is injected. The oil phase is nearly completely displaced from the fracture, with water phase saturations in the fracture,

, while the matrix has yet to be invaded by the

water phase. Only a small amount of the oil phase is retained in the fracture near the top of the core, and as is evident in the image presented in figure 4.5, the contact angle through the oil phase is quite large, close to neutral wetting conditions. The thorough displacement of the oil phase from the fracture by the water phase can be attributed to the weakly oil-wet condition. The affinity of the wetting phase for the solid is too weak for the wetting phase to be retained in the smaller apertures of the fracture. Also, as was previously mentioned, the non-wetting water phase does not invade the matrix during the initial water phase injection due to capillary forces, which can be decreased by the addition of a surfactant.

121

Figure 4.5.

Image of oil phase (white), water phase (black), glass (gray) contact in the fracture after initial water phase injection. The contact angle appears in this image to be large, near neutral wetting conditions.

At the beginning of the surfactant water phase injection, the non-wetting water phase is already established in the fracture. Thus, the fracture-matrix fluid transfer occurs in what has been previously referred to as the “instantly-filled” regime (Rangel-German and Kovscek, 2006), i.e. fracture-matrix fluid transfer does not occur until after the displacing fluid has been established in the fracture, as opposed to the “filling” regime in which fracture-matrix transfer occurs concurrently with the advancing displacement in the fracture. The regime of fracture-matrix fluid transfer is often determined by the velocity at which the displacing front advances in the fracture, with a faster moving front resulting in an “instantly-filled” fracture. The flow presented here is very slow; however, it represents an instantly-filled fracture-matrix transfer due to the delay in the injection of surfactants required for fluids to penetrate the matrix.
122

A few example images from the sequential scanning of the surfactant water phase injection of the matrix are shown in figure 4.6. During the sequential scanning there are some small amounts of the oil phase that appear in the fracture (figure 4.6). However, the total volume of the oil phase in the fracture is much smaller than the amount of the oil phase being displaced from the matrix suggesting there is very little, if any, movement of the oil phase from the matrix to the fracture. The displaced oil phase from the matrix migrates in the direction of flow through the matrix to the outlet. The small amount of oil found in the fracture appears to result from the snap-off of oil already in, or near, the fracture (figure 4.6). This suggests that the surfactant is also acting to alter the wettability of the solid to water-wet. The surfactant of the soap used is sodium dodecyl sulfate (SDDS), an anionic surfactant, which has been shown to decrease the water phase contact angle with quartz (Zdziennicka et al. 2009). We would expect that the functional surface sites of quartz would be similar to the glass used here, leading to the observed wettability alteration to water-wet.

123

Figure 4.6.

CMT images from the sequential scanning of the surfactant water phase injection showing a volume rendering (left) and a slice from a height of 1.0704 cm (right). In the volume renderings the water phase is dark orange, the solid is light orange and the oil phase is white. In the slices the water phase is light gray, the glass is dark gray, and the oil phase is white.

124

The effect this wettability alteration has on the transfer of fluids is discussed in the following pore-scale analysis. The non-wetting phase saturation of the matrix ( square root of time (

⁄

) as a function of the

) is shown in figure 4.7. At the start of the surfactant water phase

injection, the surfactant must first mix with the water phase already present, and establish a monolayer at the interface of the oil and water phase, before interfacial tension will be significantly reduced. As the interfacial tension is reduced, we observe an increasing rate of wetting-phase displacement, referred to in figure 4.7 as “early time”, which stabilizes at a near constant rate, , with rate being a linear function of

⁄

,

⁄

⁄

, where

is

the time the displacement is initiated. This stable displacement of brine-surfactant-Decane is the focus of this study, considering the majority of the displacement takes place after the surfactant has been established, and the interfacial tension significantly reduced.

The

⁄

for the

displacement is not a set number due to the slowly decreasing interfacial tension at early time.
Instead, we can determine an apparent start time for the stable displacement as determined by
⁄

fitting the rate function by linear regression. Fitting experimental measurements for

resulted in correlation coefficients of 0.9893 and 0.9801 for the right and left matrix, respectively
(figure 4.7). Including earlier times resulted in lower correlation coefficients, and excluding
⁄

did not result in significantly greater correlation coefficients. Therefore we can be

fairly confident that the stable displacement occurs in this time frame. The apparent

⁄

is 4.48

and 4.32 (min1/2) for the right and left matrix, respectively. Both sides are in agreement on the apparent ⁄

. The

is 0.0787 and 0.0551 (min-1/2) for the right and left matrix, respectively.

There is a significant difference in rate observed between the right and left matrix. Given that the only variable between the right and left matrix is the pore space geometry, the significance of
125

the effect matrix pore space geometry has on fracture-matrix fluid transfer is highlighted. To further understand the physics involved in this displacement, we compared the experimental pore-scale fluid distributions to the fluid distributions of a pore-scale LB model simulating a capillary-dominated drainage of the right and left matrix – separately but equally – in which the only variable in the simulation of the right and left matrices will be pore space geometry.

Early
Time

Figure 4.7.

Stable
Displacement

Non-wetting phase saturation of the right and left matrix as determined from sequential
CMT imaging. Early time refers to displacement before the interfacial tension reduction has stabilized. The linear fits shown represent a displacement following the rate function,
⁄

.

126

4.4.3 LB Simulation of Fracture-Matrix Fluid Transfer

For details of the LB model used refer to chapter 3. To model the surfactant water phase movement from the fracture to the matrix a 270x150x150 voxel3 (0.4256 x 0.2365 x 0.2365 cm3) sample volume (figure 4.8) of the pore space from each matrix is translated 1-to-1, voxel to lattice node, for lattice Boltzmann simulation. This translation maintains the complexity of the pore space solid boundary conditions – the image is directly translated into fluid and bounceback nodes. The size of the image volume to be simulated is limited by the prohibitive computational requirements of the LB model. Simulating the entire matrix volume scanned would be optimal, although, it should be noted, that the sample volume used is larger than the representative elementary volume determined by porosity and specific solid surface area measurements not shown here. The simulation lattices are bound by solid walls orthogonal to the direction of flow, and pressure boundaries 5 lattice layers thick at the ends. These pressure boundaries represent the non-wetting phase occupied fracture, and the pushback of the capillary force of the wetting phase occupied matrix. Images of the sample pore space used in LB simulations, and an illustration of the LB simulation set up is given figure 4.8.

127

B)

A)

Direction of Displacement

Matrix:
Initially occupied by wetting phase

Wetting phase
Non-wetting
phase Pressure Pressure Boundary
(Matrix)
Boundary
(Fracture)

Figure 4.8:

Volume images of the 270 x 150 x 150 voxel3 (0.4256 x 0.2365 x 0.2365 cm3) sample volumes of the pore space for the left and right matrix to be translated 1-to-1, voxel to lattice node, for lattice Boltzmann simulation (A). LB simulation setup showing the 5 lattice layer thick pressure boundaries, with the fracture represented as the non-wetting phase occupied pressure boundary (B).

The wettability of the core was observed previously to be very weak with a contact angle in the range of 75° to 90°. Using the equation of Huang et al. 2007 given above, a total component density,

, and a

, this contact angle can be approximated by using
.

Setting the pressure boundaries is less exact; experience and

numerical experiments were performed to determine the pressure difference (

) required for

the non-wetting phase to traverse the simulation volume, resulting in a general estimation of
0.033 (mulu-1 tu-2). Typically, the capillary pressure of the simulation is scaled ( using the physical and lattice interfacial tensions in accordance with Laplace’s law,

128

⁄

)

where

is the length scaling equal to the resolution of the image, and,

and

are the

physical and lattice interfacial tensions respectively. A lattice pressure difference of 0.033
(mulu-1 tu-2) represents a physical pressure difference of 6.33 Pa, and is less than the physical pressure difference experienced by the fluids in the portion of the core sequentially scanned.
This suggests the process of fracture-matrix fluid transfer is indeed a drainage process and not the result of wettability alteration followed by imbibition.

Had the simulation required a

pressure difference significantly greater than the pressure difference experienced in the experiment, it would suggest the drop in the interfacial tension alone could not account for the initialization of fracture-matrix fluid transfer.

The simulation runs until the convergence

criterion ( ), defined here as,

is met. Based on numerical experiments not shown here, a convergence criteria of 0.001 was chosen. The LB simulations each required near

iterations (

) to

reach equilibrium, completing the displacement.

⁄

Time scaling of LB simulations (

) is vague, typically dimensionless

measurements of the forces involved (i.e. capillary number, Reynolds number) are matched to those used in the simulation to determine a reasonable value of the time scaling (Sukop and
Thorne 2006) or chosen, as described by Hatiboglu and Babadagli 2008, “carefully” to match scaling of the interfacial tension, viscosity and density between the LB and experimental system.
Time-scaling in the manner described by Hatiboglu and Babadagli 2008 resulted here in a gross underestimation of the time-scaling on the order of 103 - 104. The displacements described here

129

occurred over a period of a few hours, whereas those described by Hatiboglu and Babadagli
2008 occurred in a matter of seconds. The LB models employed are similar, and thus the timescaling is of a similar magnitude, which was shown to be appropriate for the displacements described in Hatiboglu and Babadagli 2008, but not those presented here. We determined the time-scaling by matching the physical time of stable displacement in the experimental system to the lattice time required for the simulation to reach equilibrium. Thus, the time-scaling is determined by the measured time of displacement in the experiment, and not derived from the simulation alone. This highlights one of the weaknesses of these types of LB simulations; alone they cannot predict the rate of fracture-matrix fluid transfer for slow displacements. Unlike the physical fracture-matrix transfer, there is no delay in the start of the displacement in the LB simulations. To compensate for this, the LB simulation results are shifted by an amount equal to the mean

of the right and left matrix. The stable physical displacement required is roughly to ⁄

reach

equilibrium,

giving

a

time-scaling

of

.

4.4.4 Pore-Scale Experimental and Model Results

One of the main objectives of this study was to evaluate the accuracy with which the LB simulation captures fluid distribution at the pore-scale. As was mentioned in the Introduction, there are very few investigations that have compared LB simulation results with experimental results at the pore-scale. Also, due to the employment of a surfactant that appears to alter the wettability of the solid surface, we were interested in determining if the fracture-matrix fluid transfer can be described by a drainage process, and if the differences observed in the fracture-

130

matrix transfer rate ( ) of the right and left matrix can be attributed to differences in the pore space geometry. A comparison of

as a function of

for the experimental and LB results

is presented in figure 4.9.

Figure 4.9.

Non-wetting phase saturation of the right and left matrix as determined from sequential
CMT imaging and LB drainage simulations. The only variable between the LB simulations for the right and left matrix is the pore space geometry.

The LB simulation of the right matrix fits the experimental results very well; however, the simulation of the left matrix underestimates the rate of fracture-matrix fluid transfer and, ultimately, the equilibrium non-wetting phase saturation. This could be attributed to simulating too low of a pressure difference; at higher capillary pressures the non-wetting phase saturation is expected to be greater.

However, the right matrix simulation correctly reproduced the

131

equilibrium saturation. It is also possible that the size of the pore space translated into lattices for LB simulation was too small. Returning to figure 4.6, we observe the existence of elongated pore spaces with lengths on the order of 0.10 – 0.15 cm or 63 - 95 voxels ( ). It is possible that our simulations, which have side lengths of 150 voxels (

are too small to effectively access

these larger pore spaces during the simulated displacement, we further comment on this in the next section. Determination of appropriate lattice sizes and resolutions for LB simulations continue to be subjects of discussion, and with increasing computational accessibility, these issues can be addressed.

Although the LB model underestimates the fracture–matrix fluid

transfer for the left matrix, it does correctly predict the higher

of the left matrix over the right

matrix for the first half of the displacement. This suggests that the stable displacement of the wetting phase from the matrix by the non-wetting phase occupied fracture as a result of the introduction of an interfacial tension reducing and wettability altering surfactant can be described by a drainage process, without the consideration of wettability alteration. This is not surprising considering a surface cannot be altered until the non-wetting phase fluid makes contact with the surface, and this surface contact is controlled by a drainage-like displacement.

132

Experimental

Model

A)

B)

C)

D)

E)

F)

G)

H)

133

Figure 4.10.

Comparison of experimental and simulated vertical saturation profiles in the matrix (left and right) adjacent to the fracture for the right (A, B) and left (C, D) matrix, and horizontal saturation profiles for the right (E, F) and left (G, H) matrix.

To further investigate how effective the LB simulations were at replicating the experimentally observed pore-scale fluid distributions we compare vertical and horizontal saturation profiles.
The vertical (saturation as a function of height) and the horizontal (saturation as a function of distance from fracture) for the experimental and LB model are presented in figure 4.10.
Experimental and LB results are reported in dimensionless height (
(

⁄

) and length

) units, where the reference height and length are the respective height and lengths

of the analyzed CMT image (
(

⁄

) and LB simulation volume
). We also show some snapshots from the LB simulation

for the right and left matrix in figure 4.11.

134

Left Matrix

Figure 4.11:

Right Matrix

Sample images of the LB simulations of fracture-matrix transfer. The non-wetting phase saturation front and simulation pore space are shown, the wetting phase is occupies the pore space not occupied by the non-wetting phase.

135

The experimental vertical saturation profiles (figure 4.10A, C) do not display any dependence on the fracture aperture, even with the presence of a sizable asperity. However, if the asperity were larger we would expect the matrix adjacent to the asperity to experience some delay in the arrival of the saturation front from nearby open portions of the fracture. Similar uniformity of the vertical saturation profiles is seen in the LB simulations (figure 4.10B, D), but due to small size of the LB simulation volumes the saturation front is less uniform, a result of piston-like displacement occurring in individual pores. We can be fairly confident the capillary-dominated fracture-matrix fluid transfer presented here has very little dependence on the fracture aperture, given the fracture is occupied by the displacing phase.

The LB simulation for the left matrix (figure 4.10D) has low non-wetting phase saturations near the bottom of the simulation volume. The blockage of the saturation front in the lower portion of the simulation volume is a result of its small size. This becomes more apparent in the images of the LB simulations shown in figure 4.11.

As was previously stated, LB simulations are

computationally demanding, limiting the size of the simulation volume. The blockage of the saturation front indicates that the volume simulated may be too small, and could contribute to the underestimation of the fracture-matrix fluid transfer.

Furthermore, the higher rate of

displacement in the left matrix was captured by the LB simulation at first, until the lack of accessibility to the lower portion of the simulation volume slowed the displacement rate. The experimental horizontal saturation profiles (figure 4.10E, G) show a saturation front that is 0.30
– 0.40 cm in length, with an irreducible wetting phase saturation near 0.90-0.95. The LB model horizontal saturation profiles (figure 4.10F, H) roughly match the saturation front lengths of the experimental measurements, but more importantly they accurately predict an irreducible wetting

136

phase saturation near 0.90-0.95. Also we see here that the non-wetting phase saturation front for the left matrix reached the edge of the sample. In the LB simulation of the left matrix only the top part of saturation front transited the simulation volume, causing the overall vertical saturation profiles to underestimate experimental results.

Again the mismatch can be attributed to

simulating to small of a volume. Considering this, we do not believe any mismatches between the model and experimental results should be attributed to shortcomings of the LB method, but rather are simply a result of a rough estimation of boundary conditions. Also issues involving the size of the sample volume used to simulate fracture-matrix fluid transfer are a result of our limited computational capacity.

4.5 CONCLUSIONS

Here we have presented a pore-scale investigation of fracture-matrix fluid transfer. A surfactant water phase displacement of an oil-saturated oil-wet fractured synthetic porous medium was sequentially imaged using 3D CMT imaging. The synthetic porous medium was designed to have minimal pore space beyond the resolution of the CMT imaging used, allowing us to observe the motion of the displacement front at the pore-scale. The displacement occurred over a period of approximately 3 hours allowing us to image a 0.5043 cm length of the core 14 times to produce 4D data of fracture-matrix fluid transfer. Although the displacement was not stopped during sequential imaging, the displacement was slow enough that there was no appreciable blurring of the 3D images. A surfactant was introduced to lower the interfacial tension between the water and oil phase, and consequently the capillary entry pressure, to initiate fracture-matrix fluid transfer.

Images showed the surfactant altered the wettability of the porous medium

137

causing small amounts of the wetting phase in and near the fracture to snap-off and be transported in the fracture. But, almost all of the non-wetting phase was displaced out the top of core. There was very little, if any, counter-current migration of the wetting phase.

The

surfactant was introduced into a fracture with the non-wetting water phase already established, thus there was a delay in the stabilization of the rate of fracture-matrix fluid transfer as the surfactant established itself at the oil-water interface. After this delay the displacement stabilized and the non-wetting phase saturation of the right and left matrix correlated linearly with the square-root of time. Although the differences in the porosity and pore aperture distribution of the right and left matrix were small, the fracture-matrix transfer rate of the left matrix was significantly greater than the right matrix.

This emphasizes the importance of pore space

geometry, as this is the only variable that exists between the right and left matrix. Although the matrix pore space geometry is an important factor determining the migration of the displacing front, the movement of the displacing front showed very little dependence on the fracture aperture, and the existence of an “asperity”.

One of the objectives of this investigation was to evaluate the use of a Shan-Chen type LB models by making pore-scale comparisons of fluid distributions imaged using CMT and model results, which are lacking in the literature. We present LB drainage simulations of a sample volume of the left and right matrix to simulate the experimentally observed fracture-matrix fluid transfer. Typical methods of scaling time in accordance with the interfacial tension results in a
3-4 order of magnitude underestimate of the time scaling.

Here we scale time to match

experimentally observed physical time. This highlights one of the weaknesses of these types of
LB simulations; alone they cannot predict the rate of fracture-matrix fluid transfer for slow

138

displacements.

The LB drainage simulation of the right matrix closely matched the

experimentally observed equilibrium saturation, displacement rate, irreducible wetting phase saturation, as well as the vertical and horizontal saturation profiles. The LB drainage simulation of the left matrix closely matched the experimentally observed displacement rate for the first half of the displacement and the irreducible wetting phase saturation. However, the simulation LB simulation of the left matrix underestimated the equilibrium saturation, which was reflected in an underestimation of the displacement rate for the second half of the displacement and a pore-scale fluid distribution match that was not as convincing as the simulation of the right matrix. This underestimation is likely a result of using too small of a simulation volume, which limited access to some of the larger elongated pores of the left matrix. We do not believe any mismatches between the model and experimental results should be attributed to shortcomings of the LB method. Due to the prohibitive computational demands of multicomponent LB models the sample volumes of the image used in the LB model was significantly smaller than the ROI, considering the simulations only contained on the order of 100 grains they replicated experimental results fairly well.

This successful replication also implies that although the

surfactant altered the wettability of the porous medium to water-wet the fracture-matrix fluid transfer can be described as a drainage process.

139

4.6 REFERENCES

Al-Raoush R.I. Impact of wettability on pore-scale characteristics of residual nonaqueous phase liquids. Environmental Science and Technology 2009; 43:4796-4801. doi:10.1021/es802566s

Bakke S., and Oren P.E. 3-D pore-scale modeling of sandstones and flow simulations in the pore networks. SPE Journal 1997; 2(2):136-149. doi:10.2118/35479-PA

Boek E.S., and Venturoli M. Lattice-Boltzmann studies of fluid flow in porous media with realistic rock geometries. Computers and Mathematics with Applications 2010; 59:2305-2314. doi:10.1016/j.camwa.2009.08.063 Brusseau M.L., Peng S., Schnaar G., and Constanza-Robinson M.S. Relationships among airwater interfacial area, capillary pressure, and water saturation for a sandy porous medium.
Water Resources Research 2006; 42:W03501. doi:10.1029/2005WR004058

Brusseau M.L., Janousek H., Murao A., and Schnaar G. Synchrotron X-ray microtomography and interfacial partitioning tracer test measurements of NAPL-water interfacial areas. Water
Resources Research 2008; 44:W01411. doi:10.1029/2006WR005517

140

Buchgraber M., Al-Dossary M., and Ross C.M., and Kovscek A.R. Creation of a dual-porosity micromodel for pore-level visualization of multiphase flow. Journal of Petroleum Science and
Engineering; 86-87:27-38.

Coles M.E., Hazlett R.D., Spanne P., Soll W.E., Muegge E.L., and Jones K.W. Pore level imaging of fluid transport using synchrotron X-ray microtomography. Journal of Petroleum
Science and Engineering 1998; 19:55-63.

Dehghan A.A., Ghorbanizadeh S., and Ayatollahi Sh. Investigating the fracture network effects on sweep efficiency during WAG injection process. Transport in Porous Media 2012; 93:577595. doi: 10.1007/s11242-012-9970-7.

Er V., Babadagli T., and Xu Z. Pore-scale investigation of the matrix-fracture interaction during
CO2 injection in naturally fractured oil reservoirs. Energy Fuels 2010; 24:1421-1430. doi:
10.1021/ef901038v.

Ghassemi A., and Pak A. Numerical study of factors influencing relative permeabilities of two immiscible fluids flowing through porous media using lattice Boltzmann method. Journal of
Petroleum Science and Engineering 2011; 77:135-145. doi:10.1016/j.petrol.2011.02.007

141

Hao L., and Cheng P. Pore-scale simulations on relative permeabilities of porous media by lattice Boltzmann method. International Journal of Heat and Mass Transfer 2010; 53:1908-1913. doi:10.1016/j.ijheatmasstransfer.2009.12.066 Hatiboglu C.U. and Babadagli T. Pore-scale studies of spontaneous imbibition into oil-saturated porous media. Physical Review E 2008; 77:066311. doi: 10.1103/PhysRevE.77.066311.

Hatiboglu C.U. and Babadagli T. Experiment and visual analysis of co- and counter-current spontaneous imbibition for different viscosity ratios, interfacial tensions, and wettabilities.
Journal of Petroleum Science and Technology 2011; 70:214-228. doi:
10.1016/j.petrol2009.11.013.

Huang H., Thorne D.T., Schaap M.G., and Sukop M.C. Proposed approximation for contact angles in Shan-and-Chen-type multicomponent multiphase lattice Boltzmann models. Physical
Review E 2007; 76:066701. doi: 10.1103/PhysRevE.76.066701

Iglauer S., Paluszny A., Pentland C.H., and Blunt M.J.

Residual CO 2 imaged with X-ray

microtomography.

Letters

Geophysical

Research

2011;

38:L21403.

doi:10.1029/2011GL049680

Karadimitriou N.K. and Hassanizadeh S.M. A review of micromodels and their use in two-phase flow studies. Vadose Zone Journal 2012; 11(3): doi: 10.2136/vzj2011.0072.

142

Karpyn Z.T., Grader A.S., and Halleck P.M. Visualization of fluid occupancy in a rough fracture using microtomography. Journal of Colloid Interface Science 2007; 307(1):181-187. doi:
10.1016/j.jcis.2006.10.082

Keller, A. A., Roberts, P. V. and Blunt, M. J. Effect of fracture aperture variations on the dispersion of contaminants. Water Resources Research 1999; 35(1):55-63.

Ketcham R.A., Slottke D.T., and Sharp Jr. J.M. Three-dimensional measurement of fractures in heterogeneous materials using high-resolution X-ray computed tomography. Geosphere
6(5):499-514. doi: 10.1130/GES00552.1

Kumar M., Fogden A., Senden T., and Knackstedt M.

Investigation of pore-scale missed

wettability. SPE Journal 2012; 17(1):20-30. doi: 10.2118/129974-PA

Landry C.J., Karpyn Z.T., and Piri M. Pore-scale analysis of trapped immiscible fluid structures and fluid interfacial areas in oil-wet and water-wet bead packs. Geofluids 2011; 11:209-227. doi: 10.1111/j.1468-8123.2011.00333.x

Lebedeva E.V., and Fogden A.

Micro-CT and wettability analysis of oil recovery from

sandpacks and the effect of waterflood salinity and kaolinite. Energy and Fuels 2011; 25:56835694. doi:10.1021/ef201242s

143

Lee T.C., Kashyap R.L. and Chu C.N. Building skeleton models via 3D medial surface/axis thinning algorithms. CVGIP, Graphical Models and Image Processing 1994; 54:462-478.

Lindquist B., Lee S.M., Jones K.W., and Spanne P. Medial axis analysis of void structure in three-dimensional tomographic images of porous media. Journal of Geophysical Research 1996;
101(B4):8297-8310.

Martys N.S., and Chen H. Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method. Physical Review E 1996; 53(1):743-750.

Meakin P. and Tartakovsky A.M. Modeling and simulation of pore-scale multiphase fluid flow and reactive transport in fractured and porous media. Reviews of Geophysics 2009; 47:RG3002. doi: 10.1029/2008RG000263

Pan C., Hilpert M., and Miller C.T. Lattice-Boltzmann simulation of two-phase flow in porous media. Water Resources Research 2004; 40:W01501. doi:10.1029/2003WR002120

Porter M.L., Schaap M.G., and Wildenschild D. Lattice-boltzmann simulations of the capillary pressure-saturation-interfacial area relationship for porous media. Advances in Water
Resources; 32:1632-1640. doi:10.1016/j.advwatres.2009.08.009

144

Prodanovic M., Lindquist W.B., and Seright R.S. Porous structure and fluid partitioning in polyethylene cores from 3D x-ray microtomographic imaging. Journal of Colloid and Interface
Science 2006; 298:282-297. doi:10.1016/j.jcis.2005.11.053

Prodanovic M., Lindquist W.B., and Seright R.S. 3D image-based characterization of fluid displacement in a Berea core.

Advances in Water Resources 2007; 30:214-226.

doi:10.1016/j.advwatres.2005.05.015

Qian Y.H., d’Humieres D., and Lallemand P. Lattice BGK Models for Navier-Stokes Equations.
Europhysics Letters 1992; 17(6):479-484. doi:10.1209/0295-5075/17/6/001

Ramstad T., Oren P.E., and Bakke S. Simulation of two-phase flow in reservoir rocks using a lattice Boltzmann method. SPE Journal 2010; 15(4):917-927. doi: 10.2118/124617-PA

Ramstad T., Idowu N., Nardi C., and Oren P.E. Relative permeability calculations from twophase flow simulations directly on digital images of porous rock. Transport in Porous Media
2012; 94:487-504. doi:10.1007/s11242-011-9877-8

Rangel-German E.R. and Kovscek A.R. fracture transfer mechanisms.

A micromodel investigation of two-phase matrix-

Water Resources Research 2006; 43:W03401. doi:

10.1029/2004WR003918.

145

Renard F., Bernard D., Desrues J., and Ougier-Simonin A. 3D imaging of fracture propagation using synchrotron X-ray microtomography. Earth and Planetary Science Letters 2009: 286:285291. doi: 10.1016/j.epsl.2009.06.040

Schaap M.G., Porter M.L., Christensen B.S.B, and Wildenschild D. Comparison of pressuresaturation characteristics derived from computed tomography and lattice Boltzmann simulations.
Water Resources Research; 43:W12S06. doi:10.1029/2003WR002120

Shan X., and Doolen G. Multicomponent lattice-Boltzmann model with interparticle interaction.
Journal of Statistical Physics 1995; 81:379-393.

Silin D., Tomatsu L., Benson S.M., and Patzek T.W. Microtomography and pore-scale modeling of two-phase

fluid

distribution.

Transport

in Porous

Media

2011;

86:495-515.

doi:10.1007/s11242-010-9636-2

Sukop, M. and Thorne, D. T. (2006) Lattice Boltzmann Modeling an Introduction for
Geoscientists and Engineers, Springer-Verlag.

Sukop M.C., Huang H., Chen L.L., Deo M.D., Oh K., and Miller J.D. Distribution of multiphase fluids in porous media: Comparison between lattice Boltzmann modeling and micro-x-ray tomography. Physical Review E 2008; 77:026710. doi: 10.1103/PhysRevE.77.026710

146

Thovert J.F., Salles J., and Adler P.M. Computerized characterization of the geometry of real porous media: their discretization, analysis and interpretation. Journal of Microscopy 1993;
170(1):65-79.

Turner M.L., Knufing L., Arns C.H., Sakellariou A., Senden T.J., Sheppard A.P., Sok R.M.,
Limaye A., Pinczewski W.V., and Knackstedt M.A. Three-dimensional imaging of multiphase flow in porous media. Physica A 2004; 339:166-172. doi:10.1016/j.physa.2004.03.059

Wan J., Tokunaga T.K., Tsang C.F., and Gudmundur S. Bodvarsson. Improved glass micromodel methods for studies of flow and transport in fractured porous media. Water
Resources Research 1996; 32(7):1955-1964.

Wennberg, O. P., Rennan, L. and Basquet, R. Computed tomography scan imaging of natural open fractures in a porous rock; geometry and fluid flow. Geophysical Prospecting 2009;
57(2):239-249. doi: 10.1111/j.1365-2478.2009.00784.x

Zdziennicka A., Szymczyk K., and Janzuk B. Correlations between surface free energy of quartz and its wettability by aqueous solutions of nonionic, anionic and cationic surfactants. Journal of
Colloid and Interface Science 2009; 340:243-248. doi :10.1016/j.jcis.2009.08.040

147

CHAPTER 5: SUMMARY OF FINDINGS AND RECOMMENDATIONS FOR FUTURE
WORK

In these investigations we have highlighted some of the strengths of the LB model, and its ability to elucidate the relationship between pore space characteristics and macroscale two-phase fluid flow. Our comparisons of pore-scale fluid distributions between the model and experiment are promising. There are however many questions that can be addressed by future investigations.

In the first investigation we found that the presence of porous walls can greatly enhance the permeability of the fracture.

We also found that fracture permeability estimation using

previously proposed formulations dependent on fracture characteristics did a generally poor job of predicting the permeability of the fracture simulated.

The determination of fracture

permeability is greatly simplified when the existence of the porous walls has an insignificant effect. We can expect that as the ratio of the fracture aperture to pore aperture increases the presence of the porous walls will have less of an effect, until at a critical ratio it will have no effect. The question then arises what is this critical ratio? One should also consider that there are many different types of fractures that exist across vast length scales (10 -6 to 104 m). At what length scales and what type of fractures does the presence of a porous wall significantly affect flow? It would be optimal if a formulation could be determined that estimated fracture

permeability from the fracture characteristics, but the results seen here are not promising.
However considering there are many types of fractures, does there exist a “class” of fractures for
148

which the estimation of fracture permeability from the fracture characteristics is practical. We can speculate that fractures with a large fracture aperture to pore aperture ratio would fall into this “class”. Future investigations that collected a variety of fracture images could simulate flow using LB and begin to find the limits of fracture permeability estimation from fracture characteristics. Although it should be noted that it is unlikely that a robust fracture permeability estimation from fracture characteristics could be found for fractures like the one presented here that a small fracture aperture to pore aperture ratio – the solid boundaries of porous matrixfracture system a too complex to be distilled into a simple formulation.

In the second investigation we found that the relative permeability of both phases in homogenous-wet porous media decreased with increasing wetting strength.

It would be

interesting to compare this to experimental measurements of relative permeability. There are many techniques that can be employed to change the wettability of beads, performing conventional relative permeability measurements of homogenous-wet bead packs would help to further connect the results of the second investigation to experimental measurements.

In the third investigation we found that LB simulation matched the experimental results for the right matrix well, but not the left matrix. This was attributed to size of sample image translated into a LB lattice being likely too small. However, it is also possible that the issue was the underestimation of the pressure difference imposed in the simulations. It would be easy to determine which of these were true by simulating the fracture-matrix fluid transfer on a larger sample image. But what would be of greater interest to the community is a generalized study of lattice sizes and resolutions used in LB simulations. As in, what is the representative elementary

149

volume required to properly simulate flow? One way to answer this question would be to use a pore space image of a semi-homogenous medium, such as a bead pack, and simulate increasingly larger volumes until the macroscale flow results stabilize.

150

VITA:

Christopher James Landry

EDUCATION:
Doctor of Philosophy, Energy and Mineral Engineering
The Pennsylvania State University, University Park, Pennsylvania
May 2013
Dissertation Title: Pore-scale imaging and lattice Boltzmann modeling of single- and multiphase flow in fractured and mixed-wet permeable media.
Advisor: Dr. Zuleima Karpyn
Master of Science, Petroleum and Natural Gas Engineering
The Pennsylvania State University, University Park, Pennsylvania
December 2009
Thesis Title: Experimental pore-scale analysis of fluid interfacial areas in oil-wet and water-wet bead packs.
Advisor: Dr. Zuleima Karpyn
Bachelor of Science, Geology & Physics
Western Michigan University, Kalamazoo, Michigan
December 2006
PEER-REVIEWED PUBLICATIONS:
1) Landry C.J. and Karpyn Z. (In preparation) “Lattice Boltzmann modeling and 4D x-ray computed microtomogrpahy imaging of fracture-matrix fluid transfer”.
2) Landry C.J. and Karpyn Z. (In preparation) “Relative permeability of homogenous-wet and mixed-wet porous media as determined by pore-scale lattice Boltzmann modeling”.
3) Landry C.J. and Karpyn Z. (2012) “Single-phase lattice Boltzmann simulations of pore-scale flow in fractured permeable media” International Journal of Oil, Gas and Coal
Technology, 5(2/3), 182-206.
4) Landry C.J., Karpyn Z. & Piri M. (2011) “Pore-scale analysis of trapped immiscible fluid structures and fluid interfacial areas in oil-wet and water-wet bead packs” Geofluids,
11(2), 209-227.
5) Landry C.J., Koretsky C., Das S. & Lund T. (2009) “Surface complexation modeling of Co(II) on mixtures of hydrous ferric oxide, quartz and kaolinite” Geochimica et Cosmochica
Acta, 73 , 3723-3737
AWARDS:
1) “Second Place Society of Petroleum Engineers International Student Paper Contest – Master’s
Division” Society of Petroleum Engineers Annual Technical Conference, September
2010, Florence, Italy.
2) “First Place Society of Petroleum Engineers Regional Student Paper Contest – Master’s
Division” Eastern/Mid-Continent/Rockies SPE Regional Paper Contest, March 2010,
Pennsylvania State University, State College PA, U.S.A.
3) “Second Runner-up Graduate Student Paper Contest” Canadian International Petroleum
Conference, June 2009, Calgary Canada.
4) “2008 Marathon Alumni Centennial Graduate Fellowship” College of Earth and Mineral
Sciences, October 2008, Pennsylvania State University, State College PA, U.S.A.
5) “Outstanding Student Paper in the Division of Geochemistry” American Chemistry Society,
2006 National Meeting, Adsorption of metals to geomedia symposium, Atlanta GA
U.S.A.
6) “2005-2006 Arts & Sciences Research and Creative Activities Award” College of Arts &
Sciences, Western Michigan University, Kalamazoo MI U.S.A.