# Scale-up of reactive flow through network flow modeling

Table of Contents

LIST OF FIGURES ........................................................................................................vii LIST OF TABLES ……………………………………………….……………………...xi I.

INTRODUCTION ................................................................................................... 1

II.

NUMERICAL SCHEME ....................................................................................... 6

1.

3D image analysis .................................................................................................... 6

a.

Segmentation ......................................................................................................... 6

b.

Medial axis ............................................................................................................ 8

c.

Throat computation ............................................................................................... 9

d.

Pore-throat network construction ........................................................................ 10

2.

Pore-throat network of mineral distributions ..................................................... 11

a.

Viking sample ...................................................................................................... 11

b.

BSE and EDX analyses ....................................................................................... 12

c.

Mineral distributions of pore-throat network ...................................................... 14

3.

Single-phase network flow model ........................................................................ 19

a.

Construction of the single-phase network flow model ........................................ 19

b.

C ij based on channel shape factor ........................................................................ 21

4.

Reactive transport model ..................................................................................... 23

a.

Kinetic reactions .................................................................................................. 25

b.

Equilibrium computation of instantaneous reactions .......................................... 28

c.

Initial and boundary conditions ........................................................................... 36

d.

CO 2 solubility in aqueous NaCl solution ............................................................ 36

e.

Activity coefficients and diffusion coefficients ................................................... 39

f.

Effective and volume-averaged reaction rates [4] ............................................... 41

vi

III.

COMPUTATIONAL RESULTS .......................................................................... 43

1.

Mineral distributions ............................................................................................ 43

2.

Simulation of reactive flows ................................................................................. 60

a.

Validation of the reactive models ........................................................................ 60

b.

Application to Viking samples ............................................................................ 68

IV.

CONCLUSION & FUTURE STUDY ................................................................. 82

BIBLIOGRAPHY ........................................................................................................... 84

vii

List of Figures

Fig. 1. (a) Conceptual network model in [4] and (b) conceptual 2D pore space of a rock in [8] used for pore scale simulation of reactive flow. ..................................................... 3 Fig. 2. A schematic diagram of the simulation of the reactive flow. CO 2 dissolved water is injected from the top boundary into porous rock. The flow is similar to the case A in [4]. 5 Fig. 3. (a) The grey scale BSE image of Viking3W4. (b) The segmented image of the grey scale BSE image. Black, green, grey, and red indicate void phase, minerals of MAN < MAN quartz , quartz, and minerals of MAN > MAN quartz . ........................................... 13 Fig. 4. One pore body contains kaolinite (green) inside and it neighbors two clusters of MoI (red). Black and grey colors represent void space and quartz space respectively. S denotes surface areas between two minerals or mineral and void space. V denotes volumes of mineral clusters. ................................................................................... 17 Fig. 5. (a) A 2D example of a constructed pore-throat network. (b) Schematic of network flow model developed from the pore-throat network in (a). ........................................ 19 Fig. 6. f pH of anorthite and kaolinite as a function of pH [4]. ...................................... 27 Fig. 7. One time step in the simulation..................................................................... 28 Fig. 8. Histograms of the attenuation coefficients of three Viking samples with the threshold windows. ............................................................................................... 44 Fig. 9. (a) One slice of grey scale image of the conglomerate sandstone (Viking14W5). (b) Bi-phase segmented image of the primary pore and primary grain spaces. Black and white indicate the primary grain and primary pore respectively. (c) Segmented image of void and kaolinite in the primary pore space. Black and green indicate void and kaolinite respectively. Grey indicates the primary grain space. (d) Segmented image of quartz and

viii

MoI. Grey and red indicate quartz and MoI respectively. White indicates the primary pore space. (e) The final four-phase segmented image. Black, green, grey, and red indicate void, kaolinite, quartz, and MoI respectively. ........................................................... 46 Fig. 10. Grey scale and four-phase segmented images of the sandstone (Viking3W4). .. 47 Fig. 11. Grey scale and four-phase segmented images of the shaly sandstone (Viking10W4). ...................................................................................................... 47 Fig. 12. Distributions of (a) pore volume, (b) throat area, (c) coordination number, and (d) channel length of the pore-throats networks for the three Viking samples. The sandstone data are represented by solid lines and ‘•’. The shaly sandstone data are represented by dotted lines and ‘ x ’. The conglomerate sandstone data are represented by dashed lines and ‘ □ ’. ...................................................................................................................... 49 Fig. 13. Pore, kaolinite, and MoI volume distributions of (a) sandstone (Viking3W4), (b) shaly sandstone (Viking10W4), and (c) conglomerate sandstone (Viking14W5). ......... 50 Fig. 14. Pore, 2 nd kaolinite, and MoI surface area distributions of (a) sandstone, (b) shaly sandstone, and (c) conglomerate sandstone. ............................................................. 52 Fig. 15. Kaolinite pore volume distributions of (a) sandstone, (b) shaly sandstone, and (c) conglomerate sandstone. MoI pore volume distributions of (d) sandstone, (e) shaly sandstone, and (f) conglomerate sandstone. .............................................................. 53 Fig. 16. Plots of kaolinite pore volume vs. pore volume, (a) sandstone, (b) shaly sandstone, and (c) conglomerate sandstone. ............................................................. 54 Fig. 17. Plots of MoI pore volume vs. pore volume, (a) sandstone, (b) shaly sandstone, and (c) conglomerate sandstone. ............................................................................. 55 Fig. 18. Kaolinite pore area distributions of (a) sandstone, (b) shaly sandstone, and (c) conglomerate sandstone. MoI pore area of (d) sandstone, (e) shaly sandstone, and (f)

ix

conglomerate sandstone. ........................................................................................ 57 Fig. 19. Plots of kaolinite pore area vs. pore surface area, (a) sandstone, (b) shaly sandstone, and (c) conglomerate sandstone. ............................................................. 58 Fig. 20. Plots of MoI pore area vs. pore surface area, (a) sandstone, (b) shaly sandstone, and (c) conglomerate sandstone. ............................................................................. 59 Fig. 21. Effective reaction rates (mol/m 2 s) of anorthite and kaolinite. ......................... 66 Fig. 22. Distributions by pore: (a) pH, (b) [Ca 2+ ], (c) saturation index of anorthite, SI A , (d) saturation index of kaolinite, SI K , (e) dissolution rate of anorthite (mol/m 2 s), (f) dissolution rate of kaolinite (mol/m 2 s), and (g) precipitation rate of kaolinite (mol/m 2 s) in the steady state. ..................................................................................................... 67 Fig. 23. Effective (solid lines) and volume-averaged (dotted lines) reaction rates of Viking3W4 (red) and Viking14W5 (green) between 0 sec. and 2000 sec. (a) reaction rates of anorthite with v seep = 0.0058 cm/s (b) v seep = 0.001 cm/s (c) reaction rates of kaolinite with v seep = 0.0058 cm/s (d) v seep = 0.001 cm/s. The reaction ratio, β, is given at the end time of each graph. (e) Effective reaction rates of kaolinite of Viking14W5 with v seep = 0.0058 cm/s and v seep = 0.001 cm/s from 0 sec. to 10000 sec. Four lines are identical to the green lines in (c) and (d). .................................................................................. 72 Fig. 24. Distributions by pore of the saturation index (SI = log ). (a) case 1, (b) case 2, (c) case 3, and (d) case 4. ....................................................................................... 73 Fig. 25. Average pH vs. time for cases 1 – 4. ............................................................ 73 Fig. 26. pH distributions by pore. (a) Case 1, (b) case 2, (c) case 3, and (d) case4. ........ 75 Fig. 27. Distributions of the anorthite reaction rate by pore. (a) Case 1, (b) case 2, (c) case 3, and (d) case 4. ................................................................................................... 76 Fig. 28. Distributions of the kaolinite dissolution/precipitation rates by pore. (a) Case 1,

x

(b) case 2, (c) case 3, and (d) case 4. ....................................................................... 78 Fig. 29. Distributions of the saturation index of the kaolinite reaction by pore. (a) Case 1, (b) case 2, (c) case 3, and (d) case 4. ....................................................................... 79 Fig. 30. Effective and volume-averaged reaction rates of (a) anorthite and (b) kaolinite from t = 0 to t = 1300 second. Distributions of (c) the saturation index of the kaolinite reaction and (d) the kaolinite precipitation rate. Case 5. ............................................. 81

xi

List of Tables

Table 1. Equilibrium K eq and reaction rate k constants for anorthite and kaolinite ........ 27 Table 2. Instantaneous reactions and equilibrium constants. ....................................... 29 Table 3. 9 Derived instantaneous reactions, and corresponding equilibrium equations to compute the activities of the secondary species. ....................................................... 30 Table 4. Total concentrations of all components. ....................................................... 31 Table 5. Table for the equilibrium computation. ........................................................ 33 Table 6. Table of temperature vs. dielectric constant of pure water [41]. ...................... 40 Table 7. Threshholds, T 0 and T 1 for bi-phase segmentations. ...................................... 44 Table 8. Mineral abundances and accessibilities of three Viking samples. The abundances and accessibilities are compared with those by BSE analysis [19]. .............................. 48 Table 9. The numbers of pores which contain kaolinite or are in contact with MoI. ...... 49 Table 10. Average of volumes of pore, kaolinite, and MoI. ........................................ 51 Table 11. Summary of the generated network. .......................................................... 63 Table 12. Three reactions for the condition of the top boundary inflow. ....................... 63 Table 13. Initial state of pores. pH is 6.6. ................................................................. 64 Table 14. Composition of the injected fluid. CO 2 solubility is 2.0M and pH is 2.9. ....... 64 Table 15. Effective reaction rates,

, volume-averaged reaction rates,

, and the ratio

/

. The result is compared with Li’s. ........................................................... 65 Table 16. Concentrations in the injected saline solution. CO 2 solubility of 1.01M and pH of 3.01. ................................................................................................................ 69 Table 17. Simulation cases. .................................................................................... 69

xii

Table 18. Simulation result of the effective and volume-averaged reaction rates. The reaction rates are computed at the times given in the table. ........................................ 74 Table 19. Comparison of Ω with Ω Ω K . .............................................................. 80

xiii

Acknowledgement

First of all, I would like to thank my advisor, Professor Lindquist, for his guidance, support, and encouragement. He provided me of a great chance to start a new study on a new area. I thank committee members, Professor Joseph Mitchell, Professor Xiaolin Li, and Professor Troy Rasbury for their advices and comments on my dissertation. I appreciate Professor Peters’ work on BSE and EDX analyses which is basis of part of this work. I thank Joon Dong for teaching a lot of programming skills.

I wish to thank my wife, Moonsun for her love and support of me. I thank my parents and all my family too.

1

I. Introduction

The anthropogenic use of fossil carbon resources is generating large amounts of atmospheric CO 2 , which is the dominant contribution to global warming. There are a number of schemes proposed to reduce the increase in the atmospheric CO 2 load, one of which is storing the CO 2 in the subsurface [1-4]. Such CO 2 sequestration consists of four basic processes; capture, transport, storage, and monitoring. Capture is the process of extracting CO 2 from hydrocarbon combustion. Storage is the process of injecting the captured CO 2 in geologic formations such as deep saline aquifers, depleted oil/gas reservoirs, and unminable coal seams. Deep saline aquifers are numerous and widely distributed and, thus, represent great potential as CO 2 reservoirs. Furthermore, the choice of deep saline aquifers is economically efficient because there are greater chances of finding appropriate groundwater aquifers located near CO 2 sources (such as industrial power plants) thus minimizing transport costs. To verify that large amounts of CO 2 can safely and efficiently be stored over relevant time scales (centuries to millennia), the subsurface reactive flow in saline aquifers must be completely understood. The reaction kinetics of the CO 2 sequestration process is affected by: host rock type, secondary minerals, faulting and fracturing, as well as pore-level effects. Laboratory measurements on reactive kinetics typically involve crushed, well-mixed samples – thus excluding effects of cementation and pore network access which can result in large difference between reaction rates extracted from laboratory experiments and those observed at the field scale [5]. Reactive flow analyses

2

need to be performed at the pore scale to catch the micro scale effects, and adequate upscaling techniques are required to obtain parameters for continuum scale simulation [4, 6]. There are a number of studies of pore scale simulation of reactive subsurface flows: [4, 6-10] . In [4], pore scale reactive flow relevant to geological CO 2 sequestration is simulated using network flow modeling, and a methodology is suggested for upscaling from micro scale results to macro scale parameters. A network flow model is one tool to analyze pore scale fluid flow. It can handle single- and multi-phase flow with reactions. The method utilizes a simplified pore-channel network to capture fluid transport through porous rock; it does not include the detailed geometry of each pore or pore-to-pore connection. It can handle a rock sample of the size of a few millimeters to a centimeter. For the simulation of reactive flow, reactive mineral information is also required. The network used for the simulation in [4] is based upon a network of regular lattice connections with pore and channel properties based upon statistical models. The model considers reactions involving anorthite and kaolinite and involves a system of nine chemical reactions. The simulation in [4] shows that the network flow model captures the heterogeneities in reaction rates. However, the heterogeneities depend highly on the network data and mineral distribution; the simulation results are highly dependent on the conceptual network and mineral distributions used in the simulation. A growing body of work uses the lattice Boltzmann (LB) method to simulate pore scale reactive flows [6, 8, 10]. The LB method has several advantages. It is simple and flexible, and is appropriate for irregular shaped domains like porous media. Unlike network flow models, the LB method does not need to simplify the pore space. Segmented computed tomographic images can be used directly as the computational

3

domain with no need for surface construction and difficult grid generation which can be time consuming steps in PDE-based methods. The LB computation results capture the effect of the entire pore geometry resolving concentration gradients and velocities within a pore. Moreover, it does not need tools for analyzing the pore space to extract geometric quantities needed for network flow simulations. The LB method can handle change of pore geometry accompanying dissolution and precipitation, and mineral reactions are treated as boundary conditions on mineral surfaces [10]. The LB method is available for both single-phase and two-phase flow analyses to compute the absolute permeability and relative permeability of a porous rock [11, 12]. However, due to computational complexity, all LB simulation for reactive flow still use conceptual 2-dimensional porous media and conceptual reactions [8, 10].

(a) (b) Fig. 1. (a) Conceptual network model in [4] and (b) conceptual 2D pore space of a rock in [8] used for pore scale simulation of reactive flow.

The use of synchrotron technology to generate X-ray computed micro- tomographic (CMT) images of rock samples and 3DMA-Rock, a software package for analyses of such CMT images, enable a direct geometrical analysis of pore structure of sedimentary rocks [13-17]. 3DMA-Rock reconstructs the pore-throat network of the

4

imaged sample by segmentation, medial-axis computation, throat finding, and pore space partitioning. The resulting pore-throat network is a direct basis for network flow models. 3DMA-Rock has been used to identify and analyze images of more than 2 phases: rock, water, and oil [18]. It contains algorithms for multi-phase segmentation, which we adopt (II.2.c) to identify minerals in a CMT image. Backscatter electron (BSE) and energy dispersive X-ray (EDX) analyses [19] have been used to identify minerals of Viking samples. The BSE analysis was used to produce a four phase segmented image of a sample surface; EDX was used to identify the minerals in each of the four segmented phases. Volume fractions of minerals and surface areas fractions between void phase and minerals were estimated from this analysis. The analyses are indicators for mineral identification, but the analyses cannot directly produce a computational domain for reactive flow simulation since the analyses are 2-dimensional. The goal of this research is to develop a methodology to construct a network model from real rock, to simulate reactive flow in the network, and to upscale pore scale reaction rates to the continuum scale. Network models were constructed from analysis of three CMT images of Viking samples (a sandstone, a shaly sandstone, and a conglomerate sandstone). A reactive transport flow model is developed using the network model to simulate reactive flow in the Viking samples. In Section II.1 and 2, we describe the use of CMT to segment the three dimensional Viking samples into four phases in analogy to the BSE analysis of [19]. Mineral distributions for reactive flow modeling are computed based upon the four-phase segmented CMT images. The mineral distributions are included in the reactive pore- throat network.

5

In Section II.3, we describe the transport model to compute volume flow rate through all channels for reactive flow simulation. In Section II.4, the reactive transport model is presented. The reactive transport model computes the rate of concentration change in a pore by advection, diffusion, and two mineral reactions: anorthite and kaolinite. In Section 0.1, the three Viking samples are analyzed to construct network models of real rock. In Section 0.0, the network simulation is applied to two of the network models: the sandstone and the conglomerate sandstone. Initially, the porous rock is filled with brine at normal pH. Acidic brine with CO 2 dissolved is then injected into the rock ( Fig. 2 ). The injected brine is transported through the porous rock interacting with accessible minerals. The simulation includes transport of saline water and the aqueous species relevant to mineral reactions. It also includes dissolution and precipitation but ignores any accompanying change of geometry. We analyze the simulation results to examine pore scale effects in the reactive flow and compare the results with those of [4]. Finally, conclusion and future study is discussed in Chapter IV.

Fig. 2. A schematic diagram of the simulation of the reactive flow. CO 2 dissolved water is injected from the top boundary into porous rock. The flow is similar to the case A in [4].

6

II. Numerical scheme

1. 3D image analysis

3DMA-Rock is a software package to analyze a 3D X-ray computed micro- tomography (CMT) image of rock sample [13-16]. It uses a grey scale CMT image to identify the pore space and analyze the pore space by four steps: medial axis computation; throat computation; pore-throat network construction. Medial axis is a skeleton of the pore space, and throat computation is based on the obtained medial axis. The pore space is partitioned into pores by throats, and a pore-throat network is constructed by the connectivity of pores and throats. The pore-throat network is a basis of a network flow model for reactive flow simulation. The pore-throat network contains a network structure of the pore space, and throats and medial axis determine the channel properties.

a. Segmentation Segmentation is the term for image processing algorithm to perform phase identification on a grey-scale image. In this thesis, the kriging based algorithm [17] is used to segment the grey-scale CMT images. The method requires two threshold values, T 0 and T 1 . Voxels of intensity less than T 0 are assigned to one phase (Π 0 ) and voxels of intensity greater than T 1 are assigned to the other phase (Π 1 ). The identity of the other voxels of intensity between T 0 and T 1 are determined by the kriging algorithm, which uses indicator variables,

7

A

;P

1,i f

P

, 0,otherwise. (1)

Here,

P

is the attenuation coefficient at the position, P

, and

is a chosen threshold value. The probability that

P

for voxel, P

, is estimated by the unbiased linear estimator of the indicator variable,

P

;P

A

;P

. (2)

The covariance of the indicator variables at two different points is denoted by

;P

P

CovA

;P

,A

;P

. The

;P

are obtained by solving the kriging system,

;P

;P

P

;P

;P

P

for 1,, ,

;P

1 , (3)

where 1,, represents neighborhood voxels of P

. For each unidentified voxel, P

, the two probabilities,

P

and

P

1

P

are estimated by (2). If

, the voxel, P

is identified as phase Π 0 . Otherwise, the voxel, P

is identified as phase Π 1 .

8

In pore-grain segmentation, Π 0 is usually equivalent to the pore phase, and Π 1 to the grain phase. After pore-grain segmentation, clusters of isolated grain voxels which are physically unrealistic are removed (converted to pore phase). Clusters of isolated pore voxels are also removed (converted to grain phase) since such isolated pore space is not accessible to fluid flow.

b. Medial axis The media axis is a lower dimensional structure containing the basic network topology of the phase in which it is embedded [13]. As it retains a strict geometrical relationship to the phase in which it is embedded. It is useful for the further analyses of the phase structure due to its lower dimensional structure. A medial axis for the pore phase is obtained via a thinning algorithm [20]. Since care is taken to avoid isolated clusters of grain voxels, the medial axis for the pore phase reduces to a network of digitized curves (paths) and vertices (branches). The medial axis computation of a digitized image is sensitive to surface noise. This can result in dead-end branches or the medial axis whose presence contains little to no useful information. All branch-leaf paths of length less than the maximum burn number of the path voxels are deleted by the trimming step [13]. Intuitively, branch clusters corresponds to pores in the network and paths correspond to channels. In fact, a pore body may contain several clusters and the paths which connect them. These occurrences must be identified. We do so by defining “close” clusters. Two clusters are close if the length of a path connecting two clusters is smaller than both distances of the clusters to the closest grain. Two close clusters and the paths

9

connecting them are conceptually merged into a single cluster unit. After all such merges, the resulting clusters and paths are in 1-1 correspondence with pores and channels respectively [21].

c. Throat computation A throat is the channel cross section of minimum area corresponding to any position of medial axis path. Throats are important for fluid flow analysis as each represents a bottleneck to the flow along its channel. In 3DMA-Rock, throats are used to partition the void space into pore bodies (separated by throats). 3DMA-Rock utilizes three algorithms for throat computation, a “wedge based” algorithm, Dijkstra shortest path algorithm [13, 16], and the maximal balls algorithm [22]. We utilize the first two in this work. The wedge based algorithm constructs a local grain boundary for each medial axis path by a grass fire algorithm and constructs a rough approximation to the cross section perimeter on the local grain boundary for each voxel on the medial axis. The algorithm then compares the cross section areas determined by the rough perimeters locating the minimum cross section area. A finer reconstruction of the minimum area perimeter is computed and a triangulation performed to construct the spanning surface. This triangulated surface is defined to be the throat. The Dijkstra based algorithm dilates the medial axis path radially and sequentially records the surface grain voxels that are encountered by the dilated object. The growing set of encountered grain voxels is searched (via Dijkstra algorithm) to find a loop enclosing the medial axis. The first occurrence such a constructable loop is determined to

10

be the throat perimeter for the path. First, the wedge based algorithm is applied to all medial axis paths. To the medial axis path on which the wedge based algorithm failed to find a throat, the Dijkstra based algorithm is applied to find a throat. After this, two sets of throats computed are merged into a full set of throats [21].

d. Pore-throat network construction To construct a pore-throat network, the void space is partitioned into pore bodies. Each pore body is delimited by surface of grain voxels and by throat barriers. We compute pore volume, pore surface area, coordination number, and principal diameters for each pore [21, 23]. The principal diameters are pore widths in the three principal directions in which the products of inertia vanish. With the obtained pores, we can construct pore-throat network by using the pore-pore connectivity and pore-throat connectivity.

11

2. Pore-throat network of mineral distributions

The rates of production and consumption of a species by a chemical reaction are proportional to the reaction rate and the reactive surface area. The distributions of minerals and reactive surface areas are crucial for a reactive flow analysis. A previous research used conceptual network and mineral distributions for network flow simulation [4], and another research used conceptual pore space and mineral distribution [8] for reactive flow simulations as shown in Fig. 1 . In this section, a new methodology is introduced to compute the distributions of minerals and reactive surface areas of Viking samples for the simulation of reactive flow through the rock samples. The results of backscatter electron (BSE) and energy dispersive X-ary (EDX) analyses [19] done by Peters are referred to identify minerals in the Viking samples. BSE analysis produces a four phase segmented image of a sample surface, and EDX captures minerals at some points to identify minerals in each phase of BSE images. 3D CMT images of the Viking samples are segmented into four phases corresponding to BSE images. The mineral distributions of the Viking samples are computed using the four- phase CMT images. Furthermore, a pore-throat network is constructed for a four-phase CMT image, and the mineral distributions are added to the pore-throat network to construct a network model. The network model is used for simulation of reactive flow in Section II.4.

a. Viking sample Three Viking samples, Viking3W4, Viking10W4, and Viking 14w5 which are described as sandstone, shaly sandstone, and conglomerate sandstone respectively are

12

analyzed to construct the pore-throat networks of mineral distributions. The three Viking samples are real rock samples from the abandoned oil/gas reservoirs in the Viking formation in the Alberta basin. The samples are imaged at the X2B beam line at the National Synchrotron Light Source at Brookhaven National Laboratory. The resolution of the CMT images is 3.98μm and the size is originally 723×723×600. The images are resized to remove regions of bad image quality, and the sizes of resized Viking3W4, Viking10W4, and Viking14W5 are 723×723×356, 723×723×405, and 723×723×405 respectively. The corresponding real dimensions of the resized CMT images are 2.88×2.88×1.42 mm 3 , 2.88×2.88×1.61 mm 3 , and 2.88×2.88×1.61 mm 3 respectively.

b. BSE and EDX analyses A scanning electron microscope (SEM) is used in backscatter mode to generate images by counting backscattered electrons (BSE). The number of backscattered electrons is proportional to the mean atomic number (MAN) on the scanned surface element. The MAN surface map thus generated produces an image of mineralogical variation on the surface. However, unique identification of minerals is not, in general, possible as different minerals can bare similar MANs. Energy dispersive X-ray (EDX) analysis generates a spectrum of emitted X-ray energies for each point on a sample surface. EDX is therefore used to identify minerals on the surface by the chemical composition. EDX is much more time consuming than BSE. EDX can be used on small patches to identify species present, and BSE can be used over larger areas to infer the distribution of these species.