T H E O R E T I C A L M O D E L O F H Y D R A T E F O R M A T I O N I N N A T U R A L E N V I R O N M E N T S By Olga Zatsepina M. Sc. (Physics) Moscow State University A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF SCIENCE in T H E FACULTY OF GRADUATE STUDIES EARTH AND OCEAN SCIENCES We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA January 1997 © Olga Zatsepina, 1997 I n p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f t h e r e q u i r e m e n t s f o r a n a d v a n c e d d e g r e e at t h e U n i v e r s i t y o f B r i t i s h C o l u m b i a , I a g r e e t h a t t h e L i b r a r y s h a l l m a k e it f r e e l y a v a i l a b l e f o r r e f e r e n c e a n d s t u d y . I f u r t h e r a g r e e t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f . t h i s t h e s i s f o r s c h o l a r l y p u r p o s e s m a y b e g r a n t e d b y t h e h e a d o f m y d e p a r t m e n t o r b y h i s o r h e r r e p r e s e n t a t i v e s . I t is u n d e r s t o o d t h a t c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l n o t b e a l l o w e d w i t h o u t m y w r i t t e n p e r m i s s i o n . D e p a r t m e n t o f EOS. T h e U n i v e r s i t y o f B r i t i s h C o l u m b i a V a n c o u v e r , C a n a d a D a t e 1V /oz/$? D E - 6 ( 2 / 8 8 ) Abstract The solubility of two hydrate-forming gases, C 0 2 and C H 4 , is calculated over a range of pressure and temperature. The solubility of gas is shown to be significantly altered by the presence or absence of hydrate. In particular, gas solubility changes abruptly in hydrate presence, allowing it to crytallize from the aqueous solution without the need of any free gas. To test this prediction a set of experiments was performed. In the experiments, an aqueous solution of CO2 was cooled at a pressure of 2 MPa. A variety of methods were examined to detect the growth of hydrate. With cooling and hydrate formation, the physical characteristics of the porous medium (temperature, porosity, gas concentration) change. On the basis of known governing equations and conductivity of an aqueous solution in porous medium, the conductivity change due to hydrate formation was predicted. Conductivity was found to be partic-ularly sensitive to hydrate formation, so electrical potential measurements were used to monitor hydrate growth. These electrical measurements indicated a pronounced resis-tance increase due to a change of gas concentration in the solution, corresponding to the amount of hydrate produced. Hydrate growth in the system was also detected in temperature data, which indicated a release of latent heat. The calculated phase diagram at typical pressure and temperature conditions in ma-rine environments is applied to establish the gas concentration needed to stabilize hydrate. This information is used in simple methods of hydrate formation to estimate the vertical distribution of hydrate in marine sediments and the rate of accumulation. ii Table of Contents Abstract ii List of Tables vi List of Figures vii Acknowledgement ix 1 Introduction 1 1.1 Background 1 1.2 Scope of Study 8 2 Thermodynamics of Hydrate Formation 11 2.1 Hydrate Stability 11 2.2 Phase Equilibria 13 2.3 Thermodynamic Model for Gas Hydrate 20 2.4 Optimization with Simulated Annealing 26 2.4.1 Very Fast Simulated Annealing 31 2.5 Results 32 3 Theoretical Development 42 3.1 Mathematical Formulation of Hydrate Formation Process 42 3.2 Conductivity of Aqueous Solution 50 3.3 Conductivity of Porous Media 57 iii 4 E x p e r i m e n t a l W o r k 5 9 4.1 Experimental Apparatus 66 4.2 Gas Solubility Experiments 69 4.3 Hydrate Formation 70 5 C o n d u c t i v i t y P r o f i l e T h r o u g h E l e c t r i c a l P o t e n t i a l M e a s u r e m e n t s 7 6 5.1 Governing Equations for DC Resistivity . . 76 5.2 Fourier Transformed Equations 78 5.3 Propagator Solution 79 5.4 Inverse Hankel Transform 82 6 D i s c u s s i o n 8 4 6.1 Theoretical Model 84 6.2 Experiment 85 6.3 Geometric Factor 88 6.4 Thermal Dependence of Conductivity 88 6.5 Gas Concentration 90 6.6 Application 92 6.7 Hydrate Fraction 97 6.8 Hydrate Distribution in Marine Environments 100 7 C o n c l u s i o n s 1 0 3 R e f e r e n c e s 1 0 6 A p p e n d i c e s 1 1 1 A P e n g - R o b i n s o n a n d T r e b b l e - B i s h n o i E q u a t i o n s o f S t a t e 1 1 1 iv B Mechanism of Gas Transference in Fluid 120 B . l Equation of State 120 B.2 Compositional Convection 121 B.3 Boundary Layer Analysis 124 B.4 Application of the Boundary Layer Analysis to Compositional Convection in the Experiment 129 B .5 Gas Transport in Porous Medium of the Experiment 130 v List of Tables 1.1 Geometry of Cage 4 1.2 Hydrate Crystal Cell Structure 5 1.3 Guests and Cavity Sizes for CSI and CSII 6 1.4 Hydrate Crystal Srtucture H 8 2.1 Structure of Gas Hydrates 21 2.2 Physical and Thermodynamic Reference Properties for Gas Hydrates . . 25 3.1 Henry's Law Constant for CO2—H20 system 54 6.1 Physical Properties 98 B . l Rayleigh number parameters 122 B.2 Porous medium Rayleigh number parameters 130 vi List of Figures 1.1 Schematic presentation of hydrate types 3 2.1 Typical profiles of hydrate in permafrost (A), and in ocean sediment (B) 12 2.2 Backbone of the simulated annealing algorithm 30 2.3 Carbon dioxide solubility as a function of temperature at a pressure of 2 MPa 34 2.4 Methane solubility as a function of temperature at pressure of 20 MPa. . 35 2.5 Carbon dioxide solubility as a function of temperature at pressure of 2 MPa. 37 2.6 Phase diagram for mixture methane - water for thermodynamic conditions encounted below the seafloor 40 3.1 Boundary and initial conditions for numerical modelling 48 3.2 Distribution of gas concentration, hydrate fraction, salt concentration and temperature (in dimensionless form) over depth at different time 49 3.3 Gas concentration, hydrate fraction, and salt concentration as a function of time at Z=8 51 4.1 Conductivity for dilute solutions of salt and gas as a function of time. . 61 4.2 Three possibilities for reaching hydrate stability zone from an aqueous solution at constant pressure: AEC, ADC, and ABC 63 4.3 Schematic dependence of resistance on temperature for different gas supply schedules: AEC, ADC, and ABC 65 4.4 Circuit for resistance measurements 67 vii 4.5 Voltage, resistance and current as a function of time 68 4.6 Average R2_3 as a function of gas solubility at room temperature. . . . 70 4.7 Change in resistance over time due to increase in gas concentration at pressure 50 psi and room temperature 71 4.8 Average R2-3 as a function of time at pressure 300 psi, as temperature reduced from 6°C to 2.2°C 73 4.9 RTD temperatures around current electrodes as a function of time (pres-sure 300 psi) 74 4.10 Average R2-3 as a function of a middle temperature (pressure 300 psi). . 75 5.1 Schematic position of electrodes, denoted by circles, in a layered conductive material 77 6.1 Conductivity as a function of T at constant gas concentration 91 6.2 Carbon dioxide gas concentration as a function of temperature at pressure of 2 MPa 95 6.3 Gas concentration as a function of T in hydrate forming process 96 6.4 Temperature change during hydrate formation. 99 B . l Schematic profile of temperature distribution in convective layer 125 B.2 A comparison of theory and experiment for a convecting fluid evolving toward a stationary state 131 viii Acknowledgement There are many people whom I would like to thank for the accomplishment of my thesis. First is, certainly, Bruce. Who brought me up to the last point with gentle patience, a lot of encouragements, and time spent on discussion of the problem. As I am staying for Ph.D. with Bruce, I save the other thanks for the Acknowledgement in my next program. I would like to thank Doug for his partial financial support and help in understand-ing some aspects of the problem, Peter for providing the material, which was vital for calculations, Partha for the algorithm of VFSA, and the students from the Department who were always ready to help. I would like to thank my sister who took all my responsibilities on her shoulders, allowing me to concentrate on work; and certainly, my children for being the sense in my life. ix Chapter 1 Introduction 1.1 Background Discovered early in the last century, gas hydrates were considered mainly as a laborat-ory curiosity. At that time, studies concentrated on looking for new hydrate formers, properties, lattice structures and phase diagrams [Jeffrey and McMullan, 1967]. These days, gas hydrates have drawn attention as a new natural gas resource with huge estimated deposits and as a safe storage medium (alternative to liquid natural gas). The efficiency of the gas storage in both natural and industrial settings is impressive. For example, the amount of methane, stored at 5 MPa and 278°K in hydrate phase, is theoretically equivalent to that of a liquified gas phase at 20 MPa. Hydrates also pose a hazard for the oil and gas industry. They form from drilling fluids and plug blowout preventers, block gas pipelines, and cause wells to be abandoned. In order to prevent hydrate from blocking pipelines both thermal insulation and chemical inhibitors are often used, with costs in the range from one to tens of million dollars. Therefore, the processes of gas hydrate formation and inhibition are of research interest. A key to understanding these processes is found in the hydrate structure. The crys-tallographic studies of gas hydrate began in the 1940's and 1950's when von Stackelberg and coworkers performed fundamental x-ray diffraction experiments on hydrates. The interpretation of the experiment results showed that all gas hydrates crystallize to form 1 Chapter 1. Introduction 2 one of two cubic structures (I and II), and that they are inclusion compounds: a meta-stable lattice of water molecules is stabilized at the specific thermodynamic conditions by inclusion of the second component. The forces binding this inclusion component are sim-ilar to intermolecular forces in liquid. Therefore, gas hydrate are considered as a member of the class of compounds called clathrate (the Latin " clathratus" means to encage). The characteristic feature of these compounds is that they consist of two components: a host lattice and guest molecules. The guest molecule is trapped by vacancy inside the host lattice (without occupying a lattice position), and interact with it only through weak van der Waals' interactions. The true clathrate hydrates are a subclass of the hydrates (Figure 1.1). The host-guest hydrate class is divided into three subclasses, which include true cla-thrates, semiclathrates and ionic clathrates, depending on the interaction between guest and host. The interaction is weak for the true clathrates, and it is significant for amines and etanol as well as for ionic subclass. The study of Jeffrey and McMullan [1967] found that the two hydrate structures are composed of repetitive crystal unit; the pentagon dodecahedra 5 1 2 cavity (12 faces with 5 edges per face) is an almost spherical "cage" of hydrogen-bonded water molecules. These units, a basic building block for both hydrate structures, are linked together along their vertices to form structure I, a body centered cubic structure, or through face sharing in three dimensions to form structure II, diamond cubic structure. As a result, two types of cavity are formed. One is a small cavity formed from the basic unit 5 1 2, sometimes called 12-hedra. A second, large cavity 5 1 26 2, or 14-hedra forms in structure I hydrate and 5 1 26 4, or 16-hedra, forms in structure II hydrate. The number of cavities in the two structures are different (Table 1.1). These 12-, 14- and 16-hedra are not stable in a pure water structure. Since hydrate cavities are larger than ice cavities, hydrate can be stabilized by the presence of a gas molecule either in the cavity itself or in a large Chapter 1. Introduction 3 CRYSTAL HYDRATES H O S T - G U E S T H Y D R A T E S O T H E R H Y D R A T E S TRUECLATHRATES SEMICLATHRATES IONIC CLATHRATES Structure I Structure II amines etanol salts — r — i acids bases Structure H Figure 1.1: Schematic presentation of hydrate types [Ripmeester, 1993] percentage of neighboring cavities. It is instructive to compare these results with the structure of ice. For example, the angle between neighbouring oxygen atoms in 12—hedra differs from that in ice by only 1.2°. Similarly, the 0-0 bond lengths exceed those in ice by only 1%. Despite these Chapter 1. Introduction 4 Hydrate Crystal Structure I II Cavity Small Large Small Large Description 5 1 2 5 1 26 2 5 1 2 5 1 26 4 Number of Cavities/Structure 2 6 16 8 Average Cavity Radius, A 3.91 4.33 3.902 4.683 Variation in Radius", % 3.4 14.4 5.5 1.73 Coordination number6 20 24 20 28 Table 1.1: Geometry of Cage. a distance of oxygen atoms from center of cage, 6 number of oxygens at the periphery of each cavity small differences, the 12—hedra structure maximizes the number of hydrogen bonds (30) to molecules (20) along the surface, when compared to similar cavities in ice. It was found via computer simulation of random hydrogen-bonded networks that the 5 1 2 polyhedra arise naturally in supercooled water. It can be seen from Table 1.1 that the most spherical cavity of the three types is the 16-hedra, as variation in radius is only 1.73%. The 14-hedra deviates most from a sphere and the angle between oxygen atoms deviates from that of ice by 5.1°. Both average radii and deviations from a sphere of the 5 1 2 cavity in hydrate I and II are different. As a result, it was thought that the smallest guest molecules would occupy the 5 1 2 cavity in structure I because its effective radius is smaller (effect of distortion). However, the assumption that the same guest molecules can fit into the 5 1 2 cavity in structure II was confirmed by Proton Nuclear Magnetic Resonanse (NMR) spectroscopy, an experimental technique for determining molecular occupation of cavities. Chapter 1. Introduction 5 Structure I II Crystal System Cubic Cubic Lattice Description Body Centered Diamond Lattice Parameter, A 12 17.3 Ideal Unit Cell Formular1 2X -6Y • 46H 20 16X- 8Y • 136H20 Table 1.2: Hydrate Crystal Cell Structure. 1 - X and Y denote generic guest molecules in small and large voids, respectively. The description of hydrate crystal cell structures is given in Table 1.2. Each cavity is occupied by, at most, one molecule. Molecules which occupy small cav-ities will also occupy the large ones of the structure. The probability of occupancy is determined mainly by the size of the gas molecule, while its chemical nature and shape are of secondary importance. Therefore, the occupancy of hydrate is related to the ratio of the gas molecule diameter to that of the free cavity. Deviations from the ideal filling increases as this ratio approaches unity. Thus, the hydration number, a ratio of numbers of water to gas molecules per unit cell, is directly affected by the size of the guest mole-cule. Table 1.3 outlines structural transitions that occur as the guest molecule sizes, or the largest van der Waals' diameters increase along with the hydration numbers for some typical guest molecules. The smallest guests, such as Ar, Kr, N 2 and 0 2 form structure II hydrate, with the guests occupying both cages. Molecules, such as C H 4 ) Xe and H 2S form structure I hydrate, with all the large cages and most of the small ones occupied. With increasing guest size, the occupancy of the 5 1 2 cavity drops, until for molecules like cyclopropane and ethane (d « 0.58 nm), only the 5 1 26 2 cavities are filled. Chapter 1. Introduction Hydrat e CSI Hydrate CSII Cavity description » 5 1 2 5 1 26 2 5 1 2 5 1 26 4 Cavity size, nm 3> 0.503 0.586 0.503 0.657 Cavity/unit cell » 2 6 16 8 Guest Guest Size (nm) Guest/Cavity Size Ratio Hydration Number Ar 0.380 0.756 0.649 0.756 0.579 Kr 0.400 0.795 0.683 0.795 0.609 N 2 0.410 0.815 0.700 0.815 0.624 o 2 0.420 0.835 0.717 0.835 0.640 M- 5 | H 2 0 C H 4 0.436 0.867 0.744 0.867 0.664 M- 5 | H 2 0 Xe 0.458 0.911 0.782 0.911 0.698 H 2S 0.458 0.911 0.782 0.911 0.698 C 0 2 0.512 1.018 0.874 1.018 0.780 N 2 0 0.525 1.044 0.897 1.044 0.800 C 2 H 4 0.550 1.094 0.939 1.094 0.838 C 2 H 2 0.573 1.139 0.978 1.139 0.873 M- 7 | H 2 0 C3H8 0.628 1.249 1.072 1.249 0.957 i-GjHio 0.650 1.292 1.110 1.292 0.990 M- 17H 20 Table 1.3: Guests and Cavity Sizes for CSI and CSII [Christensen & Sloan, 1994] Chapter 1. Introduction 7 In Table 1.3 the guest/cavity size ratios in boldface indicate cavities that are occupied. If the size ratio exceeds unity, the molecule will not fit within the cavity and the structure will not form. With a ratio significantly less than 0.75, the molecule can not lend enough stability to cavity to cause formation. For example, let's consider a hydrate of ethylene, C 2 H 4 . Ethylene is too large to fit in the 5 1 2 cavity. According to the guest/cavity size ratios from Table 1.3, it can fit into either the 5 1 26 2 or 5 1 26 4 cavity. But it lends greater stability to the 5 1 26 2 cavity. In the structure I unit cell, there are just two 5 1 2 cavities for every six 5 1 26 2 cavities, while in the structure II unit cell there are sixteen smaller cavities for every eight larger ones. It means that there are fractionally fewer empty 5 1 2 cavities in the structure I than in structure II. Therefore, ethylene forms CSI hydrate. As to hydrate of methane, C H 4 , a methane molecule can fit into all cavities, with lending greater stability to the 5 1 2 cavity and very little to the 5 1 26 4. NMR spectroscopy shows that the occupancy of the larger cavities in structure I and structure II usually exceeds 95%, while that for smaller ones is on the order of 50%. So, the stabilizing effect of a guest in the 5 1 2 cavity is not as great as in the larger cavities. Hence, methane forms structure I hydrate. Gases with small radii, for example, argon, form CSII because of the fractional abundance of the 5 1 2 cavity, lending the main stability benefit according to the guest/cavity size ratios, in the unit cell. The field of clathrate hydrate structure research is thought to be relatively mature. However, a recent discovery ( by Ripmeester in 1987) indicated that hydrocarbons larger than n-butane can be incorporated into a hydrate lattice designated as structure H. This has caused researchers to re-evaluate several aspects of clathrate hydrate. Structure H hydrate was proposed to have a hexagonal crystal system, with two types of smaller cavities. Data on the cage structure and an ideal unit cell formula was developed on the basis of NMR spectroscopic and powder diffraction measurements are given in Table 1.4. Chapter 1. Introduction 8 Cavity Small, X Small, X' Large, Y Description 4 3 5 6 6 3 5 1 26 8 Number of Cavities/Structure 3 2 1 Lattice Parameter, A 12.2 & 10 Ideal Unit Cell » 3X-2X'-Y- 34H20 Hydration number >^ 34 Table 1.4: Hydrate Crystal Srtucture H How molecules fit into the cavities is essential not only for the structure but also for equilibrium properties of hydrates. The fraction of certain type cavities occupied by gas component is a parameter in hydrate thermodynamics [van der Waals and Platteeuw, 1959], as is shown in the next chapter. The size of guest molecules is also important for hydrate formation [Christensen & Sloan, 1994]. 1.2 Scope of Study Since the discovery of hydrate in pipelines [Hammerschmidt, 1934], most research has concentrated on the three-phase equilibrium, where hydrate, aqueous solution and gas exist together. However, in natural settings hydrate stability zones correspond mainly to two phase equilibrium regions. In these regions hydrate is in equilibrium with aqueous so-lution. Thermodynamics of this system has not recieved much attention in the past. One of the key questions on the formation of hydrate in marine environment is the amount of hydrate-forming gas needed in water to initiate the process. Does the solution have to be Chapter 1. Introduction 9 supersaturated or not? How much gas, dissolved in water, is in equilibrium with hydrate at the beginning and later when hydrate has already formed? The thermodynamics of hydrate formation in a marine environment is the main focus of this study. Equilibrium of multiphase multicomponent systems is a subject of chemical thermo-dynamics. If the equations of state are known, then the recipes for tackling this type of problems can be found (see, for example, Nordstrom and Munoz, 1986). Thermodynamics of a multiphase system of water and carbon dioxide is considered in Chapter 2. The gas solubility as a function of pressure and temperature is determined for a wide range of P-T values. For two-phase equilibrium, the results agree with Henry's law at high temperatures and low pressures, while in hydrate stability zone, the gas solubility decreases exponentially with departures from the three-phase equilibrium point. Global minimization of the system Gibbs free energy demonstrates that hydrate can form from aqueous solution. To prove this assertion an experiment is designed. In the experiment, electrical field measurements are used to monitor the hydrate formation, since conductivity is expected to be very sensitive to the presence of hydrate. The theoretical analysis of hydrate formation and the associated changes in conduc-tivity are given in Chapter 3. A three-component mixture of water, salt and gas in a porous medium is considered. Conductivity is modelled on the basis of Archie's law. The proposal to add salt to the solution was made to facilitate the electrical field meas-urements. This idea was based on the fact that salt is excluded from hydrate, so the concentration of ions in the fluid increases as hydrate forms resulting in conductivity increase. But after some calculations were done, it appeared that the main effect of hydrate formation on conductivity was connected with change of gas concentration. For a typical calculation, hydrate formation caused gas concentration to change by around 80%, while salt concentration changed by only 6%. Such a small change in salt concen-tration produced a conductivity change that was too small to detect. By contrast, the Chapter 1. Introduction 10 changes due to gas concentration could be easily detected. Combining the theoretical phase diagram from Chapter 2 with models for the de-pendence of gas concentration on conductivity from Chapter 3, conductivity changes due to hydrate formation were predicted. Electrical measurements which depend on con-ductivity were used to detect hydrate in the experiment. Chapter 4 is devoted to the experiment. The discussion is divided into two parts. The first part deals with solubil-ity of carbon dioxide at different pressure and temperature values, while the aim of the second part is to detect hydrate formation by monitoring resistance measurements. The interpretation of the experimental results are discussed in Chapter 6, after first present-ing a theoretical discussion of electrical conductivity in aqueous solution and forward modelling of electrical resistance measurements. The main conclusions are summarized in Chapter 7. Chapter 2 Thermodynamics of Hydrate Formation 2.1 Hydrate Stability Hydrates form at low temperature or high pressure whenever there is a sufficient supply of gas. There are two main locations on Earth where these conditions are met. One, with low temperature but modest pressure, is permafrost, and the other, with temperature above freezing, but high pressure, is deep continental margins. Figure 2.1 shows typical profiles of hydrate phase stability in permafrost and in oceans. The dashed lines represent the temperature as a function of depth in the two environ-ments. The solid lines define the boundary of the hydrate-phase equilibrium (assuming methane gas). Pressure is converted to depth assuming hydrostatic conditions in both the water and sediments. Hydrate is stable when the temperature is below the solid line in Figure 2.1, so that the intersections of the solid and dashed lines bound the depth of the hydrate stability zones. For the permafrost the hydrate stability zone is indicated by shading. For the ocean environment, hydrate accumulates only in the sediments be-cause the sediments are needed to trap the hydrate. While it is possible for hydrate to form in the water column, the buoyant solid would be carried to the surface where it is dissociated. Solid line, as noted earlier, separates the region where hydrate is stable from the region where hydrate is not present. At lower pressures or higher temperatures we expect to find only aqueous solution and free gas (e.g. two phases). Two phases are also expected 11 Chapter 2. Thermodynamics of Hydrate Formation 12 Figure 2.1: Typical profiles of hydrate in permafrost (A), and in ocean sediment (B) [Sloan, 1990]. on the other side of the solid line. One phase is solid hydrate, while the second phase Chapter 2. Thermodynamics of Hydrate Formation 13 depends on the relative abundance of gas and water. According to the temperature-composition diagram for methane and water [e.g. Sloan, 1990] there are a number of different combinations of phases, depending upon the initial concentration of components. Possible combinations include liquid water and hydrate; liquid methane and hydrate; water vapour and hydrate; methane gas and hydrate. In oceanic sediment, where water is more abundant than gas, we expect to find hydrate and sea water with dissolved gas. The solid lines in Figure 2.1 define the P-T conditions where gas, water and hydrate are stable. These conditions are known as the three-phase equilibrium, while the regions on either side represent a two-phase equilibrium. The intersection of phase boundary and geothermal gradient in Figure 2.1, for example, in oceans is a depth where all three phases exist at the same time. In permafrost, there are two possibilities for gas, liquid and hydrate to exist simultaneously (e.g., two intersections of three-phase equilibrium and temperature curves). The goal of this section is to develop a quantative model for the phase equilibrium of a gas-liquid system, including both three-phase and two^phase equilibria. 2.2 Phase Equilibria When hydrate forms from a single gas component, then there are two components and three possible phases. Other gas components or salt in the water add components to the system. We begin our discussion by considering an n-component system with given pressure, temperature and overall composition. Our problem is to determine the com-position of each existing phase in equilibrium. This problem is sometimes called a flash calculation. According to the second law of thermodynamics, the natural state for any thermo-dynamic system is the one with the lowest possible energy. The quantities playing a Chapter 2. Thermodynamics of Hydrate Formation 14 role of the system energy are called thermodynamic potentials. There are quite a few of them depending on combinations of constraints. For example, if the thermodynamical process is described by whether entropy S and volume V, or entropy and pressure P, or pressure and temperature, the system energy is conveniently described in terms of whether the internal energy U, or enthalpy, H, or the Gibbs free energy, G. In chemical thermodynamics, with the process being simulated as a set of runs at fixed P and T, the system energy is represented by G, which can be written as follows G — U + PV - TS. (2.1) The Gibbs free energy is a state function, as it is a linear combination of state functions. It means that its total differential is exact. The total differential dG is given by dG = dU + PdV + VdP - SdT - TdS. (2.2) On the base of the combined first-and-second-law equation of thermodynamics, the total differential dll can be represented as dU = TdS - PdV + £ dW', (2.3) where J2 dW' is all non-PV (chemical, gravitational, electric etc.) work in the system. Combining equations (2.2) and (2.3), we write dG = VdP + SdT + Yl dW'. (2.4) This means that at constant P and T, the non-PV work, or energy available from the chemical process (the only one we consider) can be determined by the Gibbs free energy. The second law of thermodynamics, expressed in terms of G, requires minimum for the Gibbs free energy at equilibrium state. Each chemical component i in each of the possible phases j makes a predictable contribution to G. By adjusting the proportions of Chapter 2. Thermodynamics of Hydrate Formation 15 chemical components across all possible phases, a minimum in G is sought to determine the phase equilibrium. At the minimum energy state, a useful condition exists among the contributions to G for a particular chemical component. The partial molal free energies (contribution to G of one mole) of a component i is called the chemical potentials pi (free energy per mole). In equilibrium, the chemical potential p\ for each existing phase j must be equal; otherwise an adjustment of the phase proportions could lower G. Some flash calculations make explicit use of the equality of chemical potentials p\ or, equivalently, partial fugacities (see, for example, [Peng and Robinson, 1976]). The success of methods that use the criterion of equal partial fugacities relies heavily on the initial guess as to which phases exist and their estimated compositions. In the present calculations the equality of partial fugacities (or chemical potentials) follows as a consequence of the global minimization of G using a simulated annealing algorithm [Ingber, 1989]. Gibbs free energy is an extensive thermodynamic property, so it can be written in terms of chemical potentials p\, defined by P,T,nk M. as G = £ n ^ , (2.6) where n\ is the number of moles of component i in phase j. For nonideal solutions the chemical potential pi of a component i is customarily expressed in terms of its fugacity by pi = p°i+RTln^ (2.7) J i where p° and /? are the standard state chemical potential and fugacity of the component. The standard state may be freely chosen. For the case of a COi component it may be convenient to choose pure C02 gas at T and P conditions to define the reference state Chapter 2. Thermodynamics of Hydrate Formation 16 for CO2. The chemical potential of CO2 in hydrates or in aqueous solution will then be estimated relative to this chosen reference state. The choice of the standard state is usually determined by the particular problem. The suitable reference states for the CO2-H2O problem are described below. The fugacity is denned to have some important properties. First, and perhaps most important, the fugacity must obey relationship (2.8) for the chemical potential. Second, it is convenient to define the fugacity and pressure scales to converge for the case of ideal solutions. In this case the partial pressure Pi and fugacity /; of one mole of any ideal solution are identical at any T and V . Since very dilute solutions can be approximated as ideal solutions, we require Let's rewrite equation (2.8) as lim ^ -» 1 as Pi -> 0. (2.8) Pi-p°i= RTln^. (2.9) J i According to thermodynamic expression (2.4), the volume can be expressed as « - ( S L - ' Replacing both Vj and G with their partial molal counterparts (vi and pi), we obtain or dpi = VidPi. (2.12) The difference a between the partial molal volume Vi and the hypothetical volume of an ideal solution is given by RT * „ x Chapter 2. Thermodynamics of Hydrate Formation 17 Combining last two equations (2.13) and (2.14) with expression (2.9) and (2.10), we derive an important thermodynamic relationship ( 2 1 4 ) which provides a means of calculating the fugacities from experimental P — V — T meas-urements or by assuming some appropriate equation of state (EOS). For ideal gas mixtures the fugacity is equal to the partial pressure Pi of component i, which is defined by Pi = Pxi (2.15) where P is the total pressure and Xi is the mole fraction of a component i. Most real solutions, however, are nonideal, so the fugacity /j of a component i is often expressed in terms of the fugacity coefficient fa by fi = 4>iPxi (2.16) where deviations of fa from 1 indicate departures from ideal behaviour. Going back to the Gibbs free energy, we can rewrite G using equations (2.7) and (2.8) as G = ^U + RTlnJo)- (2-17) The first term in expression (2.18) is equal to energy of all components in the standard state. Since the standard state does not vary, the only part of G that we need to minimize is G* = I X ' J n i l ( 2 - 1 8 ) We now show how the phase equilibria are calculated by minimizing G* in (2.19). Let's consider a flash calculation for the solubility of C02 gas in liquid H20. The problem may be stated as follows. There is a closed system at temperature T 0 containing Chapter 2. Thermodynamics of Hydrate Formation 18 ni moles of carbon dioxide with partial pressure Pi and n2 moles of H20. Some of the C02 may be present in the gas phase and some may be dissolved in the liquid water. Similarly, the water exists as either liquid or vapour. The total pressure is equal to the sum of the partial pressure of C02 gas and H20 vapour, although the contribution of H20 vapour is very small when liquid water is stable. This means that the partial pressure of C02 is very nearly equal to the total pressure. The goal of the flash calculation is to find the mole fractions of the components in all phases. In solubility calculations the principal interest is the number of moles of C02 in the liquid phase. The standard states for both components are chosen to make the expression for G* as simple as possible. If the standard states correspond to pure components, pure C02 gas and pure H20 liquid at P-T conditions of interest, then the fugacity of C02 gas is approximately equal to the fugacity of a pure component because only a small amount of the water vapour is present in the gas. This means that fco2 ~ fco2 (2-19) and f9 Inffi- « 0. (2.20) JCO-i Therefore, we state that G* does not contain any terms that pertain to the gas-phase. The only contribution to G* comes from the C02 which is dissolved in the water and from changes in the energy of the H20 liquid due to the presence of the dissolved C02. This means that expression (2.19) can be written as G* = nlCQ2 J n & * + nlH20 ln&&, (2.21) JC02 JH20 where n ' H 2 0 i s the number of moles of liquid water, nlco<2 is the number of carbon dioxide moles dissolved in water. Chapter 2. Thermodynamics of Hydrate Formation 19 If water vapour is accounted for in the gas phase then this change makes only a small contribution to G* compared to that of the liquid phase (results differ within 3% percent). This suggests that the neglect of the water vapour is justified, so we will abandon it in subsequent problems. Substituting expressions for the fugacities (2.17) into (2.19) and dividing G* in (2.19) by N, the whole number of moles both C 0 2 and H20 in the system, we have r>** ,„xco2<f>co2 , „ 1nXH2p(j)H2o\ G = cti \ — h XH2oln—^IFS— J , (2.22) where cti is a liquid phase fraction, xco2 a n d XH2O are component mole fractions, defined by nC02 + NH20 Ctl = - — N ri xco2 = , , (2.23) nC02 "+" NH20 NLH20 XH20 = — ——i • nC02 "+" NH20 The coefficients of fugacity both for pure components and for components in mixture can be calculated from an EOS. There are two different equations of state which were used in this work: Peng-Robinson (PR) and Trebble-Bishnoi (TB) EOS. Both are semiempirical and based on the assumptions that pressure can be expressed as a sum a repulsion PR component and attraction PA component, P = PR + P A . (2.24) The repulsion pressure PR is written according to the van der Waals hard sphere equation as PT PR = —H, (2-25) v — o Chapter 2. Thermodynamics of Hydrate Formation 20 where b is related to the size of the hard spheres and the attraction pressure PA can be represent by PA = (2.26) where g(v) is a function of the molar volume v, and a(T) is a temperature-dependent interaction parameter. The function g(v) is dependent on only one parameter b for the PR EOS whereas two additional parameters are used in TB EOS. The function a(T) depends on a single parameter. So, PR is often called a two-constant equation of state, with TB is a four-constant model. Both EOS are defined for pure substances but also include rules for dealing with mixtures. If a mixture consists of substances with similar molecule sizes then the PR EOS works well because critical both compressibility and a hard sphere size are fixed in it. On the other hand, the TB EOS allows additional adjustable parameters in the case when molecule sizes are incommensurate. Consequently, the TB EOS may be successfully applied in cases where PR fails. More information on both EOS is in Appendix A. 2.3 Thermodynamic Model for Gas Hydrate Gas hydrates are crystalline compounds formed from mixtures of water and lower mo-lecular weight gases. Gas hydrates are members of a group of solids called clathrates [Powell, 1948] which contain two components: the gas and the host. In gas hydrates the host water molecules form a lattice structure and the gas molecules occupy the intersti-tial vacancies of the lattice without actually taking a lattice position. The gas molecules can rotate and vibrate within these vacancies or cavities. The molecules are too large to move freely through the lattice and remain localized in a single cavity. For most gases the water molecules form one of two hydrate structures (a third type Chapter 2. Thermodynamics of Hydrate Formation 21 Attribute Structure I Structure II small cavity per unit cell 2 16 large cavity per unit cell 6 8 water molecules per unit cell 46 136 small cavities per water molecule, V\ I 23 2 17 large cavities per water molecule, v2 3 23 1 17 Table 2.1: Structure of Gas Hydrates [van der Waals and Platteeuw, 1959] known as structure H has recently been discovered, but it will not be discussed here). The structure formed depends primarily on the size of the guest molecule. The unit cell of each structure contains several cavities of two different sizes. The ideal composition of the hydrate corresponds to occupation of all cavities. Since the cavities are the result of a regular lattice structure, the distribution and number of cavities are uniform in each unit cell of a particular structure. Table 2.1 gives the number of water molecules and number of cavities per unit cell of each of the two hydrate structures, which according to convention are designated structure I and structure II hydrate. The thermodynamic equations for gas hydrates were derived by van der Waals and Platteeuw [1959] on the base of classical adsorption statistical me-chanics. The assumptions of their theory are listed below: • an encaged gas molecule moves in a spherical cavity; Chapter 2. Thermodynamics of Hydrate Formation 22 • each cavity contains at most one gas molecule; • there are no interactions between the encaged molecules; • the gas molecules are sufficiently small to prevent distortion of the hydrate lattice; • the internal partition functions of the encaged gas molecules are equal to those of the molecules of the ideal gas. In the model, all polar forces are assumed to be embodied in the hydrogen-bonded hydrate lattice. The chemical potential of water in the hydrate structure consists of two terms. One is the chemical potential of water in empty hydrate lattice, while the other is the energy difference between filled and empty hydrate structures. The difference, A/z^f, between the chemical potential of water in the empty hydrate lattice, ^t£, and that in the filled hydrate lattice, p,%, is found [van der Waals and Platteeuw, 1959] to be A A £ = = -RT £ vm ln(l - J2 emj) (2.27) m j where vm is the number of cavities of type m per water molecule in the lattice and 9mj is the fraction of type m cavities occupied by gas component j. The occupancy 6mj depends primarily on the fugacity fj of the gas component j according to 9mj = C m j / i / ( l + £ °mifi) (2-28) I where C m j are known as the Langmuir constant. The summation over / involves all gas components. In the absence of the gas phase the fugacity of the gas component is defined by its value in the liquid phase. This fugacity is related in the usual way to the mole fraction of the gas component, xi, and the total pressure, P, by fi = <f>iXiP (2.29) Chapter 2. Thermodynamics of Hydrate Formation 23 where <f>i is the fugacity coefficient. The Langmuir constant accounts for the gas-i^O interaction in the cavity. Using Lenard-Jones-Devonshire cell theory, van der Waals and Platteeuw [1959] showed that the Langmuir constant is roo C(T) = 4ir/kT exp[-w{r)/kT]r2 dr (2.30) Jo where T is the absolute temperature, k is Boltzmann's constant, and w(r) is the spher-ically symmetric cell potential which is a function of the cell radius, the coordination number, and the nature of the gas - in t e rac t ion . Practically, however, the Langmuir constants can be calculated in the temperature range 260 — 300K from an empirical relation Cml(T) = (Ami/T) exp[Bml/T) (2.31) where Ami and Bml are adjustable constants which are determined by experiments. In this temperature range the maximum deviation between (2.31) and (2.32) is only 0.2% or less for almost all hydrate former gases [Parrish and Prausnitz, 1972]. The chemical potential of water in the empty hydrate lattice, is defined as a sum of the chemical potential of pure water at given pressure and temperature, p,w, and a term which describes the difference between these potentials, Apw: pi = pw + &uw (2.32) The most simple and straightforward fashion for evaluating the difference between the chemical potential of the unoccupied hydrate and pure water (whether it is solid or liquid) phases is given by [Holder, Corbin and Papadopoulos, 1980]: A/i(7>, P)w _ A /x(r 0 ,0)° fTf Ahw rpAvv RTf RTQ The first term on the right of equation (2.34) is an experimentally determined chemical potential difference between the unoccupied hydrate and pure water. The experimental Chapter 2. Thermodynamics of Hydrate Formation 24 results are referred to a standard temperature To, usually 0°C, and zero absolute pressure. The second term gives the temperature dependence at constant (zero) pressure, where Ahw is the difference in enthalpy between the unoccupied hydrate and pure water. The third term corrects the pressure to the final equilibrium pressure, where Avw is the difference in volume between the unoccupied hydrate and pure water. The temperature dependence of the enthalpy difference is given by where ACPui is the heat capacity difference between the empty hydrate and pure water phases. The heat capacity difference also depends upon the temperature but the linear approximation is sufficient, so It is usually assumed that all pure water solid phases have the same heat capacity. This Substituting equation (2.36) into equation (2.35) and integrating it over temperature, we get an analytical expression for Ahw. The result goes straight into equation (2.34) to yield an estimate of Apw. In evaluating Apw it is often assumed that the difference in volume between the unoccupied hydrate and pure water, Avw, is constant and equal to that at the reference conditions of 273.15 K and zero absolute pressure, At;°. Hence, an analytical expression for the difference between pure water and empty hydrate structure chemical potentials is (2.34) ACPw = AC°pJT0) + b(T - T0) (2.35) means that AC° vanishes when the difference involves the heat capacity of ice and hydrate. AA*(T, P) (2.36) RT [Ah°w - AC°pT° + l(T°Y} A < P R \T T°J RT Chapter 2. Thermodynamics of Hydrate Formation 25 Property Structure I Structure II b {J/mol.K2)* 0.141 0.141 A C ; (J/mol.K)* -38.13 -38.13 Ah°w (J/mol)* -4860.0 -5203.5 Ah°w (J/mol)** 1151.0 807.5 Au°w (J/mol)*'** 1264.0 8830 Avl (cm3/moi)* 4.6 5.0 Av° (cm?/moi)** 3.0 3.4 T° (K) 273.15 273.15 Table 2.2: Physical and Thermodynamic Reference Properties for Gas Hydrates [Englezos & Hall, 1994]. Listed values correspond to liquid water (*) and ice (**) at T < T°. Volume Av° in ice is taken from experimental data collected by Parrish and Prausnitz [1972]. Such values as A/x°, b, AC°, A/i°, Av° are obtained using a reference hydrate (Table 2.2). For Structure I (e.g. [Holder et al., 1980]), the reference hydrate is xenon hydrate for temperature below 0°C and methane hydrate for temperatures above 0°C, while for structure II (e.g. [Holder et al., 1980]), bromochlorodifluoromethane hydrate is the reference hydrate for temperatures below 0°C, and hydrates of natural gas mixtures are the reference hydrates above 0°C. Chapter 2. Thermodynamics of Hydrate Formation 26 Therefore, the chemical potential of filled hydrates, p^, can be written as a sum H% = pw + Apw + (2.37) where the first term is the chemical potential of pure water at given P — T conditions. Both of the correction terms are expressed analytically in equations (2.37) and (2.28). If the P — T conditions correspond to the liquid water phase coexisting with the hydrate then the chemical potential of pure water, pw, can be calculated on the base of whether PR or TB EOS, according to the scheme presented in the previous section. In the case of hydrates being in equilibrium with ice, the difference between chemical potentials of water in empty and filled hydrate lattice A/i^f, equation (2.28), stays the same, while the other correction for moving from a pure ice lattice to one of the unoccupied hydrate can be written exactly as relation (2.37), but with the appropriate for ice values of parameters which are given in Table 2.2, and AC° equal zero, as it was pointed above. 2.4 Optimization with Simulated Annealing The simulated annealing algorithm is based on an analogy between the annealing of solids and the problem of solving multivariate optimization problems [Kirkpatrick et al., 1982; Cerny, 1985]. For this reason the algorithm became known as "simulated annealing". Annealing is a physical process in which a solid in a heat bath, first, is heated up by increasing the temperature of the bath so that a transition from solid into liquid phase occurs, and then, is cooled by slowly lowering the temperature of the heat bath. Two factors are very important in the annealing process. Temperature must start with a sufficiently high value, and the cooling must be carried out sufficiently slowly. In this way, all particles arrange themselves into the low energy crystalline ground state instead of being frozen into a metastable amorphous structure. Chapter 2. Thermodynamics of Hydrate Formation 27 If the starting temperature corresponds to the liquid phase, with particles in a random arrangement, the cooling part of the annealing process can be described as follows. At each temperature value T, the solid is allowed to reach thermal equilibrium. Applying the principle of equal probability for equal energy states, it can be shown [Toda et al., 1983] that the probability of a macroscopic, state i with energy Ei is given by the Gibbs or Boltzmann distribution HB is the Boltzmann constant and N is the total number of states. As the temperature decreases, the Boltzmann distribution tends to favour low energy states. To simulate the evolution of a solid toward thermal equilibrium at fixed temperature T, Metropolis et al. [1953] proposed a Monte Carlo method which generates sequences of states of the solid in the following way. The current state of the solid is characterized by the position of its particles. A small, randomly generated, perturbation is applied by a small displacement of a randomly chosen particle. If the difference in energy, AE, between the current state and the slightly perturbed one is negative, then the process is continued with the new state. If AE > 0, then the probability of acceptance of the perturbed state is given by the Boltzmann factor (2.38) where Z(T) is the partition function exp{-AE/kBT). If the probability is less than a uniformly distributed random number from the interval (0,1) than the new state is rejected, and the state is retained in the original case, helping the system to pass through local minima. Depending whether the number is less or Chapter 2. Thermodynamics of Hydrate Formation 28 more than the acceptance probability, the new state is either retained or rejected (see Figure 2.2). This acceptance rule for new states is referred to as the Metropolis criterion. Following this criterion, the system eventually evolves into thermal equilibrium, or the probability distribution of the states approaches the Boltzmann distribution given by equation (2.39). Simulated annealing can be considered as an algorithm that continuously attempts to transform the current configuration into a slightly perturbed one often called a neigh-bour. This mechanism can be described mathematically in terms of Markov chains (e.g. the sequence of trials, where the outcome of each trial depends only on the outcome of the previous one). Therefore, the process of cooling can be modelled as a sequence of homogeneous Markov chains, where the length of the chain is the number of iterations used in order to reach thermal equilibrium at fixed T, with temperature being slowly de-creasing. The connection between Markov chains and annealing is important for proving the convergence of the annealing algorithm. The Metropolis algorithm can also be used to generate sequences of configurations of a multivariate, or combinatorial problem. In this case, the configurations assume the role of the states of a solid, while the cost function C and the control parameter c take the roles of energy and temperature, respectively. Minimizing the Gibbs free energy is the multivariate problem of interest in this study. The Gibbs free energy is non-linear function of pressure P, temperature T and composition. The composition is defined by the numbers of moles n\ of each component i in each phase j. For a hydrate forming from a single component there are two components (gas and water) and three phases (vapour, liquid and solid). For a prescribed P and T, the GFE depends solely on sets of n-. They are the configurations of the multivariate problem. Initially, the control parameter c is given a high value. Given a configuration a, another configuration b can be obtained by choosing at random an element from the neighbourhood of o (that corresponds to the Chapter 2. Thermodynamics of Hydrate Formation 29 small perturbation in the Metropolis algorithm). Then the difference between cost func-tions of a and b configurations is considered. If it is negative, configuration b is accepted; if it is positive, the Metropolis criterion is applied. This process is continued until the probability of distribution of the configurations approaches the Boltzmann distribution, with the Boltzmann factor being equal to (2.39) (recall that cost function C in our application is the Gibbs free energy of a gas-water mixture). The control parameter is then lowered in steps till the final "frozen" config-uration is reached. This configuration is taken as the solution of the problem. The lower value of the control parameter c, the less the probability of accepting a higher energy state. This means that at high temperatures, or high values of the control parameter, the simulated annealing algorithm looks only at the main features of the cost function surface, jumping from one local minimum into another. While as the temperature slowly decreases, the algorithm becomes confined to the global minimum. And eventually, as the control parameter approaches zero, the algorithm converges to the bottom of one of these basins. Using mathermatical results from the theory of Markov chains, it is shown [van Laarhoven, Aarts, 1987] that the algorithm asymptotically finds the global minimum with probability 1: Urn {Pr(X(n) e Rapt} = 1 (2.40) as lim c^ = 0, (2.41) 71 —*00 where Ropt is the set of globally minimal configurations, and n is the number of itera-tions. This means that asymptotic values of parameters governing the convergence of the algorithm have to be approximated: the length of an individual Markov chain (the Chapter 2. Thermodynamics of Hydrate Formation 30 number of iterations to reach equilibrium at fixed temperature), the initial values of the control parameter, the final temperature (the number of iterations n/ to consider the control parameter c small enough to satisfy requirement (2.42)) and so on. F l o w c h a r t f o r t h e A l g o r i t h m o f S i m u l a t e d A n n e a l i n g Initialize Annealing Temperature/Control Parameter, T Random Perturb a Variable Calculate dE Generate Random number P from a Uniform Distribution (0,1) If P < exp(-dE/T) Reject Perturbation E | € i P * Reset a Variable Then Accept perturbation Reduce T f^ Stop? ^ Yes No F i n i s h Figure 2.2: Backbone of the simulated annealing algorithm. Chapter 2. Thermodynamics of Hydrate Formation 31 The strategy for changing the control parameter and for creating the neighbour con-figurations is known as the annealing schedule. As simulated annealing does not specify it, there are many variants of reducing the annealing temperature and of organizing the generation function. With respect to the control parameter decrement, two types of SA algorithm can be distinguished: homogeneous when the algorithm is represented as a set of Markov chains, with each one generated at a fixed value of c, and inhomogeneous when the whole algorithm is described by a single chain, with c decreased in between transition. Simulated Annealing Algorithm may be of a Static or Dynamic nature, de-pending on whether the values of parameters governing the algorithm convergence are defined a priori or adjusted during the cooling process, respectively. One of the ways of adjusting the parameters (Markov chain lengths, their number and etc.) is to consider simulated annealing in terms of the direct thermodynamical analogy [Vidal, 1993], with the empirical entropy estimator being written as b ~ A + R' . where A is the number of accepted configurations and R is the number of rejected ones. The entropy stays constant at equilibrium. In the next section we will consider Static Homogeneous Simulated Annealing Algo-rithm. 2.4.1 Very Fast Simulated Annealing The very fast simulated annealing (VFSA) algorithm was designed by Ingber [1989]. One of the advantages of the algorithm is the possibility for different variables to have dif-ferent annealing schedules. That is essential for the multivariate problems with different physical constraints on range of the variable values. For example, in the case of mini-mizing the GFE, if the number of moles of carbon dioxide Nco2 a n d water NH2Q in the Chapter 2. Thermodynamics of Hydrate Formation 32 system are different then the ways that is perturbed ought to be different as well. In VFSA, for a variable aj with a limited physical range o-j € [AjtBj] a neighbor generation is obtained by aj+^aj + y ^ - i l , ) , (2.42) where yj G [—1,1] is represented by V i = sgn(Uj - 1)^1(1 + l/Tjfu^ - 1], (2.43) with Uj chosen from a uniform random distribution such that Uj € [0,1]. At the fixed control parameter, or annealing temperature T, the probability density of y is given by S t M = fi 2 ( N + r , ) L ( l + l/7V) s fl ^ t o ) , (2-44) where N is the whole number of the variables. The global minimum is obtained if an annealing schedule of variable Oj for the control parameter Tj is given Tj(k) = Tfexpi-Xjk1^) (2.45) where k is the number of an iteration, T? is the initial value of Tj, and Xj is an adjustable parameter to tune the algorithm to a special problem. In our special problem of minimizing the GFE of C02 — H20 mixture, for the vari-ables nXj, the annealing temperature schedule (2.46) stays the same, while the neigbour generation, given by equation (2.43), does not [Routh & Roy, 1996]. 2.5 Results The approach developed in previous section is applied to calculate the compositional dependence of the phase equilibria for a gas-water system. If we view the P-T cross sec-tion of the phase diagram as representing the thermodynamic conditions where different Chapter 2. Thermodynamics of Hydrate Formation 33 phases are stable, then the chemical dependence (e.g. an x-T cross - section) represents the composition of the different phases needed to ensure thermodynamic equilibrium. The phase diagram is calculated at fixed pressure values, so that the system composition is presented as a function of temperature. An important component of this calculation is the gas solubility in aqueous solution or, equivalently, the equilibrium concentration of dissolved gas. The solubility is a dimensionless measure of dissolved gas concentration, usually expressed as a mole fraction. C 0 2 - H 2 0 and C H 4 - H 2 0 mixtures are considered. Carbon dioxide is convenient to use in experiment because it forms hydrate at modest pressures, whereas methane is the most common gas found in naturally occuring hydrates. The total amount of gas included in both calculations is chosen to be small compared with the amount of water, but sufficient to yield either free gas or a hydrate phase (depending on P-T conditions). Figure 2.3 represents the solubility of carbon dioxide as a function of temperature at pressure of 2 MPa, which is representative of our laboratory conditions, while Figure 2.4 shows the results for methane at a pressure of 20 MPa, which corresponds to a water depth of roughly 2000 meters. The calculated solubility corresponds to the maximum (equilibrium) concentration of dissolved gas in the systems at different pressure and temperature conditions. Any gas in excess of the solubility limit is present as bubbles at high temperature or gets consumed by hydrate formation at low temperature. The solid line in Figure 2.3 corresponds to the amount of dissolved gas in equilibrium with free gas at given P-T conditions (the same convention is used in all figures). If the amount of dissolved gas is less than value indicated by the solid line at a given P and T, then the gas remains in solution. At high temperatures, the equilibrium solubility is often described by the Henry's law (see Chapter 3). Experimental studies show that the solubility decreases with a temperature increase. The theoretical calculations shown in Figure 2.3 are in agreement with the Chapter 2. Thermodynamics of Hydrate Formation 34 1 <55 279.0 Temperature (K) Figure 2.3: Carbon dioxide solubility as a function of temperature at a pressure of 2 MPa. The dashed line represents the solubility when hydrate is present, while the solid line corresponds to the case when hydrate is absent. The arrows represent the effect of increasing and decreasing temperature from point A. Chapter 2. Thermodynamics of Hydrate Formation 35 0 . 0 0 3 0 0 . 0 0 2 0 \-o o CO 0 . 0 0 1 0 0 . 0 0 0 0 2 7 4 . 0 2 8 4 . 0 T e m p e r a t u r e (K ) 2 9 4 . 0 Figure 2.4: Methane solubility as a function of temperature at pressure of 20 MPa. observed behavior. The dashed line in Figure 2.3 represents temperature-concentration conditions of two-phase equilibrium between hydrate and aqueous solution. Hydrate forms from an Chapter 2. Thermodynamics of Hydrate Formation 36 aqueous solution as soon as concentration of dissolved gas exceeds the value represented by the dashed line at given P-T conditions. The intersection of the dashed and solid lines corresponds to three-phase equilibrium between free gas, aqueous solution, and hydrate. For a system consisting of two components and three co-existing phases, there is only one degree of freedom according to Gibbs phase rule. This means that only pressure determines three phase equilibrium in our case; there is only one temperature at a fixed pressure where the three phases can co-exist. Departures from three phase equilibrium (due to altering the temperature or Concentration) cause the system move into a two-phase equilibrium. The possible two-phase equilibria are consistent with the phase rule as well; the equilibrium state depends upon both the prescribed pressure and temperature. Alternatively, concentration can be used as a free variable. The evolution toward a new equilibrium state, after changing temperature, is shown with arrows in Figure 2.3. According to Le Chatelier's principle, the general response to an increase in temperature is to cause a change where heat is absorbed (e.g. endothermic). A temperature decrease provokes a change in the opposite direction, so heat is expelled (e.g. exothermic). This means that an increase in the temperature, for example, from point A to point B (Figure 2.3), causes hydrate to dissociate (absorbing latent heat), which transfers gas back into the aqueous solution, and establishes a higher equilibrium concentration (point C). A decrease in the temperature (point D) promotes hydrate formation (releasing latent heat), with a corresponding depletion of gas from the solution, thereby achieving the lower value of the solubility (point E). The situation is similar at higher temperatures, where a two-phase equilibrium between aqueous solution and free gas occurs. But in this case, free gas is involved instead of hydrate. If the temperature decreases, then additional free gas gets dissolved, raising the solubility to a new higher value because of the heat of mixing. On the other hand, a temperature increase lowers the equilibrium value of gas concentration, and gas is exsolved from solution, producing free gas. Chapter 2. Thermodynamics of Hydrate Formation 37 Figure 2.5: Carbon dioxide solubility as a function of temperature at pressure of 2 MPa. The arrow between points F and G represent the effect of cooling the aqueous solution while the gas concentration is fixed. Chapter 2. Thermodynamics of Hydrate Formation 38 The solubility for both CO2 and C H 4 in Figures 2.3 and 2.4 shows the same general dependence on temperature, despite the different magnitude of solubility. The solubility has its maximum value at the temperature T 3(P) corresponding to the three-phase equi-librium. The solubility decreases nearly linearly at higher temperature, and almost expo-nentially at low temperature. The abrupt decrease in solubility for T<T 3(P) compared with that for T>T 3(P) has important implications for hydrate formation. Let's consider a sample of a liquid phase which is in equilibrium with the free gas at T>T 3(P) (say point F in Figure 2.5). The sample is isolated and cooled without allowing a change in the gas concentration. This situation is represented in Figure 2.5 by a line that connects points F and G. As temperature decreases, the fluid sample becomes undersaturated because the solubility increases as temperature decreases. We do not expect hydrate to form in the sample at T 3(P) because the gas concentration is below the equilibrium value at this temperature. More cooling is needed before the amount of gas in the sample becomes equal to the equilibrium value at point (G). Hydrate is stable in point G, but there is insufficient gas to grow hydrate at this stage. Further cooling lowers the solubility and promotes hydrate formation. The amount of hydrate produced depends on the depletion of gas from the liquid phase. A low initial gas concentration in the fluid sample may require cooling well below T 3(P) before hydrate forms. The amount of gas in solutions determines the temperature at which hydrate forms, but small concentrations of gas do not preclude the possibility of hydrate formation. It is evident that chemical equilibrium does not require the presence of free gas to produce hydrate (as is often suggested). The highest temperature at which hydrate forms corresponds to T 3(P) for three-phase equi-librium. In this case free gas is present. Lower concentrations of gas in solution only lower the temperature of hydrate formation. Methane is a major hydrate former. Consequently, it is important to consider a phase diagram of the methane-water system under thermodynamic conditions appropriate for Chapter 2. Thermodynamics of Hydrate Formation 39 marine environments. A typical seafloor temperature is 2.5 °C, and a typical water depth where hydrate occurs is 2000 meters, corresponding to a pressure of 200 bars. Different depths correspond to different values of temperature and pressure; temperature is cal-culated on the basis of a geothermal gradient, while pressure increases hydrostatically. Figure 2.6 illustrates the predicted solubility as a function of depth with the geothermal gradient at 50°C/km and 25°C/km. These calculations provide an estimate of the hydrate stability zones over depth. At large depth below the seafloor, hydrate is not present. The dissolved gas concentration in equilibrium with free gas is shown by a solid line. At smaller depth, hydrate is in equilibrium with aqueous solution when the gas concentration in solution is equal to the value given by the dashed line. Thus, the dashed line represents the amount of dissolved gas which must be exceeded in order for hydrate to form at a given depth (e.g. prescribed P-T conditions). Therefore, if the amount of gas in water within the hydrate stability field is less than that of the curve, then hydrate dissociates. If the gas concentration exceeds the equilibrium value, then this gas gets consumed by hydrate formation. Note that it is impossible to grow hydrate from the solution when gas concentration is at the equilibrium value because hydrate formation would lower the gas concentration in solution below the equilibrium value. However, as soon as there is any excess of gas in the system, hydrate forms and its volume fraction is determined by the amount of gas in excess of the solubility (e.g. dashed line in Figure 2.6). When free gas and aqueous solution are stable below the hydrate zone, the equilibrium gas solubility (solid line) decreases with depth. This behavior is in agreement with experiments, since layer depths correspond to higher temperatures, and temperature plays the major role in the thermodynamics. The depth of three phase equilibrium, where free gas, aqueous solution and hydrate stable simultaneously, corresponds to the position of the Bottom Simulated Reflector Chapter 2. Thermodynamics of Hydrate Formation 40 sea f loo r 100.0 \ \ 500.0 50 C/km 300.0 \ 25 C/km 700.0 0.0000 0.0010 0.0020 Solubility (mole fractions) 0.0030 Figure 2.6: Phase diagram for mixture methane - water for thermodynamic conditions encounted below the seafloor. The seafloor temperature is 2.5°C and the geothermal gradients are 25° and 50°C/km. The pressure is assumed to be hydrostatic below a water depth of 2 km. Chapter 2. Thermodynamics of Hydrate Formation 41 (BSR). This is a prominent acoustic reflection which is often detected in marine seismic surveys when gas hydrate is present. If there is free gas in the system, then the process of hydrate formation (according to the phase diagram) starts at BSR. As soon as hy-drate forms, it consumes all available free gas, causing an abrupt change in the seismic velocities. As a result, the BSR, which marks the base of hydrate layer, is well distin-guished on vertical seismic profiles. For the phase diagram calculated with a geothermal gradient VT=50°C/km, the depth of hydrate stability boundary is 287 meters, while that of VT=25°C/km gives the depth of 565 meters (Figure 2.6). Thus, using different geothermal gradients and different water depth above the seafloor, different values for the BSR depth can be obtained. In hydrate stability zones, the equilibrium concentration of gas dissolved in water decreases nearly exponentially toward the sea floor. The equilibrium solubility ceq as a function of depth h can be described as ceq(h) = c°e/h-h°VL where c°eq corresponds to depth hG, and L is the length scale for variations in ceq. The characteristic length L for smaller geothermal gradient of 25°C/km is 170 m, while that for 50°C/km is 84 m. Chapter 3 Theoretical Development Hydrate formation in porous media produces a number of physical changes which might be detected in laboratory experiments. Previous attempts to detect temperature changes due to latent heat release have proved unsuccessful [Rempel, 1994]. The goal of this chapter is to assess whether the addition of salt into the aqueous solution can cause observable changes in electric resistance when hydrate forms. Existing models for hydrate formation must be extended to allow for the evolution of salt concentration. 3.1 Mathematical Formulation of Hydrate Formation Process For dilute solutions (in which we are interested), the presence of salt does not significantly affect the heat balance and gas conservation in the system. The governing equations for heat and gas concentration in homogeneous porous medium were previously derived by Rempel [1994]. The pore fluid consists of water and dissolved gas, and both pressure and temperature increase with depth. At certain depth, with appropriate P-T conditions and a sufficient gas concentration, hydrate begins to form. The icy solid fills pores, consumes gas from solution and changes the temperature because of latent heat release. The conservation equations were derived using a fixed control volume of porous medium. An incompressible fluid with constant dynamic viscosity TJ is used to model the gas-liquid mixture. This fluid travels through the porous medium at velocity u 42 Chapter 3. Theoretical Development 43 which is called the transport velocity and is related to the interstitial velocity v by u = — h)v, where 4> is porosity, and h is hydrate volume fraction. Heat is transported by advection through the fluid and by conduction through the bulk material, with the latent heat L of hydrate formation acting as a heat source. The energy balance equation can be written as c f + U . V ^ V . ( K V T ) + ^ f , (3.1) where L is the latent heat of hydrate formation per unit mass, re is the effective thermal diffusivity K K = P/C/' with K and C/ being the bulk thermal conductivity and the fluid heat capacity, and C is the normalized bulk heat capacity * _ PfCf(f>(l -h) + phCh<t>h + p.C.jl - <t>) JfC-f ' in which the density p and isobaric heat capacity C of the fluid, hydrate, and sediment are represented by the subscripts / , h, and s, respectively. As hydrate forms, gas is transferred from the fluid into the hydrate. The mass fraction of gas Ch in the hydrate structure is assumed to be constant, while the mass fraction of gas c in the fluid is depleted during the hydrate formation process. Gas transport in the system is due to advection and diffusion down the compositional gradient. Therefore, conservation of gas requires [Rempel, 1994] - h)jZ + • Vc = V , [(1 - h)DVc] - ^(ch - c)^, (3.2) or <p pf ot where D is the diffusive coefficient. The last term in equation (3.2) represents a sink of gas due to incorporation of this gas into the hydrate structure. Chapter 3. Theoretical Development 44 When salt, an inhibitor for hydrate formation, is present the physical picture of the process is somewhat different. Salt ionizes in solution, and its ions interact with the dipoles of water molecules. The interaction provides a bond which is much stronger than that of the van der Waals forces which bind water molecules to an apolar gas molecule. Therefore, instead of forming hydrate, the water molecules organize a new structure with the salt, inhibiting hydrate formation and decreasing the solubility of the gas in the water (a phenomenon known as "salting-out"). To overcome the structural changes due to salt and allow hydrate to form, the temperature must be lowered below the equilibrium value in the absence of salt. This reduction in temperature is called the hydrate depression temperature. The value for aqueous solutions is approximated with the expression [Hammerschmidt, 1939]: A T = : , 1 0 0 M - M W " where A T is temperature of hydrate depression in °F, M is molecular weight of an inhibitor, W is concentration ( in weight per cent) of an inhibitor in the solution, and A is an inhibitor coefficient. For sodium chloride, A= 2320 [Makogon, 1981]. The equation of salt conservation can be easily developed from the equation for gas conservation, as the transport mechanisms for both salt and gas in sediments are quite similar. The only difference is that salt is not required for hydrate formation, and no salt is incorporated into the hydrate. This means that the equation of salt conservation can be written as " h ) % + \ u • V c * = V • [(1 - h)DJ?c*} + P-c*% (3.3) at <p pf at where c* is mass fraction of salt. The rate of consumption of gas is based on an empirical first-order model ft=-K(c-ceq), (3.4) Chapter 3. Theoretical Development 45 where K is the reaction-rate constant (K « 0.06 m i n - 1 for methane [Uchida, 1996]), and ceq is the equilibrium gas concentration at the system temperature. As the gas is consumed to form hydrate, the rate of its consumption is proportional to the rate of hydrate growth, but with opposite sign. Then, we can write a kinetic law for hydrate formation as where constant K\ is proportional to K from equation (3.4). The system of equations (3.1), (3.2), (3.3) and (3.5) complete the description of the problem. The unknowns include the gas concentration c, salt concentration c*, hydrate volume fraction h and temperature; the velocity u is typically treated as a prescribed quantity. Observations collected in regions with hydrate reserves, for example, the Cas-cadia accretionary margin, show that the transport velocity is less than velocity scale for thermal diffusion ut (for details, see [Rempel, 1994]). In this case the advective transport of heat can be neglected in the governing equations. The velocity scale for diffusion is defined in terms of the characteristic length scale I (the thickness of the sediment column) and the effective diffusivity K in the absence of hydrate by As both pressure and temperature gradients in marine sediments are predominantly vertical, the solution can be restricted to one spatial coordinate (e.g. the vertical coordi-nate z). The best way to rank the relative importance of different terms in the governing equations is to reduce them to dimensionless form. Detailed description of this process is given by Rempel [1994]. As a result, the system of equations in dimensionless form can be written as — = # i ( c - C e , ) , (3.5) (3.6) Chapter 3. Theoretical Development 46 ,.9c d ,. ,sdcx Ph, ~ -xdh .„ „ . ( 1 -"'s = e a i ( ( 1 - ' ' ) s i ) - ^ - c % - ( 3 ' 7 ) ,,3c* d ,, , , 3c*, ph ~ dh . „. ( 1 - * > f l r - « « « 1 - f t ) M ) + ^ ! ' S ' ( X 8 ) § = (3.9) where £ and 5 are dimensionless counterparts to t and z (the characteristic time r is the thermal diffusion time and characteristic length I is the width of the hydrate zone). Variables f, c, c*, S, e and ch are defined by f = T ~ T ° 7^-To' c -S = e = c - c0 Coo — CQ PfCfiT^-ToY D Ch Ch - c0 Coo ~ CQ Chapter 3. Theoretical Development 47 in which T^ — TQ, CQQ — CQ, — are the chosen scales for temperature, gas and salt concentrations. Here subscripts ^ correspond to initial values, T0 is the bottom boundary temperature, CQ is the equilibrium gas concentration at T 0 , and Cg is set equal to zero. The Stefan number S measures the importance of the latent heat release relative to the heat required to change the temperature of the porous medium, and e, a reciprocal of the Lewis number, indicates the relative efficiency of gas diffusion and heat. The diffusivities of gas and salt are assumed to be comparable. The previous study of Rempel [1994] included equations (3.6) and (3.7). Equation (3.8) for salt conservation and (3.9) for phase-change kinetics are new to this study. The dimensionless rate constant K2 depends on the choice of time scale for the problem. A diffusive time scale is adopted in this work I2 T ~ «(oy where the length I « Ira corresponds to laboratory conditions. With this value of I the coefficient K2 is 0.04. The values of the other dimensionless parameters used in equations (3.6) and (3.7) are taken from Rempel [1994]. The solutions of the system of equations (3.6) — (3.9) yield the time-dependent vari-ations in temperature, hydrate fraction, gas and salt concentrations over the vertical coordinate z. These solutions are obtained numerically using the method of lines. The boundary and initial conditions for a typical calculation are shown in Figure 3.1 ( the parameters are dimensionless). The initial condition represents a system containing gas, but no hydrate. The initial amount of gas and salt are assumed to be uniform. This means that there is no gas or salt flux at the boundaries at this time. Zero flux con-ditions on the salt and gas are imposed at the boundaries for subsequent times. The system is cooled from below to form hydrate. The cooling process can be described with fixed temperature boundary conditions where the bottom temperature is low enough to produce hydrate. As the system cools, hydrate forms and lowers the gas concentration, Chapter 3. Theoretical Development 48 while salt concentration increases due to the exclusion of salt from the hydrate. q=0 C=1 q=0 q=0 C*=1 q=0 T=1 T=0 Figure 3.1: Boundary and initial conditions for numerical modelling. The solutions at three different (dimensionless) times (0.1; 1; 5) are represented in Figure 3.2. As time advances (from the solid to dashed lines), the amounts of salt and gas in the system change. The rates of change of gas and salt are quite different; the results show, for example, that with 8% of hydrate fraction, the concentration of gas decreases by almost 90%, while the concentration of salt increases by only about 7% (the bottom of the dashed line). Figure 3.3 shows how hydrate formation affects the concentrations of gas and salt at a given point in the system (this point corresponds to Z=8 in Figure 3.2). The dashed line Chapter 3. Theoretical Development 49 40.0 f. 20.0 40.0 20.0 0.0 0.0 0.4 0.8 G a s concen t ra t ion 0.0 0.00 0.04 0.08 Hydra te f rac t ion Figure 3.2: Distribution of gas concentration, hydrate fraction, salt concentration and temperature (in dimensionless form) over depth at different time. represents hydrate fraction, while the solid and point-dash lines indicate the gas and salt concentrations. This figure shows that hydrate volume and salt concentration increase by Chapter 3. Theoretical Development 50 only a few percent at Z—8 as time advances. By contract, the gas concentration changes from C=1.0 at t=0 to Cw0.2 at t=5. According to this figure, the concentration of gas can suffice as an indicator of the process of hydrate formation. 3.2 Conductivity of Aqueous Solution If an aqueous solution contains gas or salt, then the conductivity of that solution will depend upon the component concentration. The reason is that dissolved gas and salt contribute to the number of ions in solution. More charged particles yields higher con-ductivities for the solution. We restrict our consideration to dilute solutions. In this case, the contribution of salt and gas can be treated separately if both components are dissolved in water. This is important because the expressions for the dependence of conductivity upon concentration for dissolved gas and salt are different. To get the dependence of conductivity upon the gas concentration in a dilute solution, we start with the definition of conductivity. The constant of proportionality between the current density J and electrical field E in so called the point form of Ohm's law where a is the conductivity of the material. When a gas such as C 0 2 dissolves into water, some of this gas dissociates into positive and negative ions. If an electrical field is applied, cations (positive ions) will be acceler-ated towards the negative pole of the field and anions (negative ions) - to the positive pole. This acceleration is opposed by viscous drag which limits the maximum velocity to which the ions can be accelerated. The terminal velocity of the ions obtained with a unit electric field E (V/m) is defined as their mobility pa (ra2/V'/'sec): J = aE, (3.10) (3.11) Chapter 3. Theoretical Development 51 5.0 4.0 3.0 .i 2.0 1.0 0.0 0.0 0.2 0.4 0.6 Concentration (fraction) 0.8 1.0 Figure 3.3: Gas concentration, hydrate fraction, and salt concentration as a function of time at Z=8. Mobility is a function of both temperature and concentration in the solution. Decreasing the temperature decreases the viscosity of water, permitting higher terminal velocities Chapter 3. Theoretical Development 52 for the same voltage gradient. If a solution contains a high concentration of ions Ci, the motion of one ion can be inhibited by the other ions close to it, reducing terminal velocity. The current density J depends on the number density and the velocity Vdri of ions travelling through the solution. If the charge of the ions is then the current density can be written as J = niVdriqi. (3.12) Manipulating equations (3.10), (3.11) and (3.12) yields an expression for the conductivity cr of a system containing ions with different charges qi and number densities rij o- = Y^niqiPi, (3.13) i where summation is over all types of ions in the system. The expression for the charge of an ion is ft = ^ > (3-14) A where Zi is the valence of the ion, F is Faraday's constant, 9.648 • 104C/mol, and NA is Avogadro number, 6.0217 • 102 3 particles/mole. Combining (3.13) and (3.14), we get the expression for the conductivity of a dilute solution of gas, aWta. = 10s F^ZidiM, (3.15) i where Ci is the ionic concentration, in mol/L. In equation (3.15), there are two variables that depend on temperature. One is the ionic concentration Ci and the other is mobility Ui. The changes in Ci with temperature may arise from changes in the solubility of ions in the solution, while changes in pi reflect changes in the viscosity of water. For the temperature dependence of the ion mobility, we use the first two terms of the Taylor series u(T) = p° + ^(T- 2b°C) = u°[l + f39(T - 25°)] (3.16) Chapter 3. Theoretical Development 53 where p° is the ion mobility at temperature 25°C. Experimental data is available on the mobility of ions H30+ and HCO^, which arises from dissolution of C02. In a dilute solution at temperature 25°C, the mobilities of H30+ and HCO3 are equal to 4.6 • 10 - 8 and 36.2 • IO - 8 meters per second/volts per meter, respectively [Carmichael, 1987]. The coefficient {3g is found to be 0.033 using experimental data collected in this study. This estimate is based on a comparison of theoretical and experimental values of the conductivity of water with dissolved carbon dioxide as a function of pressure and temperature. The theory is presented later in this chapter and the experiment is discussed in chapter N. Ionic concentration Ci is proportional to the gas concentration in the solution. We derive here an expression for carbon dioxide. For dilute solutions of gas there is a phe-nomenological relationship between the partial pressure of the gas and its dissolved mole fraction. The relationship is called Henry's Law: PB = XBH(T), (3.17) where PB is the partial pressure of the solute B, XB is its mole fraction in the solution and H(T) is the Henry's Law Constant which depends on the particular combination of solvent A and solute B as well as the temperature in the system (see Table 3.1), The mole fraction of B is related to the mass W and molar weight M of components B and A by Wn Wn WA XB = TT/(TT + Tr) (3-18) MB MB MA Values of Henry's Constant over a wider range of P-T conditions have been obtained from experimental data [Stewart & Munjal, 1970]. It is also possible to calculate the gas concentration in solution (see Chapter 2). If the mole fraction of the dissolved gas (solubility) is known, the values of ionic Chapter 3. Theoretical Development 54 Table 3.1: Henry's Law Constant for C02—H2O system T(°K) K(MPa) 293 146 298 170 333 340 concentrations can be found by determining the equilibrium constants of the reactions that produce ions. Carbon dioxide dissolves into water to form carbonic acid, which subsequently ionizes. The set of reactions that take place in the H20—C02 system can be written as C02 + H20 ^ H2C03, H2COz + H20 ^ #30+ + HCOl. Equilibrium constants Ki for these two reactions (i = 1, 2) are given in terms of relative activities of the reagents by K I = QH2CO3 ; ( 3 1 9 ) aC02 ' aH20 K2 = aJ^l^S2L. (3.20) aH2C03 ' aH20 Water is by far the most abundant species in dilute solutions and its relative activity is almost the same as in pure water, aH2o = 1-Hence, the activity of the solvent water can be ignored in the calculation of Ki. For electrolytes in dilute aqueous solutions, the relative activities of species are assumed to be equal to their dimensionless molar concentrations: «H2CO3 ~ cH2COJc° = [H2CO,]r, (3.21) Chapter 3. Theoretical Development 55 « H 3 O + ~ cH30+/c° = [H30+}r, (3.22) aHco- « CHCO-l c° = [HCO;)r, (3.23) where c° is a reference value, equal to 1 mol/L in our case. Therefore, the equilibrium constant may be related to the concentrations of the reagents at equilibrium, e.g., * = w - ( 3 - 2 4 ) K> = [H2C03}r • ( 3 - 2 5 ) The equilibrium constants vary with the temperature for the exothermic and endothermic reactions according to LeChatelier's principle. This means that as temperature increases, the equilibrium constant decreases for exothermic reactions, and the temperature influ-ence is opposite for endothermic ones. The equilibrium constant K i for the formation of carbonic acid from dissolved C02 is small. Experiments at 25°C indicate that #i = [J^L = 2 • lO" 3 . (3.26) The second equilibrium constant K2 for the ionization reaction is K,= M°rjff°^=2.W-'. (3.27) Therefore, the relative concentration of carbonic acid at equilibrium is [H2C03]r = [C02]r • Ku (3.28) while the concentration of ions [H30+]T = [HCO$ }r are given by [H30+]r = [HC03]T = JK2 • [H2C03]r = ^ JK2 • Kx • [C02]r. (3.29) This means that if the concentration of dissolved gas and the values of the equilibrium constants are known, then the concentrations of ions can be easily calculated. Chapter 3. Theoretical Development 56 Since the conductivity of an aqueous solution of water and C02 depends on the concentration of ions C; according to (3.15), and the ion concentrations are denned by (3.29), cr has the following propotionalities o ~ C ~ VX ~ \fP, (3.30) where C is the concentration of H20+ and HCOJ ions, X is the solubility of C02 in water and P is the partial pressure of the gas. For example, if 1 atmosphere of air (with 1% C02) is replaced with 1 atmosphere of pure C02, the partial pressure of C 0 2 increases by 100 and the conductivity increases by factor of 10. Consequently, the conductivity crWgaB of a dilute aqueous solution of carbon dioxide can be estimated on the basis of (3.15) and (3.16), with the ionic concentration given by equation (3.29), as aWgas = 103 F JK2 -Kx-X £ > ° [1 + M T ~ 2 5 ° ) 1 . ( 3 - 3 1 ) i where p° = p(25°C). An empirical relation is used to calculate the conductivity of a dilute aqueous solution of salt. For concentrations c in the range 10 _ 4 M to 0.1M, the conductivity has the following form [Worthington et al., 1990] log (ac) = 0.933292 + 0.84597lZo£ c - 0.0910632(Zo£ c)2 (3.32) -0.227399(Zo(/e)3 + 0.227456(/o^c)4 - O.381O35(Zo0c)5, where c is in M, (e.g. mol./L). This expression was shown to fit the measured data to within 0.1 percent. The dominant effect of temperature upon the conductivity of the salt solutions is through the effect of the viscosity of water. There are a number of different expressions Chapter 3. Theoretical Development 57 for this temperature correction. For modest changes in temperature we use the first two terms of the Taylor series for cr(T), as we do for a solution of gas. Consequently, o(T, c) = o0(c) [1 + 0,(T - 18°)], (3.33) where 00(c) is the conductivity of a solution at 18°C, evaluated using (32), and the value of parameter (}„ for brines is 0.022 [Heiland, 1940]. 3.3 Conductivity of Porous Media The electrical conductivity of clean (clay-free) sands and sandstones can be described using Archie's law [Archie, 1942]: af = <f>mSZaw (3.34) where 07 is the bulk conductivity, aw is pore fluid conductivity, </> is porosity, Sw is the water saturation, n and m are the cementation and saturation exponents. Archie [1942] found that m was equal to 1.3 for clean, unconsolidated sands. The satuaration exponent n, for saturations ranging from 0.15 to 1.0, was determined by Archie [1942] to be approximately 2.0. We assume that the pore fluid is a dilute solution of both salt and gas. This means that the whole conductivity can be presented as a sum of conductivities of dilute solutions of salt and gas. As hydrate forms, the porosity is affected by the hydrate amount because hydrate fills the pore space and reduces the porosity. If the initial porosity prior to hydrate formation is (sand porosity), then the effective porosity after hydrate forms is (j> = 4>o(l — h), where h is the pore volume occupied by hydrates. When the porous medium is fully saturated, Sw — 1, and the equation (3.34) reduces to: <Tf = [ M 1 - J*)]1"3 K . a u + aw9as)- (3.35) Chapter 3. Theoretical Development 58 Equation (3.35) is used in both mathematical modelling of hydrate formation to predict changes in the bulk conductivity and the experimental studies to test these mathematical predictions. Chapter 4 Experimental Work An important result of the thermodynamic calculations in Chapter 2 was the chemical equilibria of hydrate in the presence of aqueous solution. This result suggested that hydrate can grow from an aqueous solution without the need of any free gas. To test this idea, an experiment was designed to monitor the growth of hydrate from aqueous solution. The objective was to use electrical measurements to detect hydrate formation. The fact that hydrate zones correspond to regions of low conductivity is known from well log analyses [Collett, 1993]. This means that electrical properties of a medium where hydrate forms are affected by the formation process. The conductivity is a general elec-trical property of materials, while the measured resistance also depends on the geometric form of the material and the electrode configurations. If C is a factor describing the form, spacing and material properties of the electrodes, then the relation between resistance R and conductivity a is R=£ (4.1) a Therefore, resistance measurements and a known geometric factor C can be used to determine the conductivity. Conductivity of dilute aqueous solutions in porous media is a function of temperature, ion concentrations and porosity. The empirical relationships describing conductivity dependence upon temperature, ion concentration and porosity are discussed in Chapters 3. As hydrate forms, temperature, ion concentration and porosity all change. The pro-cess of hydrate formation affects temperature due to the latent heat effect, while gas 59 Chapter 4. Experimental Work 60 concentration and porosity are altered by gas consumption and pore plugging. If salts are also present in the system, then their concentration is altered by the formation pro-cess as well. Salt ionizes in solution and the resulting ions interact with the dipoles of the water molecules. These ionic interactions are stronger than those due to van der Waals forces, which bind gas molecules in hydrate structure. Therefore, the temperature of hydrate formation is lowered in order to overcome structural changes produced by salt ions. But as soon as formation triggered, salt gets excluded from the structure. Any increase in the salt concentration in solution would indicate the process of formation. Temperature, ion concentrations and porosity change through time and these changes can be calculated on the basis of the governing equations discussed in Chapter 3. There-fore, the characteristics change due to hydrate formation and how these changes affect conductivity can be determined. This means that we know how the bulk conductivity changes due to hydrate formation. To facilitate the experimental detection of hydrate formation with electrical methods, salt was considered as an additive in the system. The solution of the governing equations (Chapter 3) showed that a typical hydrate fraction of 8% caused an increase in salt con-centration of 6%, and decrease in gas concentration of 80%. Combining these solutions with empirical expressions for the conductivity of dilute solutions of gas and salt based on Archie's law (see Chapter 3), the conductivity change due to hydrate formation is cal-culated. Figure 4.1 represents the predicted change in conductivity as a function of time. The dashed line represents the contribution of the changing salt concentration, while the solid line represents the contribution of the changing gas concentration. Both curves also include the effect of temperature and pore plugging on conductivity. The conduct-ivity change caused by increasing of the salt concentration is actually negative because of influences of both temperature and porosity. However, in the case of a changing gas Chapter 4. Experimental Work 61 concentration, the changes in conductivity are reinforced by the changes in temperat-ure and porosity, so the change is much more profound. In fact, the presence of salt in solution makes the conductivity change due to hydrate formation smaller that it would be without the salt. Consequently, salt was omitted from the experiments. Resistance 4.0 2.0 o.o 0.0 0.4 Conductivity (dimensionless) 0.8 Figure 4.1: Conductivity for dilute solutions of salt and gas as a function of time. Chapter 4. Experimental Work 62 measurements were used to monitor the process of hydrate formation in our experiment. As it was shown, the resistance mainly depends on the solution temperature and gas con-centration in it. The amount of gas dissolved in the solution is limited by its equilibrium value. The solubility over the range of temperatures used in our experiment is denned by the calculated phase diagram, in Figure 4.2. Also shown are three possible paths that the system might follow as it is cooled to produce hydrate. For example, suppose that the system starts in an equilibrium state at point A. As cooling, one path is to follow the equilibrium concentration at each temperature (A-E-C in Figure 4.2). This situation is obtained if the cooling is slow and there is a sufficient supply of gas (e.g. a chemically open system). The other two scenarios are for a closed system with a fixed amount of gas inside. In this case, hydrate begins to form when the temperature reaches the equilibrium value at fixed gas concentration (e.g. point D in Figure 4.2). Further cooling produces hydrate in one of two ways. If the system remains in thermodynamic equilibrium, then the gas concentration will decrease along path DC as the system is cooled. On the other hand, non-equilibrium effects may require overcooling, followed by a rapid change in gas concentration as hydrate quickly forms (e.g. path DBC). The resistance as functions of temperature for the different paths are represented in Figure 4.3. The difference in the resistance depends mainly on the amount of gas that is dissolved in the solution. At a given temperature, resistance of a closed system with fixed amount of gas (AD path), is higher than that of an open system with an equilibrium gas concentration (AE-EB path). The opposite situation occurs when the amount of gas exceeds the solubility (paths CB and DC, respectively). Line BC represents the resistance change due to a change of the dissoved gas concentration which results from hydrate formation. Therefore, resistance as a function of temperature for the three paths all differ. The dashed line AEC corre-sponds to the open system which reaches its equilibrium state at each temperature, while the solid line ABC occurs when there is a fixed amount of gas as well as kinetic barriers Chapter 4. Experimental Work 63 E Temperature Figure 4.2: Three possibilities for reaching hydrate stability zone from an aqueous solution at constant pressure: AEC, ADC, and ABC. to forming hydrate. The line composed of solid AD and dashed DC represents the closed sytem when the barriers are absent. Both AEC and ADC paths provide a smooth change in curvature due to hydrate growth in the system, while in ABC curve the moment of Chapter 4. Experimental Work 64 hydrate formation is well distinguished. The latter path could be an ideal indicator of hydrate formation if kinetic barriers are important in the experiment. In fact, the abrupt change in resistance was observed in the experiments and this change is proportional to the amount of gas depleted in the solution. Moreover, the volume fraction of produced hydrate is proportional to this change as well. The porous medium plays an important role in determining whether the system is closed or open. The sediments are confined to the bottom few centimeters of the ap-paratus, while the remainder of the overlying volume is filled with distilled water. As C 0 2 gas dissolves into the distilled water at the top of the apparatus, the fluid becomes heavier than the underlying pure water and convectively unstable. Vigorous convection rapidly mixes the gas into the water column. However, the density differences due to gas concentration are not sufficient to force convective mass transport in porous medium (see Appendix B). Hence, the gas transport is provided by diffusion. Because of the slow diffusive transport in porous medium, the region of sediments behaves like a closed system over the course of the experiment. Sediments are also needed in order to trap hydrate. Although it is possible for hydrate to form in water, the resulting solid is buoyant and it will be carried to the surface. If thermodynamic conditions at the surface are outside the hydrate stability zone, then the hydrate is dissociated there. A final function of sediments in our experiment is to lower the kinetic barriers for nucleation of a crystal. In general, there are two possibilities for crystal formation: either homogeneous or inhomogeneous nucleation. Homogeneous nucleation requires the crystal nuclei to obtain the critical size, while the inhomogeneous nucleation assumes some additive which aids crystal growth. The activation energy, an energy to overcome kinetic barriers (which is associated with overcooling in hydrate phase diagram), is usially much higher for homogeneous nucleation than for inhomogeneous nucleation. Thus, the Chapter 4. Experimental Work 6 5 Temperature Figure 4.3: Schematic dependence of resistance on temperature for different gas supply schedules: AEC, ADC, and ABC. inhomogeneities (in the form of sediments) in the system helps both to overcome kinetic barriers and diminish nonequilibrium effects. As a result, the system should stay closer to the conditions represented by the phase diagram, based on equilibrium thermodynamics. Chapter 4. Experimental Work 66 4 . 1 E x p e r i m e n t a l A p p a r a t u s The experimental apparatus came as an inheritance from an earlier experiment which was set up to test an analytical model of hydrate formation [Rempel, 1994]. The reaction chamber is a 19 mm thick, clear acrylic tube, 0.7 m long with inner diameter 0.14 m, capped at each end by an aluminun lid in which coolant is circulated to maintain a constant temperature. Carbon dioxide gas can be fed under pressure into the top of chamber. A pressure gauge situated on the gas cylinder allows different values of carbon dioxide partial pressure to be established. Therefore, a chemical mixture of H2O-CO2 can be subjected to different P-T conditions. The temperature in the system is measured with resistive thermal devices (RDTs), utilizing the temperature dependence of platinum electrical resistance. The circuit for temperature measurements is discussed in detail in Rempel [1994]. A second circuit is used to measure the electrical resistance. Figure 4.4 shows this circuit, which includes four electrodes E1-E4, a known resistance R and a voltage source V, or a wave generator. The current is denoted by I. Electrodes 1 and 4 are current electrodes, while 2 and 3 are used to measure the potential difference. Al l four electrodes are located within the porous medium. The series circuit shown in Figure 4.4 is closed, which means that the current is constant. In order to find the current, we measure a voltage drop AVR across the known resistance: 1 = ^ (4.2) Monitoring the voltage drops AV2_3 across the electrodes 2 and 3 provides a measure of the bulk resistance of the porous medium, defined by i?2-3 = ^ (4.3) Using the known current in (4.3) gives ( 4 4 ) Chapter 4. Experimental Work 67 i 4 1 — o I I V © -Figure 4.4: Circuit for resistance measure-ments. Electrical measurements are simple in principle, but sometimes difficult in practice. In a real system, there are effects such as electrode polarization, which can be caused by electrochemical action of a current. In addition, there can be background currents that are unrelated to the imposed voltage source. In order to overcome these problems, the voltage source generates square waves, with voltage as well as a current changing direction. Averaging over both polarities diminishes the effect of background currents. The result provides a better estimate of resistance over time for monitoring the process of hydrate formation. Direct measurements of R2-3 for a square wave input voltage showed that the re-sistance decreased slowly with time after the voltage changed sign. The situation is schematically shown on Figure 4.5. The phenomena is characteristic of induced polar-ization, but the time variations are probably too long. These variations in current and Chapter 4. Experimental Work 68 resistance could continue for 30 seconds or more. We interpret these variations as an electro-chemical effect. Based on these results, the square wave generation frequency was V R I Figure 4.5: Voltage, resistance and current as a function of time. set at 5Hz. This value satisfied two requirements on resistance measurements. First, it was desirable to make the resistance measurements for both polarities as constant as possible, so that averaging would eliminate unwanted electrochemical and background currents. At higher the frequencies, the resistance values were nearly constant because electrochemical effects caused changes over tens of seconds. On the other hand, at high frequencies, the number of measurements per cycle was not enough for reliable averaging. Consequently, 5Hz represented a good compromise. The known resistance R is 15 kOhm. This value was chosen because i t was comparable tO R2-3-Chapter 4. Experimental Work 69 4.2 Gas Solubility Experiments Before performing experiments on hydrate formation in porous media, some experiments were conducted to determine the dependence of resistance upon gas solubility. Both vertical and horizontal electrode locations were used. The goal of these experiments was to test the apparatus, to find out about the mechanism of the gas transport in the water, and to observe resistance as a function of temperature. On the basis of these observations, the geometrical constant in (4.1) was determined. Two sets of experiments with different electrodes locations were conducted. The first set was performed for a vertical electrode location without any porous medium. A plastic piece of equipment was constructed to hold vertically three RTDs and four electrodes, with a one-inch separation between electrodes. In these experiments, at room temperat-ure, pressure was increased from 45 psi to 200 psi, resulting in increase of gas solubility. The experimental results are represented on Figure 4.6, where resistance is plotted as a function of time. At fixed temperature, the changes in resistance are determined solely by the amount of gas dissolved in water. The change in resistance shows how fast gas is sup-plied into the current electrode area when there is no porous medium. This must occur by convection because the time scale for chemical diffusion over the entire 0.7 m length of the chamber is nearly 30 years. Figure 4.7 shows a similar result for resistance, with ther-modynamic conditions maintained at P=50 psi, T=14.5°C. This particular experiment corresponds to horizontal location of the electrodes kept 8 mm above the aluminum bot-tom. This experiment was performed over 22 hours before reaching a steady resistance value. Both figures show that gas must be transferred into pure water by convection. Convection takes place in the system of CO2-H2O, as gas starts dissolving at the top of the apparatus because the density of gas-riched water is greater than that of pure water. A more quantative treatment of this convective transport is given in Appendix B. Chapter 4. Experimental Work 70 8000.0 6000.0 J o cu o S3 m 4000.0 o H 2000.0 0.0 I ^ 4 5 p s i 0.0 7 0 p s i 2 0 0 p s i 100.0 200.0 300.0 Time (min) 400.0 Figure 4.6: Average R 2 _ 3 as a function of gas solubility at room temperature. 4.3 H y d r a t e F o r m a t i o n To carry out experiments on hydrate formation, a medium grained sand was added into the system. Because of the slow diffusive transport in porous media, the region Chapter 4. Experimental Work 71 Figure 4.7: Change in resistance over time due to increase in gas concentration at pressure 50 psi and room temperature. of sediments behaves like a closed system over the course of experiment. How far the real thermodynamic system behavior is from that of the closed one is determined by the height of the sand column. A deep layer of sediments implies a long diffusion time, which Chapter 4. Experimental Work 72 is typically longer than the duration of the experiments. The first set of the experiments was arranged for the vertical electrode distribution, with the sand column height around 0.3 m and the first current electrode located 0.01 m from the bottom. In order to diffuse gas into the sediments to the place where the resistance measurement were made, required almost two weeks. (The apparatus was kept at fixed high pressures over this interval of time). According to our phase diagram, the concentration of dissolved gas at a prescribed pressure determines the temperature at which hydrate forms, if it is possible at all. The more gas dissolved in water, the higher the formation temperature. If the amount of dissolved gas is tiny then hydrate may not form at all (see the diagram). The first set of experiments (with the vertical electrode locations) was not successful. The temperature around the current electrodes was not low enough to form hydrate, as the distance from the cooled bottom was too great. This electrode position was originally chosen because of the possibility of analytically modelling the resistance measurements (see Chapter 5). However, these resistance measurements could be used with forward modelling to determine the coefficients 0g for the temperature dependence of conductivity. The horizontal electrode array was used in the second set of experiments. Since all of the electrodes were at the distance of 8 mm above the cooled aluminum bottom, the temperature needed to form hydrate was reached more easily. In addition, the distance between the electrodes and the top of the sediments was three times smaller, so the time for gas to diffuse into the sediments was reduced to a few days (instead of a few weeks). A series of experiments were performed with the horizontal array. Individual experiments were carried out with different amounts of gas in sediments, so that different formation temperatures were obtained. Figure 4.8 shows the measured resistance as a function of time for one experiment. The whole time corresponds to approximately 15 hours. Chapter 4. Experimental Work 73 Figure 4.8: Average R2_3 as a function of time at pressure 300 psi, as temperature reduced from 6°C to 2.2°C. Figure 4.9 shows the corresponding temperature of the three RTDs as a function of time for the same experiment. The last figure (4.10) is resistance as a function of Chapter 4. Experimental Work 74 1.0 1 • 1 • 1 — 1 — 0 . 0 1 0 0 . 0 2 0 0 . 0 3 0 0 . 0 Time (min) Figure 4.9: RTD temperatures around current electrodes as a function of time (pressure 300 psi). temperature. Resistance changes at the constant temperature of roughly 2.2°C due to a decrease in the amount of gas dissolved in water. This decrease resulted from the process Chapter 4. Experimental Work of hydrate formation. The interpretation of these results is given in Chapter 6. 75 Figure 4.10: Average R2.-3 as a function of a middle temperature (pressure 300 psi). Chapter 5 Conductivity Profile Through Electrical Potential Measurements 5.1 Governing Equations for D C Resistivity An estimate of the electrical conductivity profile can be obtained by inverting electrical potetial measurements. These measurements are collected at electrodes which are placed inside the porous media. If the electrode are inserted in a way that the distance between them is much less than the cylindrical radius of the apparatus, then the boundary of the cylinder can be ignored and electrical potentials from a ID theory will suffice. In this case the conductivity a(z) varies only with vertical position z and each layer of constant a is assumed to be infinite in horizontal extent (see Figure 5.1). To reconstruct the vertical distribution of conductivity cr(z), electrodes are inserted along the 2-axis which is positive upwards, perpendicular to the layers. A steady electric current is input into one of the electrodes and the potential is measured at the others. For convenience the coordinate origin in the x — y plane is chosen to coincide with the position of the source electrode. Since the electrical potential 4> is symmetric around the source electrode the problem can be solved in the Hankel Transform Domain. A discrete inverse Hankel transform is then used to express the solution in terms of the spatial coordinates. The point form of Ohm's law, the continuity of current J equation and the represen-tation of the electrical field intensity E as a negative gradient of the potential (f> forms 76 Chapter 5. Conductivity Profile Through Electrical Potential Measurements 77 °(z3) o(z2) a(z,) Figure 5.1: Schematic position of electrodes, denoted by circles, in a layered conductive material. the system of equations for the DC resistivity problem. These equations are given by J = oE V • J = dp/dt (5.1) E = where p is net charge density. In the case of a charge source with intensity / the system —• of equations for potential 4> and current density J is V • J = -I8(f-f8) = -J/cr (5.2) Chapter 5. Conductivity Profile Through Electrical Potential Measurements 78 where f3 is the position of the point source. The explicit Cartesian components of system equations (5.2) can be written as: dJx/dx + dJy/dy + dJz/dz = —I6(x — xA)6(y — yS)8(z — zs) d(j)/dx = — Jx/a (5.3) d(j)/dy = -Jy/a d(j)/dz = — Jz/a 5.2 Fourier Transformed Equations Since the porous media extends to infinity in the horizontal directions, Fourier transforms may be used to eliminate the horizontal directions from the system of equations (5.3). Defining the Fourier transforms by roo roo Ji(kx,ky,z)= f°° f°° Ji(x,y,z)eikxXeik»ydxdy roo roo <j)(kx,ky,z) = / <f)(x,y,z)eikxXeik>ydxdy J—oo J—oo the transformed equations become -ikx Jx - ikyJy + dJjdz = -Ieik*x-eIK*Y'6(z - zs) d(f>/dz — —Jzjo —ikx(j) = —Jx/o —iky(j) = —Jy/a (5.4) Chapter 5. Conductivity Profile Through Electrical Potential Measurements 79 In deriving (5.4), the usual condition on the vanishing of both potential and current density at infinity is invoked, e.g. (j), Ji —> 0 at x, y —> ±co After substituting the expressions for horizontal components of current density from the last two equations into the first one, the system of equations can be written as follows: d<f)/dz = Jz/o dJz/dz = —I8(z — zs) — k2<t>o where k2 = k2 + k2 J — J fj-kxxa ^kyya With the source coordinates equal to the coordinate origin xs — ya = zs = 0, the system of equations reduces to d<j)/dz = Jz/cr (5.5) dJz/dz = -IS(z) - k24>o ^ 5.3 Propagator Solution For a layered system the solution is constructed using a propagator formulation. The propagator builds upon the solutions in a sequince of uniform layers. In order to prop-agate the solutions they have to be continuous through layer interfaces. The numerical problem is set up in a manner that the layer thicknesses {hi} = (hi,..., hm) are specified a priori and conductivities that correspond to the layers {oi} = (oi,...,om) are each constant. Chapter 5. Conductivity Profile Through Electrical Potential Measurements 80 For a particular wave number value k the system of equations (5.5) can be represented in a matrix form as dV(k,z) where dz V AV(k, z) + F, (5.6) is the solution vector, A = 0 -l/o -k2o 0 is a constant matrix in a uniform layer, and 0 -I8(z-z8) is the forcing vector representing the input current source. Since A is independent of z in each uniform layer, the solution can be written as V(z) = eA{Z-ZO)V0 + /"* eA(*~°F(0dC, JZQ where V0 — V(ZQ). When z < z3 the forcing term is zero and the solution is given by V(z) = eA{z-zo)V0 (5.7) Integrating across the sinqularity at zs introduces a discontinuity with a magnitude of F = [0, — I]T into the solution vector. Consequently, the solution in z > zs is given by V(z) = eA{z-Zo)V0 + F. (5.8) Chapter 5. Conductivity Prohle Through Electrical Potential Measurements 81 Since both potential <f> and vertical component of current density Jz are continuous across the interfaces that do not contain sources, the solution is continued from one layer to another by forming matrix product of the propagators for the individual layers. The propagator matrix for a layered medium is eA(*-z0) = p ( Z ) z._x)p(z._u z._2)... p/Zlt Z o ) ) (5.9) where P(ZJ, zA = eA(zi~z^ is evaluated in the jth layer. The propagator matrix P(z, z0) for a uniform layer can be given in terms of eigenvalues of the matrix A by N i-l where ^ = ( A - XJ) •••(A- Xj-il) • (A - A i + 1 I) •••(A — XNI) (Xi — Ai) • • • (Xi — Aj_i) • (Aj — A;+i) • • • (Xi — Aw) In the DC resistivity problem N = 2 and Ai = k , A 2 = —k, 1 1 -ak -(ak)'' 1 and V>2 1 (ak)'1 ak 1 Consequently, the propagator P(z, z0). = eA^z z°) is given by 1 -(ak)'1 ek{z-z0) + I 1 (ak)'1 e-k(z-z0) i -ak 1 2 ak 1 (5.10) which can also be written as Chapter 5. Conductivity Profile Through Electrical Potential Measurements 82 P{z,z0) = cosh k(z — ZQ) -sinh k(z — Zo)(ak) 1 ak sinh k(z — ZQ) cosh k(z — ZQ) The propagator solution must satisfy conditions on both upper (z = ZN) and lower (z = ZQ) boundaries as well as the jump condition (5.8) at z = zs. Denoting the boundary values at z = z0 and zN by V0 and VN, the jump condition can be written as P~(z„ z0)VQ = P+{zs, zN)VN + F, (5.11) where propagators P and P+ are matrix products of propagators for the individual layers below and above the source, e.g., P~ = P(ZS, ZQ) = P(zs, 2 m - l ) P ( ^ m - l , Zm-2) • • • P(Z\, ZQ), P+ = P(za, zN) = P(zs, Zi„i)P(zi-.i, «,_2).... P(zN-i, zN). The resulting algebraic system of equations for boundary values of potential and vertical component of current density is P11 P\2 P21 P22 4>o J. 2 0 4>. N J: + 0 (5.12) Pti Pu At P22 Hence, if two boundary values are known, the other two can be found by solving equatons (5.12) for the remaining two unknowns. Substituting both boundary values and propag-ator matrices into (5.8) and (5.9), we have the distribution of the potential and vertical component of current density in the Fourier Transform Domain. 5.4 Inverse Hankel Transform The inverse Fourier transforms for functions <j)(kx, ky, z) and Jz(kx, ky, z) are given by J roo roo cj>(x,y,z) = - f f $(kx,ky,z)e-ik*xe-ik>y 47T J-00 J-00 Chapter 5. Conductivity Profile Through Electrical Potential Measurements 83 1 roo i-oo _ Jz(x,y,z) = — / Jz(kx,ky,z)e lkxXe lk»ydkxdky_ 47T 7 - c o 7 - o o Through the appropriate change of variables the two inverse Fourier transforms can be expressed in terms of a single Hankel transform. The change of variables is x — r cos 9 y = r sin 0 kx = p cos <p ky = p sin (p So that 1 roo rl-K _ J>(x, y, z) = #r, 9,z) = -- pdp I <f>(p, zy^-^dip 1 roo r2-n _ ./,(*, y, z) = J,(r, 9,z) = — pdp I Jz(p, z)j"*"V-*)dv 47T 70 JO where 0 = 4>(k,z) = <j>(p,z) Jz = Jz(k,z) = Jz(p,z) are independent of <p. The integrals for both (j)(r, z) and Jz(r, z) reduced to 1 r°° l <f>(r,z) = —\ kJ0(kr)(j)(k,z)dk (5.13) 2 7 T 7 0 1 /-oo Jz(r,z) = — kJ0(kr)Jz(k,z)dk (5.14) 27T 7 0 where 1 r2* J0(kr) = — / * e-^'Q-rtdv Z7T 70 is the Bessel function of zero-order. The corresponding integrals in (5.13) and (5.14) define the usual inverse Hankel transform. Chapter 6 Discussion There are two main accomplishments in this thesis. The first is a theoretical calculation of the equilibrium states for a multicomponent thermodynamic system, based on the minimization of Gibbs free energy. The second is the experimental support for these theoretical predictions. 6.1 Theoretical Model The Gibbs free energy of a chemical system is evaluated on the basis of a known equation of state. For the systems considered in this study at modest pressure and temperature, two equations of state were used. One was due to Peng and Robinson [1970], and the other was given by Trebble and Bishnoi [1988]. Both equations of state were adequate for carbon dioxide - water mixtures under pressure and temperature of interest, but the methane-water mixture was only studied with the Trebble-Bishnoi equation. The Peng-Robinson equation did not work well with the methane-water mixture. One possible explanation is that the aqueous solubility of methane is 25 times less than that of carbon dioxide. In addition, the size of a methane molecule differs significantly from the size of water. Therefore, a hypothetical substance which is prescribed by an EOS for the purpose of modelling a mixture can have important properties, such as critical compressibility and covolume, that are seriously in error when the components of the mixture have unequal size or proportions. Differences are greater for methane-water mixtures than for carbon dioxide-water mixtures. Methane is usually considered to be hydrophobic substance, but 84 Chapter 6. Discussion 85 it has been shown [Haymet, 1994] that it is the water that "hates" the methane. The Peng-Robinson equation of state (PR EOS) is easy to use because it has only one empirical constant that determines the binary interaction of components, whereas the TB EOS has four constants. In order to find these coefficients, vapour-liquid equi-librium (VLE) experimental data are optimized. The values of the optimized coefficients depend upon the ranges of temperature and pressure. For CO2-H2O mixture, the inter-action coefficient was adjusted to fit the gas solubility data at pressures 1-30 atm and temperatures 5-10°C [Stewart & Munjal, 1970]. For methane - water mixture, Trebble and Bishnoi [1988] give all the four coefficients adjusted to the solubility data at higher both pressure and temperature. For the pressure and temperature of interest in this study some adjustments to these coefficients were required. Values for three of the co-efficients were adopted from [Trebble & Bishnoi, 1988], while the fourth was adjusted to satisfy the experimental solubility at lower P-T values. The resulting temperatures of three phase equilibrium, a vapour-liquid-hydrate equilibrium (VLHE), compare favour-ably with experimental data [e.g. Sloan, 1990]. For example, experimental values of the VLHE temperature for a C 0 2 - H 2 0 mixture, at 2 MPa, and for a CH4-H2O mixture, at 20 MPa, are 278°C and 291°C, respectively. The predicted values in this study are 277°C and 289°C. The differences reflect some inaccuracy in the thermodynamic description of the system, including the mixture rule and Langmuir constants. The error in tempera-ture for the C 0 2 - H 2 0 mixture is fairly constant, while the error for C H 4 - H 2 0 changes with P and T. 6.2 Experiment The second important component of this work is the experiment, which was designed to test the hydrate formation model. In the experiment, hydrate was crystallized from an Chapter 6. Discussion 86 aqueous solution of carbon dioxide in porous medium of sand at high pressure and low temperature. It was decided to test the phase diagram at a modest pressure of 300 psi; since 1 atm is equal to 14.7 psi, then this pressure is approximately 20 atm, or 2 MPa, as 1 atm=1.0M05 Pa. There are typically two degrees of freedom for a 2-component system at 2 phase equilibrium. This means that at fixed pressure, the temperature of phase equilibria is function of the gas concentration (see Chapter 2). Therefore, the temperature of hydrate formation depends on the amount of the dissolved gas in the solution. The more gas we manage to dissolve in solution, the less cooling is required to form hydrate (e.g. Figure 2.3). In a porous medium, the supply of gas is due solely to the process of diffusion (see Appendix B). The timescale for chemical diffusion is d2 where d is the distance to diffuse, and D is the chemical diffusivity. Therefore, the time to reach a desired concentration of gas is determined by the choice of d. In the experiment, hydrate growth is indicated with an abrupt change in measure-ments of the electrical field. To monitor this change, four electrodes were placed in the sediments, with two of them generating the electrical field, and the other two measur-ing the potential difference. This means that hydrate formation around the potential electrodes is most influential on the measurements. The gas concentration in the porous medium around the electrodes is affected by the distance between the electrodes and the top of the sediment column dtop, where a constant gas concentration is maintained. The temperature Te at the electrodes is controlled by the temperature at the cooled bottom dbottom- Consequently, the position of the electrodes in sediments relative to the top and bottom has an influence on the success of the experiment. Two sets of experiments were conducted with vertical and horizontal positions of Chapter 6. Discussion 87 the electrodes, respectively. The first set was not successful in monitoring the process of hydrate formation. The electrodes were located vertically in the sediments with dtop=0.05 m. At fixed P and T, three weeks were needed for the diffusing gas concentration to reach the position of the top electrode. The other electrodes were separated by 0.025 m. More than 6 additional days were required for gas to diffuse through this distance. Further difficulties were encountered in maintaining sufficiently low temperatures without freezing the water. Because the vertical array was well above the bottom, the lowest temperature we managed to achieve between electrodes 2 and 3 was only 3.56°C. The gas concentrations were typically too low at these temperatures to cross the phase boundary and produce hydrate. For the horizontal position of the electrodes, with dtop = 0.02 m, only 5 days were required for gas to diffuse to this depth. Since db0ttom = 0.008 m, temperatures as low as 2°C at the electrode depth were achieved. At this low temperature, a pronounced increase in resistance was detected in a several experiments, indicating the formation of hydrate. One example is shown in Figure 4.10. On the basis of the successful experimental data, theoretical developments (Chapter 3) and calculated phase diagram (Chapter 2), the volume fraction of hydrate can be calculated. In this sense, it is expedient to consider the following correspondences: • resistance —> conductivity; • conductivity —> gas concentration; • change of gas concentration —• fraction of hydrate formed. In order to quantify these relationships between the physical properties it was necessary to establish several constants from the experimental data. These results are summarized below. Chapter 6. Discussion 88 6.3 Geometric Factor The measured resistance R is connected with the conductivity cr by (see Chapter 4): * = (6-1) where C is a geometric factor. Therefore, to answer the first question in the chain it is necessary to find this factor. The forward problem of one dimensional DC resistivity for vertical location of the electrodes is considered in Chapter 5. Using the technique presented there, and resistance measurements collected in the experiments, the geometric factor is found to be 3.1298. Therefore, the conductivity dependence upon temperature T and gas concentration X can be presented as 3.1298 a{T>x) = imxj (6-2) Forward modelling of the resistance measurements for the horizontal location of the electrodes in our experimental setup was more difficult. Therefore, a reference value of resistance was used to interpret the measured changes in resistance in terms of changes in temperature and gas concentration. The details of this approach are discussed in sections 6.4 and 6.5. 6.4 Thermal Dependence of Conductivity To extract the dependence of conductivity upon the gas concentration, the temperature dependence needs to be accounted for. To find o = o(T) from the experiment, the concentration of dissolved gas should be constant. The difficulty is that as soon as T is lowered, more gas tends to be dissolved (see phase diagram). However, the sediments help to postpone the input of gas because of the slow chemical diffusion. Therefore, the sediments approximate a chemically closed system and allow the retrieval of the Chapter 6. Discussion 89 temperature dependence. Figure 6.1 presents the data of interest collected using the vertical array of electrodes. Over a small range of temperature values, conductivity can be represented as a linear function of T: (7 = ao ( l+&(T-r 0 ) ) , (6.3) where o0 is a conductivity value at T=T 0 , and (3g is a constant. According to this Taylor's series expansion, (3g can be written as ^ 9 cr0 ( d T ) T = T o ^ ^ The thermal coefficient (3g calculated using the data presented in Figure 6.1 is 0g = 0.03395. (6.5) This means that for a temperature varying from 5° to 14.5°C, conductivity of aqueous solution of gas as a function of temperature can be expressed as o{T) = <70(1 + 0.03395(T - 14)), (6.6) where o0 is conductivity at 14°C. This result does not depend on electrode configuration. Consequently, the same value of (3g was used to interpret resistance measurements using the horizontal array of electrodes. However, it is also possible to infer plg directly from the data using a reference value of resistance. Combining equations (6.1) and (6.3), we can write * W = ( 1 + ffi-r.))- ( 6 ' 7 ) where R(T 0) is the reference resistance which is measured at chemical equilibrium, when the gas concentration is known to be equal to the solubility at To. Note that the choice of R(T 0) is important due to its scaling function in (6.7). Using equation (6.7), the Chapter 6. Discussion 90 difference between resistance values at T i and T 2 is given by R.-R, = R(Ta)(1 + ^ T i _ T o ) - 1 + ^ n _ n ) ) . (6.8) Using equation (6.8) and the experimental measurements it is possible to find coefficient Pg. As shown in section 6.6, the value of (3g obtained with the chosen reference R(To) is in a very good agreement with the value determined using the vertical array. 6.5 Gas Concentration A phenomenological expression for conductivity of a solution with dissolved gas was given in Chapter 3 as aw(X,T) = A-y/X-f(T), (6.9) where A is a constant, f(T) describes the temperature dependence according to (6.3) and X is the gas concentration. The bulk conductivity a of a porous medium is calculated according to the Archie's law: o = (fafzaw{X)T\ (6.10) where </)Q is the porosity. Therefore, both resistance measurements and conductivity are functions of temperature, gas concentration and porosity. Let's suppose that a reference value of R is known at fixed T 0 and X 0 . In general, the equilibrium concentration Xo is solely determined by T 0 at fixed P 0 . Therefore, we view the reference resistance R(T 0) as a function of T 0 only. This known value of R(To) can be used to infer changes in gas concentration using the horizontal electrode locations. Since the geometric factor (6.1) can be written as C = R(T0)o(T0,X0), (6.11) Chapter 6. Discussion Figure 6.1: Conductivity as a function of T at constant gas concentration. Chapter 6. Discussion 92 the conductivities at other T and X are related to the resistance measurements by *(X,T,<t>0)= R { X i T M • (6-12) Substituting conductivity a and the reference value <TO into (6.10) and using then de-pendence on the solubilities X and X 0 from (6.9), we obtain after some manipulations The porosity is annihilated because the resistance measurements are taken in the same sediments. Therefore, the resistance measurements R(T) can be related to the gas con-centration through (6.13). If we define R(T,X 0 ) t / i e o ' ' at a fixed reference concentration using (6.7), then expression (6.13) can be written as where R(T ,X) e x p are the experimental measurements. 6.6 Application Figure 4.10 shows resistance as a function of temperature in the experiment with the horizontal array of electrodes. Before these data were collected, the apparatus was fed with the gas and maintained at a pressure of 100-150 psi for 2 weeks: The experiment started in the evening and continued over night (as the room temperature is the lowest at night), with coolant at temperature of -2°C. During the cooling, the T and X conditions crossed into the hydrate stability field and hydrate was produced. In the morning, the cooler was turned off and the phase boundary was crossed once again, so that the hydrate dissociated. The resistance change during the process was monitored. The slope of resistance over temperature in Figure 4.10 is fairly constant. This suggests that gas was not diffusing into the sediments during the cooling, and that the resistance change was Chapter 6. Discussion 93 due solely to the influence of temperature. As soon as the gas concentration changed due to hydrate formation, the dependence of resistance upon temperature had a very different slope. In fact, slopes are almost vertical. For the purpose of this discussion we will consider both the increase and decrease as vertical. In this case, these changes are due solely to changes in gas concentration, and they can be related to volume of hydrate formed and dissociated, respectively. The resistance increase due to hydrate formation and its decrease due to hydrate dissociation are almost the same. This means that the amounts of gas consumed and released during these rapid changes are equal. Therefore, the process is symmetric about the phase equilibrium. In Chapter 4, we discussed the dependence of resistance upon temperature under different assumptions about the corresponding changes in gas concentration (Figures 4.2 and 4.3). The curve of Figure 4.10 is similar to that of a chemically closed system with nonequilibrium hydrate formation (ABC path in Figure 4.2, 4.3). This implies that hydrate formed in the experiment from overcooled solution, where the amount of overcooling affected the volume of hydrate that was formed. Since both abrupt changes in resistance measurements are due to hydrate growth and dissociation, it follows that the overcooling, as well as overheating, is measured by the difference between the temperatures where these rapid changes occur. Since we know the equilibrium temperature of hydrate formation and how far the actual temperature of formation differed from equilibrium, we are set to quantatively interpret the experimental data on phase diagram (P=2 MPa). A highly undersaturated solution is cooled from 6°C. Hydrate begins to form after the solution is overcooled to point B in Figure 6.2. Hydrate formation consumes gas and lowers X to point A. Further cooling of a few tenths of a degree, after the abrupt increase in Figure 4.10, produces a small additional increase in resistance. Much of this increase can be attributed to the temperature dependence of R without any change in gas concentration. Since we expect a decrease in gas concentration Chapter 6. Discussion 94 with a lowering of temperature, when the solution remains in equilibrium, we interpret this small additional cooling as overcooling without much additional hydrate formation (point C in Figure 6.2). As temperature increases, the overcooling is eliminated and eventually the solution becomes overheated. Finally, hydrate dissociates at point F, releasing gas and returning the solution to its initial concentration. The gas concentration change due to hydrate formation and dissociation are roughly equal at AX «' 0.8 x 10 - 3 mole fraction. The reference concentration is determined from the equilibrium state that was es-tablished prior to the start of the experiment. As the apparatus was at a pressure of 100 psi and room temperature, which was roughly 14°C, for 2 weeks before the exper-iment started, we assume that at this P-T condition the system reached equilibrium, with a corresponding resistance measurement of 3400 Ohm. Hence, the reference value of resistance is R(T0) = 3400, where T0=14. Knowledge of the reference value for resistance and corresponding equi-librium concentration X 0 at T 0 and P=100 psi permits an estimate of Pg to be obtained from the experimental data that were collected with the horizontal array of electrodes. Using these data in equation (6.8) gives pg = 0.0330. This value differs only by 3% from that calculated for the vertical array of electrodes, according to (6.5). Figure 6.3 represents the solubility X as a function of temperature T calculated on the basis of relation (6.13). The concentration of gas in the solution is roughly consistent with a location of phase change at 2.2°C (see Figure 6.2), and the concentration change is 1.4xl0 - 3 mole fraction. Chapter 6. Discussion 95 0.020 0.005 1.0 2.0 3.0 4.0 Temperature (C) 5.0 6.0 Figure 6.2: Carbon dioxide gas concentration as a function of temperature at pressure of MPa. er 6. Discussion Figure 6.3: Gas concentration as a function of T in hydrate forming process. Chapter 6. Discussion 97 6.7 Hydrate Fraction The decrease in dissolved gas concentration resulted from the consumption of gas by hydrate. A correspondence between the hydrate fraction and the change in gas concen-tration is given by Rempel and Buffett [1996] where A X is the change in dissolved gas concentration, and X/i is the gas concentration in hydrate structure. In an ideal hydrate structure I, there are 8 moles of gas per 46 moles of water, or Xh = 0.148. But in a real situation not all vacancies are filled. In fact, the thermodynamic model of hydrate in Chapter 2 provides an estimate of the occupancy. Typical calculations yield Xh = 0.125. The concentration change inferred from the electrical resistance measurements is AX = 1.4 x I O - 3 , so the corresponding hydrate fraction is h « 1.1%. It is also possible to estimate hydrate fraction from the slope of the X-T solubility curve (Figure 6.2). It was previously assumed that the solution was overcooled by 0.2°C before hydrate began to form. When the solution return to equilibrium, the change in X inferred from the solubility curve h AX (6.15) Xh AX = 0.8 x 10 - 3 implies that h « 0.6%. Chapter 6. Discussion 98 Property Nominal Value Units Pf 1000 kg m~3 Ph 1110 kg m - 3 P* 2650 kg m - 3 C/ 4200 J kg"1 °C a 2200 J kg" 1 °C 0.33 448000 J k g - 1 Table 6.1: Physical Properties, [Rempel, 1994] There is a third way of estimating the amount of hydrate that formed, based on the temperature change due to the latent heat release. Figure 4.9 represents the temperature change during the whole experiment. The temperature records for the lower two RTDs are expanded in Figure 6.5 during the time that hydrate formed. The thermal signature of hydrate formation due to the latent heat release can be found on the basis of the equation of energy balance (Chapter 3, equation (3.1)). Since overcooling causes the hydrate to form very quickly, it is reasonable, as a first approxi-mation, to neglect the effect of thermal diffusion and advection. In this way, the latent heat release can be directly related to the temperature change, and equation (3.1) can be written as Ah = C^-^-AT. (6.16) <t>L Ph On the basis of values given in Table 6.1, the normalized bulk heat capacity (7=1.23, Chapter 6. Discussion 99 Figure 6.4: Temperature change during hydrate formation. Dashed lines represents thermal signatures of hydrate formation, a latent heat effect. Chapter 6. Discussion 100 and the coefficient of proportionality between hydrate volume Ah and A T is 0.0364. The thermal signature of hydrate formation recorded on the bottom RTD is 0.14°C, while the middle RTD measured 0.07°C. Therefore, the hydrate fraction formed near the bottom is hb « 0.5%, while the hydrate fraction near the middle, where the potential electrodes were placed, is hm « 0.25%. These values are somewhat lower than the values inferred from the electrical measure-ments. Heat losses due to diffusion could lower the temperature changes for prescribed volumes of hydrate formation. 6.8 Hydrate Distribution in Marine Environments The theoretical calculations for the phase diagram of gas hydrate can be used to determine the spatial distribution of hydrate in marine sediments. The determination is based on the idea that any gas in excess of the equilibrium concentration will be incorporated into hydrate. Conversely, when the actual concentration in the pore fluid falls below the equilibrium concentration, hydrate will dissociate and release gas to re-establish the equilibrium concentration. Therefore, the gas concentration in a hydrate zone should remain at the equilibrium concentration (e.g. gas solubility in presence of hydrate). The dashed line in Figure 2.4 represents a typical estimate of solubility in the hydrate stability zone as a function of depth below the seafloor. Now consider what happens when a parcel of fluid carries methane gas from depth into the hydrate zone. When the gas concentration exceeds the equilibrium concentration this excess gas will be removed and incorporated into hydrate. At the bottom of the hydrate Chapter 6. Discussion 101 layer rate of decrease of ceq is greatest, so the rate of hydrate formation should also be greatest close to BSR. Mathematically, the rate of hydrate formation can be determined from the conserva-tion of gas where the gas is transported mainly by fluid motions [Rempel and Buffett, 1996]. In this fluid, the concentration of gas is equal to the gas solubility X e g . Since Xeq does not depend on time, this implies that the first term in equation (3.2) (Chapter 3) vanishes. Neglecting changes in the fluid gas concentration due to diffusion, equation (3.2) becomes dh = Pfu • V X e q dt ph<t>{Xh - Xeq)' where p/ and ph are densities of fluid and hydrate, u is the transport velocity, (j) is porosity, Xh is concentration of gas in hydrate structure, and Xeq is the equilibrium concentration of gas in solution for hydrate stability zone. The concentration of gas in hydrate Xh is typically 25 times higher than Xeq under marine conditions, so the difference Xh—Xeq can be considered as a constant, even though Xeq varies across the hydrate layer. Therefore, all the parameters in (6.17) are independent of the depth, except for u • V X e q . Since u is nearly constant, the rate of accumulation over the depth is dh dXeq « a - ( 6 1 8 ) where dXeq/dz is gradient of the dashed line in Figure 2.4. This plot indicates that Xeq decreases almost exponential with distance above the BSR, so oc Xeq (6.19) and consequently h(t) « t— oc tXeq which implies that the spatial distribution of hydrate above the BSR is very similar in form to X e g in Figure 4.2. These theoretical predictions are in a good agreement with both Chapter 6. Discussion 102 geophysical and geochemical estimates of hydrate volumes [Rowe and Gettrust, 1993; Singh and Minshull, 1994; Yuan, at al., 1996], which suggests that hydrate volume decreases above the BSR. The total amount of hydrate depends on the accumulation time for a prescribed velocity. Assuming a representative vertical velocity of 4-10 - 1 1 m/s [Davis et al., 1990], we obtain a rate of hydrate formation of roughly 1% of the pore volume in 104 years. Chapter 7 Conclusions The primary goal of this study was to develop a quantitative model of hydrate formation within a porous medium and to test this model experimentally. The first step toward this goal involved extending the modelling of Rempel [1994] to include nonequilibrium effects and the influence of salt. The second step was to develop a system for using electrical field measurements to monitor hydrate growth in an experimental test of the theory. Salt was proposed to facilitate the electrical field measurements. This motivated the addition of an equation for salt conservation to the system of governing equations. A numerical program was written to solve the governing equations, and the results showed that the addition of salt did not aid the detection of hydrate. In fact, the most sensitive characteristics was the change in the concentration of hydrate-forming gas. The change in gas concentration is almost an order of magnitude more important than change in porosity, or salt, or temperature. Therefore, salt was omitted from the experiments. Measurements of electrical potential in an aqueous solution in porous medium are a function of conductivity. Conductivity of aqueous solution in porous media is described by Archie's law. It depends upon temperature, porosity and ion concentrations (dissolved gas yields two different ions). Therefore, on the basis of governing equations for the hydrate formation process, the conductivity change was determined. This value was compared with values inferred from the experiments to provide quantative estimates of hydrate formation. 103 Chapter 7. Conclusions 104 But in order to interpret the experimental data, the equilibrium concentration of gas is needed in theoretical predictions. Even more, the dependence of solubility upon the temperature is important. Therefore, the solubility of two hydrate-forming gases, CO2 and C H 4 , was calculated over a range of pressure and temperature. The calculations are based on the minimization of the Gibbs free energy for a two-component (gas and wa-ter), three-phase (gas, liquid, hydrate) system. The minimization was performed using a simulated annealing algorithm which worked accurately and efficiently for this kind of problem. It was found that the solubility of gas is significantly altered by the pres-ence or absence of the hydrate phase. When hydrate is absent at high temperature the calculations reproduce the experimentally observed increase in solubility as temperature decreases. When hydrate is present, the gas solubility drops sharply with decreasing tem-perature. Such an abrupt decrease in solubility permits hydrate to crystallize directly from an aqueous solution, without the need of any free gas. In order to test this important result an experiment was designed. An aqueous solu-tion of carbon dioxide in porous medium was cooled at a pressure of 2 MPa. According to the P-T conditions of the experiment, the solution is undersaturated until hydrate forms. Hydrate growth in the system was monitored with electrical field measurements. This choice is based on the fact that hydrate forms through depletion of gas from the solution. Therefore, resistance measurements, being sensitive to the ion concentration in solution, were assumed to use in order to monitor hydrate formation. The experiment setup approximated a compositionally closed system and showed an abrupt increase in the monitored resistance at a temperature of 2.2°C. The increase is approximately 300 Ohm. This increase at 2.2°C was monitored for three sets of the experiment which differed slightly in the starting temperature. At fixed pressure, tem-perature of hydrate formation depends solely on gas concentration (Gibbs rule), so the constant formation temperature implies that amount of gas in the solution under study Chapter 7. Conclusions 105 was also constant. Moreover, the pronounced resistance change was accompanied by a change in temperature due to latent heat release. The produced hydrate volume fraction was calculated on the basis of both the resistance and temperature change. The estim-ated volume fractions are roughly equal. Therefore, the possibility of hydrate formation from an aqueous solution was proved. This is an important result because of its implications for the formation of gas hydrate in marine environments, where the gas supply may not be sufficient to provide free gas. A resolution of this question is important to understand how hydrate forms in nature. Based on these results, a simple model of hydrate formation was posed to estimate the vertical distribution of hydrate and the rate of its accumulation in marine sediments. References Archie, G.E., 1942, The electrical resistivity log as an aid in determining some reservoir characteristics: A.I.M.E., 146, 54-67. Brooks, J.M., Kennicutt, M . C , and et al., 1984, Thermogenic gas hydrates in the Gulf of Mexico: Science, 225, 409-411. Carmichael, R.S., Ed., 1989, Practical handbook of physical properties of rocks and minerals: CRC Press. Christensen, R.L., and Sloan, E.D., Jr., 1994, Mechanism and kinetics of hydrate form-ation: Natural Gas Hydrates, editted by Sloan, E.D., Happel, Jr.J., and Hnatow, M.A., New York Academy of Sciences, New York, 715, 283-305. Co, C.C., 1994, Multicornponent vapor-liquid-liquid flash calculations through Gibbs free energy minimization using the simulated annealing algorithm: B.Sc. Thesis, University of British Columbia, Vancouver, Canada. Collett, T.S., 1993, Natural gas hydrates, of the Prudnoe Bay and Kuparuk River area, North Slope, Alaska: Am. Assoc. Pet. Geol. Bull., 77, 793-812. Davis, E.E., Hyndman, R.D., and Villinger, H., 1990, Rates of fluid expulsion across the Nortern Cascadia accretionary prism: Constraints from new heat flow and multi-channel seismic reflection data: J. Geophys. Res., 95, 8869-8889. Denbigh, K., 1981, The principles of chemical equilibrium: Cambridge Univ. Press, New York. Ellis, D.V., 1987, Well logging for earth scientists: Elsevier Science Publ. Co., Inc. Englezos P., and Hall S., 1994, Phase equilibrium data on carbon dioxide hydrate in the presence of electrolytes, water soluble polymers and montmorillonite: Can. J.Chem. 106 REFERENCES 107 Eng., 72, 887-893. Hammerschmidt, E.G., 1934, Formation of gas hydrates in natural gas transmission lines: Ind. Eng. Chem., 26, 851-857. Handa, Y.P., 1990, Effect of hydrostatic pressure and salinity on the stability of gas hydrates, J. Phys. Chem., 94, 2652-2657. Haymet, A.D.J., 1994, Water and aqueous solutions near freezing: integral equation theory and computer simulations: Natural Gas Hydrates, editted by Sloan, E.D., Happel, Jr.J., and Hnatow, M.A., New York Academy of Sciences, New York, 715, 146-158. Heiland, C.A., 1940, Exploration geophysics: Prentice-Hall, Inc. Holder, G.D., Corbin, G., Papadopoulos, K.D., 1980, Thermodynamic and molecular properties of gas hydrates from mixtures containing methane, argon, and kripton: Ind. Eng. Chem. Fundam., 19, 282-286. Holbrook, W.S., Hoskins, H., Wood, W. T., et al., 1996, Methane Hydrate and free gas on the Blake Ridge from vertical seismic profiling: Science, 273, 1840-1843. Hyndman, R.D., and Davis, E.E., 1992, A mechanism for the formation of methane hydrate and seafloor bottom-simulated reflectors by vertical fluid expulsion, J. Geo-phys. Res., 97, 7025-7041. Jeffrey, G.A., and McMullan, R.K., 1967, The clathrate hydrates: Prog. Inorg. Chem., 8, 43-108. Ingber, A.L., 1989, Very fast simulated re-annealing, J. Math. Comp. Modelling, 12, 967-973. Kirkpatrick, S., Gelatt, C D . , Jr., Vecchi, M.P., 1983, Optimization by simulated anneal-ing: Science, 220, 671-680. Makogon, Yu.F., 1981, Hydrates of natural gas: translated from Russian by W.J.Cieslewicz, PennWell, Tulsa, Okla. REFERENCES 108 Nordstrom, D.K., and Munoz, J.L., 1986, Geochemical Thermodynamics: Blackwell Sc. Publications. Parrish, W.R., and Prausnitz, J.M., 1972, Dissociation pressure of gas hydrates formed by gas mixtures: Ind. Eng. Chem. Process Des. Dev., 11, 26-35. Paull, C D . , Ussier III, W., Borowski, W.S., 1994, Source of methane to form marine gas hydrates: Natural Gas Hydrates, editted by Sloan, E.D., Happel, Jr.J., and Hnatow, M.A., New York Academy of Sciences, New York, 715, 392-409. Peng, D-Y., and Robinson, D.B., 1976, A new two-constant equation of state: Ind. Eng. • Chem., Fundam., 15, 59-64. Phillips, O.M., 1991, Flow and reactions in permeable rocks: Cambridge Univ. Press, New York. Rempel, A., 1994, Theoretical models of hydrate formation: M.Sc. Thesis, University of British Columbia, Vancouver, Canada. Rempel, A., and'Buffett, B.A., 1996, Formation and accumulation of gas hydrate in porous media: submitted to J. Geophys. res., August. Ripmeester, J.A., Ratcliffe, C.L, Klug, D.D., and Tse, J.S., 1994, Molecular perspectives on structure and dynamics in clathrate hydrates: Natural Gas Hydrates, editted by Sloan, E.D., Happel, Jr.J., and Hnatow, M.A., New York Academy of Sciences, New York, 715, 161-176. Routh, P.S., and Roy, K.K. , 1996, Crustal resistivity inversion using global optimization techniques: Deep Electromagnetic Exploration, Narosa Publications, New Delhi (to be published). Rowe, M.M., and Gettrust, J.F., 1993, Fine structure of methane hydrate-bearing sedi-ments on the blake Outer Ridge as determined from deep-tow multichannel seismic data: J. Geophys. Res., 98, 463-473. Singh, S.C., and Minshull, 1994, Velocity structure of a gas hydrate reflector at ocean REFERENCES 109 drilling program site 889 from a global seismic waveform inversion: J. Geophys. Res., 98, 221-233. Sloan, E.D., 1990, Clathrate hydrates of natural gases: Marcel Dekker, New York. Stern, L.A., Stephen H.K., and Durham, W. B., 1996, Peculiarities of methane clathrate hydrate formation and solid-state deformation, including possible superheating of water ice; Science, 273, 1843-1848. Stewart, P.B., and Munjal, P., 1970, Solubility of carbon dioxide in pure water, and synthetic sea water concentrates at -5° to 25°C and 10- to 45-atm pressure: J. Chem. Eng. Data, 15, 67-71. Stryjek, R., and Vera, J.H., 1986a, PRSV: An improved Peng-Robinson equation of state for pure compounds and mixtures, 64, 323-330. Stryjek, R., and Vera, J.H., 1986b, PRSV2: A cubic equation of state for accurate vapor-liquid equilibria calculations, 64, 820-826. Toda, M. , Kubo, R., and Saito, N., 1983, Statistical physics: Springer Verlag, Berlin. Trebble, M.A., and Bishnoi, P.R., 1986, Accuracy and consistency comparisons of ten cubic equations of state for polar and non-polar compounds: Fluid Phase Equilibria, 29, 465-474. Trebble, M.A., and Bishnoi, P.R., 1987, Development of a new four-parameter cubic equation of state: Fluid Phase Equilibria, 35, 1-18. Trebble, M.A., and Bishnoi, P.R., 1988, Extension of the Trebble-Bishnoi equation of state to fluid mixtures: Fluid Phase Equilibria, 40, 1-21. Uchida, T., 1996, Studies on formation and dissociation rates of methane hydrates in pure water - pure methane gas systems: Gas Hydrate Studies, Canada-Japan joint science and technology workshop. Van der Waals, J.H., and Platteeuw, J.C., Clathrate solutions: Adv. Chem. Phys., 2, 1-57. REFERENCES 110 Van Laarhoven, P.J.M., and Aarts, E.H.L., 1987, Simulated annealing: Theory and Applications: D.Reidel Publishing Company. Vidal, Rene V.V., Ed., 1993, Applied Simulated Annealing: Springer Verlag, Berlin. Worthington, A.E., Hedges, J.H., and Pallatt, N. , 1990, SCA guideline for sample prepa-ration and porosity measurement of electrical resistivity samplers: The Log Analyst, 31, 20-28. Yamane, K., and Aya, I., 1995, Solubility of carbon dioxide in hydrate region at 30MPa, Proceedings of the MARIENV'95 Conference, 911-917. Yuan, T., Hyndman, R.D., Spence, G.D., and Desmons, B., 1996, Seismic velocity in-crease and deep-sea hydrate concentration above a bottom-simulated reflectoron the northern Cascadia continental slope: J. Geophys. Res., 101, 655-13671. Appendix A Peng - Robinson and Trebble - Bishnoi Equations of State Equations of state describe a quantative relationship that exists between intensive vari-ables, like pressure and temperature and extensive variables, like volume. Such equations of state are the building blocks of thermodynamics, and they represent important aspects of the behavior of physicochemical systems. As with all models, however, they are not universally applicable. We consider two equations of state in this study. One due to Peng and Robinson [1976], while the other is given by Trebble and Bishnoi [1987]. Both are adjusted for real gases on the basis of the van der Waals equation of state, which quantifies the dual effects of repulsive and attractive intermolecular forces. Peng and Robinson [1976] proposed an equation of the form P _ R T • <n ( A 1 ) v-b v(v + b) + b(v - by K ' ) which can be rewritten in terms of the compressibility factor Z as Z 3 - (1 - B)Z2 + {A - 3B2 - 2B)Z - {AB - B2 - B3) = 0, (A.2) where aP (A.3) R?T2' 111 Appendix A. Peng - Robinson and Trebble - Bishnoi Equations of State 112 Equation (A.2) yields one or three roots depending upon the number of possible phases. In the two-phase region, the largest root is for the compressibility factor of the vapour phase, while the smallest positive root corresponds to that of the liquid. The compressibility factor for vapour is always larger than its value at the critical point. The critical point, characterized by a critical pressure P c , a critical temperature T c and a critical volume V c , represents the last occurence of boiling in a pure system. A parameter value at the critical point is custommarily called a critical value. The critical values of an attraction parameter a, covolume b and compressibility Z of a given substance can be found on the base of equation (A.2), noting that at the critical point with the van der Waals' condition ((dP/dV)T = (d2P/dV2)T = 0) (Z - Zcf = 0. (A.6) Equation (A.6) can be written as Z3 - 3ZCZ2 + 3Z2Z - Z\ = 0. (A.7) Considering equality of coefficients at the same power of Z in equations (A.2) and (A.7) we have three equations for unknowns A c , B c and Z c : 3ZC = 1 - Bc, (A.8) 3Z\ = AC- 3B2C - 2BC, (A.9) Z\ - ACBC - B\ - B3C, (A.10) where Ac = ^ (A.11) R2T2' Bc = (A.12) Zc = f^. (A.13) Appendix A. Peng - Robinson and Trebble - Bishnoi Equations of State 113 The solution of (A.8)-(A.10) for A c , B c and Z c is Bc = 0.07780, Ac = 0.45724, Zc = 0.307. According to relations (A. 11) and (A. 12), the critical values of parameters a and b are b{Tc) = 0 .07780^, (A.14) Pc R2T2 a(Tc) = 0.45724—^-. (A.15) At temperatures other than the critical value, the parameters a(T) and b(T) are assumed to take the form b(T) = b(TE), (A.16) a(T) = a(Tc)-a(Tr,u;), (A.17) where a(Tr,u>) is a dimensionless function of reduced temperature Tr T Tc and acentric factor u>. The dependence of a on Tr is expressed by a 1 / 2 = 1 + K(1 - Tl' 2), (A.18) where n is a constant for each substance. These constants have been correlated against the acentric factors, with a best-fitting curve given by K = 0.37464 + 1.5422a; - 0.26992a;2. (A.19) Stryjek and Vera [1986b] proposed an expression for K that depends on T r , K = K O + \ K L + K2(K3 - r , ) ( i - rr72)](i + TV*)(O.7 - r r ) , (A.20) Appendix A. Peng - Robinson and Trebble - Bishnoi Equations of State 114 where K0 = 0.378893 + 1.4897153w - 0.17131848oj2 + 0.0196544w3 (A.21) and K i , K2 and K3 are adjustable parameters for pure compounds. Values of w, K i , and the critical constants for over ninety compounds of industrial interest are published by [Stryjek and Vera, 1986a]. The values of K2 and K3 for the same compounds, which represent the temperature dependence of K, are given by [Stryjek and Vera, 1986b]. The equation of state ( A l ) is denoted PRSV2 when K has an extended form (A20) — (A21), and it is called PRSV when K2 in (A20) equals zero. In order to derive the expression for the coefficient of fugacity 0 we use the equation of state in the thermodynamic relationship [Denbigh, 1982] f rp v 1 = / ("^ - ^) d P- (A.22) P Jo RT P' v ' After some manipulations, the logarithm of the coefficient of fugacity c6 of a pure com-ponent can be written in the form lncf> = Z - 1 - ln(Z -B) + ( A . 2 3 ) where A and B are calculated according to (A3) —(A4), with parameters a and b given by (A16) and (A17) and the critical values defined by (A6) - ( A l l ) . The compressibility factor Z is a root of equation {A.2). In the case of a mixture of different components, we define a mixing rule for parameters a and b of the components and suppose that these parameter represent a hypothetical pure substance. The fugacity coefficient of component k in the mixture with mole fraction x/- is calculated from the following expression xkP bm 2 \/2Bm Z - 0.414Bm am bm Appendix A. Peng - Robinson and Trebble - Bishnoi Equations of State 115 where Am and Bm are calculated from (A.3) — (AA), with parameters a and b being equal to the mixture parameters o m and bm. These parameters are defined by a typical mixing rule a m — x i x j a i j i (A.25) i 3 bm = YJXibi, (A.26) i an = (1 - 6iAa\l2a)12, (A.27) where 8ij is an empirically determined binary interaction coefficient characterizing the binary formed by component i and component j. Many authors (e.g. [Co, 1994]) have proposed various mixing rules of the van der Waals type with composition dependent binary interaction parameters. Although these rules have been used successfully for some nonideal mixtures, they are inconsistent with the statistical mechanics requirement that the second virial coefficient obtained from the equation of state at the lower pressure limit be a quadratic function of composition. Consequently, attempts have been made to develop density-dependent mixing rules which obey this theoretical requirement. However, these mixing rules do not preserve the cubic nature of any equation of state, and they are not considered here. The Trebble-Bishnoi equation of state (TB EOS) is a four-parameter equation of state. Trebble and Bishnoi [1987] showed that if the EOS is based on the van der Waals' conditions there are four parameters required to optimize both the value of the critical compressibility Zc and the value of critical covolume bc predicted by the EOS. A third parameter is important in reproducing the experimental critical volume, while a fourth one being useful in optimizing the "hardness" of an equation of state (the slope of (dP/8V)T) at elevated pressures. Appendix A. Peng - Robinson and Trebble - Bishnoi Equations of State 116 The TB EOS can be written in a form P R T < T ) f A 2 8 ) v-b v* + (b + c)v-bc-eP' K ' } where the first term is the same as in the Peng-Robinson EOS (van der Waals repulsive pressure), and the second one contains two more parameters than that of Peng-Robinson. Both these additional parameters c and d are considered to be constant. For parameters a and b at temperatures other than the critical, the suggested dependences are b(T) = bc{3 (A.29) 13 = 1 + q2(l -Tr + lnTT) T <TC (A.30) p = l T>TC (A.31) a(T) = aca (A.32) a = exp{qx{\ - TT)) (A.33) where T r is the reduced temperature and qi, q2 are constants. Going through the usual procedure at the critical point (equations (A.6) — ( A 10)), we obtain three equations with five unknowns: Cc - 1 = -3ZC, (A.34) Ac - 2BCCC -Be-Ce- B\ - D2C = 3Z2C, (A.35) Bl + (2 - 3ZC)B2C + 3Z2BC - (D2C + Z\) = 0, (A.36) where Ac and Bc are determined by expressions (11) — (12), and c P c~ RTC' RTC In a parameter optimization procedure, Trebble and Bishnoi [1988] decided to allow Zc to be a variable along with the parameters d, q\ and q2 used in equations (A.33) and Appendix A. Peng - Robinson and Trebble - Bishnoi Equations of State 117 (A.30). They note that four optimal parameters are not always justified by the PVT data, but they succeeded in generalizing the procedure for these parameters [Trebble and Bishnoi, 1986; Trebble and Bishnoi, 1987]. The optimal values of Z c , d, qi, and q 2 for 75 pure components along with their critical properties and acentric factors are presented [Trebble and Bishnoi, 1987]. Using the TB EOS in the thermodynamic expression (A.22), the following expression for the coefficient of fugacity <f> of a pure component is derived: case 1 (r > 0): case 2 (r < 0): where ln\t\ = Z - l - ln\Z -B\ + - ^ - A i (A.37) ln\L = Z - l - ln\Z -B\ + -=^-A2 (A.38) P t>v2 u = 1 + ^ (A.39) T = l + T + v+-vr (A-40) 0i = ( r ) 0 5 (A.41) 02 - ( - r ) 0 5 (A.42) . .2Z + B(u-6) . ^ = lnl2Z + B(u + e^ ( A - 4 3 ) A2 = [ 2 i a n - 1 ( ^ ± ^ ) - 7 r ] (A.44) The coefficient of fugacity of component j in a mixture is calculated from the following expression case 1: , fj b d i\ , w r> i , ArnXx qd/n bd n9d . . In = — (Z - 1) - In2 - B m + — (A.45) X j F 6m JOm6'i a m 6m #i Appendix A. Peng - Robinson and Trebble - Bishnoi Equations of State 118 Am ZBmn9d + Q.bB2m(un9d - 9xnudQ Bmdx lZ2 + (Bm + Cm)Z - (BmCm + DI (A.46) n ° d = ^ T o - [Qb™cd - 6bdcm + 2cmcd + 8dmdd - r^-(2<4 + 8d2J) (A.47) t-V\t>m °m. case 2: In J L = £-(2 - 1) - ln\Z - Bm\ + - £ - (A.48) A m r ZBmn9d + 0.bB^(un9d - 92nud0 _^ ^ 4g^ where Bm92 l Z 2 + (23ro + C m ) Z - ( B m C m + D2m n9d = T^TJ-fibmCd - Qbdcm + 2cmcd + 8dmdd - ^-(2c 2 n + 8d2m)) (A.50) lv2om om u = 1 + ^ (A.51) 6r r2 Ad2 T = 1 + ^ + W + ^T (A-52) fcm b2m K > A a = t 2 f a n - 1 ( 2 Z p + ^ ) - T ] (A.54) Am - JfPp. (A.55) Bm = ^ (A.56) G» = ^ (A.57) A» = ^ (A.58) a m = £ £ ^ j ^ j (A.59) i j bm^Y^Yl XiXJbij (A.60) i j Cm = ^j^2XiXjCij (A.61) Y^Ylxixidv (A.62) Appendix A. Peng - Robinson and Trebble - Bishnoi Equations of State 119 on = (1 - S&faaj) 1'* (A.63) 6, = ( i -^ ) (^±M ( A 6 4 ) Ci3 = (i - ( A . 6 5 ) ad = 2^2 (A.67) bd = 2J2 Xibij - bm (A.68) i=n cd = 2 ]T XjCij - c m (A.69) i = l d d = 2 ^ 2 Xidij - dm (A.70) i=l cd bdc nud = 7 r r — (A.71) Appendix B Mechanism of Gas Transference in F lu id B . l Equation of State The density of an aqueous solution p depends on the temperature T, pressure P and the concentration of dissolved substances c». The relationship among these quantities p = p(T,P,.ci) (B.l) is called the equation of state. For the sake of simplicity, let's consider only one dissolved substance with concentration c. Variations in fluid density may result in fluid motion due to the influence on the fluid buoyancy. If the variations in temperature and concentration are independent and small, then the solution density can be written as a linear function of these variations: p = p 0 ( l - ccT + Pc), (B.2) where po = Po(-P) is the reference density at the ambient pressure and i dp i dp " = P = — -5- , Po oT Po dc are the coefficients of thermal and chemical expansion, respectively. Concentration is typically given as a mass fraction. When a thermodynamic system is at equilibrium, or, equivalently, a solution is sat-urated, the equation of state has a form p = p 0 ( l - a,T), (B.3) 120 Appendix B. Mechanism of Gas Transference in Fluid 121 which depends only on temperature bevause c is a known function of T. The coefficient aa is defined by 1 d 1 dp dp dceq Po oT Po ol dceq dl and it describes the thermal expansion coefficient of a saturated solution. It can also be represented as a. = a - / ? | p . (B.4) B.2 Compositional Convection Vigorous compositional convection was observed in our experiment as carbon dioxide was dissolved into the water. Variations in gas concentration result in changes in the solution density which affect the fluid buoyancy and cause the fluid to move. The effects of temperature variations were small, so the density varies only with composition, according to p = p 0(l+/?c), (B.5) where c is gas concentration, and /? is the coefficient of chemical expansion. Experimental data on the density of water with dissolved C 0 2 [Yamane and Aya, 1995] were used to estimate f3 for C 0 2 - H 2 0 mixture. Density variations Ap over a vertical distance d provide the driving force for convective motion. The tendency of the fluid to convect is measured by the Rayleigh number where D is chemical diffusivity for compositional variation in a fluid with viscosity p. With Ap from (B.5) the Rayleigh number can be written as fipocgd3 , . J ? a = _ ^ — , (B.7) Appendix B. Mechanism of Gas Transference in Fluid 122 Concentration of chemical expansion, ft 0.24 Concentration of C 0 2 , c 0.02 Distance for density varition, d 0.7 Chemical diffusivity, D 10~9 Viscosity, u 10 - 3 Reference density, p0 10 - 3 Gravitational constant, g 10 Table B . l : Rayleigh number parameters (SI units) Convection begins when Ra > Racr, where RaCT is the critical Rayleigh number which depends on the boundary conditions as well as geometry of the cell where this motion arises. The cell's geometry is characterized by the ratio of the horizontal width to the vertical thickness, called the aspect ratio. If the aspect ratio is equal to y/2, that of the most rapidly growing disturbances, then the critical Rayleigh numbers are 657.5 for the free boundaries, 1708 for the rigid boundaries, and 1108 for mixed boundary condition. The values given in Table B . l are appropriate for conditions in our experiment and are used for calculating the Rayleigh number for carbon dioxide-water mixture. Chemical diffusivity, viscosity and reference density are that for pure water, ft is given by [Yamane and Aya, 1995], and concentration of gas in water c is obtained from the calculated phase diagram. 0 is given in units of mass fraction - 1, and c, correspondently, is in mass fraction. Appendix B. Mechanism of Gas Transference in Fluid 123 The calculated value of the Rayleigh number using the parameters from Table B . l is 101 3. This value is almost 1010 times higher than critical. This means that the convection tends to be vigorous. This result can be compared with experimental data on gas solubility (Chapter 4). We have some useful data on resistance as a function of time at constant pressure and temperature conditions (Figure 4.7). At constant temperature, resistance solely depends on dissolved gas concentration. Therefore, this curve represents change of the gas con-centration in the fluid with time. Let's set up a physical picture of the process. Consider a fluid, with mixed mechanical boundary condition (the lower boundary is fixed, while the upper is free). This fluid has zero initial concentration of dissolved gas, a fixed gas concentration on the upper boundary, and zero flux on the lower boundary. This concentration difference affects the fluid buoyancy leading to convective motion called compositional convection. This means that gas is transported into liquid at a rate which is much greater than that provided by the process of diffusion. Convective motion ex-ists in the system as long as the density variation provides a Rayleigh numbers which is higher than critical. The density variation is due to difference in gas concentration between top and bottom. The concentration of gas is fixed on the top, but it increases on the bottom (where concentration flux is zero). Therefore, Ap decreases with time, and so does the Rayleigh number. This process continues until the gas concentration all over the container is equal to the gas solubility (which is the fixed concentration value on the upper boundary). This means that the flat part of the curve (Figure 4.7) represents a stationary state of the fluid motion when the fluid is well mixed, and concentration of dissolved gas is equal to the gas solubility ceq at given P-T condition (constant P and T). In this case, the density of the fluid is p = p0(P,T)(l + /3ceg) Appendix B. Mechanism of Gas Transference in Fluid 124 The boundary layer analysis, which is applicable for high Rayleigh numbers, provides us with a quantitative predictions about how a convective fluid reaches its stationary state. B.3 Boundary Layer Analysis For large values of the Rayleigh number, a boundary layer analysis can be used to deter-mine the convective transport. This analysis is normally done for the thermal convection. In this case, the density variation is due to heat on the lower boundary, but this is equi-valent to the input of a heavier compositional component on the upper one (see the equation of state (B.2)). Therefore, in final expressions of the analysis, we will substitute temperature with concentration, keeping in mind that, for compositional convection, the driving force of the fluid motion has opposite direction. A fluid layer of thickness d heated from below, with a temperature drop across the entire layer AT (see Figure B.l) . If we assume that the core of the convection cell is isothermal, then the temperature drop across a single boundary layer, width 6, is \AT. The heat flux across the boundary layers at the top and the bottom is _ k&T _ k AT Qconv - 2 ~ 2T"' ( ' where k is thermal conductivity. If the fluid layer were not convecting, then a linear conductive profile (dashed line, Figure B.l) develops, with a heat flux AT Qcond = k~^~- (B-9) The intensity of convection is measured by the Nusselt number which is the ratio of the total heat flux to that occuring in the system when convection is absent. At high Rayleigh numbers, the total heat flux is roughly equal to that of convection. Thus, using Appendix B. Mechanism of Gas Transference in Fluid 125 T e m p e r a t u r e Figure B . l : Schematic profile of temperature distribution in convective layer. expressions (B.8) and (B.9), the definition of the Nusselt number becomes N u = Q c f n v = d ( B 1 Q ) qCond 2o For steady-state convection, the solution of the governing equations depends on a single dimensionless number Ra. Therefore, Nu = f(Ra), (B.ll) On the basis of expressions (B.6) and (B.2), the Rayleigh number for thermal convection Appendix B. Mechanism of Gas Transference in Fluid 126 For vigorous convection, the boundary layer thickness 6 does not depend on the depth d. Since Nu is a function of Ra that must also be linearly dependent on d, according to (B.10), we require Nu « Ra* In the absence of convection, when Ra<Ra c r Nu = 1, which implies that (B.13) \ JXQJCT J Therefore, on the basis of expressions (B.9), (B.10), and (B.13) the heat flux for a vigorously convecting fluid is (RaykAT q-={Ra-) — ( R 1 4 ) Using expression for the Rayleigh number in (B.12), the convective flux can be written as 3 \ 3 « - = ( ^ - l (AT)I (B.15) \ KpRacr J It is also possible to derive an expression for characteristic velocity in vertical direction using the boundary layer analysis. The convective heat transport in fluid of density p, heat capacity cp and vertical velocity v is Qconv = pepATv (B.16) Combining this expression with equation (14), we get ( Ra \ 3 kAT , „. ^ T V = { R £ ) — ( A I 7 ) Appendix B. Mechanism of Gas Transference in Fluid 127 which leads to \RO7J d« ( R 1 8 ) where the relationship for the thermal diffusivity k was used. Expression (B.18) can be represented as _ / Ra Y CR-\Q\ where vcond (equal to n/d) is characteristic velocity for conductive heat transport. The boundary layer approach allows us to describe the cooling of a layer of convecting fluid. The layer cools as heat is removed by convection: dT •' . —pcp-g£ • Volume = qconv • Area (B.20) Combining this expression with expression for convective heat flux (B.15), the evolution for cooling layer is given by the next equation dT (apq \ 3 k / m , a / T , . sF = -c'U?) I^'' <B'21> where Cj is constant 2 4 / 3 ci Rail3 and T is the temperature of the isothermal core. Temperature T decreases with, and represents the bulk cooling of the fluid layer. At initial moment t=0, core temperature is equal to the temperature change across the top boundary, e.g. A T T0 = — • (B.22) Characteristic time for convective cooling r 0 can be calculated on the basis of dimen-sional arguments, supposing that Appendix B. Mechanism of Gas Transference in Fluid 128 Therefore, with evolution equation (B.21), the characteristic time is c\ \ up ) d At the initial time t=0, the Rayleigh number can be written in terms of To as Ra = (B.25) KU-So that, after some manipulation, r 0 becomes T0 = c2[^-jRa-s, (B.26) where * = w <B-27> As d2 jK is thermal diffusion characteristic time, the difference between convective and thermal diffusion time scales can be estimated. In order to simplify the expressions for convective cooling, let's write equation (B.21) in dimensionless form using the dimensionless time i, where 1 = l^n T^ ( B- 28) and the dimensionless temperature is defined by the temperature drop T = TAT. The resulting equation is - -3T 4/ 3, (B.29) which has the following solution: T{t) = TQ{l + t/r*0Y\ (B.30) where TQ = 3r 0 /25. Appendix B. Mechanism of Gas Transference in Fluid 129 B.4 Application of the Boundary Layer Analysis to Compositional Convec-tion in the Experiment In the case of compositional convection, composition is the source of buoyancy instead of temperature. The equation that describes the way that fluid reaches equilibrium gas concentration is c(t) = ceq(l + t/r*)-\ (B.31) where ceq is solubility, and TQ is timescale in terms of ro. The changing gas concentration will cause a change in the measured resistance of the fluid. An approximation of resistance change with the time for convective transport of gas can be represented by R(t) = ROQ + (R(0) - R(oo))(l + i / r*)" 3 . (B.32) This result can be compared to the data presented in Figure 4.7. According to data, a value of R(oo) which corresponds to resistance in solution with the reached solubility for gas concentration is R{oo) = 4740, and initial value R(0) is R(Q) = 15677. The dashed line represents an interpolation curve (B.27), where convective timescale TQ* is roughly equal to 200 minutes, which corresponds to characteristic time To T0 « 6000s. According to expression (B.24), the Rayleigh number is expressed in terms of To as „ (d2 1 \ 3 Racr Appendix B. Mechanism of Gas Transference in Fluid 130 Density coefficient for concentration, f3 0.24 Concentration of C 0 2 , c 0.004 Distance for density varition, h 0.04 Permeability, kv 2 • lO"1 2 Porosity, (j) 40% Table B.2: Porous medium Rayleigh number parameters (SI units) Since r0=6000 s, and Ra c r=10 3, Ra = 1.1 x 101 3. Therefore, the Rayleigh numbers calculated through experimental data and on the basis of the gas concentration change found from the phase diagram are in agreement. The heat transport due to convective motion is times higher than that for diffusive transport (see expression (B.19)). This means that the timescale for convection is 10 1 0/ 3 times shorter than that for diffusion. B.5 Gas Transport in Porous Medium of the Experiment The Rayleigh number for convection in porous medium is typically defined by [Phillips, where h is a thickness of a fluid layer in sediments with porosity (j) and permeability kv. The appropriate values of the parameters representing conditions in the experiment are given in Table B.2. Ra \ 3 1991]: Ra = Ppocgkyh D<j)u (B.34) Appendix B. Mechanism of Gas Transference in Fluid 131 20000.0 15000.0 E O 10000.0 5000.0 0.0 0.0 200.0 400.0 600.0 Time (min) 800.0 1000.0 Figure B.2: A comparison of theory and experiment for a convecting fluid evolving toward a stationary state. The density variation is caused by a change of gas concentration due to the process of hydrate formation. Since hydrate formation consumes gas from the fluid, it can lead Appendix B. Mechanism of Gas Transference in Fluid 132 to unstable density variation. However, it appeared that buoyancy is not enough to drive the fluid motion. The Rayleigh number for the experimental setup is 1.4, while the critical value for instabilities to grow is around 10. Therefore, in our experiment, the gas transport in the sediments is due to the process of diffusion.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Theoretical model of hydrate formation in natural environments
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Theoretical model of hydrate formation in natural environments Zatsepina, Olga 1997
pdf
Page Metadata
Item Metadata
Title | Theoretical model of hydrate formation in natural environments |
Creator |
Zatsepina, Olga |
Date Issued | 1997 |
Description | The solubility of two hydrate-forming gases, C02 and CH4 , is calculated over a range of pressure and temperature. The solubility of gas is shown to be significantly altered by the presence or absence of hydrate. In particular, gas solubility changes abruptly in hydrate presence, allowing it to crytallize from the aqueous solution without the need of any free gas. To test this prediction a set of experiments was performed. In the experiments, an aqueous solution of CO2 was cooled at a pressure of 2 MPa. A variety of methods were examined to detect the growth of hydrate. With cooling and hydrate formation, the physical characteristics of the porous medium (temperature, porosity, gas concentration) change. On the basis of known governing equations and conductivity of an aqueous solution in porous medium, the conductivity change due to hydrate formation was predicted. Conductivity was found to be particularly sensitive to hydrate formation, so electrical potential measurements were used to monitor hydrate growth. These electrical measurements indicated a pronounced resistance increase due to a change of gas concentration in the solution, corresponding to the amount of hydrate produced. Hydrate growth in the system was also detected in temperature data, which indicated a release of latent heat. The calculated phase diagram at typical pressure and temperature conditions in marine environments is applied to establish the gas concentration needed to stabilize hydrate. This information is used in simple methods of hydrate formation to estimate the vertical distribution of hydrate in marine sediments and the rate of accumulation. |
Extent | 5287647 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-03-09 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0053015 |
URI | http://hdl.handle.net/2429/5803 |
Degree |
Master of Science - MSc |
Program |
Earth and Ocean Sciences |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1997-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1997-0115.pdf [ 5.04MB ]
- Metadata
- JSON: 831-1.0053015.json
- JSON-LD: 831-1.0053015-ld.json
- RDF/XML (Pretty): 831-1.0053015-rdf.xml
- RDF/JSON: 831-1.0053015-rdf.json
- Turtle: 831-1.0053015-turtle.txt
- N-Triples: 831-1.0053015-rdf-ntriples.txt
- Original Record: 831-1.0053015-source.json
- Full Text
- 831-1.0053015-fulltext.txt
- Citation
- 831-1.0053015.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0053015/manifest