University of Southern California Dissertations and Theses

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

(USC Thesis Other)

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

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content MOLECULAR- AND CONTINUUM-SCALE SIMULATION OF SINGLE- AND TWO-PHASE FLOW OF CO2 AND H2O IN MIXED- LAYER CLAYS AND A HETEROGENEOUS SANDSTONE By Mahsa Rahromostaqim A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (CHEMICAL ENGINEERING) May 2020 ii Abstract Flow, transport, swelling, and deformation of porous media as a result of injecting CO2 in the media, and in the presence of brine, are important phenomena that have been studied over the past decade or so, because they are relevant to sequestration of CO2 in large-scale porous formations, such as depleted oil and gas reservoirs. This thesis studies these phenomena at two distinct length scales, namely, molecular and core scales, using extensive computer simulations. We present the results of computer simulations of single- and two-phase flow of water and CO2 in porous media, at both molecular and continuum scales. At molecular scale, molecular dynamics (MD) simulations have been carried out to understand swelling of mixed-layer (ML) clays. The phenomena are studied in the two most common ML clays, i.e., illite-montmorillonite (I-MMT) and chlorite-montmorillonite (CH-MMT); the hydration energy, basal spacing, density profiles and radial distribution functions are computed, and their implications are discussed. Swelling of the I-MMT is studied as a function of the water concentration and the interlayer cations (Na + and K + ). For the Na−MMT with more octahedral substitutions, weak cation−surface interactions result in fully hydrated ions and significant swelling. In the asymmetric interlayer of the ML clays, however, the illite sheets with more tetrahedral substitutions cause strong adsorption of the cations onto the clay surface. In addition, given that the hydration enthalpy of K + is smaller than that of Na + , swelling is inhibited as the ratio K + /Na + increases. We have also studied swelling of I-MMT in the presence of water and CO2, as well as the interlayer cations. At low CO2 concentrations in the MMT, weak ion-surface interactions result in fully hydrated ions and, therefore, more extensive swelling than in the I-MMT. At higher CO2 concentrations, however, since in the two- and three-layer hydration (2W and 3W, respectively) states the hydrophobicity of the MMT surface is stronger than that of illite, the disruption of the water network by CO2, reduces the difference between the extent of swelling of the ML clays and pure MMT. Swelling of CH-MMT MLCs has also been studied as a function of water and the interlayer cations concentration. In both CH-CH and CH-MMT the strong octahedral substitution of the chlorite iii layer results in, respectively, increasing both polarization and adsorption of water near, and onto the clay surface. The latter reduces hydration of the interlayer cations and, consequently, swelling of the CH-CH and CH-MMT clays. Compared with the octahedral substitutions, the interlayer cation concentration and tetrahedral substitutions are shown to have substantially weaker effect on swelling, whose pattern is also a function of the type of the interlayer cation. We also find that the higher the hydration energy and the smaller the atomic radius of the interlayer cation, the more swelling of the clay interlayers occurs. At the continuum scale, we present the results of an extensive study of characterization of the structure of the morphology and the spatial distribution of the heterogeneities of a sandstone, the Mt. Simon sandstone in which a pilot CO2 sequestration project was carried out sometime ago, and numerical simulation of single-phase and two-phase flow of CO2 and brine in the formation’s three-dimensional images. In addition to understanding the structure of the sandstone, the goal is to compare the accuracy and computational efficiency of three distinct simulation approaches, namely, the lattice-Boltzmann (LB) approach, direct numerical simulation (DNS) of the governing equations of fluid flow that uses the finite-volume method with the OpenFOAM simulator, and pore-network (PN) simulation. After validating the simulator by comparing the computed relative permeabilities that they produce for Berea sandstone, we simulate displacement of brine by CO2 at low and relatively high capillary numbers and compute the relative permeabilities and other quantities of interest. We demonstrate that all the three methods provide consistent relative permeability-saturation functions that are in close agreement with one another and, therefore, the question of which method one should use is closely linked with the available computational resources, and the computational time that one can afford. iv Acknowledgments I would like to express my deepest gratitude to my supervisor, Prof. Muhammad Sahimi for his support, encouragement and insightful directions during my Ph.D. study. It has been a privilege for me to be a member of his research group and to learn from him. I would like to extend my appreciation to Prof. Tsotsis from the Chemical Engineering Department and Prof. Felipe de Barros from the Civil Engineering Department at USC for agreeing to be my Ph.D. committee members. I would like to thank my collaborators, Amir Hossein Kohanpur and Prof. Albert Valocchi, at the University of Illinois at Urbana-Champaign for their contributions to the effort for my research. Last, but by no means least, I would like to express my appreciation to my family, especially my husband Mohammad, for their love and unconditional support, without which my work would not have been possible. Financial support for this work was provided by the Center for Geologic Storage of CO2, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award DE-SC0C12504. The calculations were carried out using 12 nodes (each with 16 processors) of the computer network at the University of Southern California High-Performance Computing Center. The support of both is gratefully acknowledged. v Table of Contents 1 Introduction ........................................................................................................................... 1 1.1 Background .................................................................................................................................................... 1 1.2 Molecular Scale .............................................................................................................................................. 2 1.2.1 Molecular Dynamics Simulation ................................................................................................................ 3 1.2.2 Illite-Smectite Clays .................................................................................................................................... 4 1.2.3 Chlorite-Smectite Clays .............................................................................................................................. 6 1.3 Continuum Scale ............................................................................................................................................ 7 2 Molecular Dynamics Simulation of Swelling of Mixed-Layer Clays in the Presence of Water ................................................................................................................................ 11 2.1 Introduction .................................................................................................................................................. 11 2.2 Molecular Dynamics Simulation ................................................................................................................ 11 2.3 Results and Discussion ................................................................................................................................. 13 2.3.1 Hydration Energies ................................................................................................................................... 14 2.3.2 Swelling .................................................................................................................................................... 15 2.3.3 Molecular Structure of the Interlayer Region ........................................................................................... 18 2.3.4 Radial Distribution Functions ................................................................................................................... 21 2.4 Summary ....................................................................................................................................................... 24 3 Molecular Dynamics Simulation of Hydration and Swelling of Mixed- Layer Clays in the Presence of Carbon Dioxide ..................................................................................... 25 3.1 Introduction .................................................................................................................................................. 25 3.2 The Force Field and Details of Simulation Procedure ............................................................................. 25 3.3 Results and Discussion ................................................................................................................................. 28 3.3.1 Swelling and Basal Spacing ...................................................................................................................... 28 3.3.2 Density Profiles ......................................................................................................................................... 32 3.3.3 Radial Distribution Functions ................................................................................................................... 35 3.3.4 Diffusion Coefficients ............................................................................................................................... 42 3.4 Summary ....................................................................................................................................................... 43 4 Molecular Dynamics Study of the Effect of Layer Charge and Interlayer Cation on Swelling of Mixed-Layer Chlorite-Montmorillonite Clays .............................................. 45 4.1 Introduction .................................................................................................................................................. 45 vi 4.2 Molecular Models and Simulation Protocol .............................................................................................. 45 4.3 Results and Discussion ................................................................................................................................. 47 4.3.1 Swelling .................................................................................................................................................... 48 4.3.2 Molecular structure of the interlayers ....................................................................................................... 51 4.3.3 Radial distribution functions ..................................................................................................................... 57 4.4 Summary ....................................................................................................................................................... 61 5 Two-phase flow of CO2-brine in a heterogeneous sandstone: Characterization of the rock and comparison of the lattice-Boltzmann, pore-network, and direct numerical simulation methods .............................................................................................................. 62 5.1 Introduction .................................................................................................................................................. 62 5.2 Characterization of Mt. Simon Sandstone ................................................................................................ 62 5.2.1 Construction of the pore-network ............................................................................................................. 64 5.2.2 Analysis of the heterogeneity .................................................................................................................... 67 5.2.3 Determining the size of the representative elementary volume ................................................................ 69 5.3 Computational Approaches ........................................................................................................................ 70 5.3.1 The lattice-Boltzmann method .................................................................................................................. 70 5.3.2 Direct numerical simulation ...................................................................................................................... 73 5.3.3 Pore-network model .................................................................................................................................. 74 5.4 Test of the Accuracy of the Numerical Approaches ................................................................................. 75 5.5 Results and Discussion ................................................................................................................................. 76 5.5.1 Single-phase flow ...................................................................................................................................... 76 5.5.2 Two-phase flow ......................................................................................................................................... 78 5.6 The importance of resolution of the computational grid ......................................................................... 83 5.7 Summary and Conclusions ......................................................................................................................... 84 Bibliography ................................................................................................................................. 85 vii List of Tables Table 2.1. Dependence of the Basal Spacing d (in Å) on the Type of Interlayer Cations in Various Hydration States ............................................................................................................... 16 Table 3.1. The CO2 parameters used in the MD calculations [100] .............................................. 27 Table 3.2. The computed diffusion coefficients (×10 −5 ) in the mixed-layers clays (in cm 2 /s). .... 43 Table 5.1. Rock source and properties. ........................................................................................ 63 Table 5.2. Porosity of eight subsamples of Mt. Simon sandstone ................................................. 64 Table 5.3. Properties of the selected Mt. Simon sandstone subsample S2. The permeabilities are in mD. ...................................................................................................................................... 65 Table 5.4. The connectivity and anisotropy of the eight subsamples ............................................ 66 Table 5.5. Berea Sandstone characteristics. .................................................................................. 67 Table 5.6. Comparison of the coefficient of variation (CV) of the properties of the Mt. Simon sandstone sample S2 and the Berea sandstone sample .................................................................. 68 Table 5.7. Properties of CO2-brine pair for the LB simulation and DNS .................................... 75 Table 5.8. Comparison of the absolute permeability of the eight subsamples computed by the PNM, LB and DNS methods. The results obtained by the LB and DNS methods are for two sample sizes of 250 3 and 500 3 that have resolutions of 5.6 and 2.8 μm, respectively. ................. 78 Table 5.9. Properties of CO2-brine two-phase flow system for LB simulation and DNS ............. 79 viii List of Figures Figure 1.1. mixed layer types .......................................................................................................... 2 Figure 2.1. Supercell structure of mixed-layer clay. Top and bottom layers are illite, while the two middle layers are MMT with the interlayer cations. .............................................................. 12 Figure 2.2. Hydration energy as a function of the water content for MMT−MMT (blue), I−MMT with Na + only (yellow), I−MMT with 72% K + and 28% Na + (black), and I−MMT with K + (red). ......................................................................................................................................... 14 Figure 2.3. Basal spacing d of (a) MMT−MMT; (b) I−MMT with Na + ; (c) I−MMT with 72% K + and 28% Na + , and (d) I−MMT with K + . .................................................................................. 16 Figure 2.4. Comparison of the absolute(a) and relative(b) basal spacings as a function of the water content for MMT−MMT (blue), I−MMT with Na + only (yellow), I−MMT with 72% K + and 28% Na + (black), and I−MMT with K + only (red). ................................................................ 17 Figure 2.5. Atomic density profile for the 1W and 2W states. (a,b) MMT−MMT; (c,d) I−MMT with Na + ; (e,f) I−MMT with 72% K + and 28% Na + , and (g,h) I−MMT with K + , for water−hydrogen (blue), water−oxygen (red), Na + (yellow), and K + (gray). ................................. 19 Figure 2.6. Radial distribution functions for Na−OW (left) and Na−OS (right) for (a,b) MMT−MMT, (c,d) I−MMT with 72% K + and 28% Na + , and (e,f) I−MMT with Na + . ................ 22 Figure 2.7. Radial distribution function of K−OW (left) and K−OS (right) for I−MMT with 72% K + and 28% Na + (a and b) and for I−MMT with K + only (c and d). ..................................... 23 Figure 3.1. Side view of the simulation cell for I-MMT. The top and bottom layers are illite, while the two middle layers are the MMT. ................................................................................... 26 Figure 3.2. Absolute basal spacings d of (a) MMT-MMT; (b) I-MMT with Na + ; (c) I-MMT with 72% K + and 28% Na + , and (d) I-MMT with K + . The number of CO2 molecules is per unit cell. ................................................................................................................................................ 29 Figure 3.3. Dependence of the relative basal spacings on the water content at a fixed number of CO2 molecules per unit cell. (a) 10 CO2; (b) 20CO2; (c) 40CO2; (d) 60CO2, and (e) 80 CO2 molecules. ...................................................................................................................................... 30 Figure 3.4. Density profiles of the 1W, 2W, and 3W hydration states, perpendicular to the clay surfaces of MMT-MMT, shown in (a), (b), and (c), and in Na+I-MMT shown in (d), (e), and (f). The results for water oxygen are represented by red, CO2 oxygen and carbon by black and ix blue, and Na + in green. In the I-MMT interlayer the left layer is illite and the right layer is MMT. ............................................................................................................................................. 33 Figure 3.5. Density profiles of 1W, 2W and 3W states perpendicular to the clay surfaces of I- MMT with 72% K + and 28% Na + , shown in (a), (b), and (c), and K+I-MMT, shown in (d), (e), and (f). The results for water oxygen are represented by red, CO2 oxygen and carbon by black and blue, Na + by green, and K + by purple. In the I-MMT interlayer the left layer is illite and the right layer is MMT. ....................................................................................................................... 34 Figure 3.6. Radial distribution functions for Na-OW (left) and Na-OS (right) at constant CO2 concentrations in (a) and (b) MMT-MMT; (c) and (d) I-MMT with Na + , and (e) and (f) I-MMT with 72% K + and 28% Na + ............................................................................................................ 36 Figure 3.7. Radial distribution functions of K-OW (left) and K-OS (right) at constant CO2 concentrations in (a) and (b) I-MMT with 72% K + and 28% Na + , and (c) and (d) I-MMT with K + only. .......................................................................................................................................... 37 Figure 3.8. Radial distribution functions for Na-OW (left) and Na-OS (right) at constant water concentrations in (a) and (d) MMT-MMT; (b) and (e) I-MMT with Na + , and (c) and (f) I-MMT with 72% K + and 28% Na + . ........................................................................................................... 38 Figure 3.9. Radial distribution functions of K-OW (left) and K-OS (right) at constant water concentrations in (a) and (b) I-MMT with 72% K + and 28% Na + , and (c) and (d) I-MMT with K + only. .......................................................................................................................................... 39 Figure 3.10. Radial distribution functions for OCO2-OS (left) and CCO2-HW (right) in (a) and (b) MMT-MMT; (c) and (d) Na + I-MMT; (e) and (f) I-MMT with 72% K + and 28% Na + , and (g) and (h) K+I-MMT. ........................................................................................................................ 41 Figure 3.11. Mean-square displacements (MSD) of water (a) and Na + (b) ................................... 42 Figure 4.1. Schematic side view of the simulation cell for the CH-MMT mixed-layer clay. 𝐌+ Represents the interlayer cation. ................................................................................................... 46 Figure 4.2. Comparison of the absolute (a)-(c); and relative (d)-(f) basal spacing of the clay interlayers as a function of interlayer cation and the water content, grouped by the cation type. ....................................................................................................................................................... 48 Figure 4.3. Comparison of the absolute basal spacing of the clay interlayers as a function of interlayer cation and the water content, grouped by the interlayer type. (a) MMT-MMT; (b) CH-MMT + cation, and (c) CH-CH, all plus the cation. ............................................................... 50 x Figure 4.4. Density profiles of (a)-(c) the 1W and (d)-(f) 2W hydration states perpendicular to the clay surface for Na plus (a, d) MMT-MMT; (b, e) CH-MMT, and (c, f) CH-CH. In the CH- MMT interlayer the left and right surfaces are CH and MMT, respectively. Water oxygen and hydrogen are represented by black and red, and Na + by blue. ...................................................... 52 Figure 4.5. Density profiles of the 1W (a)-(c) and 2W (d)-(f) hydration states, perpendicular to the clay surface for K plus (a, d) MMT-MMT; (b, e) CH-MMT, and (c, f) CH-CH. For the CH- MMT interlayer the left and right surfaces are CH and MMT, respectively. Water oxygen and hydrogen are represented by black and red, and K + by blue. ........................................................ 54 Figure 4.6. Density profile of the 1W (a)-(c) and 2W (d)-(f) hydration states perpendicular to the clay surface for Cs plus (a, d) MMT-MMT; (b, e) CH-MMT, and (c, f) CH-CH. In the CH- MMT interlayer the left and right surfaces are CH and MMT, respectively. Water oxygen and hydrogen are represented by black and red, and Cs + by blue. ....................................................... 55 Figure 4.7. Radial distribution function for 𝐍𝐚−𝐎𝐖 (a)-(c) and 𝐍𝐚−𝐎𝐒 (d)-(f) at 1W (black) and 2W (red); (a, d) Na + MMT-MMT; (b, e) Na + CH-MMT, and (c, f) Na + CH-CH. ............ 58 Figure 4.8. Radial distribution function for 𝐊−𝐎𝐖 (a)-(c) and 𝐊−𝐎𝐒 (d)-(f) at 1W (black) and 2W (red); (a, d) K + MMT-MMT; (b, e) K + CH-MMT; (c, f) K + CH-CH. ........................ 59 Figure 4.9. Radial distribution function for 𝐂𝐬−𝐎𝐖 (a)-(c) and 𝐂𝐬−𝐎𝐒 (d)-(f) at 1W (black) and 2W (red); (a, d) Cs + MMT-MMT; (b, e) Cs + CH-MMT; (c, f) Cs + CH-CH. .................... 60 Figure 5.1. Thin-section scan of Mt. Simon sandstone at a depth of 6700 feet, and a 𝟓𝟎𝟎𝟑 subsample S2 of the formation selected for the two-phase flow simulation. ................................ 63 Figure 5.2. Gray-scale (left) and segmented (right) formats of the 300 th slice (out of 500 slices) of the selected Mt. Simon sandstone subsample S2. ..................................................................... 64 Figure 5.3. The PNM of Mt. Simon sandstone and its pore-body and pore-throat size distributions. .................................................................................................................................. 66 Figure 5.4. Porosity and degree of anisotropy of 15 subsamples from the Mt. Simon sandstone sample S2 and their comparison with those of the Berea sandstone sample. ................................ 68 Figure 5.5. Statistics of the average coordination numbers of over 100 realizations of the pore- network of the Mt. Simon sandstone sample S2, and the corresponding Berea sandstone sample. ........................................................................................................................................... 69 Figure 5.6. Sample-size dependence of three properties of the subsample S2 of Mt. Simon sandstone in order to identify the size of the REV. ....................................................................... 71 xi Figure 5.7. Sample-size dependence of absolute permeability, computed by the LB simulation, of subsample S2 of Mt. Simon sandstone. .................................................................................... 72 Figure 5.8. Drainage relative permeabilities of the oil-water pair in a Berea Sandstone sample with contact angle of 𝟏𝟐𝟓 degrees, capillary number of 𝟓×𝟏𝟎−𝟓, and viscosity ratio of 𝟏.𝟑𝟐. ............................................................................................................................................. 76 Figure 5.9. The 3D pore volume of the subsample S5 (left) and the pressure distribution in the volume at the end of the simulation (right). .................................................................................. 77 Figure 5.10. CO2 invasion pattern (blue) computed by the LB simulation for capillary number 𝐂𝐚=𝟏𝟎−𝟒 at brine saturation 𝐒𝐰=𝟎.𝟓𝟎. Brine is shown in red and rock is not shown. ...... 79 Figure 5.11. Comparison of the computed CO2 saturation for capillary number 𝐂𝐚=𝟏×𝟏𝟎− 𝟒 and its dependence on the distance from the inlet at brine saturation 𝐒𝐰=𝟎.𝟓𝟎 during drainage. ........................................................................................................................................ 80 Figure 5.12. Time-dependence of the computed brine saturation obtained by the LB simulation (continuous curves) and the DNS (dashed curves) during drainage. ............................................ 80 Figure 5.13. 3D representation of the brine residual saturation (red) for 𝐂𝐚=𝟏×𝟏𝟎−𝟒 (left) and 𝐂𝐚=𝟑×𝟏𝟎−𝟓 (right) capillary numbers. The residual brine saturations are Sw = 0.30 and 0.50, respectively, reached at times 54 ms and 122 ms. ......................................................... 81 Figure 5.14. Side view (inlet on the left and outlet on the right side of each image) of the distribution of the brine (red) and CO2 in the pore space for 𝐂𝐚=𝟏×𝟏𝟎−𝟒 (left) and 𝐂𝐚= 𝟑×𝟏𝟎−𝟓 (right) at brine saturation, 𝐒𝐰=𝟎.𝟒𝟕𝟓. .................................................................. 81 Figure 5.15. CO2-brine relative permeability curves in low capillary number, 𝐂𝐚=𝟑×𝟏𝟎− 𝟓. .................................................................................................................................................... 82 Figure 5.16. CO2-brine relative permeability curves for high capillary number, 𝐂𝐚=𝟏×𝟏𝟎− 𝟒. .................................................................................................................................................... 82 Figure 5.17. Effect of resolution of lattice in LB simulation on the resulting CO2-brine relative permeability curves for 𝐂𝐚=𝟏×𝟏𝟎−𝟒 (right) and 𝐂𝐚=𝟑×𝟏𝟎−𝟓 (left). .......................... 83 1 1 Introduction 1.1 Background The huge amount of carbon dioxide produced everyday has made sequestration of CO : a vital problem for controlling greenhouse gas emissions in the atmosphere. Among the options that are currently being studied is injection and storage of CO : in large-scale porous formations, such as depleted oil and gas reservoirs. One of the main concerns in such a process is that the injection may trigger seismic events, which can sometimes be quite large [1, 2]. It is still not clear why and how injection of CO : causes the seismic events in some areas, while the same process in other regions does not produce any significant seismic activity. The fractures that can be generated in rock due to the seismicity or injection of CO : , or other types of activities, pose the risk of leakage for the supercritical CO : , stored underground, through a path to the earth surface [3]. Thus, due to its importance, the problem has been studied extensively by many research groups, investigating its various aspects [2, 4-8]. A primary reason for the complexity of the process is that a set of complex phenomena, ranging from adsorption and reaction, to fluid flow and flow-induced solid deformation may occur. In addition, a wide range of length scales needs to be considered, from the molecular scale up to the field scale. Although geological storage of CO : at field and core scales has been studied [6, 9-11], realistic modeling of the process, as well as understanding the relationship between the phenomena at distinct length scales has rarely been considered. At molecular scale, the interactions between the clay and the interlayer CO : , water and brine result in swelling of the rock that leads to significant changes in the properties of clay, such as its porosity, permeability, and wettability. It may even initiate formation of micro-fractures [4, 12]. Therefore, studying the interaction between supercritical (scCO : ) and clay is essential to understanding the long-term potential problems caused by CO : sequestration. Sequestration of CO : consists of injecting scCO : deep into the geological formation, covered by cap rocks that prevent migration of the gas back to the surface of the earth. The focus of the research in this dissertation is on studying the process at two distinct and disparate length 2 scales, namely, the molecular and continuum scales. Thus, let us discuss the phenomena at each scale. 1.2 Molecular Scale At the molecular scale, cap rocks are mostly composed of clays that are typically illite, kaolinite, chlorite or smectite [13]. In fact, more than 70% of the sedimentary rock in the United States contains combination of such clay minerals [14] in a mixed state. One important factor for the abundance of mixed-layer (ML) clays is that, under a solid-state transformation, pure clay layers are converted into mixed ones. For example, in smectite illitization the interlayer smectites are replaced by illite, resulting in gradual change of the crystal size and shape [15]. The 2:1 or tetrahedral-octahedral tetrahedral (TOT) clay minerals, namely, illite, smectite, chlorite and vermiculite, have similar stacking structures with different types and magnitude of substitutions in their internal sheets. They are composed of octahedrally-coordinated Al atoms sandwiched by the tetrahedrally-coordinated Si atoms. The TOT layers are negatively charged due to isomorphic substitutions by Mg, Al, Fe, etc. ML clays are defined as a sequence of distinct kinds of clay layers. The mixing can, in vertical direction, be regular, segregated regular, or random. It has been reported [16] that large percentages of ML clays, about 40% and higher, are randomly interstratified. Commonly found ML clays include illite-smectite, illite-vermiculite, chlorite-vermiculite, chlorite-smectite, and kaolinite-smectite [13]. Figure 1.1 presents a schematic view of the types of ML clays. Figure 1.1. mixed layer types Based on marine diagenesis, mixed illite-montmorillonite (I-MMT) and chlorite- montmorillonite are formed more than pure illite or chlorite. In pure clays only one type of cation 3 exists between the layers. However, ML clays can be formed as a result of the combination of all type of different layers, and can have a combination of the interlayer of original pure clays [14]. 1.2.1 Molecular Dynamics Simulation Molecular dynamics (MD) simulation is a powerful tool for studying swelling of ML clays. It can be used to address whether the size and charges of the interlayer’s cations play the main role in setting their spacing, or if it is the spatial distribution of the negative charges that determines the extent of the hydration and swelling. Such a combination of both sets of variables in ML clays has never been considered. Past molecular studies, both mole3cular dynamics (MD) and Monte Carlo simulations, focused on the effect of the spatial distribution of the layers’ charges, considered separate pairs of clay layers with different charge densities and distributions, and compared the results for the two. For example, McGoldrick et al. [17] studied separately the interactions of water with the MMT and beidellite, which is very similar to illite. Using only one-layer charge, Boek et al. [18] examined interlayer expansion for cations with different hydration properties. Teich- McGoldrick et al. [17] studied the interactions of water with the MMT and beidellite. They reported that at low water concentrations the basal spacing of beidellite is larger than or equal to that of the MMT, whereas at higher water contents, it is the basal spacing of the MMT that is larger. A study by Liu et al. [19] indicated that in Cs-smectites the location of the layer charge has a significant effect on the mobility of the interlayer counterions and their binding to the surface, although the difference in the swelling of the three different layer-charge distributions was reported to be negligible. Another study [20] focused on how shifting the negative charges from the octahedral to tetrahedral layers affects binding of the counterions to the surface and, thus, swelling of clays. Sun et al. [21] reported that with increasing octahedral charge fraction the swelling pressure increases, while Skipper et al. [22] claimed that increasing the tetrahedral layer charge results in the formation of the inner sphere surface complexes (ISSC) of Na. Other previous studies [23-26] also reported that since K < only perturbs the location of the water molecules, but is unable to hydrate, it prevents the clay from swelling, implying that the existence of K < in the MMT interlayer and the ratio R of K < and the overall number of cations contribute significantly to controlling the swelling, with 0.3 < R < 0.56 strongly inhibiting swelling. Several previous studies of interaction of clays’ layers with the ions focused on the smectites, rather than illite, and assumed that the clay is rigid. Very few studies [27-30] focused 4 on fully flexible clays. For example, Lammers et al. [30] used the force field flexible CLAYFF (see below) to investigate adsorption of cesium and sodium on the Illite surface and reported that, due to penetration of ions into the ditrigonal cavities, the ISSC with cation density peak is considerably closer to the clay surface. Finally, a comparison of the effect of water and CO : concentration on the hydration states of the MMT and beidellite was investigated by Makaremi et al. [31] group. They observed that in Na-beidellite, the isomorphic subsections greatly control the distribution of CO : molecules in the interlayer while in Na-MMT the hydrated sodium cations mainly affect the CO : distribution. Rao et al studied the effect of layer charge on intercalation of water and CO : . They observed that the d-spacing of the clays for monolayer and bilayer were independent of the layer charge. However, the transition between the two states happened at lower relative humidity for the case of high-charge clay than in the low-charge one [32]. Other studies focused only on investigating the effect of CO : and water in swelling of pure clays. Among them, Myshakin et al. [33] focused on the effect of CO : intercalation in swelling of hydrated MMT. They observed that the level of clay swelling caused lliteby the added CO : greatly depends on the initial water concentration in the clay interlayer. Rao et al. [34] studied intercalation of water and CO : in MMT and under geological conditions. They observed high mobility of CO : molecules, while they hardly penetrated the hydration shell of Na < that shows a low probability of a direct interaction between intercalated CO : and Na < . Using extensive MD simulations, we study in chapter swelling of two ML clays, namely, illite−MMT (I-MMT) clays for a range of water concentrations. We also study the same in the MMT−MMT clays whose interlayer contains Na < and use the results as a reference to compare with the results for the ML clats. In chapter three, in addition to water, CO : molecules will be added to the interlayer of the mixed clay structures of chapter two, and the behavior of both the interlayer species and the clay surface is studied. To our knowledge, our MD study of swelling of ML clays is the first of its kind. 1.2.2 Illite-Smectite Clays The common type of ML clays consists of expanded hydrophilic, such as the MMTs, and shrunken hydrophobic layers, such as illite or chlorite. The most abundant and interesting ML clays are illite-smectite [16]. One example is bentonite, a mixture of MMT and beidellite. Potassium-rich bentonites are dominated by illite, which have been proposed for use in addressing the problem of 5 waste disposal [35]. Their porosity is reduced by hydration, whereas the tendency of the sodium- saturated smectites to swell is the main reason for the observed instability in drilling operations, which may lead to the collapse of the wellborne. Swelling of clay minerals is of prime importance to the molecular-scale storage of CO : . Understanding the properties of swelling of ML clays in the presence of CO : will lead to its improved sequestration. Some studies have claimed that the size and charge of cations in the interlayer region of clay minerals play the main role in the swelling [18, 36-39], while others have argued that the amount and spatial location of the negative charges are the prime factors in determining the extent of hydration and, hence, the swelling. Most of the past studies, both experimental and computational, have, however, focused on pure clays of one type or another. Various parameters, such as temperature and pressure, affect swelling of clays. In addition, several other factors also influence the swelling, including [23, 40-46] the layer’s charges and their spatial distribution, their magnitude, and the ratio and type of the interlayer cations. As mentioned, the structures of both illite and the MMT is 2:1 or TOT pattern. The main difference between the structures of the two clays is in the location at which the substitutions occur. Most of the negative charges in illite are distributed in the tetrahedral sheet where substitution of Si A< by Al D< occurs, in contrast with the MMT in which the negative charges are dominant in the octahedral sheet where Al D< is substituted by Mg :< . However, overall, few substitutions occur in the tetrahedral and octahedral sheets of the two types of clays. To balance the negative charges in illite, potassium ions are distributed in the interlayer space. The K < ions do not have a strong tendency to hydrate and, thus, they lock the illite layers to each other. In contrast, the MMT is made electrically neutral by Na, < Mg :< or Ca, :< which are usually hydrated to a certain degree and, therefore, the clay swells easily. The mixed spatial distribution of Na < and K < in the interlayer region of the I−MMT is, therefore, important to balancing the layers’ negative charges. In addition, the concentration of the cations in the interlayer region and that of the water molecules and the interactions between them and the negatively charged layers are the important factors in swelling of MLCs. Using extensive MD simulations, we study in chapter two swelling of the I-MMT clays for a range of water concentrations. We also study the same in MMT−MMT whose interlayer contains Na < and use the results as a reference to compare with the results for the MCLs. 6 1.2.3 Chlorite-Smectite Clays The second most common type of ML clays is trioctahedral chlorite-smectite. Chlorites, similar to smectites and illites, have a 2:1 sandwich structure. This type of the ML clays is formed when Mg(OH) : replaces interlayer cations and water in the MMT. In the conversion of smectites to chlorites the brucite layer, generated by polymerization of the existing complex ions in sedimentary rock, locks together the two 2:1 sheets of converted smectite and forms a non- expandable 2:1:1 clay [47]. In the CH-MMT mixed clays, however, the space between the 2:1:1 chlorite and the adjacent MMT and the chlorite layers can swell during hydration. The structure of both chlorite and MMT clay layers consists of a 2:1 or tetrahedral-octahedral-tetrahedral- octahedral (TOTO) structure. In most of the mixed-layer CH-MMT layers, the chlorite clays are trioctahedral, whereas the MMT is dioctahedral. In the trioctahedral phyllosilicates structure of chlorite, the octahedral layers are similar to brucite, with Mg :< or Fe :< occupying the cation positions. On the other hand, in dioctahedral phyllosilicates of the MMT, the octahedral layers are similar to gibbsite, with Al D< occupying those positions. A chemical analysis of some of the mixed-layer CH-MMT samples are presented by Weaver and Pollard [48]. In both types of clays the negative charge of the TOT layer is balanced by the interlayer cations, such as Na < , Ca :< ,Cs < , and K < . Such cations, in the presence of water, are usually hydrated, resulting in swelling of the clays [49]. The extent of swelling depends, however, on the cation type [18, 36], the humidity, clay type, and the charge distribution [19, 50, 51] in the clay structure. The layers substations play an important role in swelling of clays. In two studies [52, 53] two MMT clays, one with both octahedral and tetrahedral substitution and a second one with only octahedral substitution, were considered. Under the same conditions, smaller extent of swelling was observed for the MMT with both octahedral and tetrahedral substitutions. Smith et al. [54] explored the effect of a range of octahedral layer charges on swelling of Na-MMT, and demonstrated that increased layer charge results in more swelling. In a separate study by Whitley and Smith [55] the increase in the layer charge affected swelling similar to the increase in the hydration energy of the interlayer cations. In an experimental study Foster [56] investigated the relationship between ionic substitution and swelling of the MMT. He found that the increase in the octahedral substitution reduces swelling of the MMT. He also considered the effect of the amount 7 of cation in the interlayer, as well as tetrahedral substitution in clay swelling, but reported that the two factors had little or no effect on swelling, whereas octahedral substitution plays a significant role in the swelling volume of the MMT. Therefore, in CH-MMT mixed layer clay with intense octahedral substitutions it is expected to see reduced amount of swelling compared to pure MMT clays. In addition to the asymmetry of the interlayer and layers substations, the interlayer cation type is also expected to have an important effect on swelling of mixed layer CH-MMT clays. In chapter four, we report the results of our MD simulations for understanding the swelling behavior of three interlayers, namely, the CH-MMT, CH-CH and MMT-MMT, as a function of water concentration. We also study the effect of interlayer cation type on swelling of the three interlayers. 1.3 Continuum Scale Flow and transport of water and CO : in porous media is a function of several forces. The flow is driven by gravity, the pressure gradient between the boundaries of the flow domain, and viscous forces that are resist the flow. The pore-scale forces and flow significantly affect the macroscopic nature of flow through porous medium, including injection and storage of CO : in depleted reservoirs and other types of large-scale porous formations [4, 57, 58]. The balance between such forces in the Navier-Stokes equation establishes the behavior at the continuum scale. Experimental studies of such phenomena are expensive and time consuming. It is also tricky and difficult to accurately measure viscous and capillary forces in experimental studies. As already mentioned, the physics of two-phase flows in porous media is of utmost importance to many problems of practical interest, such as CO2 sequestration in deep saline reservoirs, recovery of oil from hydrocarbon reservoirs, transport of non-aqueous-phase liquid contaminants in aquifers, and infiltration of rainfall into soil. Understanding and predicting the properties of such flow processes at various length scales are vital to addressing the problem of global warming, as well as managing water resources, and energy production. As emphasized earlier, there is increasing incentive for studying two-phase flow of CO2 and brine in porous media, due to societal interest in geological sequestration of CO2. The drastic increase in the amount of CO2 in the atmosphere plays the most important role in global climate change [59-61]. The main source of the atmospheric CO2 is combustion of fossil fuels for 8 producing electricity in power plants. Among the various approaches that have been suggested for mitigating the problem, CO2 capture and storage (CCS) is believed to be a viable solution [61]. The CCS is a process associated with separating CO2 from the other gases that are produced by power plants and other sources, compressing and transporting it to storage locations, and keeping it sequestered in onshore or offshore geological formations for very long times [60]. As an example, the United States can inject up to 65% of the CO2 produced by power plants deep into saline aquifers [62]. Since such a sequestration process appears to have great potential, understanding the behavior of flow of brine and CO2 in deep porous formations is vital [58, 63] to its successful implementation, and should accompany any comprehensive study of geological storage of CO2 that includes its economic feasibility, site selection, risk assessment, environmental impact, safety aspects, monitoring and verification, in addition to perspectives on retention time, physical leakage, brine displacement, and microseismicity [60]. The Illinois State Geological Survey carried out a pilot injection study to better understand the feasibility of full-scale CCS. The study site is in Decatur, Illinois, and the injection zone is the highly saline Mt. Simon sandstone [64]. Thus, one goal of the present study is to analyze two- phase flow of CO2 and brine in Mt. Simon sandstone and, more generally, in natural rock. Advances in modeling of porous media, as well as development of efficient computer simulation, together with increasing computer powers, have made it possible to model and study in detail two-phase flow of brine and CO2 at core plug and larger length scales. Compared to experiments, the computational approaches have the advantages of being generally less expensive and more flexible in implementing and changing parameters, flow and displacement mechanisms, and studying various mechanisms of displacements. At the same time, numerical simulation of flow of brine and CO2 at field scale requires such inputs as the relative permeabilities and capillary pressure as functions of the saturation at the scale of the grid blocks, which can be obtained by experiments, or by pore-scale modeling. Pore-scale modeling [65] can provide the required input data for field-scale modeling, provided that one takes into account the effect of the morphology of the core-scale porous media and the various pore-scale mechanisms of fluid displacements. This is because pore-scale flow affects significantly the characteristics of the process at larger scales, including injection and storage of CO2 in depleted reservoirs and other types of large-scale porous formations [57, 58, 63]. 9 With advances in instrumentation, it has become possible to obtain three-dimensional (3D) images of porous media with high resolution by, for example, X-ray computed tomography (CT) [66]. At pore-scale, detailed 3D geometry of rock and its void space can be captured by direct imaging using nondestructive X-ray microtomography [67, 68]. Various rock properties, such as the permeability, have been computed with the aid of micro-CT images [69-72]. Thus, one can directly simulate two-phase flows in the images, which renders as unnecessary developing models of the pore space. Image-based methods are, however, computationally expensive, although recent advances in high-performance computing, as well as a new method that smoothens the image without changing its properties using curvelet transformations (Aljasmi and Sahimi, 2019), are making use of such methods increasingly common. The methods that are used for simulating fluid flow in the images are either based on directly solving the governing equations for fluid flow, i.e., the Stokes’ equation using the finite-volume method [73-75], or based on the lattice-Boltzmann (LB) method [76-81]. The LB method has been demonstrated to be well-suited for high- performance parallel computing in the complex geometry of porous media [80, 82-84]. In addition, advances in method based on the volume of fluids, especially the new developments of surface force formulation, have made it possible to efficiently model two-phase flow by the finite-volume method at lower capillary numbers [85]. Another approach is based on pore-network models (PNMs) [66, 86] in which the pore space is simplified to a network of interconnected pore-bodies and pore-throats. While the PNMs are computationally very efficient and inexpensive, they still involve some assumptions and approximations, such as the definition of what constitutes a pore-body or pore-throat, how to assign effective sizes to them, etc. [80, 87]. Therefore, a detailed comparison between the results of the PNM computations with those obtained by the other two aforementioned methods will shed much light on its advantages and limitations. The goal of chapter five is twofold. One is computing the drainage relative permeability functions for CO2 and brine during injection of the former into a sample of Mt. Simon sandstone. This is not a trivial problem as Mt. Simon sandstone is much more heterogeneous than typical sandstones studied in the past, such as the Berea sandstone. Moreover, the viscosity ratio of the two fluids are much higher, about 10 times larger, than of that the typical oil-water systems. The two challenges lead to computational difficulties that require careful choices of the input 10 parameters, discretization, and boundary conditions. The second goal is evaluating the performance of the LB method and PNMs by comparing the relative permeabilities computed by them with those obtained through the direct numerical simulation of the fluid flow in the image of the sandstone. In addition, we also study the effect of the capillary number, i.e. the effect of flow, on the results. 11 2 Molecular Dynamics Simulation of Swelling of Mixed- Layer Clays in the Presence of Water 2.1 Introduction As described in chapter 1, one goal of this research is gaining an understanding at molecular scale of the effect of water and carbon dioxide on deformation and swelling of mixed clays, the primary structure of the porous formations that we study. In this chapter we describe our results that we have obtained, using molecular dynamics (MD) simulations, for swelling of mixed-layer (ML) clays in the presence of water. 2.2 Molecular Dynamics Simulation The MD simulations were carried out using the LAMMPS package, the standard software for MD simulations. We simulated ML clays of sodium-saturated Montmorillonite (MMT) and potassium−Illite with chemical compositions given, respectively, by Na N.OP [Si O.OP Al N.:P ](Al D.P Mg N.P )O :N (OH) A , and K N.S [(Al T.S Mg N.T )(Si D.: Al N.U )O TN (OH) : ]. The negative charges in the interlayer regions were generated by symmetrically selecting substitutions over the clays’ sheets. Half of the K + cations in illites-1.8 ions per unit cell-were in the illite−illite (I−I) interlayer, while the remaining half was equally distributed between the I−MMT interlayers. Figure 2.1 presents a snapshot of the equilibrated supercell, with the top and bottom layers being illite, surrounding the MMT layers, together with the cations. The clay layers were symmetric with respect to the central line in the middle between the two MMT sheets. As more water molecules were inserted in the interlayers, the supercell was allowed to freely expand along the z axis. We constructed a supercell of 80-unit cells that contained two layers each of MMT and illite. Each layer contained 20-unit cells, and in order to minimize the differences between the dimensions of Illite and MMT sheets along the x and y directions, they were organized as 5 × 4 unit cells; see Figure 2.1. The size of each unit cell was 5 × 9 Å 2 and 5 × 8.8 Å 2 for MMT and illite, respectively. Thus, under dry conditions the dimensions of the individual sheets or layers of MMT were approximately 26 × 36 Å 2 and 26 × 35 Å 2 for illite. Periodic boundary conditions were imposed in all directions. 12 The force field utilized to represent the interatomic interactions between the atoms in the clays was CLAYFF, based on which we generated the equilibrium molecular structures of both types of clay. The flexible SPC model was used to model the interlayer water molecules, but not the hydroxyl group of the clay layers (since CLAYFF already contains the necessary parameters for them). According to CLAYFF the total interaction energy E of the material is given by E WXW = q i q j r ij + 4ℇ ]^ _` a bc d bc e T: −` a bc d bc e f g + T : k i (r ]^ −r X ) : + 1 2 k B (θ ijk −θ o ) 2 (1) Here, k1 and k2 are force constants; θ ]^o is the angle between bonds ij and jk; e is the electron’s charge; q ] is the partial charge of atom i; ε0 is the dielectric permittivity of vacuum; and r0 and θ0 are the equilibrium values of the corresponding quantities. ℇ ] and σ ] are the usual Lennard-Jones (LJ) energy and size parameters, and the Lorentz−Berthelot mixing rules were used for the pairs ij. Figure 2.1. Supercell structure of mixed-layer clay. Top and bottom layers are illite, while the two middle layers are MMT with the interlayer cations. The simulations were initially carried out in the (NPT) ensemble at T = 348 K and P = 130 bar, which are close to the conditions for geological cap rock formation [88]. A time step of 0.01 fs was used for 100 ps, after which it was increased to 1 fs for 5 ns after which equilibrium was 13 reached. An additional 5 ns of simulation was used for the production step. Simulations longer than 10 ns were also carried out, but the difference between their results and those obtained with shorter simulations were negligible. The Nosé−Hoover thermostat and barostat were used with the coupling constants of 1000 and 100, respectively. After some preliminary simulations the long- range interactions were truncated for r > 13 Å. The electrostatic interactions were computed using the particle-mesh Ewald summation with a precision of 10 r A . To calculate the density profiles and the radial distribution functions, after equilibrium was reached via simulation in the (NPT) ensemble, we switched to the canonical (NVT) ensemble and continued the simulation for another 5 ns. A main goal of our study is gaining insight into the effect of asymmetric interlayer structure, as well as that of mixed counterions, K + and Na + , on swelling of I−MMT MLCs. Note that the illite has more negative charges in the tetrahedral sheets, whereas MMT contains more of such charges in the octahedral layer. Three series of simulations were carried out in which the counterion ratio in the interlayer was varied. In one series, the interlayer of I−MMT contained both K + and Na + , whereas in the other two series of simulations the ML clays contained only Na + or K + to compensate for the negative charges of illite and the MMT. This allowed us to monitor the effect of the spatial distribution of the layer’s charge location and that of the two types of ions on the swelling. To understand the effect of water concentration, the number of water molecules in the interlayer region was systemically increased from 1 to 10 molecules per unit cell, and for each case an independent simulation was carried out, which was enough to capture one-layer (1W) and two- layer (2W) hydration of the clays. All the simulations were initiated with inserting the cations in the midplane of each interlayer, while the initial spatial distribution of the water molecules was random. We computed the density profiles along the z axis and averaged them over the interlayer regions over a 5-ns interval in the (NVT) ensemble. 2.3 Results and Discussion We have computed several important properties of the clays that we are studying. In what follows we present and describe the results and discuss their implications. 14 2.3.1 Hydration Energies We computed the change ΔU(N) in the potential energy U(N) of the system as a function of the water content, ΔU(N) = [U(N) − U(0)]/N, where N is the number of water molecules and U(0) is the potential energy of the dry clay. The results, which represent the hydration energies, are presented in Figure 2.2 . In this and the following figures the estimated numerical uncertainty in the results, obtained by averaging over five realizations, is smaller than the size of the symbols used in the figures. Throughout the first stage of hydration, the hydration energy of clays with only Na < decreases significantly below the internal energy of bulk flexible SPC water, −41.5 kJ/mol. Due to the strong surface−ion interaction in I−MMT with only Na < , the minimum energy is relatively constant for a wide range of water content, when compared with MMT−MMT that responds more rapidly to the added water. At the same time, the broad energy minimum in the case of Na < with the I−MMT surface is due to the higher number of cations, and hence, the charge magnitude required to charge compensate the illitic surface leads to increased water content in the interlayer. This indicates that it is energetically more favorable to add water to clays with tetrahedral charge sites than to those with octahedral ones. The lack of change in the minimum hydration energy of I−MMT with Na + becomes clearer by considering the density profiles, which were presented in the main text of the paper. In the cases with K < , however, ΔU(N) decreases less rapidly. Moreover, in clays with only K < in the interlayer regions the limited interaction of the cations with water results in reduced energy change, which is expected. Thus, in the presence of K < the water molecules should be repelled by the clay. As soon as the initial energy barrier is crossed, the hydration energy of clays with K < also decreases, but to values above that of clays with Na < . Figure 2.2. Hydration energy as a function of the water content for MMT−MMT (blue), I−MMT with Na + only (yellow), I−MMT with 72% K + and 28% Na + (black), and I−MMT with K + (red). -55 -53 -51 -49 -47 -45 -43 -41 -39 0 25 50 75 100 125 150 175 200 ΔU(N)(KJ/kmol) Number of Water Molecules 15 A comparison of the four cases in Figure 2.2 reveals that it is energetically more favorable for water molecules to invade the anhydrous Na < clay. In addition, the more Na < in the interlayer region, the more the clays swell. On the other hand, keeping the interlayer cations the same, the more charges in the tetrahedral layer, the more stable and energetically favorable the material. Thus, in clays with the same structures, the more Na < are in the interlayer region, the more energetically favorable it is for the system to swell. Note, however, that the number of Na < in I−MMT is 25, compared with 14 in MMT− MMT, implying that although there are more cations in the former case with increased hydration as water molecules are inserted, MMT−MMT responds more sensitively to them due to the weaker surface−ion interactions. After the minimum energy is reached in each case, the energies rise as more water is added, close to the internal energy of bulk water. 2.3.2 Swelling The basal spacing d is defined as the sum of the clay sheet’s thickness and the interlayer distance between that sheet and its neighboring one. For I−MMT MCL all the d values were calculated based on MMT’s thickness. d was computed as a function of the number of water molecules per unit cell for I−MMT, MMT−MMT, and I−I for 10 separate cases. In the case of MMT−MMT only Na < was inserted into the interlayer region, whereas three cases of I−MMT, with K < plus Na < , with K < only, and with Na < only in the interlayer region, were also studied. In the I−I interlayer only K < was inserted in all the cases, which resulted in the layers being “locked” together, with its basal spacing d remaining fairly constant, ≈10 Å, in all the cases. As mentioned earlier, the cation content of the MMT−MMT interlayer is identical in all the simulations but different from that of the I−MMT interlayer in which we inserted, respectively, 18 and 7 K < and Na < symmetrically distributed in the y−z plane (see Figure 2.1), in order to offset the negative charges of Illite and MMT. The dependence of the absolute basal spacing d on the number of water molecules for MMT−MMT, I−MMT with 72% K < and 28% Na < , I−MMT with all K < substituted by Na < , and I−MMT with all Na < substituted by K < is shown in Figure 2.3. The overall swelling behavior of MMT−MMT clay agrees well with the existing literature on both simulation and experimental studies [36, 89, 90]. 16 Figure 2.3. Basal spacing d of (a) MMT−MMT; (b) I−MMT with Na + ; (c) I−MMT with 72% K + and 28% Na + , and (d) I−MMT with K + . Table 2.1. Dependence of the Basal Spacing d (in Å) on the Type of Interlayer Cations in Various Hydration States 0W 1W 2W MMT-MMT (𝐍𝐚 < ) 9.4 12.2 15.0 ILLITE-MMT (0.72 𝐊 < , 0.28 𝐍𝐚 < ) 9.8 12.1 14.8 ILLITE-MMT (𝐍𝐚 < ) 9.4 11.9 14.0 ILLITE-MMT (𝐊 < ) 9.9 12.2 14.9 Water adsorption in clays occurs in two steps. First, some adsorption occurs to form the hydration shell for the counterions that results in the layer expansion, whereas in the second step the increase in the water content fills up the remaining volume of the interlayer without much swelling. For cations with lower hydration enthalpy, however, due to the lower extent of hydration and the disturbed hydration shell, swelling happens more monotonically without distinguishable steps. The basal spacing d for the dry (0W), 1W, and 2W states of the various clays is presented in Error! Reference source not found. In the dry 0W state I−MMT with only K < has the largest d value, which is due to the larger ionic radius [91] of K < , 1.52 Å, than Na < , 1.2 Å. 8 10 12 14 16 18 0 25 50 75 100 125 150 175 200 d (Å) Number of Water Molecules (a) 8 10 12 14 16 18 0 25 50 75 100 125 150 175 200 d (Å) Number of Water Molecules (b) 8 10 12 14 16 18 0 25 50 75 100 125 150 175 200 d (Å) Number of Water Molecules (c) 8 10 12 14 16 18 0 25 50 75 100 125 150 175 200 d (Å) Number of Water Molecules (d) 17 Comparing the panels of Figure 2.3 reveals that MMT−MMT swells most in all the cases, and that its interlayer space manifests distinct 1W and 2W hydration states and the transition to them. Note also that there is a wider and clearer range for the 1W state than for the 2W state. The variations of d for I MMT are, however, significantly different from those of MMT−MMT with both K < and Na < . In addition, for all the water contents, swelling of I−MMT is smaller than MMT−MMT. Due to the two-step hydration of clays, the expansion due to formation of the hydration shell around the cations is smaller for both I−MMT cases, including the case with K < . Figure 2.3 also indicates that, compared with the other cases, I− MMT with only Na < has an intermediate level of expansion and that, compared with MMT−MMT, the flattened 1W state for I−MMT with Na < is postponed to higher water contents. This is due to the interaction between Na < and the more negatively charged tetrahedral layer on the illite side. We will return to this shortly. Figure 2.4. Comparison of the absolute(a) and relative(b) basal spacings as a function of the water content for MMT−MMT (blue), I−MMT with Na + only (yellow), I−MMT with 72% K + and 28% Na + (black), and I−MMT with K + only (red). Swelling in I−MMT with 72% K < and 28% Na < and one with only K < happens more smoothly in the interlayer with a smaller change of slope in the transition from the 0W to 1W state. The two MLCs that contain K < have the lowest extent of expansion. The weak transition between the hydration states in the presence of K < is due to weaker attraction between K < and the water molecules. Denis et al. [21] reported that in mixed-ion MMT K < fractions of over 33% strongly control swelling. In the present study the response of the I−MMT interlayer to both cases containing K < is very similar, presumably because the K < fraction is relatively high in both cases. To better appreciate the differences between the basal spacing of the four types of clays, both pure and mixed that we have studied, Figure 2.4 compares the relative basal spacings d as a 0 1 2 3 4 5 6 7 0 25 50 75 100 125 150 175 200 Relative d-spacing (Å) Number of Water Molecules (b) 8 10 12 14 16 18 0 25 50 75 100 125 150 175 200 d-spacing (Å) Number of Water Molecules (a) 18 function of the water content. To compute the relative basal spacing, we subtracted its value from that in the dry state. Figure 2.4 reveals that MMT−MMT has the largest amount of swelling in all cases and that its interlayer space manifests distinct 1W and 2W hydration states, as well as the transition to them. Note also that there is a wider and clearer range for the 1W state than the 2W state. The variations of the basal spacing for I−MMT are, however, significantly different from that of MMT−MMT with both K < and Na < . In addition, for all the water contents, swelling of I−MMT is smaller than that of MMT−MMT. According to the two-step hydration of clays that was described earlier, the expansion due to formation of the hydration shell around the cations is smaller for both I− MMT clays, including the case with K < . Figure 2.4 also indicates that I−MMT with only Na < has an intermediate level of expansion when compared with the other cases and that, compared with MMT−MMT, the flattened 1W state for this MLC occurs at higher water contents, which is due to the interaction between Na < and the more negatively charged tetrahedral layer on the illite side. We will return to this shortly. 2.3.3 Molecular Structure of the Interlayer Region To gain deeper insight into the interactions of K < and Na < with the clay surface, the distributions of the cations and water in the interlayer regions in the z direction, perpendicular to the sheet, were calculated. The results are shown in Figure 2.5. Our computed 1W and 2W distributions for MMT−MMT agree well with the literature [17, 36, 92-96]. We should keep in mind, however, that in this type of molecular modeling some inner sphere complexes are formed between Na < and the clay surface, which are due to the distribution of the negative charges on the tetrahedral layer as well as the overestimated surface−ion interaction for the surface oxygens incorporated in the CLAYFF. The resulting swelling is slightly smaller than the experimental data that for the basal 0 1 2 3 4 5 6 7 0 25 50 75 100 125 150 175 200 Relative d-spacing (Å) Number of Water Molecules (b) 8 10 12 14 16 18 0 25 50 75 100 125 150 175 200 d-spacing (Å) Number of Water Molecules (a) 19 spacing of the 1W state are in the range 12.3−12.6 Å and for the 2W state in the range 14.9−15.6 Å for MMT−MMT. Figure 2.5. Atomic density profile for the 1W and 2W states. (a,b) MMT−MMT; (c,d) I−MMT with Na + ; (e,f) I−MMT with 72% K + and 28% Na + , and (g,h) I−MMT with K + , for water−hydrogen (blue), water−oxygen (red), Na + (yellow), and K + (gray). 0 2 4 6 8 10 0 1 2 3 4 5 6 Total Atom Count Interlayer Distance (Å) (a) 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 Total Atom Count Interlayer Distance (Å) (b) 0 2 4 6 8 10 0 1 2 3 4 5 6 7 Total Atom Count Interlayer Distance (Å) (c) 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 Total Atom Count Interlayer Distance (Å) (d) 0 2 4 6 8 10 0 1 2 3 4 5 6 7 Total Atom Count Interlayer distance (Å) (e) 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 Total Atom Count Interlayer Distance (Å) (f) 0 2 4 6 8 10 0 1 2 3 4 5 6 7 Total Atom Count Interlayer Distance (Å) (g) 0 2 4 6 8 10 0 1 2 3 4 5 6 7 8 9 10 Total Atom Count Interlayer Distance (Å) (h) 20 As Figure 2.5 indicates, asymmetry of the atomic distributions in I−MMT occurs due to two different clay layers facing each other. Due to the low hydration energy of K < , independent of the water content, the cations bind strongly to the illite surface. In fact, they are adsorbed onto the surface as the inner sphere surface complexes (ISSCs) and are locked inside the surface’s hexagonal cavities. A previous study of K + −MMT by Boek et al. [18] indicated that adding water to the interlayer results in K < leaving the hexagonal sites and moving to the sites above the SiO A tetrahedral. Boek et al. also noted that, compared with Na < and Li < , K < is reluctant to fully hydrate and attach itself to the surface; see also Boek [97] and Chang et al. [98] In the present case there are two main reasons for K < being tied to the illite surface. One is exposure to two different clay surfaces, with illite being more negatively charged on the outer layer, which attracts immediately all the positively charged ions toward the surface. The second reason is the concentrated tetrahedral charge of illite that, compared with MMT and its more negatively charged octahedral sheet, makes the ionic bonding stronger. Interpretation of the results should be done carefully in light of the complex behavior of Na < . If we consider I−MMT, the most diverse surface complexes are formed in the 2W of I−MMT with Na < . There are two and one type of the ISSCs formed, respectively, on the illite and MMT sides of the interlayer region. Moreover, one type of outer-sphere surface complex (OSSC) is also formed in that region. Thus, with only Na < on the MMT side, the ISSCs exist in both the 1W and 2W states. With K < in the interlayer on the MMT side of I−MMT, however, the single low-density peak of Na < indicates that the ISSC forms at higher water contents. With mixed cations in the interlayer, only one type of ISSC and no OSSC of Na < are formed. The OSSC in I−MMT with Na+ explains the difference in the basal spacings of the three I−MMT studied, discussed earlier. When only Na < is present, some of them near the illite surface hydrate as the water content increases, giving rise to formation of the OSSCs similar to that in MMT−MMT that contribute to more swelling. Figure 2.5 indicates that the locations of the first density peak of both Na < and K < ISSCs are much closer to the illite surface than expected. In the case of the ISSC it is common for the location of the first peak to be the sum of the ionic radii of the cation and surface oxygen (1.26 + 1.16 = 2.42) Å for rO−Na and (1.26 + 1.52 = 2.78) Å for rO−K. Figure 2.5 indicates, however, that the location is around 0.8 Å for rO−Na and 1.5 Å for rO−K. The larger radius for K < is consistent with the ionic radius of K + that is larger than that of Na < . The large difference in the locations of 21 the density peaks illustrates the strong adsorption of the cations deep in the ditrigonal cavities of the siloxane surface. The comparably smaller penetration of K < into the illite surface implies that the ditrigonal cavities are better fits for Na < than K < Similar behavior was reported by Lammers et al. [99] In their study, however, the distance from the surface for Na < was slightly larger than what we report here, which is likely due to the difference in the chemical composition of illite that has more negative charges at its surface. Moreover, the asymmetry of the interlayer structure with the higher charge density on the Illite side contributes to its dominant role in the adsorption of the ions. The two mechanisms cause the stronger attraction of the ions to the ditrigonal cavities of the tetrahedral layer. Note that previous MD simulation of Na-smectite produced a sodium density peak that was considerably farther from the surface. Bourg et al., for example, reported that in a rigid clay the Na−ISSC is at 2.6 Å from the surface. The reason for the difference is [99] the substitutes that are distributed more in the tetrahedral layer of illite, as well as the flexible clay structure. Deep penetration of the ions into the hexagonal cavities is due to an unphysical relaxation of tetrahedral layer, expanding the ditrigonal cavities to larger hexagonal ones. The OSSCs are formed when the ions are at distances larger than 3 Å from the clay surface [99]. They are present in I−MMT only in its 2W state with Na < , while they do not form in an MLCs with only K < . This confirms our results presented above and implies that K < is primarily adsorbed as the ISSC, while Na hydrates at higher water contents. The gap in the energy density of the ISSCs and OSSCs also indicates the energy barrier of ion exchange at the clay surface. The density profile of water oxygen is also affected by the counterions. On the MMT side at lower water contents, hydrogen atoms of water bond to the oxygen on the clay’s surface and form the 1W layer. As the water content increases, the 2W state is formed, so that hydrogen preserves its distance from the surface. On the illite side, however, when K < is in the interlayer region, it only perturbs the water network. Therefore, it contributes to making the surface hydrophobic, which is the reason for not having a distinct water oxygen peak near the illite surface. Nevertheless, in terms of the water density distribution, the 2W state in the Na−I−MMT is analogous to that of MMT−MMT. 2.3.4 Radial Distribution Functions Figure 2.6 presents the calculated radial distribution functions (RDFs) g(r) of MMT−MMT and I−MMT with both Na < and K < , as well as with only Na < , while Figure 2.7 shows the results with 22 K < . In all cases the locations r t u −K and r t u −Na of the peaks corresponding to hydration are, respectively, 2.78 and 2.4 Å, which are in agreement with the previous studies. As the water content increases, the height of the peak in g(r) decreases in all cases. A comparison between Figure 2.6 and Figure 2.7 indicates that the peak corresponding to OW − K hydration is smaller than that of OW−Na, which is due to the difference in the tendency of Na < and K < to hydrate. Figure 2.6. Radial distribution functions for Na−OW (left) and Na−OS (right) for (a,b) MMT−MMT, (c,d) I−MMT with 72% K + and 28% Na + , and (e,f) I−MMT with Na + . Figure 2.6 and Figure 2.7 also show the RDFs of the counter-ions and the surface-bridging oxygen OS. In MMT−MMT with Na < with the negative charges concentrated in the octahedral sheet, the hydration shell of Na < is more stable than that of K < , as the hydration energy of the former is higher (see Figure 2.2). Therefore, as the water content increases, the cations move toward the center of the interlayer region and hydrate. Figure 2.6 (b) confirms this as it indicates that the height of the RDF peak decreases considerably in the 2W bilayer region. In the case of I−MMT with Na < , however, due to the strong adsorption of the cations onto the ditrigonal cavities, 0 5 10 15 20 25 1 2 3 4 5 6 7 8 g(r) distance (Å) monolayer bilayer (a) 0 0.3 0.6 0.9 1.2 1 2 3 4 5 6 7 8 9 g(r) distance (Å) monolayer bilayer (b) 0 0.4 0.8 1.2 1.6 1 2 3 4 5 6 7 8 9 g(r) distance (Å) monolayer (f) 0 4 8 12 16 1 2 3 4 5 6 7 8 g(r) distance (Å) monolayer bilayer (e) 0 0.4 0.8 1.2 1.6 1 2 3 4 5 6 7 8 9 g(r) distance (Å) monolayer bilayer (d) 0 6 12 18 1 2 3 4 5 6 7 8 g(r) distance(Å) monolayer bilayer (c) 23 changing the water content in the interlayer region does not change considerably the location of the ions near the surface. The same behavior is also manifested by I−MMT with K < in the interlayer region (see Figure 2.7). Comparison of the results for I−MMT with 72% K < and 28% Na < with those for I−MMT with Na < indicates that the latter swells more than the former. According to Figure 2.6 (c) and Figure 2.6 (e), when K < is in the interlayer region, Na < ions are surrounded more by water, while fewer Na < are coordinated near the surface. It should, however, be noted that in the cases with Na only more swelling occurred because there is nearly 3.6 times more Na < in the interlayer, compared with the cases with mixed cations. With regard to Figure 2.6 (d) and Figure 2.6 (f), the sharper Na−OS peak in the case with only Na is caused by two factors: First, as mentioned earlier, there are two types of the ISSC in I−MMT with Na < Second, there are more Na < in the interlayer. Figure 2.7. Radial distribution function of K−OW (left) and K−OS (right) for I−MMT with 72% K + and 28% Na + (a and b) and for I−MMT with K + only (c and d). Figure 2.7 indicates that the RDFs of the MLCs with only K < and a mixture of K < and Na < are quite similar, confirming the results for the basal spacing d and the density profiles, presented in Figure 2.3-Figure 2.5. Due to the lower hydration enthalpy of K < , however, the peaks at r t v −K, which are near the illite surface, are sharper and more distinct. In the case of I−MMT, the very slight change in the RDFs of the OS – counterion confirms their dominant distribution near the illite surface. Otherwise, the dynamic behavior of the cations near the MMT surface would have resulted in substantial changes in the RDFs with increasing water content. In contrast with 0 2 4 6 8 10 1 2 3 4 5 6 7 8 g(r) distance (Å) monolayer bilayer (a) 0 0.5 1 1.5 2 2.5 1 2 3 4 5 6 7 8 9 g(r) distance (Å) monolayer bilayer (b) 0 0.5 1 1.5 2 2.5 1 2 3 4 5 6 7 8 9 g(r) distance (Å) monolayer bilayer (d) 0 2 4 6 8 10 1 2 3 4 5 6 7 8 g(r) distance (Å) monolayer bilayer (c) 24 Na < , K < is more strongly bonded to the surface oxygens. Figure 2.7 also confirms that when the ratio K < /Na < is larger than a certain limit K < will be the dominant counterion in the interlayer, results in similar swelling, density profile, and RDF in the cases of I−MMT with 72% K < and 28% Na < . 2.4 Summary We presented the results of what we believe to be the first MD simulations of hydration energetics and swelling of mixed-layer clays, consisting of montmorillonite and illite sheets. The effect of two important cations, Na < and K < , on the hydration energies and swelling of the MLCs was studied in detail. Very significant differences were demonstrated between such properties of pure clays, such as montmorillonites and illite, and the mixed structures of the same two clays studied in this paper. As such clays constitute the majority of natural clays, it is imperative to understand their properties, particularly in the presence of CO : , which is highly relevant to the problem of its sequestration in natural porous formations. The results are discussed in the following chapter. 25 3 Molecular Dynamics Simulation of Hydration and Swelling of Mixed- Layer Clays in the Presence of Carbon Dioxide 3.1 Introduction Chapter 2 described our molecular-scale results for the effect of water on swelling and ion distribution in mixed-layer ML clays. In this chapter we present our results for the case when both water and carbon dioxide are present in the clays. We begin by describing the molecular dynamics (MD) procedure. 3.2 The Force Field and Details of Simulation Procedure The I-MMT mixed-layer representation is described comprehensively by Altaner and Ylagan [15]. We constructed molecular models for the illite-montmorillonite (I-MMT) mixed clay. The chemical composition of the illite is K0.9[(Al0.9Mg0.1)(Si3.2Al0.8)O10(OH)2], while that of the sodium-saturated montmorillonite, Na-MMT, is, Na0.75[Si7.75Al0.25](Al3.5 Mg0.5)O20(OH)4. Each MMT layer has a net negative charge of -0.75 eV per unit cell, while the corresponding net charge for illite is -1.8 eV. The substitutions were inserted symmetrically in the clays’ sheets. Figure 3.1 presents the side view of the clay layers, water, CO2, and the counterions in the supercell structure of the ML clay, which contains three distinct interlayers. The top and bottom layers are illite, whereas the two central ones are made of MMT. Therefore, the structure of the MLC is symmetric with respect to the central line located between the MMT sheets. The I-I and MMT- MMT interlayers contain, respectively, K + and Na + in order to balance the pure layers’ charge. The supercell structure contained a total of 80-unit cells. Each layer was constructed with 20-unit cells, organized as 5×4 unit cells along the x and y directions; see Figure 3.1. The dimensions of the illite and MMT layers were, respectively, 5 × 9 Å : and 5 × 8.8 Å : . The equilibrium dimensions of the individual sheets in the illites and MMTs were, respectively, 26 × 35 Å : and 26 × 36 Å : . Periodic boundary conditions were imposed on the molecular structures. 26 Figure 3.1. Side view of the simulation cell for I-MMT. The top and bottom layers are illite, while the two middle layers are the MMT. We used the force field CLAYFF [30] and the LAMMPS package in order to generate molecular models of the ML clays, as well as Na-MMT that is used as a reference in order to understand the differences between ML clays and pure clays and carried out the MD simulations. The flexible SPC model and the CLAYFF parameters for the hydroxyl group in the octahedral sheets were used to model water and the O-H groups in the supercell structure, respectively. The associated parameters, as well as those for the optimized force field parameters for CO2 were estimated by Cygan et al [100]. The fully flexible model for CO2 that we utilized in the simulations produce predictions that agree with the experimental data [101], whereas the rigid models could not capture [102] the deformation of the molecules due to vibration. As mentioned in Chapter 2, according to the CLAYFF the total interaction energy E of the material is given by, Here, k1 and k2 are force constants, θijk is the angle between bonds ij and jk, e is the electron’s charge, qi is the partial charge of atom i, ε0 is the dielectric permittivity of vacuum, and r0 and θ0 27 are the equilibrium values of the corresponding quantities. i and σi are the usual Lennard-Jones (LJ) energy and size parameters, and the Lorentz-Berthelot mixing rules were used for the pairs ij. Table 3.1 presents the optimized parameters of CO2 for both the bond and non-bond interactions that we utilized in the simulations. Table 3.1. The CO2 parameters used in the MD calculations [100] Nonbond Parameters q x +0.6512 e q t -0.3256 e ε x 0.2340 kJ/mol ε t 0.6683 kJ/mol σ x 2.800 Å σ t 3.028 Å Bond Parameters k xt 8443 kJ/mol Å 2 r °xt 1.162 Å k txt 451.9 kJ/mol rad 2 θ °txt 180.0° Before carrying out any simulation involving water and CO2, the molecular model of dry clay was generated in the NPT ensemble, which took 3 ns to reach equilibrium, after which the simulations involving water and CO2 began for each case that we studied at T = 348 K and P = 130 bar. The Nos´e-Hoover thermostat and barostat were employed to control T and P with the coupling constants of 1000 and 100, respectively. After some preliminary simulations, the cut- off distance for computing the non-bonded interactions was set at 13 Å. Therefore, long-range dispersion correction of the pressure and energy was enforced, and the long-range electrostatic interactions were computed with an accuracy of 10 −4 using the Ewald summation technique. Simulations with a time step of 0.01 fs for 1000 ps were followed by longer simulations with a time step of 1 fs. The duration of the simulations in the NPT ensemble in the presence of both water and CO2 varied between 25 and 40 ns, considerably longer than those used in our previous work [103] involving only water. Next, the simulation continued in the NVT ensemble in order to calculate the density profiles, the radial distribution functions, and the diffusion coefficients 28 (see below). After equilibrium was reached, the results were computed and averaged over the last 10 ns of the simulations. In all the simulations, the counterions were initially distributed at the midplane of each interlayer, whereas the initial spatial distributions of water and CO2 molecules were random. For each case involving the I-MMT interlayers, three series of simulations were carried out. One was with Na + only; the second case was with K + only, while the third involved a combination of two cations with 72% of the total being K + and the remaining 28% being Na + . To better understand the role of cations in the ML clays swelling and to have a basis for comparison with pure clays, two other sets of simulations were carried out in which either pure K + or Na + was present in the I-MMT interlayer, in order to neutralize the net negative charge of the MLclays. The CO2 content was varied from 0.5 to 4 molecules per unit cells, while the water content varied between one to nine molecules per unit cell. Note that since pure illite does not swell, there is only K + in the I-I interlayer, which locks the two adjacent illite layers together. Therefore, the water-CO2 mixture is only added to the MMT-MMT and I-MMT interlayers. The main purpose of our study is to understand the effect of mixtures of water and CO2, as well as the effects of the initial water concentration, the interlayer cations and the surface charges on the swelling of I-MMT ML clays and its differences with swelling of pure Na-MMT. As we describe below, since the combined number of water and CO2 molecules in the interlayer is large (up to 780 atoms of water and CO2 in each interlayer region), the 3W hydration state, in addition to the one layer (1W) and two layer (2W) states, also develops in some of the MLCs with high concentrations of water. 3.3 Results and Discussion We computed several important and characteristic properties of the MLCs in the presence of water and CO2. In what follows we present and describe the results and discuss their implications. 3.3.1 Swelling and Basal Spacing We define basal spacing d as the sum of the clay sheet’s thickness and the interlayer distance between that sheet and its neighboring ones. In the case of I-MMT all the basal spacings were calculated based on MMT’s thickness. As already pointed out, in the case of the I-I interlayer the 29 basal spacing d of pure illite remains constant about 10 Å in all the simulations that we carried out. Figure 3.2. Absolute basal spacings d of (a) MMT-MMT; (b) I-MMT with Na + ; (c) I-MMT with 72% K + and 28% Na + , and (d) I-MMT with K + . The number of CO2 molecules is per unit cell. Figure 3.2(a) presents the computed basal spacing d of pure MMT interlayer as a function of the number of intercalated water and CO2 molecules. Also shown are the results without any CO2, which were described in our previous study [103] and are included here as a reference for comparison. The computed basal spacings are very close to the experimental and computational results reported previously [31, 33].The small differences between the computed results and some of the experimental data may be due to slightly different original structure of the MMTs as a result of a slight difference in the chemical compositions. In addition, given that no force field is exact, some of the minor differences may be attributed to the CLAYFF. Water is adsorbed by clays in two steps. First, some water is adsorbed and forms the hydration shell for the counterions. This results in the expansion of the layer. In the second step, as the water content increases, it fills up the remaining volume of the interlayer without much swelling. As Figure 3.2 (b-d) demonstrate, in almost all cases the more CO2 is added to the interlayer, the larger is the extent of swelling and the expansion of the interlayer region, when they are compared with the reference cases with water only. Figure 3.2 also indicates that the rate 10 12 14 16 18 20 0 20 40 60 80 100 120 140 160 180 200 d (Å) spacing number of water molecules 0 10 20 40 60 80 (a) 10 12 14 16 18 20 0 20 40 60 80 100 120 140 160 180 200 d (Å) spacing number of water molecules 0 10 20 40 60 80 (c) 10 12 14 16 18 20 0 20 40 60 80 100 120 140 160 180 200 d (Å) spacing number of water molecules 0 10 20 40 60 80 (b) 10 12 14 16 18 20 0 20 40 60 80 100 120 140 160 180 200 d (Å) spacing number of water molecules 0 10 20 40 60 (d) 30 of the increase in the basal spacing is larger at the beginning of the transition between the hydration states. The rate reduces, however, at the end of the transition when the expanded volume is filled with more water and CO2. Figure 3.3. Dependence of the relative basal spacings on the water content at a fixed number of CO2 molecules per unit cell. (a) 10 CO2; (b) 20CO2; (c) 40CO2; (d) 60CO2, and (e) 80 CO2 molecules. Various combinations of water and CO2 give rise to three hydration states, namely, the aforementioned 1W, 2W and 3W states. In absence of CO2, the 1W-to-2W transition in the MMT- MMT interlayer begins with 80 water molecules, with the 2W state becoming fully saturated when there are 160 water molecules in the interlayer region. As shown in Figure 3.2 (a), as more CO2 was inserted into the interlayer region, the steep step signifying the transition was shifted to 4 6 8 10 0 20 40 60 80 100 120 140 160 180 200 d (Å) Number of water molecules MMT-MMT Na+I-MMT K+I-MMT (K+Na)-I-MMT (e) 2 4 6 8 10 0 20 40 60 80 100 120 140 160 180 200 d (Å) Number of water molecules MMT-MMT Na+I-MMT K+I-MMT (K+Na)-I-MMT (d) 2 4 6 8 0 20 40 60 80 100 120 140 160 180 200 d (Å) Number of water molecules MMT-MMT Na+I-MMT K+I-MMT (c) 0 2 4 6 8 0 20 40 60 80 100 120 140 160 180 200 d (Å) Number of water molecules MMT-MMT Na+I-MMT K+I-MMT (K+Na)-I-MMT (b) 0 2 4 6 8 0 20 40 60 80 100 120 140 160 180 200 d (Å) Number of water molecules MMT-MMT Na+I-MMT K+I-MMT (a) 31 lower water contents. For example, with 40 CO2 molecules the transition begins with 40 water molecules also. Since the molecular size of CO2 is larger than water’s, swelling in some of the cases with similar extents happens at lower CO2 content. It should also be noted that the transition between the hydration states happens more smoothly at higher CO2 concentrations; see the discussion of Figure 3.3 below. In general, however, and independent of the CO2 concentration, the transition between the 1W and 2W hydration states is sharper than that between 2W to 3W state at a constant CO2 concentration. This may be due to the small changes in the interactions between the interlayer species and the surface during the transition between 2W and 3W states, which also result in the 1W and 2W states being more stable states than the 3W. The computed basal spacings of the MMT-MMT in the 1W, 2W, and 3W hydration states, when both water and CO2 are present are, respectively, in the ranges 11.5−12.5 Å, 14−15 Å, and 17.5−19 Å. These are completely consistent with the experimental data [55, 104]. As Figure 3.2 indicates, for I-MMT the transition between the first two hydration states is delayed to higher water and CO2 contents, especially in the presence of K + as the counterion. The delay in the transition is caused by the asymmetry in the interlayer charges, as well as the effect of the hydration enthalpy of K + on swelling, which was alluded to earlier. To gain a deeper understanding of the differences between the basal spacing of the four types of clays, both pure and mixed that we have studied, Figure 3.3 compares the relative basal spacings d as a function of the water and CO2 content. The relative basal spacings were computed by subtracting their absolute values shown in Figure 3.2 from those of their dry states. At low CO2 content the order of the relative swelling is MMT-MMT > Na-I-MMT > K-I-MMT and (K+Na)-I-MMT. The relative basal spacings at high CO2 contents, shown in Figure 3.3(d) and Figure 3.3(e), indicate that the interlayers expand similarly, with little difference between their relative d values. Figure 3.3 also demonstrates the disruption of the hydration shells by CO2. The distinguishable sharp transition between the 1W and 2W states in (a) is gradually shifted to the more similar and continuous expansion for all the cases in Figure 3.3 (e). The reason for the difference is that high concentrations of CO2 disrupt the hydration shells. In addition, the reduced hydrogen bonding, as well as the hydrophobicity of the clay surfaces, results in similar trends in the swellings. In fact, as CO2 is added, it escapes the more charged surface and disrupts the water 32 hydration shells, which results in comparable basal spacing for all the cases that was pointed out earlier. We will return to this important point when we describe the results for the density profiles of the hydration states of the various cases in the next section. 3.3.2 Density Profiles The density profiles were calculated for all the cases that we studied in order to understand the spatial distribution of the species in the interlayer region. The results represent averages along the direction perpendicular to the clay surfaces. As mentioned earlier, each of the hydration states develops in a range of water and CO2 composition. The monolayer was developed by 2H2O + CO2, 3H2O + CO2, 4H2O + CO2, etc. However, the density profiles and the extent of the saturation were slightly different for various combinations of the number of water and carbon dioxide. The bilayer was developed for a wider range of water plus CO2, beginning from 2H2O + 3 CO2 and 3H2O + 3CO2, up to the fully saturated 2W state with 8H2O + CO2 and 9H2O + 0.5CO2 (per unit cell). The trilayer emerged at high water and CO2 concentrations, and is expected to still exist at concentrations larger than what we have considered in this study. It was first developed with 5H2O + 4CO2, until it was fully saturated with 7H2O + 3CO2 and 9H2O+4CO2. The results are presented in Figure 3.4. The results for MMT-MMT, shown in Figure 3.4 (a-c), should be compared with those reported by Myshakin et al. [33] ,who reported the 3W state with 3 peaks for CO2 and 2 peaks for water, which does not, however, represent the fully saturated 3W state. Nevertheless, as mentioned earlier, in our simulations the 3W state emerged at low water and CO2 concentrations with 5H2O + 4CO2, with the fully saturated 3W state at higher contents. The reason is that clays need to be wetted to intercalate CO2. In the work of Myshakin et al. [33] the CO2 concentration was high in the density profiles that they reported, which explains why more water was needed in order for the clay to expand to the fully saturated 3W state. In our study, however, water concentration is up to 2.5 times that of CO2 in the fully saturated 3W state. The extent of the swelling is, however, similar in both studies, if we compare the basal spacings. Our results for the density profiles for the 1W and 2W states in the MMT- MMT clays agree with the previous results [31, 33, 34, 104]. 33 Figure 3.4. Density profiles of the 1W, 2W, and 3W hydration states, perpendicular to the clay surfaces of MMT-MMT, shown in (a), (b), and (c), and in Na+I-MMT shown in (d), (e), and (f). The results for water oxygen are represented by red, CO2 oxygen and carbon by black and blue, and Na + in green. In the I-MMT interlayer the left layer is illite and the right layer is MMT. According to the density profiles in MMT-MMT and Na+I-MMT that are presented in Figure 3.4, the distributions of water and CO2 in the case of monolayers are very similar. Due to the stronger negative charge on the illite surface, the Na + cations are distributed close to the surface in I-MMT where they form inner sphere surface complexes (ISSCs), and where the basal oxygens also contribute to their first hydration shell. In MMT-MMT, however, Na + is distributed predominantly in the outer-sphere surface complexes (OSSCs) with water. The 1W hydration state for the two interlayers is very similar to that of pure water [103]. 34 Figure 3.5. Density profiles of 1W, 2W and 3W states perpendicular to the clay surfaces of I- MMT with 72% K + and 28% Na + , shown in (a), (b), and (c), and K+I-MMT, shown in (d), (e), and (f). The results for water oxygen are represented by red, CO2 oxygen and carbon by black and blue, Na + by green, and K + by purple. In the I-MMT interlayer the left layer is illite and the right layer is MMT. Although when both water and CO2 are present, the density profiles in the bilayer region of MMT-MMT are analogous to the same in the MMT-MMT but with only water [103], the same profiles for Na+I-MMT depend strongly on whether CO2 is present. When CO2 intercalates in the interlayer region, it prefers to migrate away from the clay surface substitutions. Therefore, in the MLC CO2 is at the farthest point from the more charged surface, namely, illite’s surface, and is distributed near the MMT surface. This results in more pronounced peaks in the density profiles of CO2 and water near the MMT and illite surface, respectively, with a weaker CO2 peak near the 35 illite surface. For Na + the distribution near the illite surface is similar to the case with only water. Close to the MMT surface, however, inserting CO2 in the MLCs splits the Na + density profile into two peaks of the ISSC adsorption, including one deep into the hexagonal cavities [99] of the clay and above the hexagonal ring of the tetrahedral sheet. Similar behavior occurs in the trilayer region with much stronger CO2 peaks near the MMT surface. Figure 3.5 compares the density distribution of oxygen in water, CO2 and the interlayer cations in the three hydration states of (Na+K) +I-MMT and K +I-MMT. We first note that, with the exception of the presence of Na + in the former case with a notable peak near the illite surface, the density profiles in the two clays are quite similar. Moreover, the distributions of Na + and K + are nearly independent of the amount of intercalated water and, therefore, the hydration states. A comparison of the density profiles for the two cases that contain both water and CO2 with those that contain only intercalated water in the 2W state indicates clear differences. In the bilayer states in the former cases the sharper water density peaks are near the MMT surface, whereas in the latter case they are near the illite surface. There are at least two reasons for the hydrophobicity of the surfaces in the cases with and without CO2. In the absence of CO2 the cations make the illite surface hydrophobic and contribute to the migration of water molecules to the MMT surface, where they form hydrogen bonds with the basal oxygen atoms. On the other hand, in the presence of CO2 and due to its accumulation near the MMT surface, the extent of hydrophobicity of the MMT surface is stronger than that of the illite surface. Therefore, water molecules are inclined to distribute themselves on the opposite side of the MLC interlayer. This behavior, together with the disruption of the water network by CO2, also explains the trends in the basal spacing d with no sharp jump between the states at high CO2 concentration (40 molecules and larger) in all the interlayers. 3.3.3 Radial Distribution Functions More insight into the distribution of the species in the interlayer region is gained by the radial distribution functions (RDFs) g(r) of a few pairs of the atoms. Figure 3.6 presents the results for the Na-OW (water oxygen) and Na-OS (surface oxygen) in the three interlayers containing the Na cations. In each case the water concentration was increased from 3 to 8 molecules per unit cell at constant CO2 concentration, which included the transition from the monolayer to the fully saturated bilayer in each interlayer. 36 Figure 3.6. Radial distribution functions for Na-OW (left) and Na-OS (right) at constant CO2 concentrations in (a) and (b) MMT-MMT; (c) and (d) I-MMT with Na + , and (e) and (f) I-MMT with 72% K + and 28% Na + According to Figure 3.6(a), as the water content increases, the height of the peaks decreases. Note, however, that the peaks remain identical if we increase the water content of the clays with 80 and 100 water molecules. This may be understood if we consider Figure 3.6(a-b), together with the results for the basal spacings presented earlier. The largest increase in the basal spacings in Figure 3.2(a) occurs when the water content increases from 80 to 100 molecules per unit cell. As explained earlier, the first step of hydration takes place when counterions are hydrated and transition from the ISSC- to the OSSC-dominated adsorption. Figure 3.6(a)-3.6(b) demonstrate this clearly. As Figure 3.6(b) indicates, the RDF peak with between 80 and 100 water molecules decreases significantly, hence indicating migration of more Na + from the surface 37 toward the added water. At the same time, the equal heights of the two cases in Figure 3.6(a) confirm that the extra water has formed a hydration shell around the Na + cations. Figure 3.6(c)- 3.6(f) indicate, however, that the transition between the 1W and 2W states occurs more smoothly for I-MMT than for MMT-MMT. In addition, although as Figure 3.6(c) and Figure 3.6(e) indicate, the height of the Na-OW peak reduces during the transition from fully-saturated 1W to fully saturated 2W, the RDF of Na-OS contains several peaks near the MMT surface, and indicate a small decrease in the height of the peak due to the strong adsorption of the cations onto the clay surface. In addition, the peak of the RDFs of Na-OS in I-MMT does not continuously decrease as water concentration increases. This may be due to the instability of the Na hydration shells during their transitions when CO2 is present in the interlayer. Figure 3.7. Radial distribution functions of K-OW (left) and K-OS (right) at constant CO2 concentrations in (a) and (b) I-MMT with 72% K + and 28% Na + , and (c) and (d) I-MMT with K + only. Figure 3.7 compares the RDFs of K-OS and K-OW in K+I-MMT and (Na+K) + I-MMT. Figure 3.7(a) and Figure 3.7(c) confirm the reluctance of K + to hydrate relative to Na + . As water is added to the interlayer region, the height of the peak in the RDF of K-OW decreases, whereas it increases for K-OS. This is similar to the case with only water in the interlayer space that manifested a jump from the 1W to 2W hydration state. It should also be noted that, compared 38 with Na + , the RDFs have sharper peaks with OS as one of the atoms in the pair and weaker ones with OW, which is due to the low hydration enthalpy of K + . Figure 3.8. Radial distribution functions for Na-OW (left) and Na-OS (right) at constant water concentrations in (a) and (d) MMT-MMT; (b) and (e) I-MMT with Na + , and (c) and (f) I-MMT with 72% K + and 28% Na + . In Figure 3.8 and Figure 3.9 we compare the RDFs of K + and Na + for a constant water concentration as a function of CO2 concentration, varying from 0.5 to 4 molecules per unit cell. As Figure 3.8(a-c) indicate, the RDFs of Na-OW in I-MMT are independent of CO2 concentration, whereas the height of the RDF peak for pure MMT-MMT reduces at higher CO2 concentrations. Such similarities and differences are due to the difference in the spatial distributions of the interlayer species in pure MMT and I-MMT, which were presented in Figure 3.4 and Figure 3.5. In the case of pure MMT, the locations of the peaks of density profiles of water and CO2 are at 39 similar distances from the clay surface. Therefore, as CO2 is added, they disrupt the hydration shell of the counterions, which results in reduced stable interaction of water and Na + . As discussed earlier, however, the density distribution of CO2 has a sharper peak close to the MMT surface. Recall that in the MLC interlayer most of the interlayer cations, as well as water molecules, are adsorbed onto the illite surface. Therefore, only a small portion of the total number of cations are adsorbed onto the MMT surface, where CO2 migrates, and its density distribution has a sharp peak. Thus, the added CO2 mostly increases the density peak near the MMT surface, which does not considerably affect the interaction between water and Na + with the higher density peaks near the illite surface. Figure 3.9. Radial distribution functions of K-OW (left) and K-OS (right) at constant water concentrations in (a) and (b) I-MMT with 72% K + and 28% Na + , and (c) and (d) I-MMT with K + only. Figure 3.8(e)-3.8(f) show that, contrary to the cases with constant CO2 concentration, the heights of the MLC’s RDFs increase as the amount of CO2 increases at constant water concentration. In the case of pure MMT, shown in Figure 8(d), there is a difference between its RDF with those of the MLCs in Figure 3.8(e)-3.8(f). In pure MMT with lower CO2 concentrations, where fully saturated 1W state begins transitioning to the 2W state (with 10 and 20 CO2 per unit cell, respectively), the ISSC adsorption is dominant. In fully saturated 2W state 40 (with 40 and 60 CO2 molecules per unit cell) the strongest Na peak belongs to the OSSCs. Finally, with 80 CO2 molecules per unit cell the transition from the 2W to 3W state is initiated in which the ISSCs are again more dominant than the OSSCs. The three different behaviors are better understood based on their associated density profiles, hence confirming the sensitivity of counterions in pure MMT to the intercalated species. In the case of pure MMT shown in Figure 3.8(d), however, if we consider the RDF for Na-OW shown in Figure 3.8(a), we recognize that as the interactions between Na + and water decrease, the tendency of Na + for surface adsorption increases. In the case of I-MMT with constant water concentration, as CO2 content increases, the peaks in the density profiles of the ISSCs of the counterions become consistently sharper near the illite surface, resulting in a larger height for the first peak of the ISSC larger than that the second peak at the illite clay surface. In addition, the OSSCs in Na + +I-MMT, shown in Figure 3.4(e), transition to the ISSCs at higher CO2 amounts, which also plays a role in giving rise to the sharper peaks of the density distribution of the counterions near the clay surface. Therefore, the height of the RDF for Na-OS increases considerably with the added CO2, as indicated by Figure 3.8(e)-3.8(f). K + reacts to a change in the CO2 concentration in a manner similar to its reaction to increasing water concentration in the clays. Due to the low hydration enthalpy of K + , adding water results in stronger adsorption of the cations onto the surface. Thus, as CO2 is added, the aforementioned disruption of the water network and the low hydration energy of K + conspire together to increase the chances of migration of K + to the clay surface, but this also decreases the chances of allocation of K + near water molecules. The results are shown in Figure 3.9. The RDFs for the pairs CCO2 (carbon in CO2) and HW (hydrogen in water), as well as between OS and OCO2 are presented in Figure 3.10. In the case of CCO2 - HW there is a shoulder near r ≈ 2.3 Å for both monolayer and bilayer states, which represents a hydrogen bond between water and CO2. The lifetime of the hydrogen bond decreases as water concentration increases from the 1W to 2W state, which is also consistent [33] with the average lifetimes of OCO2 and HW. For the bilayer the RDF is similar for all the interlayer regions. In the monolayer, however, the MMT-MMT interlayer and Na+I-MMT manifest more pronounced shoulders that drop sharply when the 2W state forms. The difference between the lifetimes of the hydrogen bonding in the monolayer of the MLCs with only Na + and those containing K + is linked with the 41 distributions of water and CO2 in the latter cases. In fact, a perfect single density peak for water is not observed in the 1W state when K + is present in the interlayer region that could result in a smoother change in the lifetime of the hydration states of such cases. Figure 3.10. Radial distribution functions for OCO2-OS (left) and CCO2-HW (right) in (a) and (b) MMT-MMT; (c) and (d) Na + I-MMT; (e) and (f) I-MMT with 72% K + and 28% Na + , and (g) and (h) K+I-MMT. 42 The RDFs for the OCO2 - OS, which are also shown in Figure 3.10, are very similar in all the cases at the beginning and the end of the transition from the monolayer to the bilayer states. However, similar to the transition of Na + - OS, the transition in the shape of the RDF is smoother for the cases with K + , and decays more rapidly in the cases of MMT-MMT and Na +I-MMT. The reason is that in the latter cases the emergence of the 2W state is delayed, as it is not formed until water concentrations are high. 3.3.4 Diffusion Coefficients We also computed the diffusion coefficients of water, CO2, Na + and K + in the interlayer regions by calculating the time-dependence of their mean-square displacements (MSDs) <R 2 (t)> and using the Einstein’s relation. In the case of water and CO2 the MSDs were computed for their centers of mass. Typical plots of time-dependence of <R 2 (t)> are shown in Figure 3.11. Figure 3.11. Mean-square displacements (MSD) of water (a) and Na + (b) The computed diffusion coefficients are presented in Table 3.2. 3H2O + CO2 refers to the monolayer state, while the other two cases in Table 3.2 represent bilayer hydration states. Generally speaking, the diffusion coefficients of water, CO2 and Na + increase as the transition from the 1W to 2W state takes place. K + has the lowest extent of mobility, which increases only slightly in the bilayer. This is explained by recalling the strong adsorption of K + onto the illite surface. The mobility of Na + is higher in pure MMT interlayer than in the MLCs, which is expected. In fact, the strong surface-counterion attraction in the MLCs reduces the mobility of the cations. In addition, in the two bilayer cases that we studied with a total of 8 and 9 water molecules per unit cells, the diffusivities of water and CO2 are considerably larger with a higher water-to-CO2 concentration ratio. Moreover, comparing the diffusivities in the 1W and 2W hydration states indicates that the change in them due to transition from the 1W to 2W state is 43 more pronounced for CO2 than for water, which is due to the tendency of CO2 molecules to distribute themselves parallel to the clay surface. Table 3.2. The computed diffusion coefficients (×10 −5 ) in the mixed-layers clays (in cm 2 /s). 3H 2O + CO 2 H 2O 0.17 0.11 0.09 0.16 CO 2 0.15 0.04 0.02 0.02 Na + 0.04 0.03 0.007 - K + - - 0.002 0.001 4H 2O + 4CO 2 H 2O 0.23 0.18 0.31 0.29 CO 2 0.65 0.40 0.42 0.35 Na + 0.06 0.04 0.01 - K + - - 0.01 0.002 8H 2O + CO 2 H 2O 1.72 1.19 1.94 2.03 CO 2 0.98 1.17 1.82 1.76 Na + 0.51 0.04 0.02 - K + - - 0.01 0.001 The computed diffusivities in the 3H2O + CO2 in the MMT are in good agreement with those reported previously [31, 105]. The species’ concentrations in the bilayers that we simulated are not the same as those in the previous works and, thus, no direct comparison is possible. We are also not aware of any experimental data to compare with the computed diffusivities. 3.4 Summary This paper reported on the results of extensive molecular dynamics simulation of hydration and swelling of mixed-layer clays (MLCs) consisting of intermixed stacking of illite-montmorillonite (I-MMT), the most common type of mixed clays, in the presence of water and CO2. Various combinations of the MLCs with interlayer cations Na + and K + were simulated, and their swelling was studied as a function of the water and CO2 contents. In order to have a basis for comparison and understand the differences between swelling of pure and mixed clays, we also reported on the same phenomenon in the MMT only. We found that at low CO2 concentrations in Na-MMT, which 44 has its layer charge concentrated in its octahedral sheet, weak ion-surface interactions result in fully hydrated ions and, therefore, more extensive swelling than in I-MMT under the same conditions. In the asymmetric interlayer region of I-MMT, however, the illite that has stronger surface-ion interactions gives rise to considerable cation concentration near its surface, which limits their hydration and, therefore, controls swelling of the MLCs. Further inhibition of swelling of Na-MLC may be achieved by increasing the concentration of K + in the interlayer region. At higher CO2 concentrations, however, intercalation of a mixture of water and CO2 produces a completely different behavior. In the 2W and 3W hydration states the hydrophobicity of the MMT surface is stronger than that of the illite. Therefore, the density distributions of the intercalated molecules are considerably different from those with water only. This behavior, together with the disruption of water network by CO2, reduces the difference in the extent of swelling of the MLC relative to pure MMT. 45 4 Molecular Dynamics Study of the Effect of Layer Charge and Interlayer Cation on Swelling of Mixed-Layer Chlorite-Montmorillonite Clays 4.1 Introduction The goal of this research is understanding the effect of layer substations, water concentration and interlayer cations on swelling of mixed layer clays. In this chapter we describe our molecular dynamics (MD) simulations results, for swelling of mixed-layer chlorite-montmorillonite (CH- MMT) clays in the presence of water. 4.2 Molecular Models and Simulation Protocol The chlorite structure that we simulat in this study was inspired by the work of Meunier et al. [47] who studied the smectite-to-chlorite conversion. In the conversion process saponite and chlorite interact, resulting in the formation of a structure that has the properties of saponite on one side and those of chlorite on the other side in each layer. The supercell structure that we utilized in the MD simulations contained four orthogonal clay layers consisting of two chlorite and two MMT layers with the arrangement, MMT-CH-CH-MMT. The snapshot of the dry supercell structure is presented in Figure 4.1. Each of the layers contained twenty 5 x 4 unit cells, arranged along the (xy) directions (see Figure 4.1), with the supercell containing 80 unit cells in a 5 × 4 × 4 arrangement. As mentioned in the Introduction, each chlorite layer consists of two 2:1-type layer along with a central brucite-like layer, which is why chlorite is sometimes called 2:1:1 clay mineral. Each of the 2:1 layered chlorite has an asymmetric structure, generated as a result of the interactions between saponite and chlorite with the composition M N.PD < [(Si D.DP Al N.fP )(Si :.S Al T.T )]O :N R f :< (OH) A .nH : O, where R f :< =Fe D :< +Mg D :< , and M N.PD < is the interlayer cation. The brackets enclose the atoms in the tetrahedral sheet, while the parenthesis inside the brackets indicate the asymmetric structure of the tetrahedral layer. Due to asymmetry of each 2:1 sublayer of chlorite, and in order to keep the CH-CH interlayer symmetric, the Si :.S Al T.T sides of the two CH sublayers face the brucite-like layer, while the Si D.DP Al N.fP sides face the interlayer species in both CH-CH and CH-MMT interlayers. Three cations, Na, < Cs, < and K, < were used in the MD simulations in order to understand their effect on swelling of the clays. The 46 number of water molecules n was increased monotonically in each set of simulation. The composition of the brucite layer between the two 2:1 layers of chlorite was Mg D (OH) f . The MMT was the same as what we used in our previous studies, [103, 106] with its chemical composition being Na N.OP Si O.OP Al N.:P Al D.P Mg N.P O :N (OH) A .nH : O. The lateral (x,y) dimensions (see Figure 4.1) of each individual clay layer under dry condition was approximately 26 × 36 Å : , and the substitutions in each layer was distributed symmetrically. The cations M < were distributed in the interlayer of the CH-CH, CH-MMT, and MMT-MMT, in order to balance the negative charge of the clay sheets. Figure 4.1. Schematic side view of the simulation cell for the CH-MMT mixed-layer clay. 𝐌 < Represents the interlayer cation. The MD simulations were carried out using the LAMMPS [107] simulation package. The clay structure was described by the CLAYFF force field [30] described in Chapters 2 and 3, and the interlayer water molecules were modeled by the flexible SPC model. The charges assigned to 47 each atom were according to the CLAYFF table. Periodic boundary conditions were imposed in all directions, and the system could freely expand or shrink perpendicular to the clay layers along the z axis shown in Figure 4.1. The dry system was first equilibrated employing the isobaric-isothermal ensemble (NPT) with T = 348 K and P = 130 bar, which are close to the conditions in geological cap rock. After equilibration, the dry basal spacings d were calculated. Next, for all the independent cases with various water concentrations, MD simulations in the NPT ensemble were carried out with a time step of 0.01 fs, which was gradually increased to 1 fs and were continued for at least 20 ns. Longer simulations were carried out for higher water concentrations until equilibrium was reached. Temperature was controlled by a Nose-Hoover thermostat, with the coupling constants being 100 and 1000 timesteps, respectively. The short-range interactions above 13 Å were truncated, while the long-range electrostatic interactions were calculated employing the particle-mesh Ewald summation with an accuracy of 10 r A . For each case the basal spacings were computed at the end of each series of the MD simulation. The equilibrated states obtained by the simulation in the NPT ensemble were then used as the input for the simulations in the canonical (NVT) ensemble, which were carried for up to 20 ns for the cases with higher water concentrations. The results obtained with the NVT ensemble were used to compute the density profiles, as well as the radial distribution functions (RDF). Independent simulations were also carried out for several water contents in each interlayer, which was increased from 1 to 10 molecules per unit cell, in order to monitor the one-layer (1W) and two-layer (2W) hydration of the clays. Therefore, a total of 30 cases were simulated in this study, and the results for the three different interlayers were compared. In all the simulations, the water molecules were randomly distributed in the interlayer, while the cations were initially distributed in the midplane of the interlayers. 4.3 Results and Discussion In what follows we present and discuss some important properties of the clay to better understand the reason for observed behaviors of the swelling clay under different conditions. 48 4.3.1 Swelling Figure 4.2 and Figure 4.3 present the basal spacing d for all the cases that we simulated, which, for the sake of clarity and better understanding of the factors that influence swelling, are grouped based on the interlayer cations. Since the thickness of chlorite and the MMT layers are dramatically different, the basal spacing for all the interlayers, including the CH-CH interlayer, were calculated as the sum of thickness of the MMT clay and the interlayer distance. Computed in this way, the comparison of the amount of swelling will be independent of the clay thicknesses and depends only on the extent of expansion or compression of the interlayer space. Figure 4.2. Comparison of the absolute (a)-(c); and relative (d)-(f) basal spacing of the clay interlayers as a function of interlayer cation and the water content, grouped by the cation type. 8 10 12 14 16 18 0 25 50 75 100 125 150 175 200 d-spacing (Å) Number of Water Molecules MMT-MMT+Na CH-MMT+Na CH-CH+Na (a) 0 2 4 6 8 0 25 50 75 100 125 150 175 200 Relative d-spacing (Å) Number of Water Molecules MMT-MMT+Na CH-MMT+Na CH-CH+Na (d) 0 2 4 6 8 0 25 50 75 100 125 150 175 200 Relative d-spacing (Å) Number of Water Molecules MMT-MMT+Cs CH-MMT+Cs CH-CH+Cs (f) 8 10 12 14 16 18 0 25 50 75 100 125 150 175 200 d-spacing (Å) Number of Water Molecules MMT-MMT+Cs CH-MMT+Cs CH-CH+Cs (c) 0 2 4 6 8 0 25 50 75 100 125 150 175 200 Relative d-spacing (Å) Number of Water Molecules CH-MMT+K CH-CH+K MMT-MMT+K (e) 8 10 12 14 16 18 0 25 50 75 100 125 150 175 200 d-spacing (Å) Number of Water Molecules MMT-MMT+K CH-MMT+K CH-CH+K (b) 49 Figures Figure 4.2(a) – Figure 4.2(c) show the absolute basal spacing of the three interlayers at various water concentrations. To be able to compare the results for the three interlayers, Figure 4.2(d) – Figure 4.2(f) present the basal spacing of the clay interlayers relative to their corresponding dry states. They demonstrate that swelling of the MMT-MMT interlayer and the CH-MMT mixed interlayers are, relatively speaking, similar. For all the three cations, Na, < K, < and Cs, < the CH-CH interlayer has a significantly different basal spacing and lower extent of swelling. In addition, a few points are worth pointing out. (i) For the cases in which sodium acts the cation in both CH-MMT and CH-CH interlayers, the transition to the first and second hydration states is postponed to higher water contents. The delay is very significant for the CH-CH interlayer. (ii) When potassium is the interlayers’ cation, the clay expansion occurs more monotonically, when compared with sodium. As a result, the transitions between the hydration states take place smoother and at similar water contents, with no obvious delays for the CH-CH interlayer. (iii) The cesium clay exhibits a significantly different swelling behavior with a wide first-layer hydration state that persists up to the water content of 100 molecules, which is consistent with the experimental results in the literature. According to Smith [108], the potential energy of cesium- clay clay reaches its global minimum in the monolayer hydrate. This result for the symmetric interlayers, i.e., the CH-CH and MMT-MMT interlayers, is interpreted by the differences in the hydration enthalpy of the interlayer ions and the layer substitutions in the clay. For the asymmetric interlayer - the CH-MMT mixed clay - in addition to the aforementioned factor, the different structures of the clay surfaces that are in contact with the interlayer species play an important role in the observed differences. Na < ,K < ,and Cs < have hydration enthalpies of −406,−320 and−264 kJ/mol, respectively. The hydration enthalpy affects the interaction of the cation and the clay surface and, as it decreases, so also does the amount of clay swelling. This is why K < is used as the water network breaker in order to control swelling of smectite clays [24, 25, 109]. Therefore, the larger and weakly hydrating ions have lower amount of water in their hydration shells and, therefore, lower extent of swelling. Moreover, the amount and the location of the substitutions explain the differences in the amount of swelling for each specific ion. The main difference between the structures of the MMT and CH is the full substitution of Al D< with Mg :< and Fe :< in the octahedral layer of CH. It has 50 been reported [56] that Fe :< in the octahedral sheet has a great suppressing effect on the swelling volume of the MMT, which is similar to Mg :< . Moreover, the tetrahedral substitution of the tetrahedral sheet of the CH facing the interlayer species is considerably higher than that of MMT. This can be better understood by comparing the tetrahedral structures of the two, Si D.DP Al N.fP and Si D.UOP Al N.T:P , which refer, respectively, to the CH and MMT tetrahedral sheets facing the clay interlayer. The comparison clearly indicates 5.2 times more substitution of Si A< with Al D< in the CH tetrahedral sheet, when compared with the MMT. Figure 4.3 groups the relative swelling curves by the interlayer type in order to explore independently the effect of the ion type on swelling. Figure 4.3(a) indicates that in the MMT-MMT interlayer, the overall transitions of the swelling curves are similar, except for the fact that as the hydration enthalpy of the ion decreases and the ion radius increases, the amount of swelling decreases, and the transitions are postponed to higher water contents. Figure 4.3(b) indicates behavior more similar to Figure 4.3(a) than to Figure 4.3(c), hence demonstrating that the MMT side of the CH-MMT interlayer contributes effectively to the clay swelling more than the CH side. Figure 4.3(c), on the other hand, presents an interesting behavior. The Cs-CH-CH system indicates the widest monolayer, when compared with all other cases that we have considered. Although swelling of Na-CH-CH is considerably lower than the other two interlayers, it still exhibits two distinct layers that are, however, unstable, having a short interval of hydration for monolayer with an early transition to the bilayer. Finally, the K-CH-CH depicts something in between the other two cases, which is consistent with the hydration enthalpy of K and its ion size, when compared to the other two ions. The swelling is also monotonic. Figure 4.3. Comparison of the absolute basal spacing of the clay interlayers as a function of interlayer cation and the water content, grouped by the interlayer type. (a) MMT-MMT; (b) CH- MMT + cation, and (c) CH-CH, all plus the cation. 0 2 4 6 8 0 25 50 75 100 125 150 175 200 Relative d-spacing (Å) Number of Water Molecules MMT-MMT+Na MMT-MMT+K MMT-MMT+Cs (a) 0 2 4 6 0 25 50 75 100 125 150 175 200 Relative d-spacing (Å) Number of Water Molecules CH-CH+Na CH-CH+K CH-CH+Cs (c) 0 2 4 6 8 0 25 50 75 100 125 150 175 200 Relative d-spacing (Å) Number of Water Molecules CH-MMT+Na CH-MMT+K CH-MMT+Cs (b) 51 These results show that the CH-CH interlayer is more sensitive to the interlayer cation than are the other two interlayers, hence implying that the octahedral substitutions are more responsive to the ion type than the tetrahedral substations. It should be mentioned that our results for swelling of the MMT-MMT interlayer are consistent with the existing results in the literature [25, 53, 90, 110, 111], and are used for validating the simulations. We also note that the CLAYFF overestimates the surface-ion interactions for the surface oxygen. This results in some extra inner- sphere surface complexes of the ions and the surface and, thus, a slightly smaller extent of swelling than the experimental data. 4.3.2 Molecular structure of the interlayers The results for the basal spacings can be better understood by comparing the density profiles of the various cases that we have studied. Figure 4.4-Figure 4.6 present density profiles of the 1W and 2W hydration states in the interlayer region along the z direction, perpendicular to the clay surface, in the presence of the three cations. The results are for the cases in which the water content was fixed at 60 molecules per interlayer for the monolayer and 180 molecules per interlayer for the bilayer hydration states (3 and 9 molecules per unit cell, respectively). This was done in order to make the comparison at a fixed water concentration for all the interlayers. The computed density profiles for the MMT-MMT interlayers agree well with the previous studies of 1W and 2W hydration states [17, 110, 111]. But, a few variables must be considered in order to correctly interpret the density profiles. One is the he amounts of the tetrahedral and octahedral substitutions, which are different in the CH and MMT interlayers, resulting in different total charges in the layers and, therefore, different amounts of cations in the interlayer space. The asymmetry of the clay surfaces that the CH-MM interlayer species are exposed to is another factor that must also be considered. Foster [56] investigated the effect of all such parameters on swelling of the MMTs. By considering a variety of the MMTs having different octahedral and tetrahedral substations, as well as different amounts of interlayer cations, he reported that the cation concentration, as well as the tetrahedral substitutions, have little to no effect on the swelling when compared with the octahedral substitutions, which manifested a significant correlation with the swelling volume of the MMT. In fact, both Mg and Fe substitution in the octahedral layer greatly impact clay swelling. 52 As discussed in 3.1, the same behavior is seen in the cases that we have studied. For all the cation types the CH-CH with full octahedral substitutions swell the least than the other two interlayers. This can be seen by comparing the density profiles in both 1W and 2W hydration states, as shown in Figure 4.4-Figure 4.6. Such a strong influence is due to the different polarization power of Mg and Fe, compared to Al, which affects polarization of the neighboring atoms and, thus, the entire structure of the clay, changing the interaction energies of all the atoms in the entire CH clay layer. Figure 4.4. Density profiles of (a)-(c) the 1W and (d)-(f) 2W hydration states perpendicular to the clay surface for Na plus (a, d) MMT-MMT; (b, e) CH-MMT, and (c, f) CH-CH. In the CH- MMT interlayer the left and right surfaces are CH and MMT, respectively. Water oxygen and hydrogen are represented by black and red, and Na + by blue. 0 2 4 6 8 10 0 1 2 3 4 5 Density Profile Interlayer Distance (Å) (a) 0 1 2 3 4 0 1 2 3 4 5 Density Profile Interlayer Distance (Å) (b) 0 2 4 6 8 0 1 2 3 4 5 Density Profile Interlayer Distance (Å) (c) 0 2 4 6 8 10 1 2 3 4 5 6 7 8 9 Density Profile Interlayer Distance (Å) (d) 0 1 2 3 4 5 0 2 4 6 8 Density Profile Interlayer distance (Å) (e) 0 2 4 6 8 0 2 4 6 8 Density Profile Interlayer Distance (Å) (f) 53 Comparing the MMT-MMT and CH-CH interlayers, Figure 4.4 indicates that the preference of the water molecules shifts from hydrating the interlayer Na < to forming strong hydrogen bonds with the clay surface. The denser layer substitution of the CH clay are responsible for adsorption of the more polarized water molecules onto the clay surface. The sharper peaks of the water hydrogen near the CH surface, as well as the central water oxygen, can be seen in Figure 4.4(b) and Figure 4.4(c) for the 1W hydration state. For the MMT-MMT interlayer, however, the sharper peak is at the central plane where water molecules form hydration shells around part of the Na < cations, which results in more swelling of the MMT-MMT interlayer. On the MMT side of the CH-MMT interlayer, the asymmetry of the charges attracts more of the water molecules and Na < cations to the CH- surface. Despite the significant differences in the orientations of species in the various interlayers, however, they all depict a monolayer state when 60 water molecules are present in the interlayer with different saturation levels. The same type of evolution occurs in the cases with 180 water molecules with two distinct hydration layers, 1W and 2W. For the MMT-MMT interlayer, most of the cations are transferred to the central plane to form outer sphere surface complexes (OSSC) with sharper hydrogen peaks around the central cations, and some inner sphere surface complexes (ISSC) of the cations attached to the clay surface. The OSSCs are formed when ions are distributed over distances 3-6 Å from the clay surface where the cations are hydrated, whereas in the case of the ISSCs the distances are less than 3 Å from the surface where the cation is only partially hydrated, and is directly adsorbed onto the clay surface [99]. A few points are worth mentioning here. (i) The increased polarization of water molecules can be understood by the extremely sharper hydrogen peaks at the clay surface. (ii) The central peak of Na < in the CH-CH interlayer, which is sharper than that in the MMT-MMT, does not represent more saturated Na < cations in the CH-CH. Instead, it is only a result of the “crowded” distribution of water molecules at the clay surface, which forces the cations to migrate to the central plane. To show this better, we should compare the hydrogen peaks in Figure 4.4(d) and Figure 4.4(f). At the clay surface the hydrogen peaks are at distances 6.3 and 7.4 Å, whereas around the central cations the peaks are at 8.8 and 4.7 Å for the MMT-MMT and CH-CH interlayers, respectively. Therefore, the sharper cation peaks indicate both the priority of the clay surface to adsorb more polarized water molecules, and the partial saturation of most of the cations. This 54 behavior is consistent with the basal spacings discussed in the previous section. Interestingly, the CH-MMT interlayer manifests a behavior very similar to the MMT-MMT near the MMT surface, and even more similar behavior to the CH-CH interlayer near the CH surface, due to the wider basal spacing in the 2W hydration state. This also agrees well with the basal spacings for the CH- MMT that are intermediate between the two cases. Figure 4.5. Density profiles of the 1W (a)-(c) and 2W (d)-(f) hydration states, perpendicular to the clay surface for K plus (a, d) MMT-MMT; (b, e) CH-MMT, and (c, f) CH-CH. For the CH- MMT interlayer the left and right surfaces are CH and MMT, respectively. Water oxygen and hydrogen are represented by black and red, and K + by blue. Figure 4.5 indicates similar behavior for the three interlayers when K < is inserted. In general, K < perturbs the water network and, relative to Na < , is reluctant to fully hydrate [97, 98]. 0 1 2 3 0 1 2 3 4 5 Density Profile Interlayer Distance (Å) (a) 0 1 2 3 4 0 1 2 3 4 5 Density Profile Interlayer Distance (Å) (b) 0 2 4 6 8 0 1 2 3 4 5 Density Profile Interlayer Distance (Å) (c) 0 1 2 3 4 5 0 2 4 6 8 Density Profile Interlayer Distance (Å) (d) 0 1 2 3 4 5 0 2 4 6 8 Density Profile Interlayer distance (Å) (e) 0 2 4 6 0 2 4 6 8 Density Profile Interlayer Distance (Å) (f) 55 Thus, due to the lower hydration enthalpy of K < and its reluctance to hydrate, its peaks in the central plane, shown in Figure 4.5(d), are not as sharp as those in Figure 4.4(d), which is the reason there is a more uniform distribution of K < in the interlayers, as indicated by Figure 4.5. This results directly in less swelling in the K-clay interlayers. Figure 4.6. Density profile of the 1W (a)-(c) and 2W (d)-(f) hydration states perpendicular to the clay surface for Cs plus (a, d) MMT-MMT; (b, e) CH-MMT, and (c, f) CH-CH. In the CH-MMT interlayer the left and right surfaces are CH and MMT, respectively. Water oxygen and hydrogen are represented by black and red, and Cs + by blue. Figure 4.6 manifests significant differences with the other two cases described so far. Compared with Na < and K < ,Cs < has the smallest hydration enthalpy and the largest ionic radius. 0 1 2 3 0 1 2 3 4 5 Density Profile Interlayer Distance (Å) (a) 0 2 4 6 0 1 2 3 4 5 Density Profile Interlayer Distance (Å) (b) 0 1 2 3 0 1 2 3 4 5 Density Profile Interlayer Distance (Å) (c) 0 1 2 3 4 5 0 2 4 6 8 Density Profile Interlayer Distance (Å) (d) 0 1 2 3 4 5 0 2 4 6 8 Density Profile Interlayer Distance (Å) (e) 0 1 2 3 4 5 0 2 4 6 8 Density Profile Interlayer Distance (Å) (f) 56 Thus, as discussed earlier, the global minimum energy of the Cs-clay is reached by the monolayer formation. As mentioned for Figure 4.3(c), the wide range of monolayers for swelling of Cs-CH- CH supports the results in Figure 4.6(c), where the 1W hydration state is not yet stable with 60 water molecules. In addition, due to the larger size of Cs < , the adsorption of the cations onto the MMT surface in the CH-MMT interlayer is opposite of the other two cases described so far, where at the lower water concentrations the CH surface adsorbs both cations and water molecules. Moreover, in the bilayer state of Cs-MMT-MMT, most of the Cs < cations form the ISSC and occupy hexagonal cavities adjacent to the tetrahedral substitution [108], instead of forming the OSSC, which is due to their low hydration enthalpy; see Figure 4.6(d). Overall, Figure 4.4-Figure 4.6 indicate that much stronger hydrogen peaks develop for the monolayer and bilayer of the Na-MMT-MMT, whose height is considerably reduced in the cases of K < and Cs < . Moreover, as one transitions from the interlayers with Na < to those with K < and Cs < , the cations’ hydration shells become less saturated, resulting in less swelling, due to the change in the hydration enthalpy and the size of the cations. Independent of the cation type, in both the 1W and 2W hydration states of all the CH-CH interlayers, sharper water oxygen and hydrogen peaks develop at the clay surface, which are sharper around the central cations in the MMT-MMT interlayers, where the cations are hydrated. In addition, comparison of the cations’ distributions in Figure 4.4-Figure 4.6 indicates that there are two peaks near the CH surface in the 2W state of CH-MMT and CH-CH, corresponding to K < and Cs < , which represent two types of the ISSCs formed at the CH surfaces. There is, however, only one ISSC peak corresponding to the systems with Na < . This agrees well with the study of Underwood et al. [108] that indicated that weakly hydrating cations form the ISSCs, both above the hexagonal cavities of the clay surface and the basal oxygen atoms. However, strongly hydrated cations, such as sodium, only form ISSCs above the basal oxygen atoms of the clay and, instead, prefer to form OSSCs at high water contents. Finally, comparing the density profiles of the mixed-layer CH-MMT with our previous study on the mixed-layer I-MMT [103], we see a significant difference between the distributions of the interlayer species. For illite with more tetrahedral distributions, most of the cations are adsorbed onto the surface as the ISSCs, independent of the water content in the system. However, as discussed, the dense octahedral substitutions of the CH side of mixed interlayer polarizes the water 57 molecules remarkably, with the ions distribution being a function of the water content of the system. This also explains the differences between the basal spacings of the different cation types. 4.3.3 Radial distribution functions Figure 4.7 presents the calculated radial distribution functions (RDFs) g(r) of Na < and oxygen in the water molecules, as well as the clay surface oxygen in the three interlayers, namely, Na-MMT- MMT, Na-CH-MMT and Na-CH-CH. In all the cases, independent of the intensity of the RDFs, the peak location of r t u r and r t v r is 2.4±0.1 Å, which agrees well with the previous studies [112, 113]. The location of the peak reflects the sum of the ionic radii of oxygen and sodium. In all the cases, as water content in the interlayer increases, the height of the peaks corresponding to both r t u r and r t v r decreases. The value of r t v r confirms the results for the density profiles of sodium and the clay surfaces described earlier. The height of the RDF peak for Na-MMMT-MMT decreases considerably in the bilayer state, where sodium cations migrate to the central plane to hydrate. The same significant decrease in the peak height also occurs for the Na-CH-CH. As the discussions in the last section indicated, however, this is mostly due to the replacement of the ISSC with Na < with water molecules that form strong hydrogen bonds with the clay surface. Therefore, the RDF height of Na < decreases at the clay surface but increases farther in the central plane where they are only partially hydrated, if at all. In addition, in contrast with the other two cases, in the Na- CH-MMT interlayer Na < is exposed to two different clay surfaces and, thus, to the strong charge asymmetry of the interlayer, which is why the RDF peak height does not decrease substantially, hence indicating a persistent ISSC of Na < with the surface of the CH side of the interlayer; see Figure 4.4. 58 Figure 4.7. Radial distribution function for 𝐍𝐚−𝐎 𝐖 (a)-(c) and 𝐍𝐚−𝐎 𝐒 (d)-(f) at 1W (black) and 2W (red); (a, d) Na + MMT-MMT; (b, e) Na + CH-MMT, and (c, f) Na + CH-CH. When considering the density profiles and the RDFs, the cations’ distributions should be carefully analyzed, because the dense distribution of the cations at the midplane of the CH-CH interlayer may be mistakenly viewed as hydrated cations, which, as explained earlier, is not the case. Figure 4.8 presents the RDFs for K < and the three interlayers. The peak location of r t u r and r t v r is at 2.8±0.1 Å, which agrees with the literature [112, 113]. As expected, due to the lower hydration enthalpy of K < , compared with r t u r , r t u r indicates decreased heights for both the mono-and bilayers. Value of r t v r indicates that, compared to Na-MMT-MMT, there is not significant migration of K < to the central plane in K-MMT-MMT, as one transitions from the monolayer to bilayer, hence indicating once again the reluctance of K < to hydrate even at higher water contents. Although more of the K < migrate to the central plane of the K-CH-CH interlayer, this is again due to the preference of the clay surface for adsorbing more water molecules than K < in the CH-CH interlayer, which results in less swelling than in the MMT-MMT interlayer. 0 5 10 15 20 25 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (a) 0 0.3 0.6 0.9 1.2 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (d) 0 18 36 54 72 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (c) 0 1 2 3 4 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (f) 0 16 32 48 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (b) 0 1 2 3 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (e) 59 Figure 4.8. Radial distribution function for 𝐊−𝐎 𝐖 (a)-(c) and 𝐊−𝐎 𝐒 (d)-(f) at 1W (black) and 2W (red); (a, d) K + MMT-MMT; (b, e) K + CH-MMT; (c, f) K + CH-CH. Figure 4.9 presents the RDFs for the structures containing Cs < . The peak of r t u r x and r t v r x is at 3.2±0.1 Å, which agrees with the data in the literature [114]. The distributions are similar to those for K < in the three interlayers. Due to the ionic radius and hydration enthalpy of Cs, < however, r t u r x has the lowest heights among the cations that we study in both the 1W and 2W hydration states, due to which one has a sharper RDF in the bilayer state for both Cs-MMT and Cs-CH-CH at r t u r x , when compared with the other two cations. This indicates that, even in the bilayer state, the Cs < are not hydrated and are distributed more near the clay surface than the central planes. This is also consistent with the density profiles in which there is a central peak for both the Na < and K < bilayers in the CH-CH and MMT-MMT interlayer, which confirms partial and fully-hydrated cations in the interlayer regions, respectively. On the other hand, there is no central peak in the bilayer density distribution of Cs < and the two interlayers, which again confirms 0 4 8 12 16 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (a) 0 0.5 1 1.5 2 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (d) 0 6 12 18 24 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (b) 0 1 2 3 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (e) 0 9 18 27 36 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (c) 0 1 2 3 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (f) 60 that the global minimum energy of Cs < is reached in the monolayer over a wide range of water concentrations. We point out that, due to the differences in the total negative charge of the CH and MMT clays, different number of ions exist in each interlayer. In the MMT-MMT, CH-MMT and CH-CH interlayers a total of 28, 35 and 42 cations are available for neutralizing the charge of the entire structures, which also provides the reasons for the order of r t u r xW]X in the three interlayers peaks, namely, CH-CH>CH-MMT>MMT-MMT. Figure 4.9. Radial distribution function for 𝐂𝐬−𝐎 𝐖 (a)-(c) and 𝐂𝐬−𝐎 𝐒 (d)-(f) at 1W (black) and 2W (red); (a, d) Cs + MMT-MMT; (b, e) Cs + CH-MMT; (c, f) Cs + CH-CH. 0 3 6 9 12 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (a) 0 1 2 3 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (d) 0 5 10 15 20 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (b) 0 1 2 3 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (e) 0 9 18 27 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (c) 0 1 2 3 0 1 2 3 4 5 6 7 8 g(r) distance (Å) 1W 2W (f) 61 4.4 Summary Independent of the cation type, sharper peaks for the water oxygen and hydrogen are manifested at the clay surface in both the 1W and 2W hydration states of all the CH-CH interlayers. This indicates the reduced number of water molecules bound to the cations, as well as their increased polarization, having hydrogen bonds with the clay surface. Such a behavior demonstrates how the octahedral substitutions greatly affect clay swelling, whereas the cation concentration, interlayer asymmetry and tetrahedral substitutions have considerably less influence on the expansion or compression of the clay interlayers. The hydration enthalpy and size of the cations play significant roles in the swelling of various interlayers. Compared with K < and Na < , Cs < hydrates the least, resulting in significant differences in the density profiles of interlayer species and the basal spacing, as well as the RDFs. It should also be pointed out that, due to the complexity of the mixed-layer structures along with the role of the interlayer cations, the basal spacings, density profiles and the RDFs should be analyzed simultaneously in order to avoid incorrect interpretation of the overall system’s behavior. Finally, although both illite and CH control swelling of mixed-layer clays, there are significant differences between the properties and behavior of I-MMT and CH-MMT. This means that for mixed-layer clays, combinations of a variety of parameters with unique properties play a role in their swelling and more studies and investigations should focus on MD simulations with sets of such parameters in order to gain a detailed understanding of the parameters’ contributing to swelling or shrinking of the other common mixed-layer clays. 62 5 Two-phase flow of CO 2 -brine in a heterogeneous sandstone: Characterization of the rock and comparison of the lattice-Boltzmann, pore-network, and direct numerical simulation methods 5.1 Introduction In this chapter we begin to describe our efforts for modeling of flow of carbon dioxide and brine in a heterogeneous Ssandstone formation and comparing the results of the lattice-Boltzmann, pore- network, and direct numerical simulation methods. To do so, we have been carrying out the simulation using three-dimensional (3D) images of the formation. Thus, we begin by describing the details of simulation methods. This work was carried out in collaboration with the group of Professor Albert Valoochi's of the University of Illinois at Urbana Champaign 5.2 Characterization of Mt. Simon Sandstone The sample porous medium is from Mt. Simon sandstone at a depth of 6700 feet. The formation is located at verification well number 2 of a study site in Decatur, Illinois where, as mentioned earlier, Illinois State Geological Survey carried out a pilot injection study to better understand the feasibility of full-scale CCS [64]. A core plug from the formation was scanned by micro-CT imaging technique at the National Energy Technology Laboratory (NETL) of the U.S. Department of Energy, which produced a series of grey-scale scans. Table 5.1 presents the information of a cubic sample from the core plug. A series of image processing steps were taken in Fiji [115] to filter and smooth images in order to distinguish void space from the solid grains via a thresholding algorithm. The outcome was a segmented 1200 D voxels sample with porosity of 27.1% and voxel size of 2.80 µm, i.e., a 3.36 D mm D sample. 63 Table 5.1. Rock source and properties. Rock type Mt. Simon sandstone Depth (feet) 6700 Scans resolution (µm) 2.80 Images size (voxels) 1200 D Sample size (mm D ) 3.36 D Porosity (%) 27.1 We took eight equal-size subsamples (each of size 500 D voxels) from the center of the main sample to study its geometrical, topological, and flow properties. Table 5.2 presents a list of the subsamples with their corresponding porosity. Figure 5.1 presents both a high-resolution cross section of the original sample and a 3D representation of a 500 D voxels subsample S2. One slice of this subsample in both gray-scale and segmented format is presented in Figure 5.2 and Figure 5.3. Table 5.3 presents the important properties of subsample S2. Figure 5.1. Thin-section scan of Mt. Simon sandstone at a depth of 6700 feet, and a 𝟓𝟎𝟎 𝟑 subsample S2 of the formation selected for the two-phase flow simulation. 64 Figure 5.2. Gray-scale (left) and segmented (right) formats of the 300 th slice (out of 500 slices) of the selected Mt. Simon sandstone subsample S2. Table 5.2. Porosity of eight subsamples of Mt. Simon sandstone Subsamples Name Porosity S1 0.258 S2 0.263 S3 0.272 S4 0.274 S5 0.287 S6 0.274 S7 0.283 S8 0.293 5.2.1 Construction of the pore-network To characterize the heterogeneity of Mt. Simon sandstone, as well as carrying out single- and two- phase flow simulations, we extracted from the images the equivalent PNMs of the eight sandstone subsamples. To do so, we used the maximal ball (MB) algorithm [116] that searches the entire voxelized pore space in order to identify the largest possible spheres in the porous medium. The algorithm was developed originally by Silin and Patzek [117] and was extended and improved by Al-Kharusi and Blunt [118], Dong and Blunt [116], Arand and Hesser [119], and Raeini et al. 65 [120]. In the algorithm the input is the voxelized binarized geometry of the porous medium in which the solid and pore phases are stored as 1 and 0, respectively. All the zero voxels are scanned and the largest possible voxelized sphere in the pore space for each of them is determined. The resulting voxelized sphere for each pore voxel is taken as the MB. In practice, many MBs are inside the larger MBs and, therefore, should be removed. The resulting MBs are sorted and clustered based on their size that helps identifying the ancestor MBs – the local maximums based on the size - and a chain of the MBs that are from one ancestor to another one (Dong and Blunt. 2009). Each chain is segmented such that it represents a configuration of two pore-bodies and their connecting pore-throat. Finally, by counting voxels of each element, pore-body or pore-throat, various geometrical features, such as the location, radius, volume, length, and shape factor are calculated and stored. The pore-throat length is defined as the difference between the total pore-throat length and the lengths of the neighboring pore-bodies. The pore-body length is defined as a function of its radius and the total pore-throat length is defined as center-to-center Euclidean distance between two neighboring pore-bodies (Dong and Blunt, 2009). Table 5.3. Properties of the selected Mt. Simon sandstone subsample S2. The permeabilities are in mD. Sample size (voxels) 500× 500×500 Sample size (mm D ) 1.4 × 1.4 × 1.4 Porosity (%) 26.3 Permeability by the LB simulation 4278 Permeability by the DNS 4370 Permeability by the PNM 4201 The shape factor G summarizes all the irregularities of the geometry of the pore-bodies and pore-throats by, G=VL/A : , where A is surface area, and V is volume of the voxelized element. L is defined as twice the distance between the center of the ancestor MB and the farthest voxel in that MB [116]. The shape factor is a key quantity that helps assigning familiar geometries, such as circles, squares, and triangles, to the cross sections of the pore-bodies and pore-throats in the PNM flow solvers (Patzek and Kristensen, 2001, Patzek and Silin, 2001). 66 Figure 5.3 presents the resulting PNM for sample S2 of Mt. Simon sandstone that we generated using the algorithm, as well as the size distributions of its pore-bodies and pore-throats. The ratio of the median pore-throat length and the radius is about 16.7. Table 5.4 presents the data for the connectivity of all the eight samples in terms of the average coordination number of their resulting PNMs. The connectivity density of the samples, i.e., connectivity per unit volume, is on the order of 2×10 r P [pixel r D ], while the ratio of their pore-throat length and radii is similar to that of sample S2. Figure 5.3. The PNM of Mt. Simon sandstone and its pore-body and pore-throat size distributions. Table 5.4. The connectivity and anisotropy of the eight subsamples Subsamples Name Average Coordination Number Degree of Anisotropy S1 4.77 0.23 S2 4.85 0.25 S3 4.90 0.19 S4 4.73 0.18 S5 5.24 0.19 S6 5.02 0.18 S7 4.78 0.25 S8 5.11 0.22 Also shown in Table 5.4 is the degree of anisotropy (DA) of all the eight samples. The DA is a measure of how highly oriented the pore structure of a medium is within a volume [121-124]. We used the BoneJ plugin in Fiji, the image processing package, to calculate the DA of the void 67 space of the eight samples. The method is based on the mean intercept length (MIL) method by which a large number of equal-length vectors originating from a random point within the pore space are drawn and an intercept is counted for each vector when it hits a boundary. The MIL is the vector length divided by the number of boundary hits. A cloud of points is formed where each point represents the vector times its MIL, and then an ellipsoid is fitted to the cloud. The anisotropy tensor is then constructed, and the eigenvalues and eigenvectors related to the lengths and orientations of the ellipsoid's axes are computed [123]. The DA is defined as: DA = 1 − W ] dW ] (1) The algorithm is stochastic and, therefore, new random points with the same vectors yield new MIL counts. The DA must be updated until either the minimum number of sampling points is reached or the coefficient of variation of DA becomes smaller than a threshold. As Table 5.4 indicates, the samples are slightly anisotropic. Thus, we ignore it in the flow calculations. 5.2.2 Analysis of the heterogeneity We carried extensive analysis of the morphology of the selected sample of Mt. Simon sandstone in order to characterize the severity of its heterogeneity. This was accomplished by analyzing the 3D image of sample S2, extracting its equivalent PN, and comparing the results to those of a sample Berea sandstone, which is a fairly homogeneous porous medium and whose properties are presented in Table 5.5. In the study of the heterogeneity via image analysis, we took fifteen equal- size subsamples, each 1/8 of the original sample, eight of which were from the corners, six from the sides, and one at the center. We then computed the porosity and degree of anisotropy DA of the subsamples. The results are presented in Figure 5.4. Table 5.5. Berea Sandstone characteristics. Rock type Berea sandstone Scans resolution (µm) 5.34 Images size (voxels) 400 D Sample Size (mm D ) 2.14 D Porosity (%) 19.6 68 To quantify the variability across the subsamples, we computed the coefficient of variation (CV) of each property for each sample. The results are presented in Table 5.6, according to which the CVs of the porosity and the DA from the Mt. Simon sandstone are larger than those for the Berea sandstone sample, hence indicating that the former has more variability is more heterogeneous. Table 5.6. Comparison of the coefficient of variation (CV) of the properties of the Mt. Simon sandstone sample S2 and the Berea sandstone sample Sample CV, Porosity CV, DA Berea 0.039 0.214 Mt. Simon 0.049 0.250 Figure 5.4. Porosity and degree of anisotropy of 15 subsamples from the Mt. Simon sandstone sample S2 and their comparison with those of the Berea sandstone sample. 69 In analyzing the heterogeneity of Mt. Simon sandstone using the PN, we used a stochastic approach whereby a calculation box with a size 1/8 of the PN extracted from the original sample was selected randomly. Then, the average coordination number Z of the pore-bodies inside the calculation box was calculated, since Z quantifies the connectivity of the pore space well. The results for 100 realizations of the Mt. Simon sandstone sample S2, and their comparison with those for the Berea sandstone sample are presented in Figure 5.5. According to the data, the CV of the average coordination numbers is 0.054 for the Mt. Simon PN, but 0.025 for the Berea PN. Therefore, the Mt. Simon sandstone sample S2 encompasses more variability in its pore connectivity relative to the Berea sandstone sample and, thus, it is more heterogeneous. Figure 5.5. Statistics of the average coordination numbers of over 100 realizations of the pore- network of the Mt. Simon sandstone sample S2, and the corresponding Berea sandstone sample. 5.2.3 Determining the size of the representative elementary volume A careful study was undertaken in order to determine the size of the representative elementary volume (REV) for the flow studies. It turned out that a size of 500 D voxels was large enough to represent the average flow properties. Figure 5.6 presents three computed properties of the sample versus its linear size, calculated by Fiji, while Figure 5.7 shows the variation of the absolute 70 permeability computed by the LB simulation in order to identify the size of the REV. The variation of the computed properties with the sample size becomes negligible at a size of about 300 D voxels, thus justifying the choice of simulation sizes larger than the identified size. Using the three computational approaches described in the next section, we also carried out single-phase flow simulation with the eight subsamples employing the three aforementioned approaches, namely, the LB method, direct numerical simulation of Stokes’ flow using OpenFOAM, and PN calculations, in order to evaluate the absolute permeability of each of the eight samples. The results will be discussed shortly. For now, it suffices to say that they indicate that the original sample is highly heterogeneous since different subsamples from different locations have different average properties, while they are also representative of a REV. Thus, as mentioned earlier, to study two-phase flows, we selected subsample S2 in Table 5.2, whose characteristics are presented in Table 5.3. 5.3 Computational Approaches As mentioned earlier, we have carried out extensive simulations of fluid flow using three distinct computational methods, the details of each are as follows. 5.3.1 The lattice-Boltzmann method As is well known, the LB method is based on streaming, collision and relaxation of a set of fluid particle distribution functions (PDF) on a lattice. The no-slip boundary conditions on solid surfaces are implemented by simply switching the directions of the particles on the surface nodes, the so- called bounce-back scheme. There are several LB schemes for simulating multiphase flows [125, 126]. Among them, the color-fluid model [127, 128] is capable of producing a relatively sharp interface between completely immiscible fluids, which is why it has been widely adopted [82, 129]. In addition, it can deal with high viscosity ratios due to its independent control of the surface tension and viscosity [83]. In the present work, we use a variant of the multiple relaxation time (MRT) color-fluid LB simulator [83, 84, 130]. 71 Figure 5.6. Sample-size dependence of three properties of the subsample S2 of Mt. Simon sandstone in order to identify the size of the REV. In the model, each phase has its own set of PDFs and the discrete Boltzmann’s equation is solved for each fluid phase. In our in-house LB code, we consider two sets of the D3Q19 PDFs, i.e. a 3D model with 19 velocities, representing the two fluid phases, referred to as the fluids r and b, which follow the collision-streaming procedure for the PDF: f ] (x+e ] Δt, t+Δt) =f ] (x,t)+ Ω ] (D) ¡Ω ] (T) +Ω ] (:) ¢,s=r,b (2) 72 Figure 5.7. Sample-size dependence of absolute permeability, computed by the LB simulation, of subsample S2 of Mt. Simon sandstone. where Ω ] (T) is the standard LB collision operator, Ω ] (:) is the perturbation step that generates the surface tension effect, and Ω ] (D) is the recoloring step that separates the two fluids. The collision operators Ω ] (T) and Ω ] (:) are constructed under the MRT framework that increases stability and accuracy of the model [83, 84, 131, 132] . More details of our in-house code are given by [83]. The macroscopic quantities, such as fluid velocity and pressure, are computed by calculating the moments of the PDF. Since the outputs produced by the LB simulation are defined in terms of lattice units, they must be converted to physical units [80, 133]. In the present work, we simulate a method of measuring the relative permeabilities known as the steady-state method. In this method a predefined fractional flow of both fluid phases is injected into the pore space at constant flow rates, while the pressure drop across the sample is constant. Steady state is reached when the downstream and upstream fractional flows are equal. More details about the steady-state measurement are given by [134]. In the LB model, the predefined fractional flow is implemented as an initial randomly distributed saturation that will also be the target saturation. We mirror the input geometry and impose periodic boundary conditions along the flow direction as well, in order to allow both fluid phases enter and exit the model smoothly. A body force is applied to each fluid phase to achieve the same pressure drop and avoid capillary-end effects. To ensure that steady state has been reached, various quantities, such as the flow rates, should be monitored in order to check whether they converge to steady- state values, which are then used for the relative permeability calculation at the set target 73 saturation. A similar method of measuring relative permeabilities was described by Ramstad et al. [80]. 5.3.2 Direct numerical simulation The second computational method that we employed was direct numerical simulation (DNS) using the OpenFOAM open source for both single- and two-phase fluid flow. The simulator uses finite- volume discretization and solves the continuity and momentum equations – the Stokes equations - in the pore space. We used an unstructured mesh with grid size of 2 3 and 1 voxels in the pore volumes and throat, and corner points between the solid boundaries, respectively, meaning that the grid blocks that are used in the throats and the corners are half the size in each direction of those used in the pores, employing a modified OpenFOAM mesh generator [135]. The total number of grid blocks used to discretize the image volume was 10,815,820. For each capillary number, the simulations employed 300 processors in parallel. Unlike the LB method, the finite-volume simulations are run under unsteady-state conditions, and the results are measured in a dynamic drainage simulation. More details about the method and its implementations are given in the two phase flow solver of OpenFOAM (interFOAM) [136]. As usual, Darcy’s law and its generalization to two-phase flow, are assumed. There are various methods for upscaling the microscopic dynamic and capillary pressures to the Darcy-scale pressure, but it is not clear which method averages the microscopic pressure most accurately. Raeini et al [57] analyzed five methods in a steady-state incompressible single- phase flow calculation to obtain the Darcy-scale pressure drop. Among all the methods the velocity-weighted average of the viscous forces, as well as the velocity- weighted average of the pore-scale pressure gradient, matched the experimental results most closely. The same equations with some modifications can be used to calculate the macroscopic pressure drop in two-phase flow. In the DNS simulations, the velocity-weighted average of the viscous forces was employed in order to calculate the total macroscopic pressure drop (∆P ¥ ) in each phase: ∆P ¥ =− 1 Q ¥ §u.[∇.(µ∇u)]dV ¥ , µ=βµ T +(1−β)µ : (3) (4) 74 where α may be either the wetting or non-wetting phase, Q ¥ is the flow rate of phase α, u is the velocity vector, and V ¥ is the portion of flow volume filled with phase α. In Eq. (4) µ T and µ : are viscosities of the two fluids, and β is the volume fraction of fluid 1 in each grid cell. More details on upscaling of pore-scale forces are given by Raeini et al. [57]. Various methods have been developed for including the interfacial tension in Eulerian grids. They include the continuous surface stress (CSS) [85], continuous surface force (CSF), sharp surface force (SSF) [137], and filtered surface force (FSF) [138]. The main problem with the CSF is the spurious velocity in the flow field, which, to some extent, is controlled in the SSF method. Controlling the sharpness of the capillary pressure in the SSF method is accomplished by defining a sharpened indicator function. In the FSF method the indicator is modified to have a smoother capillary force compared to the SSF approach for the interface motion. This modification compresses the transition area of the capillary pressure. With this implementation, the issue with the non-physical velocities that may arise is resolved, and the capillary pressure transition area is only one grid block. Therefore, the FSF formulation is used for the complex geometry used in this study. 5.3.3 Pore-network model The PNM simulation of drainage – displacement of brine by CO2 - was carried out under the quasi- static condition, which corresponds to low capillary number. It uses an invasion percolation (IP) algorithm [139]; for the most efficient algorithm to simulate the IP see [140, 141]. All the pore elements were initially saturated with the wetting phase. Pore elements become occupied by injecting non-wetting phase through piston-like displacement and based on the Young-Laplace equation that connects the pressure in the two phases at the interface between them. The complete procedure that we used is described by Valvatne and Blunt [142], and need not be repeated here. All the pore-bodies and pore-throats were assumed to triangular cross sections. In all three approaches the usual generalized Darcy’s law for multiphase flows, v ¥ =− o ,® (i ® )o ¯ ® ° ± ® ° ² (5) was used to compute the relative permeabilities, where k d,¥ is the relative permeability of phase α, S ¥ is its saturation, and v ¥ is its corresponding Darcy velocity, which is proportional to total flow 75 rate passing through the medium. The relative permeability is generally a function of the phase saturation, wettability, and the structure of the pore space. The competition between the capillary and viscous forces influences the displacement of one fluid by the second one, which is expressed by the capillary number Ca: Ca= ¯(³ ´µ´ /¶) a (6) where Q WXW refer to the total flowrate of phases, µ is effective viscosity, A is the cross-sectional area, and σ is the surface tension. 5.4 Test of the Accuracy of the Numerical Approaches We tested our computational methods with an oil-water system in a water-wet Berea sandstone sample. The sample, having the data listed in Table 5.5, has been extensively studied and used in the literature [57, 80, 142] and is considered a benchmark. We compare the calculated relative permeabilities of the oil-water system during drainage with the experimental data [143] in Figure 5.8. Table 5.7 presents the properties of the two-phase flow system. The LB results for the water relative permeabilities are slightly larger than the experimental data but, as we discuss below, this is due to the resolution of the lattice used in the simulation. As the resolution increases, the agreement between the LB results and the data improves significantly. The agreement between the DNS results and the data is excellent. Table 5.7. Properties of CO2-brine pair for the LB simulation and DNS Contact angle (°) 125 (water-wet) Interfacial tension (mN/m) 30.0 Water density (Kg/m D ) 1000 Water viscosity (Pa.s) 1.05×10 r D Oil density (Kg/m D ) 1000 Oil viscosity (Pa.s) 1.39×10 r D Viscosity ratio 1.32 Capillary number 5×10 r P 76 The results from the DNS and LB methods are reported up to the water saturation of about 45% (S · =0.45) where water permeability becomes negligible. We did not further continue the simulations due to the limited computational resources. Figure 5.8. Drainage relative permeabilities of the oil-water pair in a Berea Sandstone sample with contact angle of 𝟏𝟐𝟓 degrees, capillary number of 𝟓×𝟏𝟎 r 𝟓 , and viscosity ratio of 𝟏.𝟑𝟐. 5.5 Results and Discussion After validating our computational methods with the experimental data, we carried out extensive simulation of both single- and two-phase flow in the image of Mt. Simon sandstone. In what follows we present and discuss results. 5.5.1 Single-phase flow As pointed out earlier, the size of the REV for the single-phase simulations with the eight subsamples was 500 3 voxels. Using the DNS, we computed the single-phase permeability for an image size of 500 3 and resized subsamples of size 250 3 voxels (with the same physical size, i.e., 1.4 3 mm 3 , but with a coarser image resolution of 5.6 μm) with the upstream and downstream pressures of 5 Pa and 0 Pa. The fluid density and viscosity were set to be the same as those of brine, μ=0.0011 Pa.s and ρ=1100 kg/m 3 . Figure 5.9 shows the original mesh, and the pressure distribution of the subsample S5 with 500 3 voxels. 77 Figure 5.9. The 3D pore volume of the subsample S5 (left) and the pressure distribution in the volume at the end of the simulation (right). The same properties were used in the LB method to calculate the absolute permeabilities for both sizes. However, due to the low computational cost of simulation of single-phase flow with the PNM, no resizing was needed, and the simulations were performed only with the original 500 3 voxels samples. The main purpose of carrying out the LB and DNS simulations on the resized subsamples was reducing the size of the input geometry of our two-phase flow simulations that reduces significantly the computational costs for both methods. The absolute permeabilities of the eight subsamples, calculated by the DNS, LB and PNM, are reported in Table 5.8. A few features of Table 5.8 are worth pointing out. (i) The absolute permeabilities computed by the DNS and LB for the original 500 3 subsamples are quite close to that of the resized 250 3 subsamples. This indicates that the resizing process honors the connectivity of the sample well. The small increase in the permeabilities of the resized 250 3 subsamples is due to the interpolation of the voxelized geometry during downsizing that removes minor irregularities in narrow pores and, thus, reduces flow resistance. However, this does not change the overall connectivity of the sample, and the resulting permeabilities are still in good agreement. As an example, the absolute permeability of subsample S2 changes from 4278 mD to 4211 mD (computed by the LB simulation) after resizing, while its porosity remains almost unchanged. Therefore, the resized 250 3 subsamples were considered as a good-enough approximation of the original 500 3 subsamples and were used in two-phase flow simulations. (ii) The permeabilities computed by the LB simulations agree with those computed by the DNS for the 500 3 images; the difference between the two sets of results is about 10 percent or less in most cases. This is also 78 encouraging. (iii) A comparison of all the computed absolute permeabilities of all the eight subsamples indicates that the original Mt. Simon sandstone is indeed heterogeneous. Table 5.8. Comparison of the absolute permeability of the eight subsamples computed by the PNM, LB and DNS methods. The results obtained by the LB and DNS methods are for two sample sizes of 250 3 and 500 3 that have resolutions of 5.6 and 2.8 μm, respectively. Subsamples Name Porosity Permeability (mD) PNM (500 3 ) LB (250 3 ) LB (500 3 ) DNS (250 3 ) DNS (500 3 ) S1 0.258 3843 3883 3638 4147 3971 S2 0.263 4201 4211 4278 4530 4370 S3 0.272 5037 5053 4818 4555 4376 S4 0.274 5283 5746 5772 5521 5344 S5 0.287 6055 5666 5397 6721 6418 S6 0.274 5884 4830 4614 5085 4860 S7 0.283 6784 5182 5106 8186 7879 S8 0.293 8112 6645 6558 7907 7595 5.5.2 Two-phase flow We used both the LB method and DNS, in addition to the quasi-static PN computations, to carry out simulation of drainage for the CO2-brine pair in the Mt. Simon sandstone sample S2 and evaluate the drainage relative permeabilities for two capillary numbers Ca. Table 5.9 presents the properties of the CO2-brine system used in the simulations. We first carried out simulation of drainage with specified injection velocity (or flow rate) at the inlet and pressure boundary condition at the outlet to capture the invasion pattern of CO2 through the initially brine-saturated sample. In addition, initially, the first 8% of the grid blocks were filled with the non-wetting phase. At the low capillary number, the average inlet fluid velocity was 0.014 m/s, while for the high capillary number the inlet velocity was set at 0.0465 m/s. Figure 5.10 shows the results produced by the LB simulation, where the CO2 front was injected into the sample from the left at an average wetting-phase saturation of S · =0.50, indicating a fingering pattern and very heterogeneous spatial distribution of CO2 in the pore space. 79 Table 5.9. Properties of CO2-brine two-phase flow system for LB simulation and DNS Contact angle (°) 180 (brine-wet) Interfacial tension (mN/m) 30.0 Brine density (Kg/m D ) 1100 Brine kin. viscosity (m : /s) 1×10 r f CO2 density (Kg/m D ) 1100 CO2 kin. viscosity (m : /s) 1×10 r O Capillary number 1×10 r A and 3×10 r P Viscosity ratio 10 In order to have a quantitative comparison and consistency check of the results produced by the LB and DNS methods, we compare the dependence of the CO2 saturation on the distance from the inlet along the direction of macroscopic flow. The results are shown in Figure 5.11 where CO2 saturation represents an average taken over the cross-sectional area. The results produced by the two methods follow one another very closely. In addition, we computed the change in the brine saturation in the sample volume over time for the two Ca numbers. The results, shown in Figure 5.12, indicate again that the two simulation methods provide consistent and closely agreeing results. This gave us confidence that a comparison between the relative permeabilities produced by the two methods, as well as those produced by the PNM, is viable and meaningful. Figure 5.10. CO2 invasion pattern (blue) computed by the LB simulation for capillary number 𝐂𝐚=𝟏𝟎 r 𝟒 at brine saturation 𝐒𝐰=𝟎.𝟓𝟎. Brine is shown in red and rock is not shown. 80 Figure 5.11. Comparison of the computed CO2 saturation for capillary number 𝐂𝐚=𝟏×𝟏𝟎 r 𝟒 and its dependence on the distance from the inlet at brine saturation 𝐒𝐰=𝟎.𝟓𝟎 during drainage. Figure 5.12. Time-dependence of the computed brine saturation obtained by the LB simulation (continuous curves) and the DNS (dashed curves) during drainage. Figure 5.13 compares the spatial distributions of CO2, injected from left side, and brine for two capillary numbers produced by the DNS, when the latter has reached its residual saturation. Consistent with what is known for the oil-water pairs in porous media, the displacement pattern at high capillary number is more uniform, better connected and piston-like with a lower residual saturation for the brine phase, than that obtained with the lower capillary number that exhibits fingering pattern with a fractal structure. Perhaps, this can be seen better if we consider a side view of the displacement patterns, shown in Figure 5.14, which are consistent with those shown in Figure 5.13. Thus, depending on the heterogeneity of the pore space, at lower capillary numbers 81 the capillary fingering effect can be strongly dominant. In that case, there can be some fluctuations in the relative permeabilities, unlike the smoother-varying values for high capillary numbers. It may also indicate that obtaining smoothly varying relative permeabilities for low capillary numbers entails using larger REVs. Note that at the end of the simulations, the calculated brine residual saturations for the high and low capillary number are 0.30 and 0.50, reached at times 54 ms and 122 ms, respectively. Figure 5.13. 3D representation of the brine residual saturation (red) for 𝐂𝐚=𝟏×𝟏𝟎 r 𝟒 (left) and 𝐂𝐚=𝟑×𝟏𝟎 r 𝟓 (right) capillary numbers. The residual brine saturations are Sw = 0.30 and 0.50, respectively, reached at times 54 ms and 122 ms. Figure 5.14. Side view (inlet on the left and outlet on the right side of each image) of the distribution of the brine (red) and CO2 in the pore space for 𝐂𝐚=𝟏×𝟏𝟎 r 𝟒 (left) and 𝐂𝐚= 𝟑×𝟏𝟎 r 𝟓 (right) at brine saturation, 𝐒 𝐰 =𝟎.𝟒𝟕𝟓. In Figure 5.15, we compare the computed relative permeabilities of the CO2-brine pair for the lower capillary number that we simulated. All the relative permeabilities computed by the three methods are in good agreement with each other, although the methodologies are completely different. Figure 5.16 compares the relative permeabilities computed by the LB method and the 82 DNS for the higher capillary number. In this case too, the results are in agreement with each other. Since our PNM simulator is designed for quasi-static displacement by an IP-like algorithm, it could not be used for simulating an arbitrary value of the capillary number when the magnitude of the viscous forces is competitive with the capillary forces. In such case, one must use a dynamic PN simulator [87]. Figure 5.15. CO2-brine relative permeability curves in low capillary number, 𝐂𝐚=𝟑×𝟏𝟎 r 𝟓 . Figure 5.16. CO2-brine relative permeability curves for high capillary number, 𝐂𝐚=𝟏×𝟏𝟎 r 𝟒 . It should be pointed out that, to compute the relative permeabilities by the LB simulator for each brine saturation, a separate simulation should be carried out, whereas they are computed by one complete simulation when the OpenFOAM poreFOAM package [57] is used in the DNS, as it 83 simulates the entire process, from the beginning of injection of CO2 to reaching the brine residual saturation. But, if the goal is to compute the quantities of interest for a single saturation, then the LB method is more efficient. The LB method has, however, a problem with spurious velocity, which is resolved in the DNS simulator that we employ in this study by using the FSF approach, instead of the CSF method (see the earlier discussions) to account for the interfacial forces. Note also that the PNM simulation, particularly for low capillary numbers, does not require high- performance computational resources. Thus, it is still a reliable method for low Ca number systems, if the PN used is accurate representation of the pore structure. 5.6 The importance of resolution of the computational grid An important point should, however, be emphasized. If the LB simulation is used with an image of a porous medium, one must make sure that the resolution of the lattice used is high enough; that is, the results must be independent of the resolution. For example, in the calculations that we carried out, it was not enough to have one lattice unit per voxel and, thus, the lattice with higher resolutions was needed. Figure 5.17 compares the results computed by the LB simulation obtained with two lattice resolutions for two Ca numbers. In the low-resolution simulation, one lattice unit was defined for each void space voxel, whereas eight lattice units were utilized for each voxel of the input geometry in the high-resolution simulation. Figure 5.17 indicates that only when high- resolution lattices are used, do the LB results converge to those obtained by the DNS. Figure 5.17. Effect of resolution of lattice in LB simulation on the resulting CO2-brine relative permeability curves for 𝐂𝐚=𝟏×𝟏𝟎 r 𝟒 (right) and 𝐂𝐚=𝟑×𝟏𝟎 r 𝟓 (left). 84 5.7 Summary and Conclusions In order to study two-phase flow of CO2-brine in the three-dimensional image of a heterogeneous porous medium, the Mt. Simon sandstone, we first carried out a detailed quantitative study of the pore structure on several subsamples from the rock. Then, we used three distinct computational approaches, namely, the lattice-Boltzmann method, direct numerical simulation of the Navier- Stokes equations, and a pore-network model extracted from the image. The main process simulated was drainage, displacement of brine by CO2. The relative permeabilities of the two fluids for two capillary numbers were computed and compared. Provided that the computational grid in the DNS and the lattice used in the LB simulation have high-enough resolution, the computed relative permeabilities agree very closely. The DNS approach requires, however, a single drainage simulation to compute the relative permeabilities over the entire intended range of saturation, whereas the LB approach needs a separate steady-state simulation for each saturation and, therefore, it requires more computational resources. In addition, the difference in the results produced by the DNS and LB may be due to the different formulations used for the capillary forces parallel to the interfaces. The FSF formulation used in the DNS method eliminates non-physical velocities, whereas the CSF formulation employed in the LB simulation results in nonphysical currents, especially in complex geometries. The relative permeabilities computed by the PNM at a low capillary number also agree with those obtained by the LB simulation and the DNS, although the PNM does not need any high-performance computational resources. Therefore, the question of which method to use for such simulations should be addressed based mainly on the computational time that they need, and the computational resources that one has access too. In addition, one should carefully examine the effect of the resolution of the lattice used in the LB simulation. Improving the interfacial surface formulation in the LB simulation is expected to improve its accuracy, leading to much closer agreement with the results obtained by the DNS method. 85 Bibliography 1. Zoback, M.D. and S.M. Gorelick, Earthquake triggering and large-scale geologic storage of carbon dioxide. Proceedings of the National Academy of Sciences, 2012. 109(26): p. 10164-10168. 2. Arts, R., et al., Seismic monitoring at the Sleipner underground CO2 storage site (North Sea). Geological Society, London, Special Publications, 2004. 233(1): p. 181-191. 3. Verdon, J.P., Significance for secure CO2 storage of earthquakes induced by fluid injection. Environmental Research Letters, 2014. 9(6): p. 064022. 4. Juanes, R., et al., Impact of relative permeability hysteresis on geological CO2 storage. Water resources research, 2006. 42(12). 5. Arts, R., et al. -Monitoring of CO 2 Injected at Sleipner Using Time Lapse Seismic Data. in Greenhouse Gas Control Technologies-6th International Conference. 2003. Elsevier. 6. Rubino, J., D.R. Velis, and M.D. Sacchi, Numerical analysis of wave‐induced fluid flow effects on seismic data: Application to monitoring of CO2 storage at the Sleipner field. Journal of Geophysical Research: Solid Earth, 2011. 116(B3). 7. Iding, M. and P. Ringrose, Evaluating the impact of fractures on the performance of the In Salah CO2 storage site. International Journal of Greenhouse Gas Control, 2010. 4(2): p. 242-248. 8. Birkholzer, J.T., Q. Zhou, and C.-F. Tsang, Large-scale impact of CO2 storage in deep saline aquifers: A sensitivity study on pressure response in stratified systems. International Journal of Greenhouse Gas Control, 2009. 3(2): p. 181-194. 9. Gelhar, L.W., C. Welty, and K.R. Rehfeldt, A critical review of data on field‐scale dispersion in aquifers. Water resources research, 1992. 28(7): p. 1955-1974. 10. Reeves, S. and G. Koperna, Geologic Sequestration of CO2 in Deep, Unmineable Coalbeds: An Integrated Researdh and Commercial-Scale Field Demonstration Project. 2008, Advanced Resources International, Incorporated. 11. Kang, Q., et al., Pore scale modeling of reactive transport involved in geologic CO2 sequestration. Transport in porous media, 2010. 82(1): p. 197-213. 86 12. Busch, A., et al., On sorption and swelling of CO 2 in clays. Geomechanics and Geophysics for Geo-Energy and Geo-Resources, 2016. 2(2): p. 111-130. 13. Brigatti, M.F., E. Galan, and B. Theng, Structure and mineralogy of clay minerals, in Developments in clay science. 2013, Elsevier. p. 21-81. 14. Cnenr-ns, E.W., THE DISTRIBUTION AND IDENTIFICATION OF MIXED-LAYER CLAYS IN SEDIMENTARY ROCKS. 15. Altaner, S.P. and R.F. Ylagan, Comparison of structural models of mixed-layer illite/smectite and reaction mechanisms of smectite illitization. Clays and Clay Minerals, 1997. 45(4): p. 517-533. 16. Reynolds, R.C. and J. Hower, The nature of interlayering in mixed-layer illite- montmorillonites. Clays and Clay minerals, 1970. 18(1): p. 25-36. 17. Teich-McGoldrick, S.L., et al., Swelling properties of montmorillonite and beidellite clay minerals from molecular simulation: comparison of temperature, interlayer cation, and charge location effects. The Journal of Physical Chemistry C, 2015. 119(36): p. 20880- 20891. 18. Boek, E., P. Coveney, and N. Skipper, Monte Carlo molecular modeling studies of hydrated Li-, Na-, and K-smectites: Understanding the role of potassium as a clay swelling inhibitor. Journal of the American Chemical Society, 1995. 117(50): p. 12608-12617. 19. Liu, X., et al., Effects of layer-charge distribution on the thermodynamic and microscopic properties of Cs-smectite. Geochimica et Cosmochimica Acta, 2008. 72(7): p. 1837-1847. 20. Chávez-Páez, M., L. DePablo, and J. DePablo, Monte Carlo simulations of Ca– montmorillonite hydrates. The Journal of Chemical Physics, 2001. 114(24): p. 10948- 10953. 21. Denis, J., et al., Influence of potassium concentration the swelling and compaction ofmixed (Na, K) ion-exchanged montmorillonite. 1991. 22. Skipper, N., G. Sposito, and F.-R.C. Chang, Monte Carlo simulation of interlayer molecular structure in swelling clay minerals. 2. Monolayer hydrates. Clays and Clay minerals, 1995. 43(3): p. 294-303. 23. Young, D.A. and D.E. Smith, Simulations of clay mineral swelling and hydration: dependence upon interlayer ion size and charge. The Journal of Physical Chemistry B, 2000. 104(39): p. 9163-9170. 87 24. Degrève, L., S.M. Vechi, and C.Q. Junior, The hydration structure of the Na+ and K+ ions and the selectivity of their ionic channels. Biochimica et Biophysica Acta (BBA)- Bioenergetics, 1996. 1274(3): p. 149-156. 25. Hensen, E., et al., Adsorption isotherms of water in Li–, Na–, and K–montmorillonite by molecular simulation. The Journal of Chemical Physics, 2001. 115(7): p. 3322-3329. 26. Ferrage, E., et al., Hydration properties and interlayer organization of water and ions in synthetic Na-smectite with tetrahedral layer charge. Part 2. Toward a precise coupling between molecular simulations and diffraction data. The Journal of Physical Chemistry C, 2011. 115(5): p. 1867-1881. 27. Martins, D.M., et al., Toward modeling clay mineral nanoparticles: The edge surfaces of pyrophyllite and their interaction with water. The Journal of Physical Chemistry C, 2014. 118(47): p. 27308-27317. 28. Newton, A.G. and G. Sposito, Molecular dynamics simulations of pyrophyllite edge surfaces: structure, surface energies, and solvent accessibility. Clays and Clay Minerals, 2015. 63(4): p. 277-289. 29. Newton, A.G., K.D. Kwon, and D.-K. Cheong, Edge structure of montmorillonite from atomistic simulations. Minerals, 2016. 6(2): p. 25. 30. Cygan, R.T., J.-J. Liang, and A.G. Kalinichev, Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field. The Journal of Physical Chemistry B, 2004. 108(4): p. 1255-1266. 31. Makaremi, M., et al., Multiphase Monte Carlo and molecular dynamics simulations of water and CO2 intercalation in montmorillonite and beidellite. The Journal of Physical Chemistry C, 2015. 119(27): p. 15112-15124. 32. Rao, Q. and Y. Leng, Effect of layer charge on CO2 and H2O intercalations in swelling clays. Langmuir, 2016. 32(44): p. 11366-11374. 33. Myshakin, E.M., et al., Molecular dynamics simulations of carbon dioxide intercalation in hydrated Na-montmorillonite. The Journal of Physical Chemistry C, 2013. 117(21): p. 11028-11039. 34. Rao, Q. and Y. Leng, Molecular understanding of CO2 and H2O in a montmorillonite clay interlayer under CO2 geological sequestration conditions. The Journal of Physical Chemistry C, 2016. 120(5): p. 2642-2654. 88 35. Bibi, I., et al., Clay minerals: structure, chemistry and significance in contaminated environments and geological CO2 sequestration. 2016. 36. Mooney, R., A. Keenan, and L. Wood, Adsorption of water vapor by montmorillonite. II. Effect of exchangeable ions and lattice swelling as measured by X-ray diffraction. Journal of the American Chemical Society, 1952. 74(6): p. 1371-1374. 37. Brown, G. and G.W. Brindley, Crystal structures of clay minerals and their X-ray identification. Vol. 5. 1980: Mineralogical Society London. 38. Keren, R. and I. Shainberg, Water vapor isotherms and heat of immersion of na/ca- montmorillonite systems: I, homoionic clay. Clays and Clay Minerals, 1975. 23(3): p. 193- 200. 39. Suquet, H., C. De La Calle, and H. Pezerat, Swelling and structural organization of saponite. Clays and Clay Minerals, 1975. 23(1): p. 1-9. 40. Cancela, G.D., et al., Adsorption of water vapor by homoionic montmorillonites. Heats of adsorption and desorption. Journal of Colloid and Interface Science, 1997. 185(2): p. 343- 354. 41. Bourg, I.C. and G. Sposito, Molecular dynamics simulations of the electrical double layer on smectite surfaces contacting concentrated mixed electrolyte (NaCl–CaCl 2) solutions. Journal of colloid and interface science, 2011. 360(2): p. 701-715. 42. Cases, J., et al., Mechanism of adsorption and desorption of water vapor by homoionic montmorillonite; 3, The Mg (super 2+), Ca (super 2+), and Ba (super 3+) exchanged forms. Clays and Clay Minerals, 1997. 45(1): p. 8-22. 43. Laird, D.A., C. Shang, and M.L. Thompson, Hysteresis in crystalline swelling of smectites. Journal of Colloid and Interface Science, 1995. 171(1): p. 240-245. 44. Zhang, F., et al., The effect of temperature on the swelling of montmorillonite. Clay Minerals, 1993. 28: p. 25-25. 45. Viani, B.E., P.F. Low, and C.B. Roth, Direct measurement of the relation between interlayer force and interlayer distance in the swelling of montmorillonite. Journal of colloid and interface science, 1983. 96(1): p. 229-244. 46. Sato, T., T. Watanabe, and R. Otsuka, Effects of layer charge, charge location, and energy change on expansion properties of dioctahedral smectites. Clays and Clay Minerals, 1992. 40(1): p. 103-113. 89 47. Meunier, A., A. Inoue, and D. Beaufort, Chemiographic analysis of trioctahedral smectite- to-chlorite conversion series from the Ohyu Caldera, Japan. Clays and Clay Minerals, 1991. 39(4): p. 409-415. 48. Weaver, C.E. and L.D. Pollard, The chemistry of clay minerals. Vol. 15. 2011: Elsevier. 49. Fu, M., Z. Zhang, and P. Low, Changes in the properties of a montmorillonite-water system during the adsorption and desorption of water: hysteresis. Clays and Clay Minerals, 1990. 38(5): p. 485-492. 50. Slade, P., J. Quirk, and K. Norrish, Crystalline swelling of smectite samples in concentrated NaCl solutions in relation to layer charge. Clays and Clay Minerals, 1991. 39(3): p. 234- 238. 51. Greathouse, J. and G. Sposito, Monte Carlo and Molecular Dynamics Studies of Interlayer Structure in Li (H2O) 3− Smectites. The Journal of Physical Chemistry B, 1998. 102(13): p. 2406-2414. 52. Tambach, T.J., E.J. Hensen, and B. Smit, Molecular simulations of swelling clay minerals. The Journal of Physical Chemistry B, 2004. 108(23): p. 7586-7596. 53. Chávez-Páez, M., et al., Monte Carlo simulations of Wyoming sodium montmorillonite hydrates. The Journal of Chemical Physics, 2001. 114(3): p. 1405-1413. 54. Smith, D.E., Y. Wang, and H.D. Whitley, Molecular simulations of hydration and swelling in clay minerals. Fluid phase equilibria, 2004. 222: p. 189-194. 55. Whitley, H.D. and D.E. Smith, Free energy, energy, and entropy of swelling in Cs–, Na–, and Sr–montmorillonite clays. The Journal of chemical physics, 2004. 120(11): p. 5387- 5395. 56. Foster, M.D., Geochemical studies of clay minerals: II—relation between ionic substitution and swelling in montmorillonites. American Mineralogist: Journal of Earth and Planetary Materials, 1953. 38(11-12): p. 994-1006. 57. Raeini, A.Q., M.J. Blunt, and B. Bijeljic, Direct simulations of two-phase flow on micro- CT images of porous media and upscaling of pore-scale forces. Advances in Water Resources, 2014. 74: p. 116-126. 58. Pruess, K. and J. Garcia, Multiphase flow dynamics during CO 2 disposal into saline aquifers. Environmental Geology, 2002. 42(2-3): p. 282-295. 90 59. Friedlingstein, P. and S. Solomon, Contributions of past and present human generations to committed warming caused by carbon dioxide. Proceedings of the National Academy of Sciences, 2005. 102(31): p. 10832-10836. 60. Metz, B., O. Davidson, and H. De Coninck, Carbon dioxide capture and storage: special report of the intergovernmental panel on climate change. 2005: Cambridge University Press. 61. Nordbotten, J.M. and M.A. Celia, Geological storage of CO2: modeling approaches for large-scale simulation. 2011: John Wiley & Sons. 62. Shukla, R., et al., A review of studies on CO2 sequestration and caprock integrity. Fuel, 2010. 89(10): p. 2651-2664. 63. Juanes, R., et al., Impact of relative permeability hysteresis on geological CO2 storage. Water resources research, 2006. 42(12). 64. Finley, R.J., An overview of the Illinois Basin–Decatur project. Greenhouse Gases: Science and Technology, 2014. 4(5): p. 571-579. 65. Celia, M.A., P.C. Reeves, and L.A. Ferrand, Recent advances in pore scale models for multiphase flow in porous media. Reviews of Geophysics, 1995. 33(S2): p. 1049-1057. 66. Blunt, M.J., et al., Pore-scale imaging and modelling. Advances in Water Resources, 2013. 51: p. 197-216. 67. Dunsmuir, J.H., et al. X-ray microtomography: a new tool for the characterization of porous media. in SPE annual technical conference and exhibition. 1991. Society of Petroleum Engineers. 68. Flannery, B.P., et al., Three-dimensional X-ray microtomography. Science, 1987. 237(4821): p. 1439-1444. 69. Arns, C.H., et al., Computation of linear elastic properties from microtomographic images: Methodology and agreement between theory and experiment. Geophysics, 2002. 67(5): p. 1396-1405. 70. Øren, P.-E. and S. Bakke, Process based reconstruction of sandstones and prediction of transport properties. Transport in porous media, 2002. 46(2-3): p. 311-343. 71. Ramstad, T., P.-E. Øren, and S. Bakke, Simulation of two-phase flow in reservoir rocks using a lattice Boltzmann method. Spe Journal, 2010. 15(04): p. 917-927. 91 72. Krevor, S.C., et al., Relative permeability and trapping of CO2 and water in sandstone rocks at reservoir conditions. Water resources research, 2012. 48(2). 73. Ferrari, A. and I. Lunati, Direct numerical simulations of interface dynamics to link capillary pressure and total surface energy. Advances in water resources, 2013. 57: p. 19- 31. 74. Huang, H., P. Meakin, and M. Liu, Computer simulation of two‐phase immiscible fluid motion in unsaturated complex fractures using a volume of fluid method. Water resources research, 2005. 41(12). 75. Rabbani, H.S., et al., New insights on the complex dynamics of two-phase flow in porous media under intermediate-wet conditions. Scientific reports, 2017. 7(1): p. 4584. 76. Boek, E.S. and M. Venturoli, Lattice-Boltzmann studies of fluid flow in porous media with realistic rock geometries. Computers & Mathematics with Applications, 2010. 59(7): p. 2305-2314. 77. Ferreol, B. and D.H. Rothman, Lattice-Boltzmann simulations of flow through Fontainebleau sandstone, in Multiphase flow in porous media. 1995, Springer. p. 3-20. 78. Pan, C., M. Hilpert, and C. Miller, Lattice‐Boltzmann simulation of two‐phase flow in porous media. Water Resources Research, 2004. 40(1). 79. Porter, M.L., M.G. Schaap, and D. Wildenschild, Lattice-Boltzmann simulations of the capillary pressure–saturation–interfacial area relationship for porous media. Advances in Water Resources, 2009. 32(11): p. 1632-1640. 80. Ramstad, T., et al., Relative permeability calculations from two-phase flow simulations directly on digital images of porous rocks. Transport in Porous Media, 2012. 94(2): p. 487- 504. 81. Rothman, D.H., Macroscopic laws for immiscible two‐phase flow in porous media: Results From numerical experiments. Journal of Geophysical Research: Solid Earth, 1990. 95(B6): p. 8663-8674. 82. Ahrenholz, B., et al., Prediction of capillary hysteresis in a porous material using lattice- Boltzmann methods and comparison to experimental data and a morphological pore network model. Advances in Water Resources, 2008. 31(9): p. 1151-1173. 92 83. Chen, Y., et al., Lattice Boltzmann simulations of liquid CO2 displacing water in a 2D heterogeneous micromodel at reservoir pressure conditions. Journal of contaminant hydrology, 2018. 212: p. 14-27. 84. Tölke, J., S. Freudiger, and M. Krafczyk, An adaptive scheme using hierarchical grids for lattice Boltzmann multi-phase flow simulations. Computers & fluids, 2006. 35(8-9): p. 820-830. 85. Gueyffier, D., et al., Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows. Journal of Computational physics, 1999. 152(2): p. 423-456. 86. Sahimi, M., Flow and transport in porous media and fractured rock: from classical methods to modern approaches. 2011: John Wiley & Sons. 87. Joekar-Niasar, V. and S. Hassanizadeh, Analysis of fundamentals of two-phase flow in porous media using dynamic pore-network models: A review. Critical reviews in environmental science and technology, 2012. 42(18): p. 1895-1976. 88. Romanov, V.N., Evidence of irreversible CO 2 intercalation in montmorillonite. International Journal of Greenhouse Gas Control, 2013. 14: p. 220-226. 89. Meleshyn, A. and C. Bunnenberg, Swelling of Na⁄ Mg-montmorillonites and hydration of interlayer cations: A Monte Carlo study. The Journal of chemical physics, 2005. 123(7): p. 074706. 90. Marry, V., et al., Microscopic simulation of structure and dynamics of water and counterions in a monohydrated montmorillonite. The Journal of chemical physics, 2002. 117(7): p. 3454-3463. 91. BURGESS, J., Basic principles of chemical interactions. Ions in Solution, 1988. 45. 92. Michot, L., et al., Mechanism of adsorption and desorption of water vapor by homoionic montmorillonites: 2. The li § na § k § rb § and cs+-exchanged forms. Clays Clay Miner, 1995. 43: p. 324-336. 93. Chiou, C.T. and D.W. Rutherford, Effects of exchanged cation and layer charge on the sorption of water and EGME vapors on montmorillonite clays. Clays and Clay Minerals, 1997. 45: p. 867-880. 94. Rutherford, D.W., C.T. Chiou, and D.D. Eberl, Effects of exchanged cation on the microporosity of montmorillonite. Clays and Clay Minerals, 1997. 45(4): p. 534-543. 93 95. Calvet, R. Hydratation de la montmorillonite et diffusion des cations compensateurs. in Annales Agronomiques. 1973. 96. Cases, J., et al., Mechanism of adsorption and desorption of water vapor by homoionic montmorillonite. 1. The sodium-exchanged form. Langmuir, 1992. 8(11): p. 2730-2739. 97. Boek, E., Molecular dynamics simulations of interlayer structure and mobility in hydrated Li-, Na-and K-montmorillonite clays. Molecular Physics, 2014. 112(9-10): p. 1472-1483. 98. Chang, F.-R.C., N. Skipper, and G. Sposito, Monte Carlo and molecular dynamics simulations of electrical double-layer structure in potassium− montmorillonite hydrates. Langmuir, 1998. 14(5): p. 1201-1207. 99. Lammers, L.N., et al., Molecular dynamics simulations of cesium adsorption on illite nanoparticles. Journal of colloid and interface science, 2017. 490: p. 608-620. 100. Cygan, R.T., V.N. Romanov, and E.M. Myshakin, Molecular simulation of carbon dioxide capture by montmorillonite using an accurate and flexible force field. The Journal of Physical Chemistry C, 2012. 116(24): p. 13079-13091. 101. Yagi, Y., et al., Density dependence of fermi resonance of supercritical carbon dioxide. The Journal of Supercritical Fluids, 1993. 6(3): p. 139-142. 102. Qin, Y., et al., Molecular dynamics simulation of interaction between supercritical CO2 fluid and modified silica surfaces. The Journal of Physical Chemistry C, 2008. 112(33): p. 12815-12824. 103. Rahromostaqim, M. and M. Sahimi, Molecular Dynamics Simulation of Hydration and Swelling of Mixed-Layer Clays. The Journal of Physical Chemistry C, 2018. 122(26): p. 14631-14639. 104. Kadoura, A., A.K. Narayanan Nair, and S. Sun, Molecular dynamics simulations of carbon dioxide, methane, and their mixture in montmorillonite clay hydrates. The Journal of Physical Chemistry C, 2016. 120(23): p. 12517-12529. 105. Botan, A., et al., Carbon dioxide in montmorillonite clay hydrates: thermodynamics, structure, and transport from molecular simulation. The Journal of Physical Chemistry C, 2010. 114(35): p. 14962-14969. 106. Rahromostaqim, M. and M. Sahimi, Molecular Dynamics Simulation of Hydration and Swelling of Mixed-Layer Clays in the Presence of Carbon Dioxide. The Journal of Physical Chemistry C, 2019. 123(7): p. 4243-4255. 94 107. Plimpton, S., Fast parallel algorithms for short-range molecular dynamics. Journal of computational physics, 1995. 117(1): p. 1-19. 108. Smith, D.E., Molecular computer simulations of the swelling properties and interlayer structure of cesium montmorillonite. Langmuir, 1998. 14(20): p. 5959-5967. 109. Denis, J., et al., Influence of potassium concentration on the swelling and compaction of mixed (Na, K) ion-exchanged montmorillonite. Clay Minerals, 1991. 26(2): p. 255-268. 110. Tao, L., et al., Swelling of K+, Na+ and Ca2+-montmorillonites and hydration of interlayer cations: a molecular dynamics simulation. Chinese Physics B, 2010. 19(10): p. 109101. 111. Ngouana W, B.F. and A.G. Kalinichev, Structural arrangements of isomorphic substitutions in smectites: Molecular simulation of the swelling properties, interlayer structure, and dynamics of hydrated Cs–montmorillonite revisited with new clay models. The Journal of Physical Chemistry C, 2014. 118(24): p. 12758-12773. 112. Chang, F.-R.C., et al., Interlayer molecular structure and dynamics in Li-, Na-, and K- montmorillonite-water systems. 1998, ACS Publications. 113. Boulet, P., P.V. Coveney, and S. Stackhouse, Simulation of hydrated Li+-, Na+-and K+- montmorillonite/polymer nanocomposites using large-scale molecular dynamics. Chemical physics letters, 2004. 389(4-6): p. 261-267. 114. Kosakowski, G., S.V. Churakov, and T. Thoenen, Diffusion of Na and Cs in montmorillonite. Clays and Clay Minerals, 2008. 56(2): p. 190-206. 115. Schindelin, J., et al., Fiji: an open-source platform for biological-image analysis. Nature methods, 2012. 9(7): p. 676. 116. Dong, H. and M.J. Blunt, Pore-network extraction from micro-computerized-tomography images. Physical review E, 2009. 80(3): p. 036307. 117. Silin, D. and T. Patzek, Pore space morphology analysis using maximal inscribed spheres. Physica A: Statistical mechanics and its applications, 2006. 371(2): p. 336-360. 118. Al-Kharusi, A.S. and M.J. Blunt, Network extraction from sandstone and carbonate pore space images. Journal of Petroleum Science and Engineering, 2007. 56(4): p. 219-231. 119. Arand, F. and J. Hesser, Accurate and efficient maximal ball algorithm for pore network extraction. Computers & Geosciences, 2017. 101: p. 28-37. 95 120. Raeini, A.Q., B. Bijeljic, and M.J. Blunt, Generalized network modeling: Network extraction as a coarse-scale discretization of the void space of porous media. Physical Review E, 2017. 96(1): p. 013312. 121. Harrigan, T. and R. Mann, Characterization of microstructural anisotropy in orthotropic materials using a second rank tensor. Journal of Materials Science, 1984. 19(3): p. 761- 767. 122. Odgaard, A. and H. Gundersen, Quantification of connectivity in cancellous bone, with special emphasis on 3-D reconstructions. Bone, 1993. 14(2): p. 173-182. 123. Odgaard, A., Three-dimensional methods for quantification of cancellous bone architecture. Bone, 1997. 20(4): p. 315-328. 124. Toriwaki, J. and T. Yonekura, Euler number and connectivity indexes of a three dimensional digital picture. FORMA-TOKYO-, 2002. 17(3): p. 183-209. 125. Huang, H., M. Sukop, and X. Lu, Multiphase lattice Boltzmann methods: Theory and application. 2015: John Wiley & Sons. 126. Liu, H., et al., Multiphase lattice Boltzmann simulations for porous media applications. Computational Geosciences, 2016. 20(4): p. 777-805. 127. Grunau, D., S. Chen, and K. Eggert, A lattice Boltzmann model for multiphase fluid flows. Physics of Fluids A: Fluid Dynamics, 1993. 5(10): p. 2557-2562. 128. Gunstensen, A.K., et al., Lattice Boltzmann model of immiscible fluids. Physical Review A, 1991. 43(8): p. 4320. 129. Jiang, F. and T. Tsuji, Estimation of three‐phase relative permeability by simulating fluid dynamics directly on rock‐microstructure images. Water Resources Research, 2017. 53(1): p. 11-32. 130. Tölke, J., Lattice Boltzmann simulations of binary fluid flow through porous media. Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 2002. 360(1792): p. 535-545. 131. Lallemand, P. and L.-S. Luo, Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability. Physical Review E, 2000. 61(6): p. 6546. 96 132. d'Humieres, D., Multiple–relaxation–time lattice Boltzmann models in three dimensions. Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 2002. 360(1792): p. 437-451. 133. Sukop, M., DT Thorne, Jr. Lattice Boltzmann Modeling Lattice Boltzmann Modeling. 2006: Springer. 134. Honarpour, M., L. Koederitz, and A. Harvey, Relative Permeability of Petroleum Reservoirs CRC Press. Inc. Boca Raton, Florida, 1986. 135. Qaseminejad Raeini, A., Modelling multiphase flow through micro-CT images of the pore space. 2013, Imperial College London. 136. Ubbink, O., Numerical prediction of two fluid systems with sharp interfaces. 1997. 137. Francois, M.M., et al., A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework. Journal of Computational Physics, 2006. 213(1): p. 141-173. 138. Raeini, A.Q., M.J. Blunt, and B. Bijeljic, Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method. Journal of Computational Physics, 2012. 231(17): p. 5653-5668. 139. Wilkinson, D. and J.F. Willemsen, Invasion percolation: a new form of percolation theory. Journal of Physics A: Mathematical and General, 1983. 16(14): p. 3365. 140. Sheppard, A.P., et al., Invasion percolation: new algorithms and universality classes. Journal of Physics A: Mathematical and General, 1999. 32(49): p. L521. 141. Knackstedt, M.A., M. Sahimi, and A.P. Sheppard, Invasion percolation with long-range correlations: First-order phase transition and nonuniversal scaling properties. Physical Review E, 2000. 61(5): p. 4920. 142. Valvatne, P.H. and M.J. Blunt, Predictive pore‐scale modeling of two‐phase flow in mixed wet media. Water resources research, 2004. 40(7). 143. Oak, M., L. Baker, and D. Thomas, Three-phase relative permeability of Berea sandstone. Journal of Petroleum Technology, 1990. 42(08): p. 1,054-1,061.

Abstract (if available)

Abstract Flow, transport, swelling, and deformation of porous media as a result of injecting CO₂ in the media, and in the presence of brine, are important phenomena that have been studied over the past decade or so, because they are relevant to sequestration of CO₂ in large-scale porous formations, such as depleted oil and gas reservoirs. This thesis studies these phenomena at two distinct length scales, namely, molecular and core scales, using extensive computer simulations. We present the results of computer simulations of single- and two-phase flow of water and CO₂ in porous media, at both molecular and continuum scales. ❧ At molecular scale, molecular dynamics (MD) simulations have been carried out to understand swelling of mixed-layer (ML) clays. The phenomena are studied in the two most common ML clays, i.e., illite-montmorillonite (I-MMT) and chlorite-montmorillonite (CH-MMT)

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

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

PDF

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

PDF

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

PDF

Molecular dynamics studies of protein aggregation in unbounded and confined media

PDF

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

PDF

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

PDF

Machine-learning approaches for modeling of complex materials and media

PDF

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

PDF

Laser manipulation of atomic and molecular flows

PDF

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

PDF

A process-based molecular model of nano-porous silicon carbide membranes

PDF

Modeling and simulation of complex recovery processes

PDF

A study of dispersive mixing and flow based lumping/delumping in gas injection processes

PDF

Hybrid physics-based and data-driven computational approaches for multi-scale hydrologic systems under heterogeneity

PDF

Molecular-scale studies of mechanical phenomena at the interface between two solid surfaces: from high performance friction to superlubricity and flash heating

PDF

Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials

PDF

Asphaltene deposition during co₂ injection in enhanced oil recovery applications

PDF

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

PDF

In situ studies of the thermal evolution of the structure and sorption properties of Mg-Al-CO3 layered double hydroxide

PDF

Supersonic flow study of a capsule geometry using large-eddy simulation

Asset Metadata

Creator Rahromostaqim, Mahsa (author)

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

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Chemical Engineering

Publication Date 01/24/2020

Defense Date 12/19/2019

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

Tag continuum scale,micro CT image,mixed layer clay,molecular dynamics,molecular scale,OAI-PMH Harvest,simulation,two-phase flow

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Sahimi, Muhammad (committee chair), De Barros, Felipe (committee member), Tsotsis, Theo (committee member)

Creator Email mahsa.rahro@gmail.com,rahromos@usc.edu

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

Unique identifier UC11673094

Identifier etd-Rahromosta-8128.pdf (filename),usctheses-c89-261962 (legacy record id)

Legacy Identifier etd-Rahromosta-8128.pdf

Dmrecord 261962

Document Type Dissertation

Rights Rahromostaqim, Mahsa

Type texts

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

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

Repository Name University of Southern California Digital Library

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