Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Multiscale and multiresolution approach to characterization and modeling of porous media: From pore to field scale
(USC Thesis Other)
Multiscale and multiresolution approach to characterization and modeling of porous media: From pore to field scale
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
Multiscale and Multiresolution Approach to Characterization and Modeling of Porous Media: From Pore to Field Scale by Hassan Dashtian A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (CHEMICAL ENGINEERING) August 2018 Copyright 2018 Hassan Dashtian Dedication To Sahar and my parents. ii Acknowledgements I would like to express my deepest appreciation to all the people who have assisted me during my Ph.D. studies with their support, help and encouragement. Above all, I want to thank my M.S. and Ph.D. dissertation advisor, Professor Muhammad Sahimi. His unique and comprehensive understanding of porous materials has been an inspiration for me to learn and apply new and novel approaches to study small and large-scale porous materials. I also appreciate him for giving me the freedom to explore, learn and apply many dierent research topics, techniques and simulation approaches, in addition to allowing me to taking courses related to high performance computing, simulation and visualization. I have learned many things from him beyond the boarders of research and academic. I would also like to thank my qualifying exam committee members, Professor Aiichiro Nakano, Professor Katherine S.-F. Shing, Professor Kristian Jessen and Professor Felipe de Barros. I appreciate their useful feedback, guidance, support, comments, and questions that helped me conduct a better research. I thoroughly enjoyed the many courses that I took during my graduate studies with Professors Sahimi and Nakano. Professor Nakano was always available to discuss a wide variety of topics that I learned so much from; I will be always be grateful to him. I would also like to thank Professor Nima Shokri of the University of Manchester. It was a pleasure to work with him. I have learned many things related to evaporation and salt precipitation from his useful comments and feedback. iii I also would like to thank the administrative assistants and sta of our department, Andy Chen and the late Martin Olekszyk. I want to thank my friends that have always been available to help me. Without their support I would not be able to complete my Ph.D. thesis. So, thank you: Soheil Negahbani, Maryam Nakhjiri, Majid Monji, Mohammad Evazi, Shahram Farhadinia, Sasan Dabir, Ashkan Garshasbi and Nariman Piroozan. I would also like to express my most heartfelt appreciation and gratitude to my wife, Sahar, who has directly and indirectly contributed to this dissertation in so many ways. Finally, but most of all, I want to thank my family who helped me and supported me all the way. In particular, I would like to thank my parents who have sacriced so much for my benet, and my sister in-law, Sonia, for her unbounded and sel ess support during the roughest days of my life at USC. This work was supported as part of the Center for Geologic Storage of CO 2 , an Energy Frontier Research Center, funded by the U.S. Department of Energy, Oce of Science, Basic Energy Sciences, under Award Number DE-SC0C12504. iv Table of Contents Dedication ii Acknowledgements iii List of Tables vii List of Figures viii Abstract xii Chapter 1: Introduction 1 1.1 Background and Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Pore-Network Scale versus Field Scale . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Pore-Network Scale: Salt Crystallization and Precipitation . . . . . . . . . . . . . . 3 1.4 Field-Scale: Multifractal and Multiresolution Analysis of Data . . . . . . . . . . . 5 1.5 Salt Precipitation in Porous Media . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.6 Salt Precipitation Due to CO 2 Injection . . . . . . . . . . . . . . . . . . . . . . . . 7 1.7 Drying of Porous Media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Chapter 2: Nucleation of Salt Crystals in Clay Minerals: Molecular Dynamics Simulation 9 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Materials and Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Chapter 3: Ecient Simulation of Fluid Flow and Transport in Heterogeneous Media Using Graphics Processing Units (GPUs) 22 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2 Pore and Random Resistor Networks . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3 CPU and GPUs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.4 Solving for the Pressure (Voltage) Distribution on CUDA . . . . . . . . . . . . . . 29 3.5 Generation of the Fractional Brownian Motion . . . . . . . . . . . . . . . . . . . . 33 3.6 Solving for the Pressure/Voltage Distribution with a GPU . . . . . . . . . . . . . . 34 3.7 Eciency of the Computations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.8 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.8.1 Isotropic networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.8.2 Anisotropic networks and permeability anisotropy . . . . . . . . . . . . . . 41 3.9 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.10 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 v Chapter 4: Drying of Layered Porous Medium: The Eect of Pore-Size Correla- tion and Disorder 45 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.2 Generating Correlated Pore Network . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.3 Pore Network Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.4 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.4.1 Long-range correlations and short-range disorders . . . . . . . . . . . . . . 58 4.4.2 Stratication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 Chapter 5: Pore-network model of evaporation-induced salt precipitation in porous media: the eect of correlations and heterogeneity 64 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.2 Pore-network model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.3 Pore-network model and the governing equations . . . . . . . . . . . . . . . . . . . 73 5.4 Computational algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 5.5 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 Chapter 6: Non-universality of the Archie exponent due to multifractality of re- sistivity well logs 94 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6.2 The Reservoirs and the Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 6.3 Analysis of the Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.3.1 Identifying the Clay Sectors . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 6.3.2 Analysis of Correlations in the Data . . . . . . . . . . . . . . . . . . . . . . 100 6.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 6.5 Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 Bibliography 110 vi List of Tables 3.1 Eciency and speed-up of the MPI implementation of FBM algorithm on four processes. The times are in seconds. . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.1 Pore size distribution properties of dierent cases. . . . . . . . . . . . . . . . . . . 57 5.1 Pore size distribution properties of dierent cases. . . . . . . . . . . . . . . . . . . 79 vii List of Figures 2.1 (a) and (c) show the kaolinite molecular structure, while (b) and (d) depict that of Na-MMT clays. Water molecules in (a) and (b) are shown by sticks. Colors represent oxygen (red), hydrogen (white), Si (yellow), Mg [green, only in (b)], Al [pink in (a) and (b)] and [green in (c) and (d)], Na (purple), and Cl (green; very small). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Top view of the positions of Na + (pink) and Cl (yellow) within the interlayer space of kaolinite and Na-MMT. Time increases from (a) to (d). (a) Early positions (after 0.1 ns) of the ions after the system reached equilibrium. (b) The initial spatial ordering of the ions after 7 ns for kaolinite and 6 ns for Na-MMT. (c) and (d) show the nucleation steps and spatial ordering after 30 and 70 ns for, respectively, kaolinite and Na-MMT. Colors show Na + (pink) and Cl (yellow). . . . . . . . . . 15 2.3 Evolution of the radial distribution function (RDF), indicating nucleation and spa- tial ordering of the ions in the interlayer space of kaolinite (a)-(d) and Na-MMT (e)-(h). (a) and (e) show the early-stage (< 1 ns) RDF of the ions. . . . . . . . . . 16 2.4 Initial density proles of the salt ions and oxygen in water along thez axis (perpen- dicular to the clay planes) for salt concentration of 30.9% weight. z< 6 represents locations inside the lower clay sheet. . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.5 Evolution with time of the mean-square displacements (MSD) of the ions in the interlayer space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.1 GPU architecture hierarchy with four streaming multiprocessors (SM). . . . . . . . 28 3.2 The CUDA memory model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.3 Mixed-precision algorithm for the conjugate-gradient method that was used in the computations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.4 Communications between four nodes (left) and point-to-point communication model. 34 3.5 Examples of the networks with correlated conductances. The top and bottom rows show, respectively, the networks with the Hurst exponent of H = 0:75 and 0.35, and x = y of (a) and (d) 1; (b) and (e) 1/5, and (c) and (f) 5. . . . . . . . . . . . 36 viii 3.6 a) Dependence of the speed-up on the size of isotropic networks withH = 0:75 and 0.35 b) Speed up for one iteration in CG algorithm. . . . . . . . . . . . . . . . . . 37 3.7 a) Dependence of the speed-up on the size of isotropic networks and the anisotropy parameter x = y . The Hurst exponent is H = 0:35 b) Speed up for one iteration in CG algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.8 The speed-up of the computations with a random resistor network with 3:6 10 6 nodes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.9 a) Dependence of the speed-up of the computations for network with 3:6 10 6 nodes on the fraction of the removed bonds, the Hurst exponent and the anisotropy parameter x = y . b) Speed up for one iteration in CG algorithm. . . . . . . . . . . 40 3.10 Pressure distribution in a 50 50 network with the Hurst exponent of H = 0:75 and 0.35, and x = y of (a) and (d) 1; (b) and (e) 1/5, and (c) and (f) 5. . . . . . . 41 3.11 K x =K y ratio for networks with various sizes: (a) H = 0:75 and (b) H = 0:35. . . . 42 4.1 Examples of correlated FBM elds. The top and bottom rows show, respectively, the PNs with the Hurst exponents H = 0:75 and 0.35, and anisotropy x = z = (a) and (d) 1; (b) and (e) 1/10, and (c) and (f) 10. . . . . . . . . . . . . . . . . . . . . 48 4.2 Size distribution of the throats same as Figure 4.1. . . . . . . . . . . . . . . . . . . 48 4.3 The phase distributions in an isotropic pore network with the Hurst exponent H = 0:75 ( x = z = 1), at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.9. The top row of the network is open. . . . . . . . . . . . . . . . . . . . . . . . . 49 4.4 The phase distributions in an isotropic pore network with the Hurst exponent H = 0:35 ( x = z = 1), at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.5 The phase distributions in an isotropic pore network with the Hurst exponent H = 0:75 ( x = z = 1=10), at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.6 The phase distributions in an isotropic pore network with the Hurst exponent H = 0:35 ( x = z = 1=10), at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.7 The phase distributions in an isotropic pore network with the Hurst exponent H = 0:75 ( x = z = 10), at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.8 The phase distributions in an isotropic pore network with the Hurst exponent H = 0:35 ( x = z = 10), at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 ix 4.9 Dynamic evolution of the transversely-averaged brine saturations in six types of pore networks. Black curves show the average saturations; cyan, blue and pink shadings represent, respectively, the times at which vapor saturation is 0.1, 0.2, and 0.6, while red indicates the time at which vapor saturation is 0.9. The colored areas indicate the variations over multiple realizations. . . . . . . . . . . . . . . . . 59 4.10 Dynamic evolution of the average drying rates (black curves), averaged over 10 realizations (shadings). The gray areas indicate variations over multiple realizations. 60 5.1 Examples of pore networks with correlated throat sizes. The top and bottom rows show, respectively, the PNs with the Hurst exponents H = 0:7 and 0.3, and anisotropy x = y = (a) and (d) 1; (b) and (e) 1/4, and (c) and (f) 4. . . . . . . . . 70 5.2 Illustration of the assumption of equilibrium vapor pressure in the central vapor pore, together with a neighboring throat containing liquid. . . . . . . . . . . . . . . 71 5.3 The phase distributions in an isotropic pore network with the Hurst exponent H = 0:7 at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.8, and (d) 0.95. Blue, white, and yellow indicate, respectively, the throats lled with brine, the emptied throats, and those that are partially or completely blocked. The top row of the network is open. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.4 The phase distributions in a pore network with the Hurst exponent H = 0:7 and anisotropy x = y = 1=4, at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.8, and (d) 0.95. Colors indicate the same as in Fig. 5.3. The top row of the network is open. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.5 The phase distributions in a pore network with the Hurst exponent H = 0:7 and anistropy x = y = 4, at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.8, and (d) 0.95. Colors indicate the same as in Fig. 5.3. The top row of the network is open. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.6 The phase distributions in an isotropic pore network with the Hurst exponent H = 0:3 ( x = y = 1), at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.8, and (d) 0.95. Colors indicate the same as in Fig. 5.3. The top row of the network is open. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.7 The phase distributions in a pore network with the Hurst exponent H = 0:3 and anisotropy x = y = 1=4, at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.8, and (d) 0.95. Colors indicate the same as in Fig. 5.3. The top row of the network is open. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.8 The phase distributions in a pore network with the Hurst exponent H = 0:3 and anisotropy x = y = 4, at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.8, and (d) 0.95. Colors indicate the same as in Fig. 5.3. The top row of the network is open. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.9 The phase distributions in an uncorrelated PN at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.8, and (d) 0.95. Colors indicate the same as in Fig. 5.3. The top row of the network is open. . . . . . . . . . . . . . . . . . . . . . . . . . . 86 x 5.10 Dynamic evolution of the transversely-averaged brine saturations in six types of pore networks. Black curves show the average saturations; cyan, blue and pink shadings represent, respectively, the times at which vapor saturation is 0.2, 0.4, and 0.8, while red indicates the time at which vapor saturation is 0.95. The colored areas indicate the variations over multiple realizations. . . . . . . . . . . . . . . . . 88 5.11 Dynamic evolution of the transversely-averaged saturations in the uncorrelated pore network. Black curves show the average saturations; cyan, blue and pink shadings represent, respectively, the times at which vapor saturation is 0.2, 0.4, and 0.8, while red indicates the time at which vapor saturation is 0.95. The colored areas indicated the variations over multiple realizations. . . . . . . . . . . . . . . . . . . 89 5.12 Dynamic evolution of the average drying rates (black curves), averaged over 10 re- alizations (shadings). The colored areas indicate variations over multiple realizations. 90 5.13 Dynamic evolution of the average length of the drying front. Shaded areas show the extent of the variations over ten realizations. . . . . . . . . . . . . . . . . . . . 90 5.14 Size distribution of the partially or completely plugged throats before and after evaporation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.15 Size distribution of the partially or completely plugged throats before and after evaporation in the uncorrelated pore network. . . . . . . . . . . . . . . . . . . . . . 92 6.1 A portion of a sample well log used in this study. The data belong to the TBK- 1. Top: typical neutron porosity (NPHI), gamma-ray (GR) and resistivity logs. Bottom: wavelet transform coherence of gamma ray and resistivity, and gamma ray and neutron porosity. The black curve indicates 5% signicance level. The coherence ranges from dark red (nearly 1) to dark blue (nearly 0). The red colored regions inside the thick dark counters show the aected data (porosity or resistivity) by clay minerals (indicated by gamma ray) at various scales (periods). . . . . . . . 96 6.2 The formation factors as a function of eective porosity ( t ). Solid lines are the predictions based on percolation-eective-medium approximation with vari- ous crossover porosity x . The data are from three wells used, two from Maroon (MN296 and MN297) and one from Ahvaz (AZ273) oil elds . . . . . . . . . . . . . 98 6.3 Logarithmic plot of M(q = 2;s) versus scale s for four porosity logs from four dierent reservoirs. The slope of the line is the Hurst exponent H . . . . . . . . . . 102 6.4 Logarithmic plot of M(q = 2;s) versus scale s for four resistivity logs from four dierent reservoirs. The slope of the line is the Hurst exponent H R . No data from the logs were removed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.5 Logarithmic plot of M(q = 2;s) versus scale s for four resistivity logs of clean for- mations (those in which the data indicating high correlations with the clay contents removed) from four dierent reservoirs. The slope of the line is the Hurst exponent H R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 xi Abstract Salt transport and precipitation in porous media constitute a set of complex and fascinating phenomena that are of fundamental interest to several important problems, ranging from storage of CO 2 in geological formations, to soil fertility, and protection of pavements and roads, as well as historical monuments. The phenomena occur at the pore scale and are greatly in uenced by the heterogeneity of the pore space morphology. In this thesis salt precipitation is studied at two scales: molecular and pore scales. First, we use molecular dynamic simulation to study evaporation and salt precipitation in a single pore. Salt crystals precipitate on the pores' surface, modify the pore space morphology, and reduce its ow and transport properties. Despite numerous eorts to understand the mechanisms of nucleation of salt crystals in conned media, the eect of the rock's chemistry, which is mostly clay, on the growth, distribution and properties of the crystals is not well understood. We report the results of extensive molecular dynamics simulation of nucleation and growth of NaCl crystals in a clay pore using molecular models of two types of clay minerals, the Na-montmorillonite and kaolinite. Clear evidence is presented for the nucleation of the salt crystals that indicates that the molecular structure of clay minerals aects their spatial distribution, although the nucleation mechanism is the same in both types of clays. We then develop a pore-network model of ow of brine, and salt transport and precipitation in porous medis. To carry our ecient simulations with large pore networks, one must develop new computational algorithm. The computation time and required memory are two limiting factors that hinder the scalability of the computations to very large network sizes. We present a dual approach, based on the use of a combination of the central processing units (CPUs) and xii graphics processing units (GPUs), to simulation of ow, transport, and similar problems using pore networks. A mixed-precision algorithm, together with the conjugate-gradient method is implemented on a single GPU solver. The eciency of the method is tested with a variety of cases, including pore- and random-resistor network models in which the conductances are long- range correlated, and also contain percolation disorder. Both isotropic and anisotropic networks are considered. To put the method to a stringent test, the long-range correlations are generated by a fractional Brownian motion (FBM), which we generate by a message-passing interface method. For all the cases studied an overall speed-up factor of about one order of magnitude or better is obtained, which increases with the size of the network that one can use. Even the critical slow- down in networks near the percolation threshold does not decrease the speed-up signicantly. We also obtain approximate but accurate bounds for the permeability anisotropyK x =K y for stratied porous media. Before developing a pore-network model for salt transport and precipitation, one must also develop an model for evaporation in porous media. At the pore-scale, a threshold capillary pressure P =a t is needed for the invading (vapor) phase at the throat entrance in order to expel the wetting uid from a throat/pore in a porous medium. The threshold pressure depends on the interfacial tension between the wetting and non-wetting uids and the entrance wetting radius r t of the throat/pore. Isotropic and homogeneous porous media consist of throats/pores of a sole average radius and, hence, the threshold capillary pressure usually outweighs the ow- induced viscous pressure in a pore. Thus, a minimal variation in r t changes the uid phase distribution and transport properties of the porous medium during drainage invasion/drainage process. Natural porous media, such as geological formations, are stratied, consisting of more or less parallel layers, with each layer having a dierent morphological properties. The extra variations in the pore space strongly aects ow and transport in porous media. Thus, we develop a three-dimensional pore-network model of evaporation in layered porous media, in which we use a fractional Brownian motion (FBM) to include the impact of the pore-size disorder, heterogeneity, xiii and correlations. Our simulations indicate that diusion of the evaporating phase, drying patterns of liquid, and saturation proles are all aected strongly by the correlations and disorder in pore space, and that, while the increase in the stochasticity of the pore sizes stabilizes the drying front, the increase in the anisotropy parameter results in a more disordered drying front. We then present a pore-network model to study salt transport and precipitation in pores media. Vapor diusion, capillary eect at the brine-vapor interface, ow of brine, and transport of salt and its precipitation in the pores that plug the pores partially or completely are all accounted for. The drying process is modeled by the invasion percolation, while transport of salt in brine is accounted for by the convective-diusion equation. We demonstrate that the drying patterns, the clustering and connectivity of the pore throats in which salt precipitation occurs, the saturation distribution, and the drying rate are all strongly dependent upon the pore-size distribution, the correlations among the pore sizes, and the anisotropy of the pore space caused by stratication that most natural porous media contain. In particular, if the strata are more or less parallel to the direction of injection of the gas that dries out the pore space (air, for example) and/or causes salt precipitation (CO 2 , for example), the drying rate increases signicantly. Moreover, salt tends to precipitate in clusters of neighboring pores that are parallel to the open surface of the porous medium. In order to demonstrate the existence of the heterogeneity, correlations and long-range cor- relations in large-scale porous medium, we also study the well logs of hydrocarbon reservoirs. Archie's law expresses a relation between the formation factorF of porous media and their poros- ity , F/ m , where m is the Archie or the cementation exponent. Despite widespread use of Archie's law, the value of m and whether it is universal and independent of the type of reservoir have remained controversial. We analyze various porosity and resistivity logs along 36 wells in six Iranian oil and gas reservoirs using wavelet transform coherence and multifractal detrended uctuation analysis. m is estimated for two sets of data: one set contains the resistivity data that include those segments of the well that contain signicant clay content, and one without. xiv The analysis indicates that the well logs are multifractal, and that due to the multifractality the exponent m is non-universal. Thus, analysis of the resistivity of laboratory or outcrop samples that are not multifractal yields estimates of m that are not applicable to well logs in oil or gas reservoirs. xv Chapter 1 Introduction 1.1 Background and Motivation Flow, transport, and reaction in porous media occur in a wide variety of contexts and have been studied for decades. Examples of porous media and materials that are of fundamental scientic and industrial importance vary from membranes, catalysts, adsorbents, skin and other biological tissues, wood, printing paper, thin dielectric lms, polymeric composites, all of which are relatively \small" length scales and are referred to as the laboratory scale, to large-scale porous formations that include oil, natural gas, and geothermal reservoirs and groundwater aquifers that are at the largest length scales, and are usually called eld-scale porous media. Such porous media are heterogeneous at many distinct and disparate length scales, ranging from pore to core to eld scales. The heterogeneities manifest themselves in a variety of manners. At the laboratory scale, a porous medium is characterized by its pore-size distribution, pore connectivity, and roughness of its pore surface. But, even at this scale, other complications may arise. For example, a sample of a porous medium taken from a eld-scale formation, contains complex variations in its porosity distribution, a feature that is not expected to occur at such scale, but has been shown to occur (Sahimi, 2011; Dashtian et al., 2015). 1 At the eld scale, the heterogeneities are in the form of the spatial distributions of the perme- ability and porosity, as well as discrete fractures or an interconnected network of fractures, as well as faults, severe channeling, etc. In addition, chemical composition of eld-scale porous media is highly complex and heterogeneous, resulting from chemical processes that have occurred over geological time scales, spanning hundreds of millions of years. The complex chemical composition aects multiphase ow, reaction and adsorption in eld-scale porous media. Any meaningful study of a given phenomenon in porous media must, as the rst step, charac- terize the morphology of the media and develop appropriate models for them. Due to signicant advancements in instrumentation, as well as the dramatic increase in the computational power, characterization of laboratory-scale porous media has made great progress over the past two decades. The progress has been so advanced that one can now obtain three-dimensional images of porous media and carry out ow and transport simulations directly with the images. Such advances have opened up the door to the study of highly important, but little studied, problems in porous media, one of which will be studied in this thesis. On the other hand, characterization of eld-scale porous media is fraught with diculties and complexities. Lack of adequate data is one major problem, but analyzing such data when available is just as dicult, if not more. For example, well logs provide valuable information and insights into the structure of eld-scale porous media. But, the analysis of such logs is still a dicult problem, as the data are collected over hundreds of meters in depth, with a well passing through multiple layers of the porous formations that have widely varying chemical compositions. At the same time, such logs are often noisy and, therefore, separating noise from actual physical attributes of a porous formation, interpreting the data, and using them meaningfully in models of eld-scale porous media are still unsolved problems of interest. Finally, a highly important subject has not received the attention that it deserves. This is the problem of cross-correlations between the various types of logs, and the interpretation of such correlations. 2 1.2 Pore-Network Scale versus Field Scale This thesis studies two sets of important and inter-related problems at two distinct length scales. But, while one set of phenomena that we study is not understood even at the pore or laboratory scales, the second set of problems have been studied at such scales for decades, and a reasonable understanding of them has emerged. But, they must still be studied at eld scale, because not only they are not well-understood at such length scales, they are also important to a wide variety of problems at the largest scales. Thus, the problems that we study are at two dierent length scales, one at the laboratory or pore-network scale, and a second one at the eld scale. What follows is a brief description of each. 1.3 Pore-Network Scale: Salt Crystallization and Precipitation Salt crystallization in porous media or at their surface is a phenomenon of interest in rela- tion with several important applications, such as soil physics, underground storage of CO 2 , civil engineering and the protection of our cultural heritage, to name but a few. Therefore, there is an increasing importance to understanding the processes behind salinization. Salt deposition in porous media can be classied into three main dierent categories: drying, evaporation/wicking with the medium fully saturated and evaporation/wicking with an internal evaporation front. The three related but dierent situations can occur in any natural porous medium, depending on the available saline source. Although the three problems are quiet similar, we focus on the drying of salt water in porous media in which the volume of salt water is limited. The phenomenon of salinization causes several changes in the physical properties of porous medium, and in uences the ow and transport processes of the components of a mixture within the system. Given the negative impact that salt deposition has on various types of porous media, 3 this research purposes to study the process of drying-salinization in porous media and its impact on the static and dynamic properties of porous media. But, even though salinization is highly important to a wide variety of problems, ranging from irrigation, agriculture and drinking water from groundwater aquifers, to oil and gas reservoirs and sequestration of CO 2 in large-scale porous media, it is still not well understood even at the laboratory scale, and in fact it still lacks a comprehensive and detailed model at this scale. Thus, one goal of this thesis is to develop such understanding at the laboratory scale. For this purpose we will develop a pore-network simulator to model salt deposition in porous media during transport of water. The transport mechanisms that will be studied in drying porous media are vapor diusion in multicomponent systems, capillary ow in the liquid phase, ow in liquid lms, menisci moving in the throats and evaporation at the menisci. At the same time, evaporation of salt water leaves the salt crystals on the throats' wall, leading to reduction of the permeability and porosity and, hence, the uid velocity eld. In our approach, we will include all the aforementioned phenomena in the simulation. Pore network modeling is the standard method of studying scalar and vector transport pro- cesses in heterogeneous materials and media, such as uid ow in porous media, and conduction, deformations, and electric and dielectric breakdown in heterogeneous solids. In any pore network simulation of porous medium, there are two two challenging problems: (1) how to represent and generate the pore and throat network, and (2) how to solve the uid ow and transport equations. In order to assist such study of the phenomena in porous media, any modeling technique must be capable of controlling various parameters of the pore space, such as the particle size distribution and the extended correlations that typically exist in porous media. Thus, in addition to studying the salinization problem in porous media, we study two related problems, namely, developing an algorithm that accounts for the heterogeneity and correlations in the pore space, and developing an algorithm for using the computational power of graphics processing units (GPUs) to speed up the calculations with pore networks. 4 1.4 Field-Scale: Multifractal and Multiresolution Analysis of Data An important task in geophysics is to delineate as much information as possible about the structure of the interior of rock. Since direct exploration by penetrating rock is impractical, various types of data, such as seismological, electromagnetic, and gravitational measurements, collectively referred to as geophysical data, are used. Such data provide much insight into the structure and properties of rock and eld-scale porous media, such as oil, gas, and geothermal reservoirs and groundwater aquifers. Thus, a major task is making inferences from such data. The method of analysis depends, however, on the volume of the data. Many properties of eld-scale porous media are represented by well log data. Such data consist of depth proles of the media's important properties measured along wells, such as the spontaneous potential of production, sonic transient time (STT), natural gamma ray (GR) emission, resistivity, neutron porosity (NP), bulk density (BD), and photoelectric absorption factor, and are referred to as well logs. They provide information on the spatial heterogeneity of the properties of the eld-scale porous media, such as porosity, density, and the lithology at distinct length scales. The production potential and economic feasibility of investment in oil and natural gas reservoirs and their development are assessed based on the well logs. Moreover, any realistic model of such porous media is based in part on the existing well log data and, therefore, their correct analysis is critical to the development of the model. During sedimentation and evolution of geological strata, a collection of several stochastic vari- ables controls the evolution of the porosity and other important properties of porous formations. Such stochastic variability may give rise to extended correlations in the data, and make their anal- ysis a complex task. If the correlation functions are of power-law type, then important properties 5 of eld-scale porous media may follow fractal and multifractal distributions. Thus, understanding the degree of fractality in well log data and their eect on each other are very important. In this part of our thesis we use multifractal methods to explore multifractality of well logsin eld-scale porous media. The correlations and cross-correlations in the logs are also very impor- tant to characterization of porous media, as is the interpretation of the data. We are developing multiscale multiresolution methods, such as wavelet and curvelet transformations, to analyze geo- physical data, such as well logs and seismic data. The last two decades have seen tremendous activity in the development of new mathematical and computational tools based on multiscale ideas. Today, multiscale/multiresolution concepts permeate many elds of contemporary science and technology. We will develop the application of such concepts in order to extract more infor- mation from geological data. 1.5 Salt Precipitation in Porous Media Salt deposition and salinization over years was seen as a common sight causing problems in various elds of interest. Though the problem is quiet general, in this study we focus on tow main problems. One is the salt precipitation in aquifers and underground reservoirs due to the injection of CO 2 , and the other is the salinization of soil due to free air ow. Salinization causes signicant negative impacts in various sectors. To name a few: Energy/Oil and Gas: problems related to the production of natural gas, gas or CO 2 injection or in geo-thermal wells; Agriculture: where fertile land degradation is seen to cause reduction in crop yield; Infrastructure: damage to pipelines, boreholes, roads, building, historical heritages and etc; 6 Environment/Water: turning fresh inltrating water and ground water sources saline, thereby spreading the saline water boundaries. 1.6 Salt Precipitation Due to CO 2 Injection Injection of large quantities of CO2 for underground geological storage in shallow or deep saline aquifers will lead to a strong water desaturation of the near well-bore because of drying mechanisms. Even the capillary trapped water can be removed by thermodynamical transfer of water vapor in the CO2 phase. The soil-physics-related issues concern the impact of salt on evaporation of water from soils and the associated impacts on many aspects of water management, global water cycle, soil and groundwater salinization, agriculture (e.g. U. et al. (2011), and references therein). As discussed in Peysson et al. (2014), the evaporation process occurring during the injection of CO2 in a saline aquifer can induce the dissolved salt precipitation in the porous structure. The result is changing of the pores' surface area by the crystallized salts which decreases the permeability. Permeability reduction has a serious risk of CO 2 injectivity decline. As discussed for example by (Scherer, 1990a; Shokri, 2014), the crystallization of salt in the pore space induces a pressure which may damage to the porous material. 1.7 Drying of Porous Media Drying and evaporation from porous materials are important problems for numerous tech- nological and engineering applications ranging from agriculture (Ben-Noah and Friedman, 2018; Aydin et al., 2008) to recovery and remediation of soil contaminated with hydrocarbon (Nadim et al., 2000; Shokri et al., 2010; Hosseini and Al, 2016; Soltanian et al., 2016; Lu et al., 2017), soil technology, recovery of volatile hydrocarbons from oil reservoirs, cosmetics, building restoration and material processing, such as the production of food, wood, paper, textiles, pharmaceuticals 7 and washing powders. As an example, drying form soil (Scanlon et al., 2012) is a vital process in many environmental phenomena such as soil salinization (Qazi et al., 2017). The rate of evapo- ration results in nutrient disturbances by aecting the availability, transport, and partitioning of nutrients as it controls the transfer of heat, air and water between the atmosphere and subsurface. Fluid ow and transport in heterogeneous media in general, and evaporation in particular, closely controlled by the grain and pore size distribution, pore connectivity, wettability and several other factors. In general, uid ow is characterized by heterogeneous patterns arising form spatial heterogeneity of pore size distribution in porous medium (Sahimi, 2011; Bernab e and Bruderer, 1998). Drying rate from porous medium depends on both uid and porous structure properties. Al- though, evaporation is very complex, and to some extend, unpredictable phenomena, it generally consists of three main stages: (i) the initial constant evaporation; (ii) the relatively fast inter- mediate falling rate, and (iii) the residual slow rate of evaporation. The rst stage is essentially independent of the morphology of porous media and is controlled by the external conditions. It appears as a very short period in the initial time steps. Evaporation in the intermediate stage is dictated by the ability of the porous medium to supply the evaporation zone with brine, which in turn depends strongly on the morphology of the pore space, and in particular the PSD, the correlations, and the pore connectivity, and lasts much longer than the initial stage (Shokri et al., 2008a,b; Lehmann and Or, 2009; Chauvet et al., 2009). Similar to other transport process in porous media, there exist a complex interplay between capillary, viscous, and gravitational forces (Or, 2008; Yiotis et al., 2015; Pegler et al., 2016), wettability eects (Bergstad and Shokri, 2016), and the underlying heterogeneous pore geometry during drying which leads to ramied, preferen- tial ow paths or uid distribution patterns. 8 Chapter 2 Nucleation of Salt Crystals in Clay Minerals: Molecular Dynamics Simulation 2.1 Introduction Nucleation and growth of salt crystals from supersaturated solutions by evaporation of water is a common phenomenon in many natural systems, both on the surface (Desarnaud et al., 2015, 2016) of Earth and in porous geological formations (Baumann et al., 2014a). For example, salt weathering damages porous rock and the soil prole. While geologic sequestration of CO 2 is a promising way of mitigating global warming, salt precipitation, both near and far from the injection wells, hampers ow and transport of CO 2 in the porous formations, as it leads to pressure build up and reduction in the permeability and porosity, and ultimately the capacity of the formations for holding CO 2 . Deep insight into crystal nucleation and growth in the pores, and more generally in conned media, is critical to controlling such processes and avoiding salt weathering. The problem is more severe when the geological formation contains signicant amount of clay minerals. Surface chemistry at the clay-brine interface plays a fundamental role in a wide range of phenomena that occur in geological formations, including solute transport and salt precipitation. Clay minerals consist of parallel tetrahedral and octahedral structures, and are further divided 9 into subgroups that are distinguished based on the type of the isomorphic cation substitution that they contain. Kaolinite, illite, vermiculite, and smectite are the four major clay groups (Freedman, 2015). Clay particles are phyllosilicates that consist of numerous combinations of bonded tetrahedral and octahedral sheets stacked together. Silicon-oxygen tetrahedrons, SiO 4 4 are linked together by sharing the basal oxygen atoms and form a tetrahedral sheet, referred to as the T layer. The layers do not exist by themselves. Instead, they are always bonded to octahedral sheets, the O layers that are formed from aluminum-oxygen octahedrons, AlO 5 4 , with shared apical and basal oxygen. Hydrogen bonding bridges the T and/or O sheets. Thus, one may classify clay minerals into subgroups that are distinguished based on the way the T and O sheets are stacked together. The compositions of the T and O layers dene the overall charge of the clay layers, which can be zero or negative. In the latter case the charge is balanced by introducing such cations as K + and Na + into the interlayer space (clay pore) between the clay particles. The 1:1 clay particles, such as kaolinite, are formed by binding of one T and one O sheet. Their structure does not contain additional cations and, thus, they are electrically neutral. The 2:1 clay particles are formed by bonding of two T and one O sheets. Due to isomorphic substitu- tion in their structure, they have a net negative charge on their surfaces. Such clays, including the well-known montmorillonites (MMT), also have high surface area and adjustable spacing between their layers that contribute to their swelling and/or shrinkage and reactions. The dierences in the chemical compositions and structures of clay minerals determine they diverse behavior under similar conditions. For example, clay-water interaction is completely dierent for kaolinite and the MMT. Whereas kaolinite exhibits little or no swelling on hydration, the MMTs swell considerably when water is added to their pore space. They also dier in the degree of isomorphous replace- ments in their structure, as well as in the amount and nature of their associated exchangeable cations. 10 Clay minerals, water and sodium chloride, the most abundant salt on Earth, co-exist every- where. Although recent theoretical and experimental studies (DeCarlo and Shokri, 2014; Shokri et al., 2015) have provided valuable insights into the mechanisms of nucleation and growth of NaCl crystals out of brine in the presence of clay minerals, many questions remain unresolved. Given the signicant role that salt nucleation on clay particles plays in nature, gaining deep understanding of the phenomenon is of great importance. Experimental studies always face di- culties associated with the small size of the clay particles and their mixing with other minerals. On the other hand, fundamental, molecular-scale understanding of clays and their properties can be obtained by molecular dynamics (MD) simulation, if accurate force elds (FFs) are available. The utility of MD simulation is manifested by the fact that water near the surface - the bound water - and its properties are completely dierent from those under bulk conditions, leading to the so-called dual water model. Macroscopic theories cannot dierentiate between the two types of water and, hence, are incapable of describing accurately the process that occur there. It is, of course, imperative to include the interactions of all the atoms and molecules in any model of clays, if accurate MD simulation is to be carried out (Cygan, 2001). This requires very accurate FFs. The clay force eld, CLAYFF (Cygan et al., 2004), which was developed based on the metal- oxygen interactions connected with hydrated phases, is one such accurate FF for MD simulation that we utilize in the present chapter. In this chapter we report on extensive MD simulation of evaporation of brine in clay minerals, and nucleation of salt crystals. We demonstrate that Na + and Cl ions prefer to be absorbed, and then precipitate, on specic sites of clay surfaces, hence aecting the spatial ordering of the nucleated crystals. Macroscopic crystals may form from the initial nuclei with only a few atoms. The growth of the nuclei aects the solution's thermodynamics. Although capturing such mechanisms by experiment is dicult, MD studies provide clear indications for them. In addition, MD simulations allow one to study pure systems, whereas experimental systems may contain impurities. 11 Z X Y (a) (b) (c) (d) Figure 2.1: (a) and (c) show the kaolinite molecular structure, while (b) and (d) depict that of Na-MMT clays. Water molecules in (a) and (b) are shown by sticks. Colors represent oxygen (red), hydrogen (white), Si (yellow), Mg [green, only in (b)], Al [pink in (a) and (b)] and [green in (c) and (d)], Na (purple), and Cl (green; very small). 2.2 Materials and Method As the model clay particles we use kaolinite and the Na-MMTs, for molecular modeling of which we have extensive experience (Ghassemzadeh and Sahimi, 2004b; Kim et al., 2005, 2007). The Na-MMT model is based on a triclinic pyrophillite structure, which has been used widely in clay-related studies. The chemical structure of the MMT is characterized by random isomorphic substitution of Al by Mg atoms. In the present study the chemical composition of the MMT unit cell was set to be (Kim et al., 2007) Na 3 [Si 31 Al 1 ][Al 14 Mg 2 ]O 80 (OH) 16 . The size of the unit cell was 20:92 18:12 12:50 A 3 . The unit cell was replicated six times along the x and y axes and twice along the z axis (see Figure 2.1), in order to obtain a super cell with a size of 125:52 108:72 25:00 A 3 . The same approach was used to construct a model kaolinite (Freedman, 2015; Molina-Montes et al., 2008), for which the chemical formula is Al 2 Si 2 O 5 (OH) 4 and the size of its super cell was 123:68 107:30 14:78 A 3 . The structures of 12 the MMT and kaolinite are shown in Figure 2.1. Note the dierences in the thickness of two types of clay. Clay minerals swell and/or shrink when they come into contact with water, either pure or saline, which may in uence the nucleation and growth of salt crystals. Here, we focus on the eect of the clay types on the nucleation of salt crystals. As mentioned earlier, CLAYFF was used for generating the molecular structures of both clays. The FF contains a single bond-stretching parameter for the structural O-H groups. As a smectite particle, the MMTs hold a net negative structural charge, due to isomorphic substitution by elements of lower valences in the tetrahedral and octahedral sheets (Ou et al., 2014). The charge deciency caused by the isomorphic substitution of Mg is compensated by inserting four Na + ions in the structure. According to CLAYFF the total energy E of the material is given by (Cygan, 2001; Ou et al., 2014), E = X bonds k 1 (r ij r 0 ) 2 + X angle k 2 ( ijk 0 ) 2 + X ij ij " ij r ij 12 2 ij r ij 6 # + e 2 4" 0 X i;j q i q j r ij : (2.1) Here,k 1 andk 2 are two force constants, ijk is the angle between bondsij andjk,e is the electron's charge, q i is the partial charge of atom i, " 0 is the dielectric permittivity of vacuum, and r 0 and 0 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, ij = ( i + j )=2 and ij = p i j were used for the pairs ij. Sodium chloride was used to represent salt, with its molecular structure and parameters given by CLAYFF. After some preliminary simulations, the cuto distance for the LJ interactions was set at 12 A. Water was represented by the SPC/E model. The electrostatic interactions were computed using the particle mesh Ewald summation method. The simulations were carried out in the (NPT ) ensemble at 1 atmosphere and 298 K. Pressure and temperature were held constant using, respectively, the Andersen (Andersen, 1980) and Nos e- Hoover algorithms (Nos, 1984; Hoover, 1985). The masslike damping parameter for the Andersen 13 barostat was 100 amu. For the Nos e-Hoover thermostat the parameter that we used to scale the ctitious mass was 1. The time step was always 0.5 fs. In the Supplementary Information all the parameters of the FF are given. We consider a slit pore that consists of two parallel kaolinite or MMT sheets. In both cases 16128 water molecules were used in the simulation cell. The number of NaCl molecules were 1576 in the Na-MMT to 1584 in koalininte, which is equivalent to 30.9 NaCl weight percent in the solution, representing the saturation limit of salt in water at 298 K. This concentration is slightly above the crystallization limit for NaCl at 298 K and 1 atm under which the previous experiments had been carried out. The salt concentration was set at a relatively high value because it facilitates nucleation of salt crystals during the MD simulations. Note that depending on the precise values of the molecular parameters of water molecules and the interaction parameters, the crystallization condition slightly varies, as reported in some MD simulation studies (Aragones et al., 2012; Vega and Abascal, 2011). Alejandre and Hansen (2007) and Mendoza and Alejandre (2013) investigated the sensitivity of the solubility limit of NaCl in water to the NaCl model parameters. The MD simulations were carried out using the LAMMPS package (Plimpton, 1995). The initial spatial distribution of the water molecules and the salt ions was obtained by minimizing the total energy of the system that consisted of the water molecules, the salt ions and the clay pore and plates. The system was then equillibrated for 2 ns, followed by a production run of up to 70 ns. Water molecules were removed from the simulation cell at a very slow rate (Mucha and Jungwirth, 2003) to simulate evaporation, as the nucleation and crystallization occur at higher salt concentrations than the initial value. Each time a water molecule was removed, the system was allowed to reach equilibrium. Our simulations indicated that, so long as the evaporation rate is small, which is the case in our simulations, it has a minor eect on the nucleation mechanism and growth of salt crystals. We used three slightly dierent rate of evaporation, namely, 168, 161 and 155 water molecules/ns, which are very small. Of course, if evaporation is rapid, then non-equilibrium eects become important, but this limit is beyond the scope of our present work. 14 Kaolinite (a) (b) (c) (d) (a) (b) (c) (d) Na-MMT Figure 2.2: Top view of the positions of Na + (pink) and Cl (yellow) within the interlayer space of kaolinite and Na-MMT. Time increases from (a) to (d). (a) Early positions (after 0.1 ns) of the ions after the system reached equilibrium. (b) The initial spatial ordering of the ions after 7 ns for kaolinite and 6 ns for Na-MMT. (c) and (d) show the nucleation steps and spatial ordering after 30 and 70 ns for, respectively, kaolinite and Na-MMT. Colors show Na + (pink) and Cl (yellow). 2.3 Results and Discussion Figure 2.2 shows the top view of the nucleation and growth of NaCl crystals in the interlayer space of kaolinite and Na-MMT, as they evolve with time. For clarity the water molecules and clay sheets are not shown. In the case of Na-MMT, there is no distinct order among the ions at early times, less than 3 ns. After about 5 ns, the local ion density around some particular sites increases, but there is still no distinguishable order amongst Na + and Cl . In addition, by this time some water molecules have left the interlayer region. After about 7 ns, some initial small nuclei of crystals have formed, but they are still not stable. After about 10 ns larger crystals have formed and the central part of the system appears to be stable. At the same time, the outer surface of the nuclei indicates some sort of disorder, indicating that the structure there is not stable yet, but the inner part includes parallel arrangements of Cl and Na + . The nal stage is marked by stable crystals lling some portions of the interlayer space of the clay, on the one hand, and some empty regions with essentially zero ion density, on the other hand. Note that since the 15 Na-MMT has extra Na + in its structure, they are attached to the clay or the salt crystal lattice. Qualitatively, the same sort of phenomena is seen in the kaolinite system, except that the time scales and the spatial distributions of the crystals are dierent; we will return to the kaolinite case shortly. 5 10 15 0 0.1 0.2 0.3 0.4 0.5 r (nm) g(r) (a) 5 10 15 0 0.1 0.2 0.3 0.4 0.5 r (nm) g(r) (b) 5 10 15 0 0.1 0.2 0.3 0.4 0.5 r (nm) g(r) (c) 5 10 15 0 0.1 0.2 0.3 0.4 0.5 r (nm) g(r) (d) 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 r (nm) g(r) (e) 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 r (nm) g(r) (f) 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 r (nm) g(r) (g) 5 10 15 0 0.1 0.2 0.3 0.4 0.5 0.6 r (nm) g(r) (h) Figure 2.3: Evolution of the radial distribution function (RDF), indicating nucleation and spatial ordering of the ions in the interlayer space of kaolinite (a)-(d) and Na-MMT (e)-(h). (a) and (e) show the early-stage (< 1 ns) RDF of the ions. Our simulations indicate that crystallization in both the MMT and kaolinite follow quali- tatively the well-known two step mechanism (An et al., 2016; Lupi et al., 2016; Vekilov, 2005; Chakraborty and Patey, 2013). First, uctuations in the concentration form several regions with local ion concentrations higher than the rest, followed by a second step that includes nucleation and development of spatial ordering of the ions. The salt concentration uctuations in the inter- layer space lead to spatial inhomogeneity of brine, with nucleation of ions expected to be in the regions with higher concentrations. The regions that are formed in the rst step are unstable, so that only a few of them will eventually produce stable crystals. 16 One may be tempted to describe the crystal nucleation based on the classical nucleation theory (CNT) that consist of a single step. The CNT has, however, been challenged in recent experimental (Billinge, 2009a; Lutsko and Nicolis, 2006; Chung et al., 2009; Zhang and Liu, 2007) and simulation (Chakraborty and Patey, 2013) studies. According to the CNT crystallization involves a single free energy barrier, and the initial state of nucleation is a simplied version of the nal, stable crystal structure. Recent studies have shown, however, that the phase diagrams of systems experiencing nucleation and crystallization involve a phase transition region between two stable disordered states with distinct concentration densities, implying that before formation of stable crystal, the system transitions to an intermediate state in which the ions form disordered aggregates (Lupi et al., 2016). Such observations are usually described based on the Ostwald rule (Ostwald, 1897; Threlfall, 2003), according to which the free energy of the initial emerging nucleation in a solution is closer to that of the solution, rather than that of the salt crystals, and that the initial nuclei eventually represent the nal stable crystals after going through the stages. Chakraborty and Patey (2013) used large-scale MD simulation and showed that nucleation and crystallization of NaCl solution also follows the two-step mechanism. As Figure 2.2 indicates, our simulation too indicates that, initially, small size aggregates or clusters of ions form, and their number increases with time. Initially, the solution tends to form several small nuclei because (ten Wolde et al., 1995) it is entropically more favorable for such systems as the NaCl solution plus the clay minerals. Then, the small clusters combine gradually and form larger aggregates that are referred to as the critical nuclei. The Gibbs free energy of a LJ system approaches its maximum when the the rst critical nucleus appears in the system (ten Wolde et al., 1995). For kaolinite, the birth time of the crystal nucleation is longer than the MMT clay. The birth time and location of the rst nucleation site have important eect on the subsequent crystal growth (Zimmermann et al., 2015). Moreover, as Figure 2.2 indicates, the structure of the two systems are dierent. One obvious reason is the extra Na + in the MMT system that increase 17 the chances of Cl and Na + interactions, on the one hand, and can act as seeds to initialize nucleation, on the other hand. In fact, seeding the solution with a crystal cluster has been used in several MD simulations of nucleation and crystallization (Bai and Li, 2006; Knott et al., 2012; Pereyra et al., 2011; Sanz et al., 2013; Espinosa et al., 2016), including those by Espinosa et al. (2015) and Zimmermann et al. (2015) for nucleation and crystallization of NaCl solution. In Figure 2.3 we show the dynamic evolution of the radial distribution function (RDF) g(r). They represent averaged distributions, with the average taken over all the ions in the interlayer space. The initial RDFs for kaolinite and Na-MMT are, of course, similar, because there is no cluster in the two systems and the salt concentrations in the two are very close. As the water leaves the system and nucleation proceeds, the RDF peaks become more evident. This is due to the fact that the range of ion-ion correlation increases as salt crystals grow, and a larger lattice of Cl and Na + is formed. For example, after 30 ns sharper peaks are evident that correspond to the small patches of ordered ions. At still longer times the long-ranged structural peaks become even sharper. Similar RDFs of NaCl systems, computed by MD simuation, have been reported elsewhere (Galamba et al., 2005; Wang et al., 2014). The dierence in the shape of the nal RDFs for kaolinite and Na-MMT arises from the distinct distribution of the nucleation centers in the two systems. The spatial arrangement of the ions and water in the interlayer space is an important factor in nucleation and growth of the salt crystals. The ions and water molecules bind to the surface of the clay and form two major sorption sites in the interlayer space between the two clay sheets. The area (site) close to the surface of clay is called the inner sphere (at a distance of 2.5-3.0 A from the surface), while the one farther away from the clay surface is called the outer sphere (at distances of about 2.0-2.5 A beyond those of the inner spheres). In the former the electrostatic forces hold the ions and water molecules attached to the clay surface. The water molecules surround the ions in the interlayer space and form hydration shells. This prevents direct contact of ions in the outer sphere with the clay. 18 2 4 6 8 10 12 14 16 0 0.01 0.02 0.03 0.04 0.05 0.06 Z (A) Density (A −3 ) Kaolinite Cl − Na + O Cl − + Na + 2 4 6 8 10 12 14 16 18 0 0.02 0.04 0.06 0.08 Z (A) Density (A −3 ) Na−MMT Cl − Na + O Cl − + Na + Figure 2.4: Initial density proles of the salt ions and oxygen in water along thez axis (perpendic- ular to the clay planes) for salt concentration of 30.9% weight. z < 6 represents locations inside the lower clay sheet. As Figure 2.1 indicates, the interlayer space of kaolinite contains two dierent surfaces. The Na + ions migrate close to the AlO 6 octahedral surface and are paired with Cl , leading to adsorption of Na + on the octahedral basal surface. The number of such ions are, however, small when compared with the total number of ions in the interlayer space. In addition, Na + forms the outer sphere adsorption complexes on the tetrahedral SiO 4 surface. As discussed by Vasconcelos et al. (2007), the mobility of Na + parallel to the clay surface is higher than the other ions, such as Cs 2+ or Pb 2+ (Wang et al., 2014). In addition, octahedral surface (AlO 6 ) attracts Cl ions and forms inner sphere complexes. This is clear in the narrow peak of the Cl density prole, shown in Figure 2.4(a). That side of kaolinite is slightly positively charged and, hence, attracts Cl ions. However, the opposite tetrahedral SiO 4 surface has a small negative charge, because Si-bridging oxygen atoms relax outward slightly with respect to the Si atoms. Therefore, the surface repels Cl ions. 19 On the other hand, as Figure 2.1 indicates, the interlayer space of the Na-MMT clay includes two surfaces, but with identical structures. Each layer is formed by two tatrahedral (T) silicate sheets surrounding one octahedral (O) aluminum sheet, giving rise to the aforementioned T-O-T structure. Therefore, both surfaces in the interlayer space are of the tetrahedral (SiO 4 ) type. As Figure 2.4(b) shows, the density proles of the species in the interlayer are symmetrical. As discussed earlier, the SiO 4 surface repels the Cl ion. Hence, they form peaks at the midplane of the interlayer. Similar to the kaolinite tetrahedral surface, Na + ions form the outer sphere adsorption complexes near each clay particle surface that includes two maxima. It is due to such diverse arrangements and adsorption of ions in kaolinite and the Na-MMT interlayer space that distinct nucleation and crystallization patterns emerge in the two types of clay particles. Indeed, as discussed earlier, the initial step of nucleation is the emergence of high ion concentration region. The interlayer space of kaolinite contains one region in which the total concentration of Na + and Cl exceeds signicantly that of the bulk, whereas the Na-MMT has two of such regions with highest ion concentration. The density proles of Na + and Cl together, shown in Figure 2.4, clearly indicate this. Time (ns) MSD (A) Kaolinite Na−MMT 0 10 20 30 40 50 60 70 0 1 2 3 4 5 6 7 8 9 Figure 2.5: Evolution with time of the mean-square displacements (MSD) of the ions in the interlayer space. 20 Another way of looking at the evolution of the system is through the time-dependence of the mean-square displacements (MSDs) of the ions, as nucleation of the crystals is approached. They were calculated and averaged over all the ions with respect to a reference position. In Figure 2.5 we show the resulting MSDs for both clay types. In both cases the MSDs increase over time, albeit with dierent rates for two two types of clay, but they approach saturation, hence indicating that the ions eventually stop moving since nucleation has commenced and salt crystals have begun to form. The rate of approach to the saturation limit is higher in the Na-MMT case, indicating that the spatial ordering of the nucleation in kaolinite under the conditions that we study is not as stable as that in the Na-MMT. As discussed earlier, the competition for the adsorption sites in kaolinite and the tendency of ions to form nucleation centers are responsible for the phenomenon. 2.4 Summary Summarizing, we carried out extensive MD simulation to study nucleation and growth of salt crystals in the nanopore space of clay minerals. Two types of clay minerals with distinct structures were studied, namely, kaolinite and Na-montmorillonite. We provided clear evidence that crystal nucleation occurs in the interlayer space of the clay particles that is initially lled by oversaturated brine, and that it follows a two-step mechanism. The rate of nucleation in the Na-MMT particles is higher than in kaolinite, which is due to isomorphic substitution in the structures of the Na-MMT clay particles pointed out earlier. 21 Chapter 3 Ecient Simulation of Fluid Flow and Transport in Heterogeneous Media Using Graphics Processing Units (GPUs) 3.1 Introduction Fluid ow and transport in heterogeneous media is an important problem (Torquato, 2002), in view of its relevance to a wide variety of phenomena in natural and industrial processes. A partial list of such phenomena include ow through porous media (Blunt, 2017; Sahimi, 2011), conduction and hopping transport in composite solids (Shklovskii and Efros, 1984), mechanical properties of disordered materials (Sahimi, 2011), fracture of heterogeneous solids (Chakrabarti and Benguigui, 1997), and many more. For over four decades the standard model of a heterogeneous medium has been a random resistor or pore network (RRN and PN, respectively) in the case of ow and scalar transport, and a network of springs or beams for studying such vector transport processes as deformation and fracture propagation. Often, disorder is introduced in the model by a percolation process (Sahimi, 2011) by which a fraction of the bonds of the network - pores, resistors, or springs - do not allow uid ow, or scalar or vector transport to occur, while the rest of the bonds represent the uid ow or transport paths. Even with such a relatively simple form of disorder a wide variety 22 of interesting phenomena occur, a detailed understanding of which requires simulating very large networks and, thus, solving very large set of ow and transport equations. Numerous approaches, ranging from continuum to networks models, integrate (Chu et al., 2012) the complexity of heterogeneity of porous media with uid ow. For example, the PN models discretize the void space of a porous structure into a network consisting of pore bodies (nodes or sites) connected via pore throats (bonds) - hereafter referred to as, respectively, pores and throats - and have been widely used to investigate a wide variety of phenomena related to multiphase ow (Ovaysi and Piri, 2010; Oren et al., 1998; Piri and Blunt, 2005; Blunt et al., 2002). The applications includes drying of porous media (Borgman et al., 2017; Yiotis et al., 2015), reactive transport (Mehmani and Tchelepi, 2017; Nogues et al., 2013; Li et al., 2006), dissolution (Dillard and Blunt, 2000), unstable miscible displacements and ngering phenomenon (Hekmatzadeh et al., 2016), grain boundary wetting (Ghanbarian et al., 2014) and gas ow in shaly formations. The computational limitations of the PN modeling restrict, however, its application to relatively small networks. Thus, attention has also been focused on developing ecient methods for simulating very large PNs or RRNs and solving the associated large sets of ow or transport equations. The eorts has led to the development of the transfer-matrix method (Derrida and Vannimenus, 1982; Zabolitzky et al., 1986; Normand and Herrmann, 1990), and coarsening the disordered PNs or RRNs based on wavelet transformations and then solving a signicantly-reduced number of equations (Mehrabi et al., 1997). With the exception of the wavelet coarsening method, which is ecient both near (Sahimi and Tajer, 2005; Pazhoohesh et al., 2006) and away from the percolation threshold, all the aforementioned methods are ecient only if the PN or RRN is near the percolation threshold. Aghaei and Piri (2015) developed a computational strategy for simulating multiphase ows in porous media that is capable of simulating large PNs. Their approach is, however, a massively parallel scheme. 23 In practice, many heterogeneous media do not contain percolation-type disorder or, if they do, their disorder is not random but highly correlated. Moreover, with the advent of sophisticated experimental techniques, such as x-ray computed tomography (Jasti et al., 1993; Knackstedt et al., 2001), it is now possible to obtain detailed data for the morphology of heterogeneous media. Taking proper account of the correlations and utilizing the detailed data in the model entails employing a high-resolution PN or RRN with several million nodes, a still daunting task. Moreover, the problem is much more dicult when one must solve unsteady-state problems over a period of time, which entails solving a very large number of equations for thousands of time steps. In this thesis we present an ecient approach to PN or RRN modeling, which is based on graphic processing units (GPUs) that have opened up a new approach for high performance computing, particularly for the researchers who do not have access to massively-parallel machines with a large number of nodes, or to vector supercomputers. The GPUs were originally designed to analyze three-dimensional graphic images. To achieve extremely high performance with geometric data, the GPUs have been designed with simple and tiny processors. More modules are used for data processing, but not for the data cache, nor for the ow control. Hence, by design the GPUs are very dierent from typical computations with central-processing units (CPU), and are ideal for highly intensive and parallelized calculations, if the computational algorithm can take advantage of their special features. In facilitating nongraphical calculations on its graphic units, Nvidia Corporation released Compute Unied Device Architecture (CUDA) (NvidiaCorporation, 2015). Several groups have tested the performance of CUDA parallel computing in various calculations (Markall et al., 2010, 2013; Liu et al., 2008; Umtsev and Martinez, 2009; Olivares-Amaya et al., 2010; Zheng et al., 2013). 24 3.2 Pore and Random Resistor Networks The rst step in the PN/RRN modeling is to generate the networks. It was suggested (Sahimi, 2011; Molz et al., 1997) a long time ago that many natural porous media exhibit long-range correlations. Such correlations, which in uence the percolation (Sahimi and Mukhopadhyay, 1996) and ow and transport properties of porous media (Dashtian et al., 2015), follow the statistics of fractional Browning motion (FBM) (Mandelbrot and Ness, 1968), a nonstationary stochastic process. Thus, in addition to a random distribution of the pores' or resistors' conductances, we also use the FBM to generate correlated conductances, or throats' size, in the PN or the RRN. We use a square network in which each bond represents a throat with an eective radius r, or a resistor, selected from a statistical distribution. The bonds' length ` is assumed to be constant, although it poses no diculty to select it from a statistical distribution as well. The bonds' radii (or conductances) are selected from a FBM. An ecient method for generating a FBM array is through its power spectrum S(!), which for 2D systems is given by S(!) = a (! 2 x +! 2 y ) H+1 ; (3.1) wherea is a constant, while! x and! y are the Fourier components in thex andy directions. Here, H is the Hurst exponent such that forH > 0:5 (H < 0:5) one has positive (negative) correlations, while H = 0:5 represents the case in which the successive increments of the FBM are completely random. Due to stratication, many porous media are anisotropic. To introduce the anisotropy, we follow Ansari-Rad et al. (2012) and generalize the power spectrum to S(!) = b ( x ! 2 x + y ! 2 y ) H+1 : (3.2) 25 Here, x and y are the anisotropy parameters, and b is a constant. To generate anisotropy induced by layering with the layers being parallel to the x direction, we set x = y > 1. For large PNs that we utilize in the present study, however, generating a large FBM array is computationally expensive, because as the size of the PN increases, the computation and required memory increase exponentially, if sequential programming is used for generating the FBM. Thus, to make the running time manageable, we used message-passing interface (MPI) strategy and generated the FBM algorithm as a parallel scheme in the GPU solver (see below). Assuming steady-state and slow (laminar) ow of an incompressible uid, the volume ow rate in a bond ij that connects nodes i and j is given by, q ij =r 4 P ij =(8`)g ij P ij , where is the uid's viscosity, P ij = P j P i is the pressure drop between nodes i and j, and g ij the conductance of bond ij. Writing a mass balance for every node i, one obtains, P j2i q ij = 0, where the sum is over all the nodesj that are connected to nodei. Thus, substituting forq ij and writing the mass balance for every interior node of the network results in a set of linear equation of the following form, GP = b; (3.3) where G is the conductance matrix that depends only on the geometric properties of the network, P is the nodal pressure vector (or voltage vector in the RRN) to be calculated, and b is a vector related to the external boundary conditions. We impose a xed pressure gradient on the network in one direction, and used periodic boundary conditions in the second direction. Once the pressure (voltage) distribution in the PN is calculated, the eective permeability K (conductivity) of the network is computed based on the Darcy's law, K = Q S L P ; (3.4) whereQ is the total ow rate passing through the cross sectionS of the PN,L is the PN's length, and P is the macroscopic pressure drop imposed on the network. A similar procedure is used 26 to compute the eective conductivity of the RRN. The set of the equations must be solved by the conjugate-gradient (CG) method, which is the standard method in the computations with the CPU, but must be implemented in the GPU. The matrix-vector multiplication is the most dominant part in any CG iteration and, thus, it is this part that is accelerated by the GPU. To implement the solution of the pressure/voltage distribution in a PN/RRN on GPUs, one must use CUDA programming that makes it possible to implement the CG algorithm on a single GPU. CUDA programming, as well as the memory allocation and performance from the perspec- tive of our GPU solver kernels, which is by far the most time-consuming and resource-intensive computational step in the simulations, are described brie y in next section. 3.3 CPU and GPUs Memory allocation in, and performance of the implementation of the algorithm are by far the most time-consuming and resource-intensive computational step. On a GPU, parallel tasks, called threads, are scheduled and executed simultaneously in groups referred to as warps. One warp contains 32 threads that processed in parallel by one Compute Unied Device Architecture (CUDA) streaming multiprocessor (SM). The GPUs have many SMs that run in parallel to in- crease the parallelism. The threads are organized also into larger structures called blocks, which are themselves organized into grids (NvidiaCorporation, 2015). Figure 3.1 shows an example of the GPU architecture hierarchy with four SM. The CUDA memory contains various segments with dierent scopes and properties. All the threads have access to the same global memory. The shared memory is visible to all the threads within a block with the same lifetime as the block. Each thread has its own registers and local memory. The global memory is located in the GPU memory, which is the largest, but also the slow- est. A programmer can allocate free global memory. In order to access the GPU memory, 32, 64, or 128 byte memory transactions are needed (from the host to the device). These \transactions" 27 Figure 3.1: GPU architecture hierarchy with four streaming multiprocessors (SM). must have the same sizes, but because they are costly in terms of their computational perfor- mance, they must be aligned and minimized. In our implementation of the conjugate-gradient (CG) algorithm for solving the pressure/voltage distribution in a pore or random-resistor network (PN and RRN, respectively), we use one-dimensional (1D) arrays to store the elements in a vector. The local memory is allocated by the compiler for large structures that do not t into the register space, and are part of the global memory that can be used to avoid costly memory transfers. Despite local and global memories, shared memory is an on-chip memory and, hence, has much lower latency with higher bandwidth. The threads in a block can use shared memory to work together. For example, in the case that multiple threads in a block need to use the same data, shared memory is called upon to access the data from the global memory only once. The GPU that we use has a Kepler architecture, one of the most ecient high performance computing architecture in which the shared memory is 48 KB, and is organized into 32 banks. In 28 Figure 3.2: The CUDA memory model. our CUDA kernel - functions that are executed on the GPU - we store temporal variables in the shared memory. In Fig. 3.2 we show the concept of memory in a typical GPU. 3.4 Solving for the Pressure (Voltage) Distribution on CUDA As described earlier, the pressure/voltage distribution in a PN/RRN is calculated by the CG method that solves the set of equations, GP = b, where G is the conductance matrix, P is the (unknown) vector of the nodal pressures/voltages, and b is a vector that is zero everywhere, except on the boundaries of the network at which an external pressure or voltage is imposed. G is, of course, symmetric and positive denite. Once G and b, the initial guess for P, and the acceptable error for the CG iterations are specied, the algorithm solves for the pressure/voltage distribution in at most M iterations, where M is the number of equations to be solved. It is evident that the CG algorithm is composed mainly of matrix-vector multiplication, dot products, and scalar and vector addition and multiplication. We implemented the CG algorithm 29 in both C++ and CUDA. In CUDA the central processing unit (CPU) is a host for the GPU, often referred to as a device. Computer programs that are executed on the CPU (host) can manage memory on both the the CPU and the GPU. The CPU lunches the kernels, which are then executed in a parallel scheme by the GPU threads. The most computationally-intensive kernel in the CG (and any other iterative solver) is the matrix-vector multiplication (MVM) that consumes most of the execution time of the main loop of a solver. One important issue is the oating-point accuracy of the computations. In order to ensure con- sistent computations across platforms and to exchange oating-point data, the IEEE-754 denes basic and interchange formats to convert decimal oating point too 32 bit and 64 bit hexadeci- mal representations, along with their binary equivalents. The 32 and 64 bit basic binary oating point formats correspond to the C data types oat and double. The GPU devices with high com- pute capability support both single and double precision (SP and DP, respectively) oating point computation. On CUDA the compute capability is the \feature set," for both the hardware and software, of the device. The higher the compute capability, the more powerful the device. The GPU device that we use in this study supports both the SP and DP. There is, however, a signicant gap between the SP and DP computations with the GPUs, as the SP calculations run much faster than the DP computations. The acceleration of the computations depends on the ratio of arithmetic operations and data input/output, which is also important in the CG calculations and represents its main bottleneck. Thus, by using the SP data we reduce their \trac" by a factor of two between the GPU processor and its device memory. Then, to further improve the performance of the CG calculations, we used a mixed-precision method, implemented by an iterative renement algorithm. The main idea is using two types of iterative loops to obtain the true solution of the equations. First, by using the SP iteration we approach rapidly a rough solution within the inner loop. Next, iterations with the DP are used to converge to the nal solution within the outer loop allowed error. Since the iteration loop of the 30 Figure 3.3: Mixed-precision algorithm for the conjugate-gradient method that was used in the computations. CG algorithm is responsible for most of the computation time, we explain in detail the parallel implementation of this part on the GPU; see Fig. 3.3. First kernel: A 1D block of size 256 is used. For each part, the size of the network is set pro- portional to the total number of rows and blocksize, such that there exists a thread corresponding to each row. Second kernel: This is the sparse matrix-vector multiplication part of the computation, for which a thread is assigned to each row of matrix G, which is responsible for calculating the elements of Z, the vector that acts as temporary b during the iteration process. Consecutive rows of G have redundant access to P through the calculations associated with the three main diagonals of G. Thus, P is cached in the shared memory for improving the access pattern. After the calculation of each entry of Z in a thread, it is multiplied by its corresponding P element and the result is held in the shared memory. Then, by performing a reduction operation over each 31 block, the sums, which arise when one multiplies one row of a matrix by the column of another matrix or a vector, are calculated and stored in the global memory. Third Kernel: For this part a thread is assigned for the calculation of each element of U, R, and Z and, similar to the second kernel, partial result of dot products for R and Z are calculated and saved in the global memory for each block. Here, R is the residual vector for each iteration, while U represents the temporary approximate solution after each iteration. Merging the calculations of the vectors' entries into a single kernel leads to a more ecient performance by decreasing the number of times that the the global memory must be accessd. We optimized the number of threads and blocks based on the computational capability of the GPU device that we used. The number of threads per block must be a multiple of the warp size, which is 32 in all the current hardwares. We used 256 threads per block and 64 blocks in a xed 1D grid arrangement. A loop was placed inside the threads to calculate multiple outputs per thread, so that the rst iteration through all the threads is responsible for the rst 2 14 row calculations; the next iteration for the next 2 14 row calculations, and so on. After completion of the iteration loop, the results are saved in the shared memory. By perform- ing a reduction operation in each block, the sum is calculated and saved in the global memory. As there are only 2 6 blocks, there will be 2 6 partial sums in the global memory, to be added together to calculate the nal solution. By copying the sums' values to the main memory and performing the nal sum on the CPU, a more ecient performance is achieved. Thus, a multipro- cess algorithm was designed in an attempt to use the global memory throughput more eciently. Since most of the summation operations are performed in registers or the shared memory, there is little eect on the global memory bandwidth, and the ecacy of the mixed-precision design is not aected. Task latency - the elapsed time between initiation and completion of a task - and throughput - the rate at which the system can process tasks - are two fundamental measures of processor performance. Improving throughput or reducing latency results in a better speed-up. The CPUs 32 perform better for latency-oriented tasks, whereas the GPUs are designed for throughput-oriented computing systems. There are some specic measures that are used to decide whether a certain application (or task) is suitable for the GPU or the CPU implementation. The memory footprint of a task, the amount of main memory that a program uses or references while running, is the primary measure. For applications that require large memory, the CPUs can be equipped with more random access memory (RAM) to execute the application, whereas the GPUs are very limited in this regard. Other measures are used for evaluating parallel computations and their optimization. The fast memory access of the GPUs and their massively-parallel units are well suited for ordered data patterns. Vector operations perform better on the CPUs, whereas matrix operations are more ecient on the GPUs. It is for such reasons that we used a dual combination of CPU + GPU to carry out the computations. 3.5 Generation of the Fractional Brownian Motion The most time-consuming part of the generation of a PN/RRN is generation of a 2D FBM by a fast Fourier transform (FFT) that computes the spectral density, Eq. 3.2, with a time complexity ofO(N log 2 N). Another main issue associated with the sequential implementation of the FBM algorithm is that, for large values of network size N, huge memory is required. To address the problem, we parallelize the program by using the message-passing interface (MPI) strategy on four processes numbered (0, 1, 2, 3); see Fig. 3.4. Each processor executes a part of the FBM generation algorithm, and then sends the results to the root process, number 0 in Fig. 3.4. 33 Figure 3.4: Communications between four nodes (left) and point-to-point communication model. 3.6 Solving for the Pressure/Voltage Distribution with a GPU We used the CG algorithm in both C++ and CUDA. In CUDA, the CPU is a host for GPU, with the latter called a device. Programs that are run on the CPU (host) can manage the memory on both the CPU and GPU. The functions that are executed on the GPU are called kernels. The CPU lunches kernels, which and are then executed in a parallel scheme by GPU threads. As far as the CG algorithm (as well as other iterative solvers) is concerned, the most computationally intensive kernel is the matrix-vector multiplication (MVM) kernel (Barrett et al., 1994). Most of the execution time of the main loop of a solver is spent inside the MVM kernel. One important issue related to the GPU and CPU is the oating point accuracy of compu- tations. The main point is that the GPU devices with high computational capabilities support both single and double precision (SP and DP, respectively) oating point computation. However, since the SP calculations are executed much faster, a signicant speed gap exists between them and the DP computations, when they both are carried out with the GPUs. The main bottleneck 34 of the CG iterations is the ratio of the arithmetic operations and the data input/output steps, and the acceleration of the computations depends critically on this ratio. Thus, by using the SP we reduce by a factor of two the \trac" between the GPU processor and its device memory. We also use a mixed-precision method to increase the speed-up in the computations (see the SI for further details). To do so, two types of iterative loops are utilized to obtain the true solution of Eq. 3.4. We rst use the SP in the iterations within the inner loop to approach rapidly a rough solution of Eq. 3.4, followed by DP iterations within the outer loop for computing the precise solution of Eq. 3.4 within the pre-set error. The GPU that we used in this study supports both SP and DP computations. The iteration loops of the CG algorithm consume most of the computation time. We emphasize that we use a combination of CPU and GPU to carry out the PN/RRN simu- lation. Using such a combination of two distinct paradigms and programming models for solving a given problem is completely nontrivial and has, in fact, been discussed extensively in the past. In our work we rst generate the PN/RRN using four CPU processors, which are then used in a GPU to solve for the pressure/voltage distribution in the PN/RRN. The percentage of run times spent on the CPU host vary with the size of the network, as well as the number of connected resistors or pores in the network. As the network's size increases, the run time on the CPU also increases due to the limited memory of the GPU. In the problems that we study the GPU code is memory-bounded and most of the time is spent on memory and communication between the nodes. As we show below, in this problem a considerable time is spent for the communications and memory transfer. What this means is that the eciency will increase dramatically over what present in this study, if the GPU memory does so. Our preliminary tests indicated that such a combination of the CPU and GPU yields the most ecient performance because of the following. The distribution of the data in a dual computing system, such as CPU and GPU, is an important factor that can restrict the performance. Al- though both CPUs and GPUs are capable of performing similar types of calculations, there are 35 a) b) c) d) e) f) 50 100 150 200 Figure 3.5: Examples of the networks with correlated conductances. The top and bottom rows show, respectively, the networks with the Hurst exponent of H = 0:75 and 0.35, and x = y of (a) and (d) 1; (b) and (e) 1/5, and (c) and (f) 5. several important dierences in their execution models and architectures. The number of cores on a CPU is limited, but they are large and fast. GPUs, on the other hand, have hundreds of cores that are slower and, compared with the CPU cores, possess limited capability. It is expensive, in terms of the computer time, to transfer data from one processing unit from the CPU (GPU) to a unit in GPU (CPU), and vice versa. Generating a PN/RRN and solving for the pressure/voltage distribution, are, however, two independent problems that are solved sequentially. The fast mem- ory access of the GPUs and their massively-parallel units are ideal for ordered data patterns. Vector operations are carried out more eciently on the CPUs, whereas matrix operations are executed faster on GPUs. Thus, it is ecient to use several CPUs to parallelize the PN/RRN generation, because it requires very large computer memory. 36 10 2 10 3 10 4 10 5 10 6 0 2 4 6 8 10 12 14 Network Size Speed−up H a) 0.75 0.35 10 2 10 3 10 4 10 5 10 6 0 10 20 30 40 50 60 70 Network Size Speed−up H b) 0.75 0.35 Figure 3.6: a) Dependence of the speed-up on the size of isotropic networks with H = 0:75 and 0.35 b) Speed up for one iteration in CG algorithm. The computations were carried out with various PN and RRN sizes in both CPU alone and in GPU, in order to measure the speed-up of the computations. The computing system that we employed consisted of one Intel core i7-4710 CPU (2.5 GHz) with 16 GB of main memory, and the NVIDIA Geforce GTX 880 graphic card with 1536 cores. All speed-up comparisons have been made with respect to runtime of sequential PN generation and simulation algorithm in C++. Applications were written in CUDA version 6.5 and C++ using Visual Studio 2012. 3.7 Eciency of the Computations Since one goal of this chapter is demonstrating the eciency of the computations with standard PN or RRN models, we put the algorithm to use in a most stringent environment. In the former case, the PN network contains long-range correlations between the bonds' permeabilities, as well as anisotropy, as described earlier. In the case of the RRN, we delete a fractionq of the resistors and compute the voltage distribution. As q approaches the percolation threshold of the network, q = 1=2, there is usually critical slow-down because the structure of the sample-spanning percolation 37 10 2 10 3 10 4 10 5 10 6 0 2 4 6 8 10 12 14 Network Size Speed−up β x /β y a) 5 1/5 10 2 10 3 10 4 10 5 10 6 0 10 20 30 40 50 60 70 Network Size Speed−up β x /β y b) 5 1/5 Figure 3.7: a) Dependence of the speed-up on the size of isotropic networks and the anisotropy parameter x = y . The Hurst exponent isH = 0:35 b) Speed up for one iteration in CG algorithm. cluster becomes increasingly tortuous. Thus, a key test of the method is whether the gained eciency is lost as the percolation threshold is approached. For a 2D PN/RRN of size NN the speed-up of the computations is dened as the ratio of the sequential execution time t(1;N) and the corresponding parallel execution time t(P;N) with P processors, S(P) = t(1;N) t(P;N) : (3.5) In the PN/RRN simulation, the parallel part consists of both the CPUs (for four processors) and the GPUs. Thus, the parallel execution time is the sum of running time for both generating the PN/RRN on the CPU and solving the governing equations on the GPU. Another way of measuring the performance of the parallel implementation is the eciencyE(P), dened by E(P) = S(P) P = t(1;N) Pt(P;N) : (3.6) The eciency is only used for the parallel implementation of PN/RRN generation. In Table I we report the performance of the MPI implementation of the generation of the 2D FBM eld and the PN/RRN of various sizes, which indicate that the eciency and speed-up both increase as the size increases. 38 0 0.1 0.2 0.3 0.4 0.5 0 2 4 6 8 10 12 14 Fraction of Insulating Bonds Speed−up a) 0 0.1 0.2 0.3 0.4 0.5 0 10 20 30 40 50 60 70 Fraction of Insulating Bonds Speed−up b) Figure 3.8: The speed-up of the computations with a random resistor network with 3:6 10 6 nodes. Table 3.1: Eciency and speed-up of the MPI implementation of FBM algorithm on four pro- cesses. The times are in seconds. Array Size Required Memory (GB) Time (Sequential) Time (MPI) Speed-up Eciency 2 24 0.128 1.2 0.4 2.8 0.71 2 28 2 8.0 2.7 2.9 0.73 2 30 8 71.6 20.6 3.5 0.87 We then use the generated PN/RRN as the input to the GPU-based solver. Solving for the pressure/voltage distribution consists of the CG algorithm, which is matrix operation, well suited for the GPUs. Several PNs/RRN were generated with various values of the Hurst exponent H, the anisotropy parameter x = y , and network sizes N. In Fig. 3.5 we present samples of the generated PNs/RRN in which the conductances vary over two orders of magnitude. 3.8 Results We carried out extensive computations with a variety of PN and RRN in order to test the speed-up in the computations. In what follows we present and discuss the results. 39 0 0.1 0.2 0.3 0.4 0.5 0 2 4 6 8 10 12 Fraction of Impermeable Bonds Speed−up H β x /β y a) 0.75 1/5 0.75 5 0.75 1 0.35 1/5 0.35 5 0.35 1 0 0.1 0.2 0.3 0.4 0.5 0 10 20 30 40 50 60 70 80 90 Fraction of Impermeable Bonds Speed−up H β x /β y b) 0.75 1/5 0.75 5 0.75 1 0.35 1/5 0.35 5 0.35 1 Figure 3.9: a) Dependence of the speed-up of the computations for network with 3:6 10 6 nodes on the fraction of the removed bonds, the Hurst exponent and the anisotropy parameter x = y . b) Speed up for one iteration in CG algorithm. 3.8.1 Isotropic networks Figure 3.6 presents the eect of the Hurst exponent H on the speed-up of computing the pressure/voltage distribution. The network is isotropic [ x = y in Eq. 3.2], and each data point represents the average of multiple realizations. The speed-up increases with increasing the size of the PN/RRN. For a network withH = 0:75 and 3.6 million nodes, the speed-up is 12, more than an order of magnitude improvement in the computations. More importantly, the same speed-up is obtained with a FBM distribution with H = 0:35, the case in which the correlations between the permeabilities of conductances are negative, which is typically the case in natural porous media (Sahimi, 2011). Thus, the value of the Hurst exponent H has only a minor eect, if at all, on the speed-up of the GPU-based solver, which is strong evidence for the eciency of the algorithm. Next, we consider the case of the classical RRN model in which a randomly-selected fraction p = 1q of the bonds have a unit conductance, while the rest are insulating. The speed-ups are shown in Fig. 3.7. As the percolation threshold q c = p c = 1=2 is approached, the speed-up decreases from about 9.6 forp = 1 to about 7.5 very close top c . Thus, even though, similar to any 40 a) b) c) d) e) f) Figure 3.10: Pressure distribution in a 50 50 network with the Hurst exponent ofH = 0:75 and 0.35, and x = y of (a) and (d) 1; (b) and (e) 1/5, and (c) and (f) 5. critical phenomenon, there is critical slow down in a RRN model as p c is approach, the speed-up obtained by using the GPU decreases by only about 20 percent, hence indicating the eciency of the GPU-based computations. 3.8.2 Anisotropic networks and permeability anisotropy As mentioned earlier, natural porous media, as well as many solid materials, are anisotropic. Thus, we studied the eect of the anisotropy parameter x = y on the performance of CPU+GPU parallel computations. The results are shown in Fig. 3.9, where we present the speed-up for two values of of x = y and a Hurst exponent ofH = 0:35. The speed-up is essentially the same that of the isotropic media, and also independent of the value of x = y . Thus, although in the traditional approaches the computations for anisotropic media take longer times, in the proposed algorithm both the Hurst exponent H and the anisotropy parameter x = y have a minor eect, if any, on the speed-up gained by the dual CPU-GPU. 41 10 2 10 4 10 6 Network Size 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 K x /K y x / y a) 5 1/5 10 2 10 4 10 6 Network Size 0 0.5 1 1.5 2 2.5 3 3.5 K x /K y x / y b) 5 1/5 Figure 3.11: K x =K y ratio for networks with various sizes: (a) H = 0:75 and (b) H = 0:35. Next, we consider the eect of percolation, i.e., attributing zero conductance to the bonds in PNs or RRNs in which the conductances are correlated, with the correlations described by the FBM. The motivation for considering such networks is that natural porous media, in addition to be mostly anisotropic, are highly heterogeneous and contain very low-permeability zones that contribute very little, if any, to uid ow. Thus, their permeability or hydraulic conductance may be set to zero. Composite materials that consist of conducting and insulating phases behave the same. But, since the conductance eld is correlated, one cannot remove bonds randomly, which would destroy the correlations. Instead, the bonds' conductances are sorted from the smallest to the largest, and a given fractionq of the bonds with the smallest conductances are removed. The results are presented in Fig. 3.10. In all the case the network contained 3:6 10 6 nodes. The speed-up for all the cases has a trend similar to Fig. 3.9. Note that as the percolation threshold q c = 1=2 is approached, the speed-up decreases. As pointed out earlier, this represents critical slow down as the percolation threshold is approached. However, even very close toq c the speed-up is nearly one order of magnitude. An important problem in modeling of anisotropic porous media is the ratio of the permeabilities K x =K y , where x represents the direction that is more or less parallel to the strata that make 42 such porous media anisotropic. There are many empirical estimates of the ratio K x =K y in the literature (Sahimi, 2011) that vary between 1 and 7. Using our ecient solver, and generating the permeability distribution according to Eq. 3.2, we obtain approximate, but accurate estimates for the permeability anisotropy. The results are presented in Fig. 3.11. For H = 0:75 the anisotropy varies between 0.6 and 1.6, whereas the corresponding numbers for H = 0:35, i.e., the much more heterogeneous PNs, are 0.5 and about 2.5. 3.9 Discussion Several important points must be mentioned here. (i) Generally speaking, unlike CPU-based computations, the denser the matrix G is, the more ecient the GPU-based computations are. Thus, similar calculations for 3D systems and with networks with higher connectivities should yield higher speed-ups. (ii) The calculations with the lattice models of vector transport, such as brittle fracture propagation and similar phenomena will be even more ecient than the scalar transport considered in this study, because the corresponding matrix G will be denser, even for 2D systems. (iii) Clearly, any of the past methods of speeding-up the calculations will be more ecient if carried out with the GPUs. (iv) We utilized a single GPU-based solver. Use of several GPUs in parallel should result in very high speed-up factors, a matter that we are currently studying. 3.10 Summary We proposed an algorithm for simulating the pore-network and random-resistor network mod- els on GPUs. The single GPU-based solver improves the eciency and running time of the calculations by at least one order of magnitude, when compared with the same calculations with CPU-based simulators. The speed-up is even larger for larger networks, and for calculations in 43 which the matrix of the coecients is denser. Thus, the GPU-based computations make simu- lation of very large PNs/RRNs possible, regardless of the correlations, anisotropy, or any other pertinent parameter. 44 Chapter 4 Drying of Layered Porous Medium: The Eect of Pore-Size Correlation and Disorder 4.1 Introduction The importance of randomness and correlation between pore and throat size in both natural and engineered porous materials has been proved and studied (Holtzman, 2016; Fantinel et al., 2017). For example, in many engineering applications such as membranes and lters, a low disorder in the size of porous and grains are desired to produce a uniform ow pattern through the porous structure. One of the inherent properties of natural porous medium and geological formations if the correlation between pore and grain sizes (Sahimi, 2011). Numeruos studies proved that (Molz et al., 1997) many natural porous media exhibit long-range correlations. Such correlations, which in uence the percolation (Sahimi and Mukhopadhyay, 1996) and ow and transport properties of porous media (Knackstedt et al., 2001, 1998; Dashtian et al., 2015), follow the statistics of fractional Browning motion (FBM) (Mandelbrot and Ness, 1968), a nonstationary stochastic process. Muljadi et al. (2016) investigated the eect of pore-scale heterogeneity on the onset of non- Darcy ow behaviour. They used direct ow simulations through 3D images of a beadpack, 45 Bentheimer sand- stone and Estaillades carbonate. Alim et al. (2017), studied the relationship between the microstructure of porous medium and ow distribution in the pores. They used both analytical models and numerical simulations to show the transition of ow distribution in 2D porous medium with increasing disorder. Datta and Weitz (2013), used confocal microscopy to study the spatial uctuations in uid ow through a three dimensional porous structures. They indicated that the pore-scale correlations in the ow are determined by the geometry of the medium which suggest that despite the complexity of the pore structures, uid ow through it is not completely random. Volpe and Volpe (2017) showed that the structure of the porous environment signicantly alters the optimal search strategy of a particle (searcher) performing a blind cruise search for uniformly distributed non-regenerating sparse targets. Fantinel et al. (2017), showed through experimental and simulation studies that reducing the pore-size disorder in porous 2D media increases the displacement eciency while decreasing the uid- uid interfacial area. Liu et al. (2012) and Liu et al. (2017) showed that the heterogeneity of porous medium is one of the main factors in which in uence the formation of wormhole in carbonate reservoir during acidizing operation. Viscous ngering and uid trapping were accounted as the mechanisms that in uence the sweep eciency. There are several porous medium related processes in which preferential pathways and heterogeneous uid distributions is preferred over a uniform uid phase distribution (Das et al., 2004). These applications include acidizing of carbonate reservoirs where uid mixing and reactions are needed (Liu et al., 2017). On the other hand, there are several other processes in which a heterogeneous uid distribution is undesirable. These processes involves enhanced oil recovery when ecient sweep of hydrocar- bons are needed (Lake et al., 2014). Drying and salt precipitation in porous medium, similar to the mentioned, porous medium related phenomena, are strongly aected by the heterogene- ity and morphology of the medium. The spatio-temporal dynamics of a phase distribution and drying front in a porous structures is controlled by two categories of factors. The rst category 46 includes the static properties of porous medium such as available pore space, pore and throat size distribution, morphology and topology of fractures, pores and grains (Prat, 2002; Fantinel et al., 2017; Borgman et al., 2017). The second category of factors concerns the pore scale phenomena, such as capillary forces-driven liquid transport (Scherer, 1990b; Shokri et al., 2008a; Le et al., 2009), in-pore evaporation and vapour diusion (Miri et al., 2015a; Lehmann and Or, 2009) and formation of liquid lms on the pore surface (Chauvet et al., 2009; Lehmann and Or, 2009; Yiotis et al., 2001). The eects of spatial correlations and disorder between pore size of porous medium has not been studied in three dimensional systems. The goal of this study is to implement a 3D PN model of drying and evaporation in porous medium which includes heterogeneity and anisotropy. In a recent study, we have shown that correlations aects the uid phase distribution and salt precipitation in a drying 2D porous medium (Dashtian et al., 2018). Borgman et al. (2017), also included such eects in their experimental and modeling studies and increasing the correlation length in porous medium prolongs the liquid connectivity to the surface and hence increasing the drying rates. Such models have been developed for 2D systems. The transport phenomena in 3D porous structures with heterogeneity is much more complex than 2D systems. Also the computational complexity of implementing evaporation in 3D is a bottleneck which limits the size of such systems. The literature of evaporation and drying from porous media is inundated with experimental and theoretical studies. The evaporation from porous space is also a displacement process, in which the liquid phase in the pore-space recedes as a result of vaporization into gaseous phase. Thus, drying can be viewed as a drainage process, in which the non-wetting gas phase displaces evaporating wetting liquid. Therefore, in the absence of viscous eects, the pattern of the receding interface is described by invasion percolation. Pore network modeling is a discrete model and powerful approach to study the link between the discrete transport pore level phenomena and the macroscopic properties at the sample scale (Sahimi, 2011; Mohammadmoradi and Kantzas, 2016; 47 Figure 4.1: Examples of correlated FBM elds. The top and bottom rows show, respectively, the PNs with the Hurst exponentsH = 0:75 and 0.35, and anisotropy x = z = (a) and (d) 1; (b) and (e) 1/10, and (c) and (f) 10. Figure 4.2: Size distribution of the throats same as Figure 4.1. 48 Mehmani and Tchelepi, 2017; Balho et al., 2012; Mehmani et al., 2013). Simplied and more sophisticated pore network models have been developed to understand the pore-scale dynamics of drying. Most of the current PN models are applicable to 2D porous structures (Yiotis et al., 2001; Le et al., 2017; Rahimi et al., 2016). Three dimensional models have also been proposed by (Yiotis et al., 2015; Bray and Prat, 1999; Moghaddam et al., 2017). Recently, Vorhauer et al. (2017) proposed a 3D PN model for drying of thin porous medium in cylindrical coordinate. Using the IP process, many PNs have been developed for simulating drying (Prat, 1995; Yiotis et al., 2001; Surasani et al., 2008). The performance of PN models of evaporation in porous media have been tested in the previous studies. In our modeling, we use the same approach to simulate the evaporation process. As we mentioned before, one important aspect of porous structures is the existence and extend of heterogeneity and disorder in their microscopic and macroscopic properties (Dashtian et al., 2015) at various length scales. Such eects has not been implemented in existing 3D PN models of drying and evaporation. The importance of the correlation and disorder in pore/throat size of porous medium and its eect on the uid phase distribution, ngering and uid displacement has been proved (Zhao et al., 2016). z x y (a) (b) (c) Figure 4.3: The phase distributions in an isotropic pore network with the Hurst exponentH = 0:75 ( x = z = 1), at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.9. The top row of the network is open. 49 Figure 4.4: The phase distributions in an isotropic pore network with the Hurst exponentH = 0:35 ( x = z = 1), at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.9. 4.2 Generating Correlated Pore Network The rst step in the PN/RRN modeling is to generate the networks. As we discussed above, the main focus of this study is to include the spatial correlations of pore sizes in the 3D PN model of porous medium. Thus we use the FBM to generate correlated pore/throat size distributions. We use a cubic connected graph. Each edge in this graph represents throat with an eective radius r. We assume that the length of the edges (bonds) ` are constant. The throats' radii are selected from a FBM. The existence of fractal correlations and long-range memory in the porosity data of large-scale porous structures was rst proposed by Hewett (1986b), followed by similar work by (Neuman, 1994) for the permeability distribution. Such properties have been investigated in numerous studies (Sahimi and Mukhopadhyay, 1996; Dashtian et al., 2015). FBM is a non-stationary process B H (r), and should have the following properties: hB H (r)B H (r 0 )i = 0; (4.1) 50 and h[B H (r)B H (r 0 )] 2 ij r r 0 j 2H ; (4.2) In which, H is the Hurst exponent and r is the position in which for a 3D system we have r = (x;y;z). The two-point correlation functionC(r) of a FBM is given by,C(r) =C 1 r 2H , where C 1 =C(r = 1). An ecient method for generating a FBM array is through its power spectrum S(!), which for d dimensional systems is given by S(!) = a(d) P d i=1 i ! 2 i H+d=2 ; (4.3) where a(d) is a constant, while ! = (! x ;! y ;! z ) are the Fourier components in the x and y and z directions and i = ( x ; y ; z ) are the anisotropy parameters for 3D systems. Hurst exponent, H, for H > 0:5 (H < 0:5) one has positive, long-term memory (negative, short-term memory) correlations. H = 0:5 represents a random walk process. Due to strati- cation, many porous media are anisotropic. To generate anisotropy induced by layering with the layers being parallel to the x direction, we set x = y > 1, and similarly for the layers being parallel to the y direction. Figure 4.1 presents examples of the heterogeneities generated by the model. 4.3 Pore Network Model Pore network models have been used to simulate drying and evaporation in porous medium over two decades (Prat, 2002). We use a cubic network in which each pore is connected to 6 neighboring pores through cylindrical throats. Pressure build up and saturation of gas and liquid in pore bodies are ignored. We also ignore the eect of liquid lm in the throat and pore bodies. We consider the viscous eects in the liquid phase. In drying and evaporation in porous media, 51 z x y (a) (c) (b) Figure 4.5: The phase distributions in an isotropic pore network with the Hurst exponentH = 0:75 ( x = z = 1=10), at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.9. z x y (a) (b) (c) Figure 4.6: The phase distributions in an isotropic pore network with the Hurst exponentH = 0:35 ( x = z = 1=10), at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.9. the practice is to assume that mass transfer occurs in gas phase by diusion (Prat, 2011; Yiotis et al., 2001). The driving force for evaporation, thus, is the concentration gradient alongside the drying front. Therefore, our model must include diusion in gas phase. Assuming steady-state and slow (laminar) ow of an incompressible uid, the volume ow rate in a bond ij that connects nodes i and j is given by, q ij =r 4 P ij =(8`) =g ij P ij , where is the uid's viscosity, P ij = P j P i is the pressure drop between nodes i and j, and g ij the conductance of bond ij. Writing a mass balance for every node i, one obtains, P j2i q ij = 0, 52 where the sum is over all the nodesj that are connected to nodei. Thus, substituting forq ij and writing the mass balance for every interior node of the network results in a set of linear equation of the following form, GP = b; (4.4) where G is the conductance matrix that depends only on the geometric properties of the network, P is the nodal pressure vector to be calculated, and b is a vector related to the external boundary conditions. The PN's top row is assumed to be a layer whose throats are initially empty and through which a gas (such as air or CO 2 ) ows into the pore space, while the other sides are sealed. Liquid ows in the pore space's throats, while the vapor diuses to the boundary layer and induces further evaporation. In the IP model a gas invades the brine-saturated porous medium throat-by-throat according to the lowest capillary pressure needed for entering a throat at the interface between the vapor and liquid. The capillary pressure for entering a throat of radius r ij between pores i and j at the liquid-vapor interface is given by P c;ij = 2 r ij ; (4.5) where is the surface tension. At every step of the simulation we identify the state of the pores, which is determined by the lling state of the pore throats connected to them: if at least one neighboring throat does not contain brine, it is a vapor pore; otherwise, it is a liquid pore. A throat can be lled, fully or partially, with brine, or it can contain vapor that is leaving it through a path of pores and throats that are also lled by the vapor. A vapor pore can be either at equilibrium or at an unknown vapor pressure (to be determined). In reality, equilibrium vapor pressure prevails only at the menisci between brine and the vapor, as shown in Fig. 4.2, but we assume that a pore is at equilibrium vapor pressure if at least one of its neighboring throats still contains brine, even if 53 it is only partially so. Therefore, diusion of vapor between menisci of partially-lled throats and the adjoining vapor pore happens without any resistance. Our preliminary numerical simulation indicated that the assumption leads to slightly, only about 5 percent, overestimates of the drying rate, but simplies the computations greatly. Once the states of pores and throats are identied, the second stage of the simulation, namely, modeling of transport of brine throughout the PN and the resulting salt precipitation, begins. After one or a few steps of the IP that advance the drying front, we solve the governing transport equations. Assuming that the ow is slow enough, the liquid ow rate q (l) ij in a throat ij that connects pores i and j is given by the Hagen-Poiseuille equation q (l) ij = r 4 ij 8 (l) ` ij h P (l) i P (l) j i ; (4.6) where the liquid viscosity (l) is assumed to be constant, as we simulate an isothermal process. Here, P (l) i is the liquid pressure in pore i, and ` ij is the length of the throat ij. Yiotis et al. (2001) showed that both advection and viscous eects aect the drying patterns and rates. Since there is no liquid accumulation in the pores, we must have X fijg (l) q (l) ij = 0; (4.7) with (l) being the brine's mass density. The sum is over the set of the throatsfijg that are connected to pore i. The vapor is transported by diusion and, thus, we use the equation for evaporation in the so-called Stefan tube (Bird et al., 2007) to calculate its mass ow rate in the throats. Ignoring the pressure in the vapor-lled throats, vapor mass conservation at the pores must hold, implying that, X fijg (v) q (v) ij = X fijg S ij D (v) ` ij PM (v) RT ln " PP (v) i PP (v) j # = 0; (4.8) 54 where D (v) is the vapor diusion coecient, M (v) its molecular weight, P is the total pressure in the vapor phase, and (v) is the vapor mass density. We used, (v) = 8:12 10 4 gr/cm 3 and D (v) 2:6 10 5 m 2 s 1 , both corresponding to 298 K. The vapor concentration at the vapor- liquid interface is equal to the equilibrium concentration, which is zero outside the boundary layer of the PN. Writing down Eqs. 4.6 - 4.8 for all the pores in the liquid and vapor phases results in two sets of equations that govern the pressures in the pores in the two phases, which are solved by the conjugate-gradient method. The pressure at the liquid-vapor interface is assumed equal to the equilibrium vapor pressure, P = 2339 Pa. At a meniscus between the brine and vapor the Kelvin equation describes the dependence of parameters on the vapor pressure: ln P (v) P = 2M (v) RTr ij (l) = P c M (v) RT (l) ; (4.9) where R is the gas constant. z x y (a) (b) (c) Figure 4.7: The phase distributions in an isotropic pore network with the Hurst exponentH = 0:75 ( x = z = 10), at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.9. 55 z x y (a) (b) (c) Figure 4.8: The phase distributions in an isotropic pore network with the Hurst exponentH = 0:35 ( x = z = 10), at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.9. 4.4 Results and discussion Despite numerous experimental and theoretical studies related to the drying and evaporation from porous media, the eects of correlations, disorder and heterogeneity has not been fully understood. Such properties has enormous practical importance in natural porous structures such as underground reservoirs in which the interplay between capillary and viscous forces determines the uid phase distribution in a 3D layered porous structure. This requires incorporating the pore-scale variations of pore/throat radii in the model. Almost all the existing evaporation and drying models uses a random normal distribution to generate the sample porous medium. As we discussed before, the path taken by a wetting and non-wetting phase is altered by a slight variation in a single pore. Thus, similar to other uid displacement processes in porous media, drying is strongly aected by the correlation of throat sizes. The common practice in investigation of transport processes in stratied porous medium is to assume several parallel layers each of them having a unique average pore size distribution, porosity or permeability. Thus, in each layer, the properties of porous medium is random. In this study, we used a modied fBm stochastic process to generate the throat/pore radii. The H exponent controls the correlations and disorder in 56 each layer and the entire 3D porous medium while the i parameter controls the stratication in each desired direction. Thus, our approach has Superior advantages over current methods for generating porous media. In Figure 4.1 we show a sample of the generated fBm elds. The initial input to the algorithm which generates fBm is a PSD to generate random numbers, H exponent and the ( x ; y ; z ). Figure 4.1 shows some sample outputs. As can be seen from the gures, the correlations and stratication are evident and function ofH and i . Also, the corresponding PSD of the resulted PN is shown in Figure 4.2. We used similar PSD as input to the fBm generating algorithm but as can be seen in this gure, it produces fBm elds with dierent PSD. In Table. 1 we show the properties of the sample PN that we used to simulate the drying. Table 4.1: Pore size distribution properties of dierent cases. H x = y x = z y = z Mean throat radius (m) 0.35 1 1 1 120.0 0.35 1 1/10 1/10 119.3 0.35 1 10 10 116.8 0.75 1 1 1 112.2 0.75 1 1/10 1/10 123.4 0.75 1 10 10 117.8 After generating the sample PN, we apply the PN drying algorithm (Dashtian et al., 2018; Prat, 1995; Yiotis et al., 2015). We assume that the top row of the PN is emptied and exposed to the atmosphere. Such assumptions are compatible with actual porous medium. For example, when injecting natural gas or CO 2 into fractured reservoirs(Hamzehpour and Khazaei, 2016; Amooie et al., 2017), the liquid phase, which is usually trapped in the matrix of the rock due to the capillary forces, may evaporate and ow through the fractures as a gas phase. In most drying PN simulations, the porous medium is assumed to be hydrophilic and the wetting phase is assumed to be the liquid phase and non-wetting is gas phase. We also have the same type of assumptions but the algorithm is extendable to hydrophobic porous materials. 57 As we pointed out earlier, our PN model for drying is indeed a Invasion Percolation (IP) Vorhauer et al. (2017), in which the next throat to be invaded by the non-wetting phase has the smallest capillary threshold. After some time that the simulations proceed, there will be more than one liquid cluster in the system, thus, each cluster will follow its own cluster rule. 4.4.1 Long-range correlations and short-range disorders First, we start with the isotropic cases in which only the long and short term correlations changes. Thus we set the x = y = z but apply several values of H exponent to generate PN. Shown in Figure 4.3 is the evolution of the wetting and non-wetting phase distribution in an isotropic PN ( x = y = z ) with the H = 0:75, implying long-rage correlations in the throat size distribution. Figure 4.3 (a) represent the early stage of drying in which the drying rate is essentially constant. Close to the open end, the throats empty (dry out) at a faster rate than those deeper in the pore space. Figure 4.4 shows the case in which H = 0:35, a porous medium with short-range correlation. The drying front invades from the open end surface toward the bottom. The invasion leads to the disconnection of the wetting phase which interrupts capillary mass transfer. As the drying proceeds evaporated liquid/moisture is transferred in the gas phase. The gas concentrating build up decreases the drying rates. This also aects the rate at which clusters dry out. Clusters which are formed near the open boundary of porous medium are emptied faster because are subjected to a faster evaporation. Figure 4.3 (b) and (c) shows the phase distributions at the later stages of drying. As can be seen the sizes of the liquid clusters closer to the open boundary are smaller than those are formed at the bottom of the porous medium. These variations show the eects of competition between mass transfer and viscous eects (Shokri et al., 2008b) in which in this case, mass transfer overcomes the viscous forces. The drying patterns for the case with H = 0:35 is totally dierent for what we have seen in gure 4.3 for H = 0:75. The sizes of the liquid clusters are smaller for the case with smaller H. Positive correlation (H > 0:5) augments the probability 58 of connectivity of throats with similar sizes and hence invading phase has accessibility to throats with minimum capillary pressure (larger throats) through the PN. As described by Borgman et al. (2017) and Dashtian et al. (2018) the larger throats will dry out rst while the smaller pores remain saturated with wetting phase. This further indicates the role of the underlying geometry of porous medium on the displacement of wetting phase. Figure 4.9: Dynamic evolution of the transversely-averaged brine saturations in six types of pore networks. Black curves show the average saturations; cyan, blue and pink shadings represent, respectively, the times at which vapor saturation is 0.1, 0.2, and 0.6, while red indicates the time at which vapor saturation is 0.9. The colored areas indicate the variations over multiple realizations. 4.4.2 Stratication In geological formations, sedimentary rocks are formed via a variety of processes such as weath- ering of rocks, deposition and transport of the weathering rocks, cementation and compaction. Stratication is a key feature of sedimentary rock and can occur on the scale of micrometer up 59 Figure 4.10: Dynamic evolution of the average drying rates (black curves), averaged over 10 realizations (shadings). The gray areas indicate variations over multiple realizations. to hundreds of meters. Small scale stratication is due to alternately operating depositional pro- cesses in the same environment. Through experimental and numerical studies, Zhang et al. (2013), such small scale stratications have a dominant eect on the capillary pressure eld, CO 2 /brine saturation patterns, ow regime and apparent relative permeability of sandstone. (Shokri et al., 2010), developed a model to predict the drying front depth and veries the model with experi- mental data. They showed that presence of interfaces in stratied porous medium aects the uid phase distribution, drying rates and other dynamics of evaporation form porous structures. Vorhauer et al. (2017) showed that the impact of the nonuniform distribution of evaporation rate becomes increasingly less important with increasing number of porous layers. The transition from thin porous medium, where the drying process is highly dependent on the local structure of the external mass transfer, to thick porous medium, where the drying process is much less sensitive to the detail of the external mass transfer at the surface, is progressive under isothermal conditions. As we discussed in the previous section, we used the x = y and x = z to introduce anisotropy and layering in the porous medium. We used a x = y = 1 and x = z = 10 to create layers parallel to the open end of the porous medium. Both short (H = 0:35) and long (H = 0:75) 60 range correlations were considered. In Figure 4.5 and 4.6 we show the phase distributions for the case in that stratication is perpendicular to the open end of the porous medium, for H = 0:75 and H = 0:35 respectively. There is an evident dierence between the uid distributions of case with no layering (Figure 4.3 and 4.4) and with layering (Figure 4.5 and 3.6). The presence of stratication and anisotropy, modies the liquid phase distribution through preferential drying of throats in a certain direction. Also the drying curves associated with the porous medium are dierent. Figure 4.7 and 4.8 show the dynamic evolution of phase distribution for cases in which x = z = 1=10 for H = 0:75 and H = 0:35. In these cases stratication is parallel to the open end to the porous medium. As can be seen, the drying front prefers to invade throats in parallel to the open end, rst, then expand through the depth of the porous medium. We further analyze the behavior of the drying 3D porous medium, with long-range correlations, by calculating the evaporation rates and liquid saturation distribution through the PN. Figure 4.9 shows the transversely-averaged saturations in the same PNs of Figures 4.3-4.8. Four dierent time steps were selected to plot the saturations (averaged over the z direction). We used blue color to show the range of variations over multiple realizations and the black curves indicating the averages over all the realizations. Comparing Figure 4.9 (a)-(c) with Figure 4.9 (d)-(f) the stronger eect of positive correlations with H = 0:75 on the saturation curves is evident. PNs with short-range disorder (H = 0:35) show a narrower range of variations in the saturations over multiple realizations. The saturation proles for the case in which stratication is perpendicular to the open end of the porous medium ( x = z = 1=10), the saturation are dierent. Over the entire times, there is no region in which the saturation is zero and the liquid phase stays connected to the surface. In order to quantify our discussion on saturation curves, we note that there exist three distinct regions in saturation curves. The rst region represents the emptied PN with zero saturation. The second region is associated with a transition zone in which 0 <S < 1, and the third region 61 presents the PN regions with, S 1. Generally, at each time step, the slope of the saturation cure, @S @z , for case with x = z 6= 1 is smaller than the case with x = z = 1. The invasion of the drying front into the PN increases the resistance imposed by the emptied larger throats in the dried region of the PN. This resistance, slows down the evaporation and vapor transport dominates drying mechanism. As we discussed before, the evaporation process is divided into three main stages: (i) the constant rate period; (ii) the intermediate falling rate, and (iii) the residual and nal drying rate. Some studies assume a two stage process for drying of porous medium (Or et al., 2013; Shokri et al., 2010). The initial stage of drying is independent of the heterogeneity and disorder of the morphology of porous media and is controlled by conditions such as the ow velocity of the gas over the open end of the porous medium. The rst stage is characterized by as a very short period in the initial time steps of simulations (see Figure 4.10). Drying rate in the second stage is a direct function of the ability of the porous medium to supply the emptied zone with wetting phase. Thus, this stage depends strongly on the morphology of the pore space, and hence the short and long range correlations. Comparing two the cases, the intermediate stage with H = 0:75 and x = z = 1=10 has the highest evaporation rate, as shown in Fig. 4.9 (b). In this case the liquid stays connected to the surface of the pore space and, therefore, evaporation can happen in throats near the surface. To better understand the results, consider ow of brine through the PN between two points separated by several throats. If the throats' sizes between the two points are positively correlated, there will be smaller resistant against the ow. In the PN in which the strata are parallel to the y direction - Figs. 4.9(c) and 4.9(f) - the falling rate is slower and, thus, this stage lasts until almost the end of the process. The trend of evaporation curves and saturation proles in our study are similar to those reported by experimental studies (Shokri et al., 2010, 2008a,b). Deviation from experimental data is reasonable since we do not include parameters such as condensation, adsorption (Bakhshian and Sahimi, 2016, 2017) or temperature in our study. 62 4.5 Summary In this study we investigated the eects of correlation and stratication of the drying and evaporation from porous medium. We use 3D pore network modeling and IP to simulate the evaporation. We incorporated the stratication and correlation into pore network modeling using the FBM. Our results indicate that stratication and correlation aect the evaporation from porous media. The stratication perpendicular to the open end of porous medium presents the highest evaporation rate. In this case the liquid stays connected to the surface of the pore space and, therefore, evaporation can happen in throats near the surface. This leads to higher drying rates and longer constant drying-rate periods. 63 Chapter 5 Pore-network model of evaporation-induced salt precipitation in porous media: the eect of correlations and heterogeneity 5.1 Introduction Salt precipitation in porous media occurs in numerous systems, including in soil that represents a global issue aecting salinization of agricultural land (Nachshon et al., 2001; U. et al., 2011; Schoups et al., 2005) and in geological formations that are being tested for CO 2 sequestration (Baumann et al., 2014b). Several mechanisms contribute to salt precipitation in porous media, the most important of which is drying and evaporation/wicking of saline water (Veran-Tissoires and Prat, 2014), as well as mineral replacement and reactions that contribute to salt crystallization (Dashtian et al., 2017) and growth in soil and rock (Noiriel et al., 2010). Salt precipitation damages porous media by reducing their porosity and permeability (Bacci et al., 2011), causing mechanical failure (Schiro et al., 2012). Thus, numerous theoretical and experimental investigations have been undertaken to understand the eect of various factors on drying of brine and salt distribution (Shokri et al., 2008b; Shokri and Or, 2011; Norouzi Rad and Shokri, 2012; Norouzi Rad et al., 2013; Bergstad and Shokri, 2016) in order to develop methods for limiting the damage. In particular, 64 Shahidzadeh-Bonn et al. (2010) investigated the eect of the kinetics of salt dissolution and crystallization on salt-induced damage to sandstones, while Flatt (2002) studied the eect of salt concentration on the extent of damage to porous media. Salt precipitation is preceded by evaporation of saline water. Thus, many experimental studies of brine evaporation (Shahraeeni and Or, 2010; Shahraeeni et al., 2012; Shokri and Or, 2011) and precipitation of salt (Shokri et al., 2010; Bergstad and Shokri, 2016) at the pore scale have been reported, elucidating the eect on the drying rates of several factors, such as the pore-size distribution (PSD), connectivity of the pore space - the way pore bodies and pore throats are connected to one another - and the wettability. For example, experiments with various PSDs indicate that two types of eorescence, crusty and patchy, may form on the surface of porous media (Veran-Tissoires and Prat, 2014) that aect the evaporation rates dierently. Thus, the PSD is a major component of porous medium that should be considered in modeling of drying- induced salt precipitation. But, despite the large number of theoretical and experimental studies of salt crystallization in porous media, there is still no comprehensive model and ecient numerical simulation method that can predict salt precipitation patterns in porous media at the pore scale, usually referred to as sub orescence. Salt crystallization in a porous medium is a complex molecular-scale phe- nomenon (Dashtian et al., 2017) that is in uenced by the thermodynamic state of the system - temperature, pressure, and humidity - and other factors, such as the salt type and the morphol- ogy of the medium. Even repeating a given experiment in a sample porous medium undergoing salt crystallization does not guarantee the same precipitation patterns. Such complexities in the physics of the salt precipitation in porous media exacerbate the inherent diculties of its model- ing. The goal of our study is to propose a pore-network (PN) modeling approach that is capable of predicting the spatial distributions of brine, vapor and the precipitated salt. To develop a physically-based PN model, one must take into account the mechanism by which the precipitated salt accumulates. Under the right conditions, Na + and Cl ions in a solution 65 give rise to nucleation and growth of salt crystals. Experimental studies by Shahidzadeh-Bonn et al. (2010) demonstrated the eect of the salt type on the extent of damage to the pore space of porous media. They showed that the formation of the hydrated sodium sulfate crystals enhances the spreading power of the salt solution on the PN network of sandstone, giving rise to very rapid growth of the hydrated phase of sulfate and forming clusters. Sodium chloride, on the other hand, only forms anhydrous crystals, which explains why less eorescence is observed compared to sodium sulfate. Shahidzadeh-Bonn et al. (2010) did not, however, consider the eect of the chemistry of the pore surface of porous media, or the co-existence of several types of salt. In our PN model, we consider only the formation of sub orescence and, thus, we work with NaCl as the salt in the study. Some experimental and theoretical studies indicate that salt accumulates in a very narrow layer at the soil-atmosphere interface, where the initial evaporation occurs, after which precipi- tation within the pore space is limited. On the other hand, experimental studies of Roels et al. (2014) with core samples indicated that capillary-driven ow does not have substantial eect on salt precipitation, because of which accumulation of salt near the core surface was low. Their computed-tomography images after drying indicated uniform salt distribution inside the core. B onhorst et al. (2016) reported that high evaporation rates limit transport of solids on the exter- nal surface, leading to a more uniform precipitation distribution, whereas low drying rates result in strong solid accumulation on the porous medium's external surface, but that small throats reduce the solute transport and restrict solid accumulation on the surface. One must also rec- ognize the self-enhancing of salt growth and salt transport in water lms (Miri et al., 2015a). Aggregated salt crystals form a microporous medium (Malmir et al., 2016b,a, 2017) with a high degree of capillarity in each throat that enhances solute transport through imbibing brine over long distances. But, since such phenomena as the growth of salt crystals at the drying interface occur at the molecular scale (Dashtian et al., 2017), it is not possible to include them in a PN model. Instead, one must resort to molecular dynamics simulation (Dashtian et al., 2017). 66 Kim et al. (2013) reported on two distinct salt crystallization in synthetic porous media. One was large bulky crystals that form away from the vapor-liquid interface, while aggregated, randomly-oriented crystals that form near the drying interface constituted the second type. In addition, experiments in circular capillary tubes have shown that salt precipitation is also aected by the boundary conditions imposed on the system. When the two ends of capillary tubes are open, salt crystals tend to migrate even outside the tube and form a patchy packing. In closed- end capillaries, however, salt crystallization occurs inside the tubes (Miri et al., 2015b; Miri and Hellevang, 2016). Therefore, the question of whether salt precipitates mainly on the external surface of a porous medium, or more uniformly in its pore space, depends on the medium's morphology, the competition between convection and capillarity, and the boundary conditions. In the PN model that we describe below, salt precipitate occurs inside the pore throats and, after precipitation, its crystals do not migrate from one throat to other throats. At the eld scale salt precipitation is also in uenced by the initial spatial distributions of the porosity and permeability of a porous formation (Zeidouni et al., 2008; Celia et al., 2011). Several experimental studies of salt precipitation in geological formations during CO 2 injection have also been reported (Muller et al., 2009; Wang et al., 2014; Bacci et al., 2011; Ott et al., 2011; Peysson et al., 2014), indicating that salt precipitation reduces permeability and porosity in Berea sandstone. This supports the view that in such formations crystallization occurs inside the pore throats in the pore space, rather than migration to the outside surface. Given such insights, the eect of such distributions on the eorescence formation has also been studied (Veran-Tissoires and Prat, 2012, 2014). In particular, Sghaier et al. (2014) reported on the in uence of evaporation ux from the boundary on the sub orescence formation and growth, while the eect of grain angularity was studied by Norouzi Rad and Shokri (2014). The porous media used by Sghaier et al. (2014) and Norouzi Rad and Shokri (2014) represented, however, isotropic medium, and did not also include the correlations in the grain or pore sizes that are prevalent in natural 67 porous media (Sahimi, 2011). One major purpose of our study is to show the eect of pore-size correlations on drying and the precipitation patterns. Models of evaporation and salt transport and precipitation in porous formations can be divided into two groups. In one group are the thermodynamic models of salt deposition in Earth's crust; see, for example, Toner et al. (2015). The second groups of models contains two subgroups. Continuum models of uid ow and transport (Huinink et al., 2002; Guglielmini et al., 2008) - those that are based on the classical averaged equations of hydrodynamics and transport - belong to the rst subgroup. They are capable of modeling evaporation and salt precipitation at the eld scale. They include models for the evolution of salt crystals (Le et al., 2009) and salt concentration in liquids (Mahadevan et al., 2006) in porous media. They only predict, however, the location of the maximum salt concentration inside porous medium. Salt precipitation and invasion of the pore space by a gas (e.g. air) are not included in such models. Though continuum models provide valuable insights, they cannot take into account the eect of the pore-scale heterogeneity of a pore space, which is of utmost importance because, as already pointed out, drying and salt precipitation are pore-scale phenomena that are in uenced strongly by pore-scale heterogeneity. Thus, one needs to develop PN models in order to study the eect of pore-scale disorder on evaporation and salt precipitation. Once the phenomena are understood at the PN scale, one can up-scale the models to much larger porous media. The PN models have been used extensively for decades to study various phenomena in porous media (for reviews see Sahimi (2011); Blunt (2017)). Recent applications of the PNs include coating of surface of paper (Ghassemzadeh et al., 2001; Ghassemzadeh and Sahimi, 2004a), transport and separation of gaseous mixture in nanoporous membranes (Chen et al., 2008; Rajabbeigi et al., 2009; Mourhatch et al., 2010), ow permporometry for measuring the PSD of a mesoporous medium (Mourhatch et al., 2010), and simulation of multiphase uid ow in reservoir rock (Blunt et al., 2002; Piri and Blunt, 2005; Blunt, 2017). 68 Development of the PN model of evaporation-induced salt precipitation involves simulation of the drying process. Slow drying of porous media is an invasion percolation (IP) process (for a simple introduction to the IP see Ebrahimi (2010)), whereby the drying front invades a porous medium (Shaw, 1987) according to the least capillary pressure that a nonwetting uid (air) requires to enter a pore throat lled by saline water (see below). Using the IP process, many PNs have been developed for simulating drying (Daian and Saliba, 1991; Prat, 1995; Laurindo and Prat, 1998; Yiotis et al., 2001, 2006; ?). The performance of PN models of evaporation in porous media have been tested in the previous studies. In our modeling, we use the same approach to simulate the evaporation part of the process. Thus, the goal of the present study is developing a PN model for evaporation, brine transport and salt precipitation in porous media, and presenting the results of extensive PN simulations by which we have computed the various important properties of the phenomena. To our knowledge only B onhorst et al. (2016) presented a PN model of precipitation of solids in porous media, which is presumably applicable to salt precipitation. Their model diers, however, from what we present in this paper. In their work transport of the solids is by diusion, rather than by the convection- diusion process. It is known, however (Imdakm and Sahimi, 1987; Sahimi and Imdakm, 1991), that if the particles are relatively large, diusion is not an eective mode of transport, and the range of the particle sizes that their model is applicable for is not clear to us. At the same time, even if diusion is an eective mechanism of the particles, it could be so only for porous media whose pore throats are very small, so that convection may be ignored, or that the ow eld is very slow. Moreover, the important eect of the heterogeneity of porous media, which is due to the PSD and the correlations between the pore sizes is taken into account by our PN model, whereas this important eect was mostly ignored in the model by B onhorst et al. (2016). Indeed, a major aspect of salt precipitation at the pore scale is the eect of the heterogeneity. The rest of this paper is organized as follows. In the next section we describe the PN model. The details of the models of drying and salt transport and precipitation are presented in Sec. 69 3, while the computational algorithm is described in Sec. 4, followed by the presentation and discussion of the results in Sec. 5. The paper is summarized in Sec. 6. (a) (b) (c) (d) (e) (f) 50 100 150 200 Figure 5.1: Examples of pore networks with correlated throat sizes. The top and bottom rows show, respectively, the PNs with the Hurst exponents H = 0:7 and 0.3, and anisotropy x = y = (a) and (d) 1; (b) and (e) 1/4, and (c) and (f) 4. 5.2 Pore-network model This part was already explained in Chapter 4, but the sake of completeness we provide the description. We use a square network with a pore connectivity of 4. Extension of the model to 3D is straightforward but, relative to 2D networks, computationally expensive and will be reported elsewhere. Pressure and concentration gradients in pore bodies are ignored, the standard practice in many PN modelings. Thus, each pore body is characterized by a single value of liquid (or vapor) concentration and pressure, the calculation of which is described below. As a result, when we refer to the PSD of the PN, we mean the size distribution of the pore throats. The pore throats 70 are assumed to be cylindrical. So long as the eect of thin liquid lms on the internal surface of the pore throats is insignicant, the assumption of cylindrical throat is accurate, and has been used many times in the past. In the present paper in which we are interested in precipitation of salt, the thin lms are indeed unimportant. Figure 5.2: Illustration of the assumption of equilibrium vapor pressure in the central vapor pore, together with a neighboring throat containing liquid. An important aspect of the morphology of porous media that has not been considered in the simulation of drying is the correlations between the size of the pore throats. The existence of such correlations in the porosity distribution of eld-scale porous media was rst demonstrated by Hewett (1986a), followed by similar work by Neuman (1994) for the spatial distribution of the permeability of the same type of porous media. In fact, such correlations manifest themselves over 71 multitudes of scales (Sahimi and Mukhopadhyay, 1996; Mukhopadhyay and Sahimi, 2000; Dash- tian et al., 2015). Using experiments and PN simulations, Knackstedt et al. (2001) demonstrated that, even at the core scale, the computed properties of multiphase ow in porous media agree quantitatively with the experiments when such correlations are taken into account. To include the eect of the correlations in the size of the throats, we used a fractional Brownian motion (FBM) that has been shown to not only represent well the distribution of the permeabilities in eld-scale porous media (for a review see Sahimi (2011)), but also the correlated distribution of the throat sizes in cores (Knackstedt et al., 1998). The two-point correlation function C(r) of a FBM is given by, C(r) = C 1 r 2H , where C 1 = C(r = 1), and H is the Hurst exponent. The correlations are positive or persistent if H > 1=2, and negative or anti-persistent if H < 1=2. In the limit H = 1=2 the successive increments of a FBM are completely uncorrelated. An ecient method for generating a FBM array is through its power spectrum S(!), which for 2D systems is given by S(!) = a (! 2 x +! 2 y ) H+1 ; (5.1) where a is a constant, while ! x and ! y are the Fourier components in the x and y directions. To introduce anisotropy (layering) in the distribution, we follow Ansari-Rad et al. (2012) and generalize the power spectrum to S(!) = b ( x ! 2 x + y ! 2 y ) H+1 : (5.2) Here, x and y are the anisotropy parameters, and b is a constant. To generate anisotropy induced by layering with the layers being parallel to the x direction, we set x = y > 1, and similarly for the layers being parallel to the y direction. Figure 5.1 presents examples of the heterogeneities generated by the model. 72 5.3 Pore-network model and the governing equations As the evaporation proceeds, the vapor-brine interface advances in the pore space, and salt begins to precipitate, which we assume to be mainly NaCl. Thus, one must simulate both the drying and salt transport and precipitation. The computations consist of two distinct parts. One is simulation of the drying process by the IP algorithm. The PN's top row is assumed to be a layer whose throats are initially empty and through which a gas (such as air or CO 2 ) ows into the pore space, while the other sides are sealed. Brine ows in the pore space's throats, while the vapor diuses to the boundary layer and induces further evaporation. In the IP model a gas invades the brine-saturated porous medium throat-by-throat according to the lowest capillary pressure needed for entering a throat at the interface between the vapor and liquid. The capillary pressure for entering a throat of radius r ij between pores i and j at the liquid-vapor interface is given by P c;ij = 2 r ij ; (5.3) where is the surface tension. Each time one or a few throats are dried out by the IP algorithm, we switch to the second part of the simulation, namely, transport of brine and precipitation of salt, described below. At every step of the simulation we identify the state of the pores, which is determined by the lling state of the pore throats connected to them: if at least one neighboring throat does not contain brine, it is a vapor pore; otherwise, it is a liquid pore. A throat can be lled, fully or partially, with brine, or it can contain vapor that is leaving it through a path of pores and throats that are also lled by the vapor. A vapor pore can be either at equilibrium or at an unknown vapor pressure (to be determined). In reality, equilibrium vapor pressure prevails only at the menisci between brine and the vapor, as shown in Fig. 5.2, but we assume that a pore is at equilibrium vapor pressure if at least one of its neighboring throats still contains brine, even if 73 it is only partially so. Therefore, diusion of vapor between menisci of partially-lled throats and the adjoining vapor pore happens without any resistance. Our preliminary numerical simulation indicated that the assumption leads to slightly, only about 5 percent, overestimates of the drying rate, but simplies the computations greatly. Once the states of pores and throats are identied, the second stage of the simulation, namely, modeling of transport of brine throughout the PN and the resulting salt precipitation, begins. After one or a few steps of the IP that advance the drying front, we solve the governing transport equations. Assuming that the ow is slow enough, the liquid ow rate q (l) ij in a throat ij that connects pores i and j is given by the Hagen-Poiseuille equation q (l) ij = r 4 ij 8 (l) ` ij h P (l) i P (l) j i ; (5.4) where the liquid viscosity (l) is assumed to be constant, as we simulate an isothermal process. Here,P (l) i is the liquid pressure in porei, and` ij is the length of the throatij. Yiotis et al. (2001) showed that both advection and viscous eects aect the drying patterns and rates. Since there is no liquid accumulation in the pores, we must have X fijg (l) q (l) ij = 0; (5.5) with (l) being the brine's mass density. The sum is over the set of the throatsfijg that are connected to pore i. The vapor is transported by diusion and, thus, we use the equation for evaporation in the so-called Stefan tube (Bird et al., 2007) to calculate its mass ow rate in the throats. Ignoring the pressure in the vapor-lled throats, vapor mass conservation at the pores must hold, implying that, X fijg (v) q (v) ij = X fijg S ij D (v) ` ij PM (v) RT ln " PP (v) i PP (v) j # = 0; (5.6) 74 where D (v) is the vapor diusion coecient, M (v) its molecular weight, P is the total pressure in the vapor phase, and (v) is the vapor mass density. We used, (v) = 8:12 10 4 gr/cm 3 and D (v) 2:6 10 5 m 2 s 1 , both corresponding to 298 K. The vapor concentration at the vapor- liquid interface is equal to the equilibrium concentration, which is zero outside the boundary layer of the PN. Writing down Eqs. 5.4 - 5.6 for all the pores in the liquid and vapor phases results in two sets of equations that govern the pressures in the pores in the two phases, which are solved by the conjugate-gradient method. The pressure at the liquid-vapor interface is assumed equal to the equilibrium vapor pressure, P = 2339 Pa. At a meniscus between the brine and vapor the Kelvin equation describes the dependence of parameters on the vapor pressure: ln P (v) P = 2M (v) RTr ij (l) = P c M (v) RT (l) ; (5.7) where R is the gas constant, and Eq. 5.3 was used. Note that the dissolved salt increases the osmotic potential , which in turn decreases the saturation vapor pressure of the saline water (Battistelli et al., 1997; Jambhekar et al., 2015). In addition, capillary pressure might also lower the vapor pressure. can be easily computed via, =RT (l) ln(x (s) ), wherex (s) is the solubility- limit mole fraction of salt in water. One can include the eect of osmotic pressure by replacing P c in Eq. 5.7 by (P c + ) (Jambhekar et al., 2015), but the eect is believed to be small and, thus, we ignore it. Transport of salt through a porous medium is by convection and diusion (Shokri, 2014). The electroneutrality of brine implies that the concentration C of the salt at time t is governed by convection-diusion equation (CDE): @C @t =D e @ 2 C @x 2 v ij @C @z ; (5.8) 75 wherev ij is the mean liquid velocity in the throatij, andD e is the eective diusion coecient in the solution. Numerous studies have proposed various methods for including the eect of diusion in ow through PNs (Sahimi and Imdakm, 1991; Koplik et al., 1988; Sorbie and Cliord, 1991; Billinge, 2009b). One can simply take D e to be the molecular diusivity D m , which means that the throats must be narrow enough that the concentration at any given axial position z within the throats is single-valued. Doing so simplies the computations. Alternatively, one may take D e = D L , where D L is the Taylor-Aris dispersion coecient (Taylor, 1953; Aris, 1956). The choice D e =D L , which we make in this paper, entails imposing some constraints. According to Mehmani and Balho (2015) D L can be used in a PN simulation if the aspect ratio A = ` ij =r ij of the throats is 10 or larger. Note also that it has often been assumed and there is experimental evidence (Sahimi, 2011) that, ` ij / 1=r ij . Given that the throats' radii are typically small, ` ij are typically large, which is what we also use. Thus, the condition for using the Taylor-Aris dispersion is satised, at least approximately. In a throat ij D L is given by D L =D m + r 2 ij v 2 ij 48D m : (5.9) One can show that the salt concentration near the liquid-vapor interface is maximum (Huinink et al., 2002; Guglielmini et al., 2008). Thus, salt precipitates mostly at the interface throats. In the liquid-lled pores i the net accumulation of salt at time step n is given by X fijg S ij J ij = X fijg v ij C n+1 i C n i t ; (5.10) whereJ ij =v ij C i D L @C i =@z =v ij C i D L (C j C i )=` ij is the mass ux in the throatij over a time step t. If a throat contains partially the vapor that can escape (i.e., it is connected to other pores and throats that are lled with vapor), then, mass balance will take the mass ow rate of the vapor in that throat into account. Writing the brine mass balance for each pore i yields a 76 set of equations for the concentrations C i of the salt that we also solve by the conjugate-gradient method. (a) (b) (c) (d) Figure 5.3: The phase distributions in an isotropic pore network with the Hurst exponentH = 0:7 at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.8, and (d) 0.95. Blue, white, and yellow indicate, respectively, the throats lled with brine, the emptied throats, and those that are partially or completely blocked. The top row of the network is open. 77 5.4 Computational algorithm Given the above formulations, the computational algorithm is as follows. (i) The pressure distribution in the vapor phase is computed. (ii) The mass evaporation rate in each throat and that from the boundary layer are calculated. (iii) The pressure distribution in the liquid phase and, therefore, the mean ow velocity and the local dispersion coecient D L in each throat are computed. (iv) The time step for advancing the vapor-liquid interface by the IP algorithm is selected, following Hekmatzadeh et al. (2016). There are two time scales in the problem. One is associated with evaporation, while the second one is the time scale over which one solves the CDE. We use the smaller of the two. Thus, the evaporated mass from all the throats in the drying front containing liquid throats, which leaves the PN, is computed. The throat at the interface with the lowest capillary pressure is one with the largest radius. The time step for emptying that throat is given by t = (l) V ij P fijg2L D M ij ; (5.11) whereV ij is the volume of the throatij,M ij is the mass rate of evaporation, and the sum is over all the throatsfijg that are at the drying frontL D . (v) The CDE is solved in the brine part of the PN in order to compute the salt concentration in the throats. (vi) As drying continues, salt concentration in the liquid phase increases. As already pointed out, the solution of the CDE for salt transport in porous media indicates that salt concentration at the drying front is maximum (Huinink et al., 2002; Guglielmini et al., 2008). Thus, salt concen- trations in the throats next to the drying frontL D are compared with a threshold concentration for its precipitation. If they exceed the threshold, which we take it to be the solubility limit of salt in water at room temperature (36 percent in mass), then, salt precipitates in such throats. 78 Note also that preliminary simulations indicated that the salt content of practically no throat away from the interface exceeds the saturation limit. Then, the amount of precipitated salt in the throats is calculated as the dierence between the current salt concentration in the throats and the saturation concentration. The precipitated salt forms a porous packing of its crystals and, therefore, it would take some time to plug a throat completely. Thus, the volume of the precip- itated salt in each throat is calculated and subtracted from the throat's volume, and its radius is updated, assuming that salt precipitates uniformly. Depending on its size and the amount of precipitated salt, the throat can either become partially or completely plugged. (vii) The IP algorithm advances further the drying front in the PN by one or a few throats, and the procedure is repeated until the evaporation rate from the boundary layer is negligible. We note that recent work (Desarnaud et al., 2014) indicated that the precipitation threshold may be as high as 1.5 times the solubility limit. But, this would not make any dierence to the PN model, as the threshold may be set at any value. The PN is initially saturated with brine at 30% w. NaCl. All the results presented below represent averages over at least 10 realizations of the PN. As discussed in Section 2, a FBM algorithm was used to generate 2D PNs with two values of the Hurst exponents H and three values of x = y . Table 1 presents the numerical values of H, x = y , and the mean throat size (radius) for each case. Table 5.1: Pore size distribution properties of dierent cases. H x = y Initial mean throat size (m) Final throat radius (m) 0.3 1 70.5 53.2 0.3 4 74.1 55.9 0.3 1/4 67.7 51.1 0.7 1 104.5 74.4 0.7 4 110.1 88.4 0.7 1/4 101.1 66.7 79 (a) (b) (c) (d) Figure 5.4: The phase distributions in a pore network with the Hurst exponent H = 0:7 and anisotropy x = y = 1=4, at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.8, and (d) 0.95. Colors indicate the same as in Fig. 5.3. The top row of the network is open. 5.5 Results and discussion Let us rst inspect the eect of the PSD distribution and the correlations on the spatial distribution of the pore throats that are either lled with brine, or with the vapor, or are partially or completely plugged by the precipitated salt. Shown in Figures 5.3 is the evolution of the phase distribution in an isotropic PN ( x = y ) with the Hurst exponent H = 0:7, implying positive 80 correlations between the sizes of the throats. In particular, Fig. 5.3(a) corresponds to the stage in which the drying rate is essentially constant (this will be discussed quantitatively shortly). Close to the open end, the throats empty (dry out) at a faster rate than those deeper in the pore space. Since the correlations in the PSD are positive, the throats with similar sizes are clustered together, and the drying patterns mimic the underlying connectivity of the throats. This is manifested by the large clusters that are clearly visible in Figure 5.3. Moreover, due to the isotropy of the pore space, the drying front advances in both directions, shown in Figures 5.3(b) and 5.3(c). Next, consider anisotropic (stratied) porous media. Figures 5.4 present the evolution of the spatial phase distribution with x = y = 1=4, so that the throats with similar sizes tend to extend orthogonal to the open end of PN, as also indicated by Figure 5.1(b). Thus, preferential evaporation in the throats with similar sizes that are in they direction enhances the advancement of the drying front deep into large throats, while the liquid phase stays connected to the surface. AS a result, the drying front stretches out to the bottom of the PN faster than the isotropic porous media shown in Fig. 5.3; see Figs. 5.4(b) and 4.4(c). In the case of x = y = 4, the strata are parallel to the open end of the PN. Thus, the throats with similar sizes cluster together and extend parallel to the open end. As depicted in Figs. 5.5(b) and 4.5(c), the drying patterns mimic the direction of the strata. Altogether, Figs. 5.3-5.5 demonstrate a universal eect in that they show the strong eect of positive and extended correlations on the drying patterns, and the fact the clusters of dried out throats form based on the underlying correlations. Salt precipitation does not, however, follow a universal pattern in all the three cases. But, looking at the throats that are partially or completely plugged in Figs. 5.3 - 5.5, it can be seen that salt tends to precipitate in neighboring throats that are parallel to the open end of the PN. As reviewed comprehensively by Sahimi (2011), extensive data indicate that if the correlations in porous media are of the FBM type, then the corresponding Hurst exponent H is often less than 0.5, implying negative correlations between the throats' sizes and, therefore, a much more heterogeneous pore space than when H > 0:5. Therefore, we also studied the cases in which 81 (a) (b) (c) (d) Figure 5.5: The phase distributions in a pore network with the Hurst exponent H = 0:7 and anistropy x = y = 4, at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.8, and (d) 0.95. Colors indicate the same as in Fig. 5.3. The top row of the network is open. H = 0:3, with the corresponding means of the PSD presented in Table 1. When the correlations are negative, a larger throat size is more likely to have small throats as its neighbors, and vise versa. The resulting evolution of the phase distributions are shown in Figs. 5.6 - 5.8 that correspond, respectively, to x = y = 1, 1/4 and 4. Comparison of these gures with their counterparts when H = 0:7 demonstrates the eect of correlations and PSD. Consider, for example, Figs. 5.6, and compare them with those in Fig. 5.3. Whereas clustering of the throats type - liquid lled, vapor 82 (a) (b) (c) (d) Figure 5.6: The phase distributions in an isotropic pore network with the Hurst exponentH = 0:3 ( x = y = 1), at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.8, and (d) 0.95. Colors indicate the same as in Fig. 5.3. The top row of the network is open. lled, or partially or completely plugged - is clearly visible in Figs. 5.3 (and likewise in Figs. 5.4 and 5.5), the negative correlations essentially destroy such clustering eects. Figure 5.9 presents the phase distribution in a PN in which there is no correlation between the size of the throats, but with the same range of throat sizes as before. The dierences between the evolution of the phase distribution in Fig. 5.9 with those shown in Figs. 5.3-5.8 for H = 0:7 are evident. In particular, the number of clustered brine-lled throats, as well as isolated liquid 83 (a) (b) (c) (d) Figure 5.7: The phase distributions in a pore network with the Hurst exponent H = 0:3 and anisotropy x = y = 1=4, at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.8, and (d) 0.95. Colors indicate the same as in Fig. 5.3. The top row of the network is open. lled throats, is much higher in the case of random PN than when the morphology of the pore space is correlated. We now make a quantitative comparison between the various cases whose phase distribution evolution has been described and discussed so far. An important property is the brine saturation S and its evolution with the time. Figures 5.10 and 5.11 show the transversely-averaged (averaged over the x direction in Figs. 5.3-5.9) saturations in the same six PNs of Figs. 5.3-5.9, computed 84 at four times with the colors showing the range of variations over multiple realizations and the black curves indicating the averages over all the realizations. As one may expect based on Figs. 5.3-5.8, the correlations aect the saturation distribution strongly. (a) (b) (c) (d) Figure 5.8: The phase distributions in a pore network with the Hurst exponent H = 0:3 and anisotropy x = y = 4, at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.8, and (d) 0.95. Colors indicate the same as in Fig. 5.3. The top row of the network is open. At any given time there exist three distinct regions in Fig. 5.10: one is associated with zero saturation over a range of the distance from the open end of the porous medium. The second region represents a transition zone in which 0 < S < 1, while in the third region the medium 85 (a) (b) (c) (d) Figure 5.9: The phase distributions in an uncorrelated PN at time steps when vapor saturation is (a) 0.2; (b) 0.4; (c) 0.8, and (d) 0.95. Colors indicate the same as in Fig. 5.3. The top row of the network is open. is fully saturated, S = 1. Note the stronger eect of positive correlations with H = 0:7 on the saturation proles, indicated by the fact that for negative correlations (H = 0:3) the range of variations in the saturations over multiple realizations is narrower. The eect is most pronounced when the porous medium is isotropic, as the advancing drying interface does not have any preferred direction to move. Another notable dierence is the depth at which the transition zone ends, i.e. where S = 1. For example, consider the time at which the liquid saturation is 40% (dark blue 86 color). The depth at which the transition zone ends is 40 (in units of the throats' length) for H = 0:7 and x = y = 1=4, whereas it is 30 for x = y = 4. Similar patterns exist with negative correlations, as shown in Fig. 5.10. The ranking of the depth at which the transition zone ends is x = y = 4< x = y = 1< x = y = 1=4. Under isothermal condition the evaporation rate depends on such driving forces as the capillary and viscous forces (Shokri and Or, 2011; Lehmann and Or, 2009). The capillary pressure depends, of course, on the PSD. Scherer (1990b) discussed the dierence between the capillary pressures of large and small throats that result in ow of the liquid to the surface of a porous medium, keeping the evaporation close to the early stages, which is called capillary pumping. The evaporation rate at the open side of the porous medium connected to the boundary pores is given by _ m (v) = X i2 boundary X j2fig S ij D (v) L PM RT ln " PP (v) i PP (v) j # : (5.12) Figure 5.12 compares the evolution of the evaporation rates, usually called the drying curve, which were computed for the six types of the PNs of Figs. 5.3-5.8, averaged over multiple realizations of the PNs and rescaled by the number of pores in the open boundary of the PNs and the initial drying rate. The qualitative trends of Figure 5.12 are in complete agreement with the recent experimental data of Shokri-Kuehni et al. (2017). The results in Figure 5.12 manifest the three main stages of drying and precipitation: (i) the initial constant evaporation; (ii) the relatively fast intermediate falling rate, and (iii) the residual slow rate of evaporation. The rst stage is essentially independent of the morphology of porous media and is controlled by the external conditions. It is visible in Fig. 5.12 as a very short period in the initial time steps. Evaporation in the intermediate stage is dictated by the ability of the porous medium to supply the evaporation zone with brine, which in turn depends strongly on the morphology of the pore space, and in particular the PSD, the correlations, and the pore connectivity, and lasts much longer than the initial stage. Comparing all the cases, the 87 intermediate stage with x = y = 4 and H = 0:7 has the highest evaporation rate, as shown in Fig. 5.12. As we discussed earlier when we described Figure 5.4, in this case the liquid stays connected to the surface of the pore space and, therefore, evaporation can happen in throats near the surface. Evaporation from the open end of the pore space dominates drying from the throats deeper into the pore space. Figure 5.10: Dynamic evolution of the transversely-averaged brine saturations in six types of pore networks. Black curves show the average saturations; cyan, blue and pink shadings represent, respectively, the times at which vapor saturation is 0.2, 0.4, and 0.8, while red indicates the time at which vapor saturation is 0.95. The colored areas indicate the variations over multiple realizations. To better understand the results, consider ow of brine through the PN between two points separated by several throats. If the throats' sizes between the two points are positively correlated, there will be smaller resistant against the ow. In the PN in which the strata are parallel to the y direction - Figs. 5.12(b) and 5.12(e) - the falling rate is slower and, thus, this stage lasts until almost the end of the process, followed by fast declining. The third stage emerges when the surface zone has completely dried out, so that further ow of liquid through it eectively ceases. There are some studies (e.g., Or et al. (2013)) that dene the drying in two stages. Depending on the porous media and the evaporation condition, there might be other phenomena, such as 88 Figure 5.11: Dynamic evolution of the transversely-averaged saturations in the uncorrelated pore network. Black curves show the average saturations; cyan, blue and pink shadings represent, respectively, the times at which vapor saturation is 0.2, 0.4, and 0.8, while red indicates the time at which vapor saturation is 0.95. The colored areas indicated the variations over multiple realizations. sorption of water on the throats' surface (Bakhshian and Sahimi, 2017). We have not included such secondary eects in the simulation, because our primary interest is salt precipitation and the eect of the pore space heterogeneity on the phenomenon. Depending on the morphology porous medium and the chemical composition of its pores' surface, sorption of water might have some eect on the last stage of drying. On the other hand, through comprehensive experiments, Bergstad and Shokri (2016), who studied evaporation of brine in sand grains, reported drying curves and the images of the phase distributions that are comparable with what we are reporting in the present paper, even though they utilized porous media with various structures. Another important characteristic of drying in porous media is the length L D and shape of its front. Figure 5.13 presents the time evolution of L D (t) in the six types of porous media of Figs. 5.3-5.8. It is clear that the morphology of the porous medium, and in particular its correlation structure, have a signicant eect onL D (t). For example, the average drying length is larger when the strata are parallel to the y direction. Results shown in Fig. 5.13 indicates that at practically 89 Figure 5.12: Dynamic evolution of the average drying rates (black curves), averaged over 10 realizations (shadings). The colored areas indicate variations over multiple realizations. Figure 5.13: Dynamic evolution of the average length of the drying front. Shaded areas show the extent of the variations over ten realizations. all times, L D (t) is larger for x = y = 1=4; that is, when the strata are parallel to the y direction. The physical eect of this type of morphology was already discussed. 90 One may also study the eect of the pore space morphology on salt precipitation by considering the distribution of the throats' sizes that are partially or completely plugged by salt precipitation. Figure 5.14 depicts such distributions in the same six types of porous media that are shown in Figs. 5.3-5.8. Although the shapes of the distributions are more or less similar before and after the precipitation, there are also quantitative dierences between the two. In Table 1 we also present the average throat sizes before and after salt precipitation. The simulations indicate that the porosity reduction due to precipitation is around 15%, which is close to what have been reported for injection of CO 2 in oil elds, 20 percent. The results shown in Fig. 5.14 should also be compared with those shown in Fig. 5.15, which presents the PSDs before and after salt precipitations in a completely random (uncorrelated) porous medium. Miri et al. (2015a) studied salt precipitation during CO 2 injection, using a fabricated 2D PN, initially saturated with brine. Similar to our PN the connectivity of their PN was 4, but their boundary conditions were dierent. They injected CO 2 from one side and produced brine from the opposite side. The rate of injection was very low, with most of the water leaving the PN through evaporation. Despite the dierences in the boundary conditions, Miri et al. (2015a) reported uniform salt precipitation in their experiments, similar to what we report here. Experimental results reported by Kim et al. (2013) also share similarities with our results. They used fabricated 3D PNs to study drying and salt precipitation, and reported formation of uniform large salt crystals away from the vapor-liquid interface inside the porous medium. 5.6 Summary This paper presented a pore-network model of isothermal evaporation of brine owing in a porous medium and the resulting salt precipitation. Vapor diusion in the gas phase, ow of brine, transport of salt and its precipitation in the pore space were all accounted for by the model. The dynamic evolution of the brine saturation, the evaporation rate, the length of the drying front, 91 Figure 5.14: Size distribution of the partially or completely plugged throats before and after evaporation. Figure 5.15: Size distribution of the partially or completely plugged throats before and after evaporation in the uncorrelated pore network. and the size distribution of the plugged or partially-plugged throats were computed and studied. Furthermore, we demonstrated the strong eect of the heterogeneity and correlations between the size of the throats on the evaporation rate and the various drying stages. Despite extensive experimentation, the knowledge about factors that contribute to the phe- nomenon is still incomplete (Miri and Hellevang, 2016). The PN model provides a framework to 92 study salt precipitation under more complex conditions than what we have considered in this pa- per. In particular, the lithology of geological formations possesses complex chemistry that aects the wettability, which in turn in uences not only salt precipitation, but also the drying behavior (Or et al., 2013). Such eects will be included in the PN model, which we hope to report in the near future. 93 Chapter 6 Non-universality of the Archie exponent due to multifractality of resistivity well logs 6.1 Introduction A fundamental backbone of petrophysics is Archie's law (Archie, 1942), which states that the saturationS w of water (brine) in a porous formation with porosity and tortuosity is given by S n w = R w R t m : (6.1) Here, m is called the Archie or cementation exponent, n is a constant, R w is the resistivity of brine, and R t is the overall resistivity of the brine-saturated rock. Originally proposed as an empirical correlation, Archie's law has been derived theoretically using various models and approximations (see, for example, Sen1981 et al. (1981); Yonezawa and Cohen (1983)). Although it is widely used [see Sahimi (2011), for a comprehensive review], it does have some limitations, such as its inability to account for the rock matrix geometry (Haro, 2010). Many granular porous media with well-connected pore space have an Archie exponent, 1:5< m < 2:5 [(Glover, 2009). In carbonate reservoirs, on the other hand, the pores are not well connected. While, in general, m increases as the interconnectivity of the pores diminishes, the 94 precise value of m, or the range of its variations, for carbonate reservoirs has remained contro- versial. In fact, the range of variations of m in such reservoirs is considerable and, hence, the use of an average value, even if it is based on extensive data, may lead to considerable inaccuracy in the saturations (Focke and Munn, 1987). Several studies have presented evidence that a detailed description of rock type is needed for obtaining accurate estimates ofm for carbonate formations. Typical rock types in such formations are intergranular, intercrystalline, vuggy, microporous (dis- connected pores), and fractured zones. Another factor that aects estimates of m is the presence of shale and clay minerals in the formation that lead to a higher and unrealistic estimate of the porosity with a lower electrical resistivity and the corresponding exponent m. Despite extensive use of Archie's law and signicant research, the question of the numerical value of m and its de- pendence on the rock morphology has remained unresolved. The issue is important as the water, oil, and/or gas saturations of porous formations calculated by Eq. 6.1 are sensitive to changes in m. The Archie exponent m may be estimated by the analysis of core data (Bassiouni, 1993) or well logs but, as we show in this chapter, they lead to completely dierent estimates. Moreover, despite evidence for non-universality of m, several studies have claimed that m is universal, at least for certain classes of rock. For example, the solution of the problem of the resistivity of a porous medium in which the pore space is modeled as a packing of spheres yields m = 3=2 (Sen1981 et al., 1981). A recent paper (Ghanbarian et al., 2014); see also Ewing and Hunt (2006), suggested that percolation theory can accurately predict experimental data for the resistivity of porous media and, thus, the Archie exponent m is linked to percolation theory. Another paper (Wei et al., 2015) used a model of fractal porous media to develop an expression for the resistivity. In this Letter we present a careful analysis of extensive porosity and resistivity logs for six oil and natural gas reservoirs in Iran, based on which we estimate the Archie exponent m. We show that m is non-universal, and explain the cause of non-universality, as well as the apparent universality of m reported recently by Ghanbarian et al. (2014). 95 0 0.2 0.4 NPHI 0 50 100 150 GR 10 0 10 2 10 4 Resistivity Period 2670.81 2708.91 2747.01 2785.11 2823.21 4 8 16 32 64 128 256 0 0.5 1 Depth (m) Period 2670.81 2708.91 2747.01 2785.11 2823.21 4 8 16 32 64 128 256 0 0.5 1 NPHI−GR Resistivity−GR Figure 6.1: A portion of a sample well log used in this study. The data belong to the TBK-1. Top: typical neutron porosity (NPHI), gamma-ray (GR) and resistivity logs. Bottom: wavelet transform coherence of gamma ray and resistivity, and gamma ray and neutron porosity. The black curve indicates 5% signicance level. The coherence ranges from dark red (nearly 1) to dark blue (nearly 0). The red colored regions inside the thick dark counters show the aected data (porosity or resistivity) by clay minerals (indicated by gamma ray) at various scales (periods). 6.2 The Reservoirs and the Data We used 36 well logs belonging to six reservoirs in Iran, each of which has 2000 - 12000 data points, far more data than any of the previous studies. Sample data are presented in Figure 6.1. In the Supplementary Information we provide extensive details on the data and the calculations described below. Typically, the data analyzed in the past belonged to outcrops, which may not have the same properties as deeply buried formations. The well logs belong to Tabnak (TBK), 96 Ahwaz (AZ) and Maroon elds in Iran. The TBK eld produces natural gas, the AZ eld is a major oil reservoir, while the Maroon eld produces both natural gas and oil. The uid content and structure and morphology of the three elds are dierent. Well logs from three other reservoirs that have the same characteristics were also analyzed. Most of the data belong to the Asmari formation, which has an average thickness of about 400 m, and represents the most prevalent rock type among the Iranian oil and natural gas elds. The reservoirs we study may be divided into two main groups. One group includes the AZ and Maroon elds, the bulk of which consists of the Asmari formation. The Oligocene-Miocene Asmari formation consists predominantly of carbonates that are interbedded with sandstones, and forms one of the principal reservoir types in the Iranian oil and natural gas elds. Geological studies indicated that carbonate deposition was initiated in a shallow-marine environment and continued through shallowing upward conditions, which led to a more restricted lagoonal environment. The geological interpretation and spatial distribution of the sandstone layers indicate that they may be of deltaic origin, and provenanced from the west and southwest. The limestone facies range from wackestone to bioclastic, pelletoidal, and in-part oolithic packstonegrainstone that have been more or less dolomitized. The porosity types include interparticle, intercrystalline, moldic and vuggy. Their permeability is low, but is enhanced by fractures. The second group of the reservoirs include the TBK, Aghar and two other reservoirs. Their hydrocarbon-containing layers are Dashtak, Kangan and upper Dalan. The formations are car- bonates of Permian and Triassic age type. In this group, from exposed outcrops to the target, the formations consist mainly of limestone and dolomite, as well as secondary shale with interbeds of gypsum. The upper formations are interbedded with chert and yellow iron ore. Laanh and Kazhidumi are consisted primarily of shale. The formations below Gadvan are dominated by lime- stone and dolomite with beds of blackish gray, grayish-green and mauve shale. Some limestone is much argillaceous. (All the names refer to local regions where the reservoirs are located.) 97 10 −3 10 −2 10 −1 10 0 10 0 10 1 10 2 10 3 10 4 10 5 10 6 10 7 Effective Porosity (φ − φ t ) Formation Factor (F) MN−296 MN−297 AZ−273 φ x = 1.00 φ x = 0.40 φ x = 0.75 Figure 6.2: The formation factors as a function of eective porosity ( t ). Solid lines are the predictions based on percolation-eective-medium approximation with various crossover porosity x . The data are from three wells used, two from Maroon (MN296 and MN297) and one from Ahvaz (AZ273) oil elds 6.3 Analysis of the Data Each well passes through several layers. Typically, the data near the boundaries between the layers were ignored. One goal of this study is understanding the eect of the clay contents on the Archie exponent m, because porosity logs are strongly in uenced by the clay content, and not correcting for the eects of clays would provide unrealistic estimates of the porosity (Waxman et al., 1968). Thus, we rst identify the "clean" sectors along the wells and their corresponding depth, resistivities and porosities, where the clay content is not high. Since in geological formations 98 clays represent a source of radiation as they contain radioactive elements, they exhibit, relative to the clean formations, more intense radioactivity in the gamma-ray (GR) logs. Thus, one can identify the data associated with the clay contents by studying the mutual or cross correlations between the GR and porosity logs. The zone in which such correlations are strong are associated with the clay contents. Well logs are also typically nonstationary series, giving rise to spurious correlations (Mehrabi et al., 1997; Jafari et al., 2011), if they are not analyzed properly. Thus, calculation of the correlation functions through the logs' spectral density and their associated properties is not accurate. Therefore, we use the cross-wavelet transform (XWT) (Grinsted et al., 2004), the accuracy of which has been tested on various nonstationary logs, and shown to provide reliable results. 6.3.1 Identifying the Clay Sectors The discrete WT of a well log (or series) x(i) given by W x (s;k) = z p s N1 X i=0 x(i) [z(ik)=s]; (6.2) wheres andk are, respectively, the scale expansion and translational parameters. The distance between two neighboring points in the well log is z, while (z) is the mother wavelet. We use the Morlet wavelet, (z) = 1=4 exp(i! 0 z 2 =2): (6.3) Given two well logs (or series) x(i) and y(i), their XWT is dened by, W xy (s;n) =W x (s;n)W y (s;n); (6.4) 99 with indicating the complex conjugate of W y . The normalized XWT, which is called the wavelet transform coherence (WTC), is used in this study to dierentiate the clay contents from the rest of the data (Grinsted et al., 2004; Dashtian and Sahimi, 2014). It is calculated in a manner similar to the traditional correlation coecient: R xy (s;n) = jS [W xy (s;n)=s]j S [jW x (s;n)j=s] 1=2 S [jW y (s;n)j=s] 1=2 : (6.5) Thus, R xy can be thought of as a localized correlation coecient in the wavelet space. Since the wavelets' support is compact and nite, R(x;y) is localized. The operatorS is given by S(W ) =S ws S ss [W x (s;n)]: (6.6) Here,S ws is the smoothing operator along the wavelet scales, whereasS ss represents the same for the spatial scales. For the wavelet that we use, S ws (W )j s = h W x (s;n)a t 2 =2s 2 1 i s ; (6.7) S ss (W )j s = [W x (s;n)a 2 (0:6s)] s : (6.8) Here, is the gate or rectangular function, while a 1 and a 2 are two constants. 6.3.2 Analysis of Correlations in the Data It has been demonstrated that the porosity (Hewett, 1986b) and other types of well logs (Sahimi and Tajer, 2005; Sahimi, 2011) exhibit correlations whose extent is as large as the system under study, and follow nonstationary stochastic processes, such as the fractional Brownian motion (FBM). Thus, we analyze the well logs in order to understand the nature of such correlations and 100 their eect on the Archie exponent m by using the multifractal detrended uctuation analysis (MF-DFA). Consider a nonstationary log y(j) with j = 1;;N. We rst set up a new series whose mean and variance are, respectively, zero and one, and is given by y(j) =y(j)hyi: (6.9) and then integrate it to determine a new series Y (j): Y (j) = j k=1 y(k): (6.10) Next, the series Y (j) is partitioned into L = [N=s] equal nonoverlapping segments, where s is the segments' size. Each segment contains some local trends, such as unphysical noise. In the MF-DFA method the trends are usually tted to a linear funtion, but to ensure better accuracy we t the local trends of Y (j) in segment , Y (j), to quadratic polynomials. The variance of the segments is given by, (s;) = 1 s s k=1 Y [( 1)s +k]Y (k) 2 : (6.11) The detrended variance of order q is then given by M(q;s) = 1 2L 2L =1 (s;) q=2 1=q ; (6.12) if q6= 0, and M(0;s) = exp 1 4L 2L =1 ln [(s;)] (6.13) q> 0 and q< 0 weight, respectively, large and small uctuations. The scaling of the uctua- tions is given by Kantelhardt et al. (2002, 2003), 101 3 3.5 4 4.5 5 5.5 6 −2 −1.5 −1 −0.5 0 0.5 1 1.5 ln(s) ln[M(2,s)] AZ−273 D MN−297 TBK−260 Figure 6.3: Logarithmic plot ofM(q = 2;s) versus scales for four porosity logs from four dierent reservoirs. The slope of the line is the Hurst exponent H . M(q;s)s h(q) : (6.14) As is well-known, the nature of long-range correlation is characterized by the Hurst exponent H. Positive or persistent correlations correspond to 1=2 < H < 1, whereas 0 < H < 1=2 is a characteristic of negative or anti-persistent correlations. In the detrended uctuation analysis H is equal to h(q = 2). 6.4 Results Figure 6.1 presents an example of the WTC between a gamma ray log and neutron porosity log, as well as between the former and a resistivity log. Regions with high WTC are marked by dark red. Thus, we also analyzed the data that correspond to low WTC (blue in Fig. 6.1) in order to estimate the cementation exponent m, and refer to them as the \clean" data. We do so because the WTC is similar to the covariance of two series, but is much more accurate 102 3 3.5 4 4.5 5 5.5 6 6 7 8 9 10 ln(s) ln[M(2,s)] AZ−273 D MN−297 TBK−260 Figure 6.4: Logarithmic plot of M(q = 2;s) versus scale s for four resistivity logs from four dierent reservoirs. The slope of the line is the Hurst exponent H R . No data from the logs were removed and robust as a descriptor of the cross correlations and, thus, much better suited for the type of noisy, non-stationary data that we analyze. Assuming that for a given porosity the tortuosity is constant, we rst rewrite Eq. 6.1 as S w =I 1=n ; (6.15) I = R t FR w ; (6.16) F/ m ; (6.17) where I is resistivity index, and F is the formation (resistivity) factor. Therefore, log(R t ) =m log() + log(R w ) + log(I): (6.18) 103 Equation 6.18 indicates that if logR t is plotted versus log, one should obtain a straight line with slopem. It is not necessary to know R w to estimate m. While Eqs. 6.15 - 6.18 can be used for both core and well-log data, estimating m from the latter type of data is dicult, even if we set aside the eect of rock type and the presence of clay minerals. It has been suggested (Ghanbarian et al., 2014) that F can be predicted by a combination of percolation theory (Sahimi, 2011) and the eective-medium approximation (EMA): F = 8 > > < > > : (1 c )( x c ) ( c ) 2 c << x (1 c ) ( c ) x << 1; (6.19) where c is the critical porosity so that for < c the pores do not form a sample-spanning cluster, and x is the porosity at whichF crosses over from a percolation description to one by the EMA. The rst equation in 6.19 is the power law provided by percolation theory, hence predicting that the Archie exponent is m = 2, while the second equation is the EMA. Core data were used by (Ghanbarian et al., 2014) to test their theory. They reported that x ranges from 0.45 to 1. We used three well logs to check the accuracy of Eqs. 6.19. For each well we plotted F versus and varied x to identify its precise value. Following Hunt (2004), we set the critical (threshold) porosity t to be 10% of total porosity. Sample results are shown in Fig. 6.2. Clearly, Eqs. 6.19 do not accurately predict or t our data. Note that, generally speaking, we do not dispute the applicability of the EMA or percolation theory. What we are studying is the non-universality of the Archie exponent m and its cause. Note also that Hunt (2000) showed that any relevance of the statistics of percolation cluster to the variability of hydraulic conductivity would automatically generate a proportionality of the range of a semivariogram at the scale of the system. Hunt's analysis did not depend on the specic form of the conductivity and, thus it is equally applicable to the electrical conductivity that we study here. 104 We then used the the MF-DFA method to analyze the extended correlations in the logs, analyzing both the porosity logs and the corresponding resistivity logs. Two dierent cases were considered for the latter. In one case the entire resistivity log for each well was analyzed, while in the second case the clean data, as dened earlier, were studied. If the scales s are very large or small, the computed M(q;s) is unreliable. Thus, in accordance with the standard practice in using the MF-DFA method, the computed M(q;s) for both regions was neglected. We then computed the Hurst exponent H = h(q = 2) for both the resistivity and porosity logs of all the wells. Figure 6.3 presents the typical results for M(q;s) for the porosity logs with q = 2, indicating that power law (14) is followed accurately. In the Supplementary Information we present such plots for all the data that we analyzed. Analyzing all the porosity logs, we found that the Hurst exponent H =h (q = 2) associated with the porosity logs is 0:73H 0:94; porosity logs Such values ofH represent strong long-range persistent correlations in the porosities along the wells. Figure 6.4 presents the typical results for M(2;s) for the resistivity logs, when all the data were analyzed. Power law (14) is roughly followed. When all the resistivity logs for all the well were tted to Archie's law, we found that, 0:21m 1:75; all the data which are much lower than the previous reported values of m for other types of reservoirs. The reason, we believe, is the presence of clays that aect the resistivities signicantly. Archie 105 exponents less than 2 were also reported by Porter and Carothers (1971), although they did not discuss the reason for such values of m. Figure 6.5 presents the typical results forM(2;s) for the resistivity logs in which the data that correspond to the clay contents, identied by the WTC, were set aside. Once again, power law (14) is followed more accurately than in Fig. 6.4. Fitting the clean resistivity data to Archie's law, we nd that 1:6m 2:5; clean data: Moreover, we found for all the resistivity data that the Hurst exponent H R is, 0:52H R 0:90; resistivity logs Such values ofH R represent persistent long-range correlations in the distribution of the resis- tivities, and are completely consistent with the porosity logs. It is therefore clear that the Archie exponent m is nonuniversal. 6.5 Discussions Well logs represent stochastic series. Their stochasticity is due to the sedimentation, com- paction, tectonic motion, etc. As rst demonstrated by Hewett (1986b), stochastic variations of the porosity along wells can be described by self-ane distributions, e.g. fractional Brown- ian motion. This implies that the spatial distribution of the porosity and other rock properties contain long-range correlatns. It has been shown (Sahimi and Mukhopadhyay, 1996; Knackstedt et al., 2000) that if the spatial distributions of the porosity, permeability and other properties of large-scale porous media follow self-ane stochastic functions, then, their ow, transport and 106 3 3.5 4 4.5 5 5.5 6 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 ln(s) ln[M(2,s)] AZ−273 D MN−297 TBK−260 Figure 6.5: Logarithmic plot of M(q = 2;s) versus scale s for four resistivity logs of clean forma- tions (those in which the data indicating high correlations with the clay contents removed) from four dierent reservoirs. The slope of the line is the Hurst exponent H R . morphological properties follow power laws that are characterized by non-universal exponents. Our analysis demonstrates this convincingly for the Archie exponent that is associated with the resistivity. Ghanbarian et al. (2015) discussed theoretically the possible causes of the non-universality of the exponent m in Eq. 6.19, although the type of long-range correlation and multifractality described here was not discussed. The apparent universality ofm reported by (Ghanbarian et al., 2014) is due to their use of the resistivities of samples of porous media at laboratory scales. Only when the resistivity logs are analyzed over several hundred meters, do long-range correlations manifest themselves. Indeed, extensive computer simulations (Knackstedt et al., 2000) indicated that if a cuto length scale ` c is introduced such that for length scales `<` c the local properties - porosity, resistivity, etc. - are correlated, but become completely uncorrelated for length scales `` c , m and other exponents associated with various properties of porous media are universal. 107 Wei et al. (2015) proposed an equation linking the Archie exponent m with the morphological characteriostics of porous media. They expressed m in terms of a "tortuosity fractal dimension" D T , the pore surface fractal dimension D f , and the Euclidean dimension D E of the pore space: m = 1 + D T 1 D E D f : (6.20) The tortuosity fractal dimension resembles closely the fractal dimension of the minimum path D min in percolation theory (Sahimi, 2011), which is dened by, L min L Dmin , where L min is the length of the shortest path between two widely separated points, and L is the linear size of the sample. If there are no extended correlations, D min 1:13 and 1.34 for D E = 2 and 3. In the presence of such correlations, however, D min 1 for all 0 H 1 (Knackstedt et al., 2000). That is, such correlations make the percolation clusters smoother and more homogeneous, a fact that was rst pointed out by Sahimi and Mukhopadhyay (1996). Thus, if D T D min , then Eq. 6.20 would predict an approximate universal value m 1, which does not agree with our data. 6. Summary We analyzed vary extensive porosity, resistivity, and gamma-ray logs along 36 wells in six oil and gas reservoirs in Iran with dierent lithology and uid content, using the wavelet transform coherence and multifractal detrended uctuation analysis. The former was used to dierentiate the resistivity data associated with the clay content of the reservoirs from the rest of the logs, while the latter was utilized to analyze long-range correlations in the well logs. The cementation exponent m was estimated for two sets of data, one with clay content, and one without. All the logs, including the resistivity logs, exhibited long-range positive correlations with a Hurst exponent H > 0:5. The cementation exponent m is nonuniversal and varies signicantly. We propose that the nonuniversality of m is due to the presence of the long-range correlations, and that the analysis of the resistivity data for laboratory or outcrop samples, which do not contain 108 such correlations, yields estimates of m that are not applicable to oil and gas reservoirs and their well logs. 109 Bibliography Aghaei, A. and Piri, M. (2015). Direct pore-to-core up-scaling of displacement processes: Dynamic pore network modeling and experimentation. Journal of Hydrology, 522:488 { 509. Alejandre, J. and Hansen, J.-P. (2007). Ions in water: From ion clustering to crystal nucleation. Phys. Rev. E, 76:061505. Alim, K., Parsa, S., Weitz, D. A., and Brenner, M. P. (2017). Local pore size correlations determine ow distributions in porous media. Phys. Rev. Lett., 119:144501. Amooie, M. A., Soltanian, M. R., and Moortgat, J. (2017). Hydrothermodynamic mixing of uids across phases in porous media. Geophysical Research Letters, 44(8):3624{3634. 2016GL072491. An, S., Li, J., Li, Y., Li, S., Wang, Q., and Liu, B. (2016). Two-step crystal growth mechanism during crystallization of an undercooled ni50al50 alloy. Sci. Rep., 6:31062. Andersen, H. C. (1980). Molecular dynamics simulations at constant pressure and/or temperature. J. Chem. Phys., 72:2384. Ansari-Rad, M., Vaez Allaei, S., and Sahimi, M. (2012). Nonuniversality of roughness exponent of quasistatic fracture surfaces. Phys. Rev. E, 85:021121. Aragones, J. L., Sanz, E., and Vega, C. (2012). Solubility of nacl in water by molecular simulation revisited. J. Chem. Phys., 136:244508. Archie, G. (1942). The electrical resistivity log as an aid in determining some reservoir charac- teristics. Trans. AIME, 146:54. Aris, R. (1956). On the dispersion of a solute in a uid owing through a tube. Proc. R. Soc. Lond. A, 235:67. Aydin, M., Yano, T., Evrendilek, F., and Uygur, V. (2008). Implications of climate change for evaporation from bare soils in a mediterranean environment. Environmental Monitoring and Assessment, 140(1):123{130. Bacci, G., Korre, A., and Durucan, S. (2011). Experimental investigation into salt precipitation during co 2 injection in saline aquifers. Energy Procedia, 4:4450. Bai, X. M. and Li, M. (2006). Calculation of solid-liquid interfacial free energy: A classical nucleation theory based approach. J. Chem. Phys., 124:124707. Bakhshian, S. and Sahimi, M. (2016). Computer simulation of the eect of deformation on the morphology and ow properties of porous media. Physical Review E, 94:042903. Bakhshian, S. and Sahimi, M. (2017). Adsorption-induced swelling of porous media. Int. J. Greenhouse Gas Control, 57:pp. 1{13. 110 Balho, M., Sanchez-Rivera, D., Kwok, A., Mehmani, Y., and Prodanovi c, M. (2012). Numerical algorithms for network modeling of yield stress and other non-newtonian uids in porous media. Transport in Porous Media, 93(3):363{379. Barrett, R., Berry, M., Chan, T., Demmel, J., Donato, J., Dongarra, J., Eijkhout, V., Pozo, R., Romine, C., and van der Vorst, H. (1994). Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. Society for Industrial and Applied Mathematics. Bassiouni, Z. (1993). Measurement and Interpretation of Well Logs. SPE Textbook Series, Volume 4. Battistelli, A., Calore, C., and Pruess, K. (1997). The simulator tough2/ewasg for modelling geothermal reservoirs with brines and non-condensible gas. Geothermics, 26:437. Baumann, G., Henninges, J., and De Lucia, D. (2014a). Monitoring of saturation changes and salt precipitation during co2 injection using pulsed neutron-gamma logging at the ketzin pilot site. Int. J. Greenhouse Gas Control, 28:134. Baumann, G., Henninges, J., and De Lucia, D. (2014b). Monitoring of saturation changes and salt precipitation during co2 injection using pulsed neutron-gamma logging at the ketzin pilot site. Int. J. Greenhouse Gas Control, 28:134. Ben-Noah, I. and Friedman, S. P. (2018). Review and evaluation of root respiration and of natural and agricultural processes of soil aeration. Vadose Zone Journal, 17. 1. Bergstad, M. and Shokri, N. (2016). Evaporation of nacl solution from porous media with mixed wettability. Geophysical Research Letters, 43(9):4426{4432. 2016GL068665. Bernab e, Y. and Bruderer, C. (1998). Eect of the variance of pore size distribution on the transport properties of heterogeneous networks. Journal of Geophysical Research: Solid Earth, 103(B1):513{525. Billinge, S. J. L. (2009a). How do your crystals grow. Nat. Phys., 5:13. Billinge, S. J. L. (2009b). How do your crystals grow? Nat. Phys., 5:13. Bird, R., Stewart, W., and Lightfoot, E. (2007). Transport Phenomena. Wiley, New York, 2nd edition. Blunt, M. (2017). Multiphase Flow in Permeable Media. Camridsge University Press, London. Blunt, M., Jackson, J., Piri, M., and Valvatne, P. (2002). Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase ow. Adv. Water Resour., 25:1069. B onhorst, M., Walzel, P.and Rahimi, A., Kharaghani, A., Tsotsas, E., Nestle, N., Besser, A., Kleine J ager, F., and Metzger, T. (20162016). In uence of pore structure and impregnation drying conditions on the solid distribution in porous support materials. Drying Technol, 34:1964. Borgman, O., Fantinel, P., Lhder, W., Goehring, L., and Holtzman, R. (2017). Impact of spa- tially correlated pore-scale heterogeneity on drying porous media. Water Resources Research, 53(7):5645{5658. Bray, Y. L. and Prat, M. (1999). Three-dimensional pore network simulation of drying in capillary porous media. International Journal of Heat and Mass Transfer, 42(22):4207 { 4224. 111 Celia, M., Nordbotten, J., Court, B., Dobossy, M., and Bachu, S. (2011). Field-scale application of a semi-analytical model for estimation of co 2 and brine leakage along old wells. Int. J. Greenhouse Gas Control, 5:257. Chakrabarti, B. K. and Benguigui, L. (1997). Statistical Physics of Fracture and Breakdown in Disordered Systems. Oxford University Press, London, 2nd edition. Chakraborty, D. and Patey, G. N. (2013). How crystals nucleate and grow in aqueous nacl solution. J. Phys. Chem. Lett., 4:573. Chauvet, F., Duru, P., Georoy, S., and Prat, M. (2009). Three periods of drying of a single square capillary tube. Phys. Rev. Lett., 103:124502. Chen, F., Mourhatch, R., Tsotsis, T., and Sahimi, M. (2008). Pore network model of transport and separation of binary gas mixtures in nanoporous membranes. J. Membr. Sci., 315:48. Chu, J., Engquist, B., Prodanovi, M., and Tsai, R. (2012). A multiscale method coupling network and continuum models in porous media i: Steady-state single phase ow. Multiscale Modeling & Simulation, 10(2):515{549. Chung, S. Y., Kim, Y. M., Kim, J. G., and Kim, Y. J. (2009). Multiphase transformation and ostwalds rule of stages during crystallization of a metal phosphate. Nat. Phys., 5:68. Cygan, R. T. (2001). Molecular modeling in mineralogy and geochemistry. Rev. Mineral. Geochem., 42:1. Cygan, R. T., Liang, J. J., and Kalinichev, A. G. (2004). Molecular models of hydroxide, oxy- hydroxide, and clay phases and the development of a general force eld. J. Phys. Chem. B, 108:1255. Daian, J. and Saliba, J. (1991). Determining a representative random pore-network for moisture sorption and migration in cement mortar. Int. J. Heat Mass Transfer, 34:2081. Das, D. B., Hassanizadeh, S. M., Rotter, B. E., and Ataie-Ashtiani, B. (2004). A numerical study of micro-heterogeneity eects on upscaled properties of two-phase ow in porous media. Transport in Porous Media, 56(3):329{350. Dashtian, H. and Sahimi, M. (2014). Coherence index and curvelet transformation for denoising geophysical data. Phys. Rev. E, 90:042810. Dashtian, H., Shokri, N., and Sahimi, M. (2018). Pore-network model of evaporation-induced salt precipitation in porous media: The eect of correlations and heterogeneity. Advances in Water Resources, 112:59 { 71. Dashtian, H., Wang, H., and Sahimi, M. (2017). Nucleation of salt crystals in clay minerals: molecular dynamics simulation. J. Phys. Chem. Lett., 8:3166. Dashtian, H., Yang, Y., and Sahimi, M. (2015). Nonuniversality of the archie exponent due to multifractality of resistivity well logs. Geophysical Research Letters, 42(24):10,655{10,662. Datta, S. S. and Weitz, D. A. (2013). Drainage in a model stratied porous medium. EPL (Europhysics Letters), 101(1):14002. DeCarlo, K. F. and Shokri, N. (2014). Salinity eects on cracking morphology and dynamics in 3-d desiccating clays. Water Resour. Res., 50:3052. 112 Derrida, B. and Vannimenus, J. (1982). A transfer-matrix approach to random resistor networks. Journal of Physics A: Mathematical and General, 15(10):L557. Desarnaud, J., Bonn, D., and Shahidzadeh, N. (2016). The pressure induced by salt crystallization in connement. Sci. Rep., 6:30856. Desarnaud, J., Derluyn, H., Carmeliet, J., Bonn, D., and Shahidzadeh, N. (2014). Metastability limit for the nucleation of nacl crystals in connement. J. Phys. Chem. Lett., 5:890. Desarnaud, J., Derluyn, H., Molari, L., de Miranda, S., Cnudde, V., and Shahidzadeh, N. (2015). Drying of salt contaminated porous media; eect of primary and secondary nucleation. J. Appl. Phys., 118:114901. Dillard, L. A. and Blunt, M. J. (2000). Development of a pore network simulation model to study nonaqueous phase liquid dissolution. Water Resources Research, 36(2):439{454. Ebrahimi, F. (2010). Invasion percolation: A computational algorithm for complex phenomena. Comput. Sci. Eng., 12:84. Espinosa, J. R., Vega, C., Valeriani, C., and Sanz, E. (2015). The crystal- uid interfacial free en- ergy and nucleation rate of nacl from dierent simulation methods. J. Chem. Phys., 142:194709. Espinosa, J. R., Vega, C., Valeriani, C., and Sanz, E. (2016). Seeding approach to crystal nucle- ation. J. Chem. Phys., 144:034501. Ewing, R. and Hunt, A. (2006). Dependence of the electrical conductivity on saturation in real porous media. Vadose Zone, 5:731{741. Fantinel, P., Borgman, O., Holtzman, R., and Goehring, L. (2017). Drying in a micro uidic chip: experiments and simulations. Scientic Reports, 7(1):15572. Flatt, R. (2002). Salt damage in porous materials: how high supersaturations are generated. J. Crystal Growth, 242:435. Focke, J. and Munn, D. (1987). Cementation exponents in middle eastern carbonate reservoirs. SPE Form. Eval., 2:155{167. Freedman, M. A. (2015). Potential sites for ice nucleation on aluminosilicate clay minerals and related materials. J. Phys. Chem. Lett., 6:3850. Galamba, N., Nieto de Castro, C. A., and Ely, J. F. (2005). Shear viscosity of molten alkali halides from equilibrium and non-equilibrium molecular-dynamics simulations. J. Chem. Phys., 122:224501. Ghanbarian, B., A.G. Hunt, R. E., and Skinner, T. (2014). Universal scaling of the formation fac- tor in porous media derived by combining percolation and eective medium theories. Geophys. Res. Lett., 41:3884{3890. Ghanbarian, B., A.G. Hunt, T. S., and Ewing, R. (2015). Saturation dependence of transport in porous media predicted by percolation and eective medium theories. Fractals, 23:1540004. Ghassemzadeh, J., Hashemi, M., Sartor, L., and Sahimi, M. (2001). Pore network simulation of uid imbibition into paper during coating processes: I. model development. AIChE J., 47:519535. Ghassemzadeh, J. and Sahimi, M. (2004a). Molecular modelling of adsorption of gas mixtures in montmorillonites intercalated with al13-complex pillars. Mol. Phys., 102:1447. 113 Ghassemzadeh, J. and Sahimi, M. (2004b). Pore network simulation of uid imbibition into paper during coatingiii: Modeling of two-phase ow. Chem. Eng. Sci., 59:2281. Glover, P. (2009). What is the cementation exponent? a new interpretation. The Leading Edge, 28:82. Grinsted, A., Moore, J., and Jevrejeva, S. (2004). Application of the cross wavelet transform and wavelet coherence to geophysical time series. Nonlinear Process Geophys., 11:561. Guglielmini, L., Gontcharov, A., Jr., A., A.J., and Stone, A. (2008). Drying of salt solutions in porous materials: intermediate-time dynamics and eorescence. Phys. Fluids, 20:077101. Hamzehpour, H. and Khazaei, M. (2016). Eective permeability of heterogeneous fractured porous media. Transport in Porous Media, 113(2):329{344. Haro, C. (2010). The equations archie forgot: anisotropy of the rocks. SPE Reser. Eval. Eng., 13:823{836. Hekmatzadeh, M., Dadvar, M., and Sahimi, M. (2016). Pore-network simulation of unstable miscible displacements in porous media. Transp. Porous Med., 113:511. Hewett, T. (1986a). Fractal distributions of reservoir heterogeneity and their in uence on uid transport. SPE Paper 15386, New Orleans:Louisiana. Hewett, T. A. (1986b). Fractal Distributions of Reservoir Heterogeneity and Their In uence on Fluid Transport. Society of Petroleum Engineers, New Orleans, Louisiana. Holtzman, R. (2016). Eects of pore-scale disorder on uid displacement in partially-wettable porous media. Scientic Reports, 6:36221. Hoover, W. G. (1985). Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys., 31:1695. Hosseini, S. A. and Al, M. (2016). Time-lapse application of pressure transient analysis for monitoring compressible uid leakage. Greenhouse Gases: Science and Technology, 6(3):352{ 369. Huinink, H., Pel, L., and Michels, M. (2002). How ions distribute in a drying porous medium: A simple model. Phys. Fluids, 14:1389. Hunt, A. G. (2000). Percolation cluster statistics and conductivity semi-variograms. Transport in Porous Media, 39:131{141. Hunt, A. G. (2004). Percolative transport in fractal porous media. Chaos Solitons Fractals, 19:309{325. Imdakm, A. and Sahimi, M. (1987). Transport of large particles in ow through porous media. Phys. Rev. A, 36:5304. Jafari, G., M. Sahimi, M. R., and Tabar, M. R. (2011). Analysis of porosity distribution of large- scale porous media and their reconstruction by langevin equation. Phys. Rev. E, 83:026309. Jambhekar, V. A., Helmig, R., Schr oder, N., and Shokri, N. (2015). Free- ow{porous-media coupling for evaporation-driven transport and precipitation of salt in soil. Transport in Porous Media, 110(2):251{280. Jasti, J. K., Jesion, G., and Feldkamp, L. (1993). Microscopic imaging of porous media with x-ray computer tomography. 114 Kantelhardt, J., D. Rybski, S.A. Zschiegner, P. B. E. B. V. L. S. H., and Bunde, A. (2003). Multifractality of river runs and precipitation: comparison of uctuation analysis and wavelet methods. Physica A, 330:240{245. Kantelhardt, J., S.A. Zschiegner, E.K. Bunde, S. H. A. B., and Stanley, H. (2002). Multifractal detrended uctuation analysis of nonstationary time series. Physica A,, 316:87{114. Kim, M., Sell, A., and Sinton, D. (2013). Aquifer-on-a-chip: understanding pore-scale salt pre- cipitation dynamics during co 2 sequestration. Lab. Chip, 13:2508. Kim, N., Harale, A., Tsotsis, T. T., and Sahimi, M. (2007). Atomistic simulation of nanoporous layered double hydroxide materials and their properties: Ii. adsorption and diusion. J. Chem. Phys., 127:224701. Kim, N., Kim, Y., Tsotsis, T. T., and Sahimi, M. (2005). Atomistic simulations of nanoporous layered double hydroxide materials and their properties i. structural modeling. J. Chem. Phys., 122:214713. Knackstedt, M., Sahimi, M., and Sheppard, A. (2000). Invasion percolation with long-range correlations: First-order phase transition and nonuniversal scaling properties. Phys. Rev. E, 61:4920. Knackstedt, M. A., Sheppard, A. P., and Pinczewski, W. V. (1998). Simulation of mercury porosimetry on correlated grids: Evidence for extended correlated heterogeneity at the pore scale in rocks. Phys. Rev. E, 58:R6923{R6926. Knackstedt, M. A., Sheppard, A. P., and Sahimi, M. (2001). Pore network modelling of two- phase ow in porous rock: the eect of correlated heterogeneity. Advances in Water Resources, 24(3):257 { 277. Pore Scale Modeling. Knott, B. C., Molinero, V., Doherty, M. F., and Peters, B. (2012). Homogeneous nucleation of methane hydrates: Unrealistic under realistic conditions. J. Am. Chem. Soc., 134:19544. Koplik, J., Redner, S., and Wilkinson, D. (1988). Transport and dispersion in random networks with percolation disorder. Phys. Rev. A, 37:2619. Lake, L. W., Johns, R., Rossen, B., and Pope, G. (2014). Fundamentals of Enhanced Oil Recovery. Society of Petroleum Engineers, 2nd edition. Laurindo, J. and Prat, M. (1998). Modelling of drying in capillary porous media: a discrete approach. Drying Technol., 16:1769. Le, D., Hoang, H., and Mahadevan, J. (2009). Impact of capillary-driven liquid lms on salt crystallization. Transport in Porous Media, 80(2):229. Le, K. H., Kharaghani, A., Kirsch, C., and Tsotsas, E. (2017). Discrete pore network modeling of superheated steam drying. Drying Technology, 35(13):1584{1601. Lehmann, P. and Or, D. (2009). Evaporation and capillary coupling across vertical textural contrasts in porous media. Phys. Rev. E, 80:046318. Li, L., Peters, C. A., and Celia, M. A. (2006). Upscaling geochemical reaction rates using pore- scale network modeling. Advances in Water Resources, 29(9):1351 { 1370. Liu, M., Zhang, S., and Mou, J. (2012). Eect of normally distributed porosities on dissolution pattern in carbonate acidizing. Journal of Petroleum Science and Engineering, 94-95:28 { 39. 115 Liu, P., Yao, J., Couples, G. D., Ma, J., and Iliev, O. (2017). 3-d modelling and experimental comparison of reactive ow in carbonates under radial ow conditions. Scientic Reports, 7(1):17711. Liu, W., Schmidt, B., Voss, G., and Mller-Wittig, W. (2008). Accelerating molecular dynamics simulations using graphics processing units with cuda. Computer Physics Communications, 179(9):634 { 641. Lu, J., Darvari, R., Nicot, J.-P., Mickler, P., and Hosseini, S. A. (2017). Geochemical impact of injection of eagle ford brine on hosston sandstone formationobservations of autoclave waterrock interaction experiments. Applied Geochemistry, 84:26 { 40. Lupi, L., Peters, B., and Molinero, V. (2016). Pre-ordering of interfacial water in the pathway of heterogeneous ice nucleation does not lead to a two-step crystallization mechanism. J. Chem. Phys., 145:211910. Lutsko, J. F. and Nicolis, G. (2006). Theoretical evidence for a dense uid precursor to crystal- lization. Phys. Rev. Lett., 96:046102. Mahadevan, J., Sharma, M., and Yortsos, Y. (2006). Flow-through drying of porous media. AIChE J., 52:2367. Malmir, H., Sahimi, M., and Rahimi Tabar, M. (2016a). Microstructural characterization of random packings of cubic particles. Sci. Rep., 6:paper 35024. Malmir, H., Sahimi, M., and Rahimi Tabar, M. (2016b). Packing of nonoverlapping cubic particles: computational algorithms and microstructural characteristics. Phys. Rev. E, 94:062901. Malmir, H., Sahimi, M., and Rahimi Tabar, M. (2017). Statistical characterization of microstruc- ture of packings of polydisperse hard cubes. Phys. Rev. E, 95:052902. Mandelbrot, B. B. and Ness, J. W. V. (1968). Fractional brownian motions, fractional noises and applications. SIAM Review, 10(4):422{437. Markall, G. R., Ham, D. A., and Kelly, P. H. J. (2010). Generating optimised nite element solvers for gpu architectures. AIP Conference Proceedings, 1281(1):787{790. Markall, G. R., Slemmer, A., Ham, D. A., Kelly, P. H. J., Cantwell, C. D., and Sherwin, S. J. (2013). Finite element assembly strategies on multi-core and many-core architectures. Interna- tional Journal for Numerical Methods in Fluids, 71(1):80{97. Mehmani, A., Prodanovi c, M., and Javadpour, F. (2013). Multiscale, multiphysics network mod- eling of shale matrix gas ows. Transport in Porous Media, 99(2):377{390. Mehmani, Y. and Balho, M. (2015). Eulerian network modeling of longitudinal dispersion. Water Resour. Res., 51:8586. Mehmani, Y. and Tchelepi, H. A. (2017). Minimum requirements for predictive pore-network modeling of solute transport in micromodels. Advances in Water Resources, 108:83 { 98. Mehrabi, A., Rassamdana, H., and Sahimi, M. (1997). Characterization of long-range correlations in complex distributions and proles. Phys. Rev. E, 56:712{722. Mendoza, F. N. and Alejandre, J. (2013). The role of ion-water interactions in the solubility of ionic solutions. J. Mol. Liq., 185:50. 116 Miri, R. and Hellevang, H. (2016). Salt precipitation during co 2 storage: a review. Int. J. Greenhouse Gas Control, 51:136. Miri, R., van Noort, R., Aagaard, P., and Hellevang, H. (2015a). New insights on the physics of salt precipitation during injection of co2 into saline aquifers. International Journal of Greenhouse Gas Control, 43:10 { 21. Miri, R., van Noort, R., Aagaard, P., and Hellevang, H. (2015b). New insights on the physics of salt precipitation during injection of co2 into saline aquifers. Int. J. Greenhouse Gas Control, 43:10. Moghaddam, A. A., Kharaghani, A., Tsotsas, E., and Prat, M. (2017). Kinematics in a slowly drying porous medium: Reconciliation of pore network simulations and continuum modeling. Physics of Fluids, 29(2):022102. Mohammadmoradi, P. and Kantzas, A. (2016). Petrophysical characterization of porous media starting from micro-tomographic images. Advances in Water Resources, 94:200 { 216. Molina-Montes, E., Timon, V., Hernandez-laguna, A., and Sainzdiaz, C. I. (2008). Dehydroxy- lation mechanisms in al3+/fe3+ dioctahedral phyllosilicates by quantum mechanical methods with cluster models. Geochim. Cosmochim. Acta, 72:3929. Molz, F. J., Liu, H. H., and Szulga, J. (1997). Fractional brownian motion and fractional gaussian noise in subsurface hydrology: A review, presentation of fundamental properties, and extensions. Water Resources Research, 33(10):2273{2286. Mourhatch, R., Tsotsis, T., and Sahimi, M. (2010). Network model for the evolution of the pore structure of silicon-carbide membranes during their fabrication. J. Membr. Sci., 356:138. Mucha, M. and Jungwirth, P. (2003). Salt crystallization from an evaporating aqueous solution by molecular dynamics simulations. J. Phys. Chem. B, 107:8271. Mukhopadhyay, S. and Sahimi, M. (2000). Calculation of the eective permeabilities of eld-scale porous media. Chem. Eng. Sci., 55:4495. Muljadi, B. P., Blunt, M. J., Raeini, A. Q., and Bijeljic, B. (2016). The impact of porous media heterogeneity on non-darcy ow behaviour from pore-scale simulation. Advances in Water Resources, 95:329 { 340. Pore scale modeling and experiments. Muller, N., Qi, R., Mackie, E., Pruess, K., and Blunt, M. (2009). Co 2 injection impairment due to halite precipitation. Energy Procedia, 1:3507. Nachshon, U., Weisbrod, N., Dragila, M., and Grader, A. (2001). Combined evaporation and salt precipitation inhomogeneous and heterogeneous porous media. Water Resour. Res., 47:W03513. Nadim, F., Hoag, G. E., Liu, S., Carley, R. J., and Zack, P. (2000). Detection and remediation of soil and aquifer systems contaminated with petroleum products: an overview. Journal of Petroleum Science and Engineering, 26(1):169 { 178. Neuman, S. P. (1994). Generalized scaling of permeabilities: Validation and eect of support scale. Geophysical Research Letters, 21(5):349{352. Nogues, J. P., Fitts, J. P., Celia, M. A., and Peters, C. A. (2013). Permeability evolution due to dissolution and precipitation of carbonates using reactive transport modeling in pore networks. Water Resources Research, 49(9):6006{6021. 117 Noiriel, C., Renard, F., Doan, M., and Gratier, J. (2010). Intense fracturing and fracture sealing induced by mineral growth in porous rocks. Chem. Geol., 269:197. Normand, J. M. and Herrmann, H. J. (1990). Precise numerical determination of the supercon- ducting exponent of percolation in three dimensions. International Journal of Modern Physics C, 01(02n03):207{214. Norouzi Rad, M. and Shokri, N. (2012). Nonlinear eects of salt concentrations on evaporation from porous media. Geophys. Res. Lett., 39:L04403. Norouzi Rad, M. and Shokri, N. (2014). Eects of grain angularity on nacl precipitation in porous media during evaporation. Water Resour. Res., 50:9020. Norouzi Rad, M., Shokri, N., and Sahimi, M. (2013). Pore-scale dynamics of salt precipitation in drying porous media. Phys. Rev. E, 88:032404. Nos, S. (1984). A unied formulation of the constant temperature molecular dynamics methods. J. Chem. Phys., 81:511. NvidiaCorporation (2015). Nvidia cuda programming guide. http://docs.nvidia.com/cuda/cuda- c-programming-guide/abstract. Olivares-Amaya, R., Watson, M. A., Edgar, R. G., Vogt, L., Shao, Y., and Aspuru-Guzik, A. (2010). Accelerating correlated quantum chemistry calculations using graphical processing units and a mixed precision matrix multiplication library. Journal of Chemical Theory and Computation, 6(1):135{144. PMID: 26614326. Or, D. (2008). Scaling of capillary, gravity and viscous forces aecting ow morphology in unsat- urated porous media. Advances in Water Resources, 31(9):1129 { 1136. Or, D., Lehmann, P., Shahraeeni, E., and Shokri, N. (2013). Advances in soil evaporation physics{ a review. Vadose Zone Journal, 12. 4. Oren, P.-E., Bakke, S., and Arntzen, O. J. (1998). Extending predictive capabilities to network models. Ostwald, W. (1897). Studien uber die bildung und umwandlung fester korper. Z. Phys. Chem., 22:289. Ott, H., de Kloe, K., Marcelis, F., and Makurat, A. (2011). Injection of supercritical co 2 in brine saturated sandstone: pattern formation during salt precipitation. Energy Procedia, 4:4425. Ou, X., Li, J., and Lin, Z. (2014). Dynamic behavior of interfacial water on mg(oh)2 (001) surface: A molecular dynamics simulation work. J. Phys. Chem. C, 118:29887. Ovaysi, S. and Piri, M. (2010). Direct pore-level modeling of incompressible uid ow in porous media. Journal of Computational Physics, 229(19):7456 { 7476. Pazhoohesh, E., Hamzehpour, H., and Sahimi, M. (2006). Numerical simulation of ac conduction in three-dimensional heterogeneous materials. Phys. Rev. B, 73:174206. Pegler, S. S., Huppert, H. E., and Neufeld, J. A. (2016). Stratied gravity currents in porous media. Journal of Fluid Mechanics, 791:329357. Pereyra, R. G., Szleifer, I., and Carignano, M. A. (2011). Temperature dependence of ice critical nucleus size. J. Chem. Phys., 135:034508. 118 Peysson, Y., Andre, L., and Azaroual, M. (2014). Well injectivity during co 2 storage operations in deep saline aquiferspart 1: experimental investigation of drying eects, salt precipitation and capillary forces. Int. J. Greenhouse Gas Control, 22:291. Piri, M. and Blunt, M. (2005). Three-dimensional mixed-wet random pore-scale network modeling of two- and three-phase ow in porous media. i. model description. Phys. Rev. E, 71:026301. Plimpton, S. J. (1995). Fast parallel algorithms for short-range molecular. J. Comput. Phys., 117:1. Porter, C. and Carothers, J. (1971). Formation factor-porosity relation derived from well log data. The Log Analyst, 12:16{26. Prat, M. (1995). Isothermal drying on non-hygroscopic capillary-porous materials as an invasion percolation process. International Journal of Multiphase Flow, 21(5):875 { 892. Prat, M. (2002). Recent advances in pore-scale models for drying of porous media. Chemical Engineering Journal, 86(1):153 { 164. Drying 2000: Selected Papers from the 12th International Symposium IDS2000. Prat, M. (2011). Pore network models of drying, contact angle, and lm ows. Chemical Engi- neering Technology, 34(7):1029{1038. Qazi, M. J., Lieerink, R. W., Schlegel, S. J., Backus, E. H. G., Bonn, D., and Shahidzadeh, N. (2017). In uence of surfactants on sodium chloride crystallization in connement. Langmuir, 33(17):4260{4268. Rahimi, A., Metzger, T., Kharaghani, A., and Tsotsas, E. (2016). Interaction of droplets with porous structures: Pore network simulation of wetting and drying. Drying Technology, 34(9):1129{1140. Rajabbeigi, N., Tsotsis, T., and Sahimi, M. (2009). Molecular pore-network model for nanoporous. J. Membr. Sci., 345:323. Roels, S., Ott, H., and Zitha, P. (2014). -ct analysis and numerical simulation of drying eects of co2 injection into brine-saturated porous media. Phys. Rev. Lett., 66:1169. Sahimi, M. (2011). Flow and Transport in Porous Media and Fractured Rock. Wiley-VCH, Weinheim, 2nd edition. Sahimi, M. and Imdakm, A. (1991). Hydrodynamics of particulate motion in porous media. Sahimi, M. and Mukhopadhyay, S. (1996). Scaling properties of a percolation model with long- range correlations. Phys. Rev. E, 54:3870{3880. Sahimi, M. and Tajer, S. (2005). Self-ane distributions of the bulk density, elastic moduli, and seismic wave velocities of rock. Phys. Rev. E, 71:046301. Sanz, E., Vega, C., Espinosa, J. R., Caballero-Bernal, R., Abascal, J. L. F., and Valeriani, C. (2013). Homogeneous ice nucleation at moderate supercooling from molecular simulation. J. Am. Chem. Soc., 135:15008. Scanlon, B. R., Faunt, C. C., Longuevergne, L., Reedy, R. C., Alley, W. M., McGuire, V. L., and McMahon, P. B. (2012). Groundwater depletion and sustainability of irrigation in the us high plains and central valley. Proceedings of the National Academy of Sciences, 109(24):9320{9325. Scherer, G. (1990a). Theory of drying. J. Am. Ceram. Soc., 73:3. 119 Scherer, G. W. (1990b). Theory of drying. Journal of the American Ceramic Society, 73(1):3{14. Schiro, M., Ruiz-Agudo, E., and Rodriguez-Navarro, C. (2012). Damage mechanisms of porous materials due to in-pore salt crystallization. Phys. Rev. Lett., 109:265503. Schoups, G., Hopmans, J., Young, C., Vrugt, J., Wallender, W., Tanji, K., and Panday, S. (2005). Sustainability of irrigated agriculture in the san joaquin valley, california. Proc. Natl. Acad. Sci. USA, 102:15352. Sen1981, P., Scala, C., and Cohen, M. (1981). A self-similar model for sedimentary rocks with application to dielectric constant of fused glass beads. Geophysics, 46:781{795. Sghaier, N., Georoy, N., Prat, M., Eloukabi, H., and Ben Nasrallah, S. (2014). Evaporation- driven growth of large crystallized salt structures in a porous medium. Phys. Rev. E, 90:042402. Shahidzadeh-Bonn, N., Desarnaud, J., Bertrand, F., Chateau, X., and Bonn, D. (2010). Damage in porous media due to salt crystallization. Phys. Rev. E, 81:066110. Shahraeeni, E., Lehmann, P., and Or, D. (2012). Coupling of evaporative uxes from drying porous surfaces with air boundary layer: Characteristics of evaporation from discrete pores. Water Resour. Res., 48:W09525. Shahraeeni, E. and Or, D. (2010). Thermo-evaporative uxes from heterogeneous porous surfaces resolved by infrared thermography. Water Resour. Res., 46:W09511. Shaw, M. (1987). Drying as an immiscible displacement process with uid counter ow. Phys. Rev. Lett., 59:1671. Shklovskii, B. I. and Efros, A. (1984). Electronic Properties of Doped Semiconductors. Springer, Berlin, 2nd edition. Shokri, N. (2014). Pore-scale dynamics of salt transport and distribution in drying porous media. Phys. Fluids, 26:012106. Shokri, N., Lehmann, P., and Or, D. (2008a). Eects of hydrophobic layers on evaporation from porous media. Geophysical Research Letters, 35(19):n/a{n/a. L19407. Shokri, N., Lehmann, P., and Or, D. (2010). Evaporation from layered porous media. Journal of Geophysical Research: Solid Earth, 115(B6). B06204. Shokri, N., Lehmann, P., Vontobel, P., and Or, D. (2008b). Drying front and water content dynamics during evaporation from sand delineated by neutron radiography. Water Resources Research, 44(6):n/a{n/a. W06418. Shokri, N. and Or, D. (2011). What determines drying rates at the onset of diusion controlled stage-2 evaporation from porous media? Water Resour. Res., 47:W09513. Shokri, N., Zhou, P., and Keshmiri, A. (2015). Patterns of desiccating cracks in saline bentonite layers. Transp. Porous Media, 110:333. Shokri-Kuehni, S., Vetter, T., Webb, C., and Shokri, N. (2017). New insights into saline water evaporation from porous media: Complex interaction between evaporation rates, precipitation and surface temperature. Geophys. Res. Lett., 44:5504. Soltanian, M. R., Amooie, M. A., Cole, D. R., Graham, D. E., Hosseini, S. A., Hovorka, S., Pner, S. M., Phelps, T. J., and Moortgat, J. (2016). Simulating the craneld geological carbon sequestration project with high-resolution static models and an accurate equation of state. International Journal of Greenhouse Gas Control, 54:282 { 296. 120 Sorbie, K. and Cliord, P. (1991). The inclusion of molecular diusion eects in the network modelling of hydrodynamic dispersion in porous media. Chem. Eng. Sci., 46:2525. Surasani, V., Metzger, T., and Tsotsas, E. (2008). Consideration of heat transfer in pore network modelling of convective drying. International Journal of Heat and Mass Transfer, 51(9):2506 { 2518. Taylor, G. (1953). Dispersion of soluble matter in solvent owing slowly through a tube. Proc. Roy. Soc. Lond. A, 219:186. ten Wolde, P. R., Ruiz-Montero, M. J., and Frenkel, D. (1995). Numerical evidence for bcc ordering at the surface of a critical fcc nucleus. Phys. Rev. Lett., 75:2714. Threlfall, T. (2003). Structural and thermodynamic explanations of ostwalds rule. Org. Process Res. Dev., 7:1017. Toner, J., Catling, D., and Light, B. (2015). Modeling salt precipitation from brines on mars: Evaporation versus freezing origin for soil salts. Icarus, 250:451. Torquato, S. (2002). Random Heterogeneous Materials. Springer, New York, 2nd edition. U., N., Shahraeeni, A., Or, D., Dragila, M., and Weisbrod, N. (2011). Infrared thermography of evaporative uxes and dynamics of salt deposition on heterogeneous porous surfaces. Water Resour. Res., 47:W12519. Umtsev, I. S. and Martinez, T. J. (2009). Quantum chemistry on graphical processing units. 3. analytical energy gradients, geometry optimization, and rst principles molecular dynamics. Journal of Chemical Theory and Computation, 5(10):2619{2628. PMID: 26631777. Vasconcelos, I. F., Bunker, B. A., and Cygan, R. T. (2007). Molecular dynamics modeling of ion adsorption to the basal surfaces of kaolinite. J. Phys. Chem. C, 111:6753. Vega, C. and Abascal, J. L. F. (2011). Simulating water with rigid non-polarizable models: A general perspective. Phys. Chem. Chem. Phys., 13:19663. Vekilov, P. G. (2005). Two-step mechanism for the nucleation of crystals from solution. J. Cryst. Growth, 275:65. Veran-Tissoires, S. and Prat, M. (2012). Discrete salt crystallization at the surface of a porous medium. Phys. Rev. Lett., 108:054502. Veran-Tissoires, S. and Prat, M. (2014). Evaporation of a sodium chloride solution from a satu- rated porous medium with eorescence formation. J. Fluid Mech., 749:701. Volpe, G. and Volpe, G. (2017). The topography of the environment alters the optimal search strategy for active particles. Proceedings of the National Academy of Sciences, 114(43):11350{ 11355. Vorhauer, N., Tsotsas, E., and Prat, M. (2017). Drying of thin porous disks from pore network simulations. Drying Technology, 0(0):1{13. Wang, J., Sun, Z., Lu, G., and Yu, J. (2014). Molecular dynamics simulations of the local structures and transport coecients of molten alkali chlorides. J. Phys. Chem. B, 118:10196. Waxman, M., , and Smits, L. (1968). Electrical conductivities in oil-bearing shaly sands. SPE J., 8:107{122. 121 Wei, W., J. Cai, X. H., and Han, Q. (2015). An electrical conductivity model for fractal porous media. Geophys. Res. Lett., 42:4833{4840. Yiotis, A., Stubos, A., Boudouvis, A., and Yortsos, Y. (2001). A 2-d pore-network model of the drying of single-component liquids in porous media. Advances in Water Resources, 24(3):439 { 460. Pore Scale Modeling. Yiotis, A., Tsimpanogiannis, I., Stubos, A., and Yortsos, Y. (2006). Pore-network study of the characteristic periods in the drying of porous materials. J. Colloid Interface Sci., 297:738. Yiotis, A. G., Salin, D., and Yortsos, Y. C. (2015). Pore network modeling of drying processes in macroporous materials: Eects of gravity, mass boundary layer and pore microstructure. Transport in Porous Media, 110(2):175{196. Yonezawa, F. and Cohen, M. (1983). Granular eective medium approximation. J. Appl. Phys., 54:2895. Zabolitzky, J. G., Bergman, D. J., and Stauer, D. (1986). Precision calculation of elasticity for percolation. Journal of Statistical Physics, 44(1):211{223. Zeidouni, M., Pooladi-Darvish, M., and Keith, D. (2008). Analytical solution to evaluate salt precipitation during co 2 injection in saline aquifers. Energy Procedia, 1:1775. Zhang, T. H. and Liu, X. Y. (2007). Multistep crystal nucleation: A kinetic study based on colloidal crystallization. J. Phys. Chem. B, 111:14001. Zhang, Y., Kogure, T., Chiyonobu, S., Lei, X., and Xue, Z. (2013). In uence of heterogeneity on relative permeability for co2/brine: Ct observations and numerical modeling. Energy Procedia, 37:4647 { 4654. Zhao, B., MacMinn, C. W., and Juanes, R. (2016). Wettability control on multiphase ow in patterned micro uidics. Proceedings of the National Academy of Sciences, 113(37):10251{10256. Zheng, M., Li, X., and Guo, L. (2013). Algorithms of gpu-enabled reactive force eld (reax) molecular dynamics. Journal of Molecular Graphics and Modelling, 41:1 { 11. Zimmermann, N. E. R., Vorselaars, B., Quigley, D., and Peters, B. (2015). Nucleation of nacl from aqueous solution: Critical sizes, ion-attachment kinetics, and rates. J. Am. Chem. Soc., 137:13352. 122
Abstract (if available)
Abstract
Salt transport and precipitation in porous media constitute a set of complex and fascinating phenomena that are of fundamental interest to several important problems, ranging from storage of CO₂ in geological formations, to soil fertility, and protection of pavements and roads, as well as historicalmonuments. The phenomena occur at the pore scale and are greatly influenced by the heterogeneity of the pore space morphology. In this thesis salt precipitation is studied at two scales: molecular and pore scales. First, we use molecular dynamic simulation to study evaporation and salt precipitation in a single pore. Salt crystals precipitate on the pores' surface, modify the pore space morphology, and reduce its flow and transport properties. Despite numerous efforts to understand the mechanisms of nucleation of salt crystals in confined media, the effect of the rock's chemistry, which is mostly clay,on the growth, distribution and properties of the crystals is not well understood. We report the results of extensive molecular dynamics simulation of nucleation and growth of NaCl crystals in a clay pore usingmolecular models of two types of clay minerals, the Na-montmorillonite and kaolinite. Clear evidence is presented for the nucleation of the salt crystals that indicates that the molecular structure of clay minerals affects their spatial distribution, although the nucleation mechanism is the same in both types of clays. ❧ We then develop a pore-network model of flow of brine, and salt transport and precipitation in porous medis. To carry our efficient simulations with large pore networks, one must develop new computational algorithm. The computation time and required memory are two limiting factors that hinder the scalability of the computations to very large network sizes. We present a dual approach, based on the use of a combination of the central processing units (CPUs) and graphics processing units (GPUs), to simulation of flow, transport, and similar problems using pore networks. A mixed-precision algorithm, together with the conjugate-gradient method is implemented on a single GPU solver. The efficiency of the method is tested with a variety of cases, including pore- and random-resistor network models in which the conductances are long-range correlated, and also contain percolation disorder. Both isotropic and anisotropic networks are considered. To put the method to a stringent test, the long-range correlations are generated by a fractional Brownian motion (FBM), which we generate by a message-passing interface method. For all the cases studied an overall speed-up factor of about one order of magnitude or better is obtained, which increases with the size of the network that one can use. Even the critical slow-down in networks near the percolation threshold does not decrease the speed-up significantly. We also obtain approximate but accurate bounds for the permeability anisotropy Kₓ/Kᵧ [K sub x over K sub y] for stratified porous media. ❧ Before developing a pore-network model for salt transport and precipitation, one must also develop anmodel for evaporation in porous media. At the pore-scale, a threshold capillary pressure Pᵧ ∽ ɣ/aₜ is needed for the invading (vapor) phase at the throat entrance in order to expel the wetting fluid from a throat/pore in a porous medium. The threshold pressure depends on the interfacial tension ɣ between the wetting and non-wetting fluids and the entrance wetting radius rₜ of the throat/pore. Isotropic and homogeneous porous media consist of throats/pores of a sole average radius and, hence, the threshold capillary pressure usually outweighs the flow-induced viscous pressure in apore. Thus, a minimal variation in rₜ changes the fluid phase distribution and transport properties of the porous medium during drainage invasion/drainage process. Natural porous media, such as geological formations, are stratified, consisting of more or less parallel layers, with each layer having a different morphological properties. The extra variations in the pore space strongly affects flow andtransport in porous media. Thus, we develop a three-dimensional pore-network model of evaporation inlayered porous media, in which we use a fractional Brownian motion (FBM) to include the impact of the pore-size disorder, heterogeneity, and correlations. Our simulations indicate that diffusion of theevaporating phase, drying patterns of liquid, and saturation profiles are all affected strongly by the correlations and disorder in pore space, and that, while the increase in the stochasticity of thepore sizes stabilizes the drying front, the increase in the anisotropy parameter results in a more disordered drying front. ❧ We then present a pore-network model to study salt transport and precipitation in pores media. Vapor diffusion, capillary effect at the brine-vapor interface, flow of brine, and transport of salt and its precipitation in the pores that plug the pores partially or completely are all accounted for. The drying process is modeled by the invasion percolation, while transport of salt in brine is accounted for by the convective-diffusion equation. We demonstrate that the drying patterns, the clustering and connectivity of the pore throats in which salt precipitation occurs, the saturation distribution, and the drying rate are all strongly dependent upon the pore-size distribution, the correlations among the pore sizes, and the anisotropy of the pore space caused by stratification that most natural porous media contain. In particular, if the strata are more or less parallel to the direction of injection of the gas that dries out the pore space (air, for example) and/or causes salt precipitation (CO₂, for example), the drying rate increases significantly. Moreover, salt tends to precipitate in clusters of neighboring pores that are parallel to the open surface of the porous medium. ❧ In order to demonstrate the existence of the heterogeneity, correlations and long-range correlations in large-scale porous medium, we also study the well logs of hydrocarbon reservoirs. Archie's law expresses a relation between the formation factor F of porous media and their porosity ϕ, F∝ϕ⁻ᵐ, where m is the Archie or the cementation exponent. Despite widespread use of Archie's law, the value of m and whether it is universal and independent of the type of reservoir have remained controversial. We analyze various porosity and resistivity logs along 36 wells in six Iranian oil and gas reservoirs using wavelet transform coherence and multifractal detrended fluctuation analysis. m is estimated for two sets of data: one set contains the resistivity data that include those segments of the well that contain significant clay content, and one without. The analysis indicates that the well logs are multifractal, and that due to the multifractality the exponent m is non-universal. Thus, analysis of the resistivity of laboratory or outcrop samples that are not multifractal yields estimates of m that are not applicable to well logs in oil or gas reservoirs.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Efficient simulation of flow and transport in complex images of porous materials and media using curvelet transformation
PDF
Machine-learning approaches for modeling of complex materials and media
PDF
Molecular- and continuum-scale simulation of single- and two-phase flow of CO₂ and H₂O in mixed-layer clays and a heterogeneous sandstone
PDF
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
Effective flow and transport properties of deforming porous media and materials: theoretical modeling and comparison with experimental data
PDF
The effects of gravity and thermal gradient on the drying processes in porous media
PDF
Pore-scale characterization and fluid-flow simulation: application of the Lattice Boltzmann method to conventional and unconventional rocks
PDF
Stability and folding rate of proteins and identification of their inhibitors
PDF
Inverse modeling and uncertainty quantification of nonlinear flow in porous media models
PDF
The role of counter-current flow in the modeling and simulation of multi-phase flow in porous media
PDF
On the transport dynamics of miscible fluids in porous media under different sources of disorder
PDF
Flow and thermal transport at porous interfaces
PDF
Synergistic coupling between geomechanics, flow, and transport in fractured porous media: applications in hydraulic fracturing and fluid mixing
PDF
Exploring properties of silicon-carbide nanotubes and their composites with polymers
PDF
Hybrid physics-based and data-driven computational approaches for multi-scale hydrologic systems under heterogeneity
PDF
Oscillatory and streaming flow due to small-amplitude vibrations in spherical geometry
PDF
Deep learning architectures for characterization and forecasting of fluid flow in subsurface systems
PDF
Multiscale modeling of cilia mechanics and function
PDF
Waterflood induced formation particle transport and evolution of thief zones in unconsolidated geologic layers
Asset Metadata
Creator
Dashtian, Hassan
(author)
Core Title
Multiscale and multiresolution approach to characterization and modeling of porous media: From pore to field scale
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Chemical Engineering
Publication Date
07/16/2018
Defense Date
07/16/2018
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
evaporation,GPU,Modeling,molecular dynamics,OAI-PMH Harvest,pore network,porous media,salt precipitation,simulation
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Sahimi, Muhammad (
committee chair
), De Barros, Felipe (
committee member
), Shing, Katherine S-F (
committee member
)
Creator Email
dashtian@usc.edu,khayamdashti@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c89-17521
Unique identifier
UC11670727
Identifier
etd-DashtianHa-6405.pdf (filename),usctheses-c89-17521 (legacy record id)
Legacy Identifier
etd-DashtianHa-6405.pdf
Dmrecord
17521
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Dashtian, Hassan
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
Tags
evaporation
GPU
molecular dynamics
pore network
porous media
salt precipitation
simulation