THEORETICAL AND EXPERIMENTAL INVESTIGATIONS INTO T H E FORMATION AND ACCUMULATION OF GAS HYDRATES By Alan Rempel B. A. Sc. (Engineering Physics) University of British Columbia , 1 9 9 2 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L 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 T H E D E G R E E O F M A S T E R O F S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S Department of G E O P H Y S I C S We accept this thesis as conforming to the required standard 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 December 1994 © Alan Rempel, 1994 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. The University of. British Columbia Vancouver, Canada DE-6 (2/88) Abstract The substantial volumes of gas hydrates found in the Arctic and in marine sediments are both a possible source of global climate change, and a potential future energy resource. The rate at which a hydrate layer forms, and the spatial distribution of hydrate in the layer are controlled by the physical conditions of the formation environment. To better understand the physical conditions that affect hydrate layer characteristics, I present a quantitative model for the formation of hydrates in a porous medium. The theory is tested using the the results of laboratory simulations of the modelled conditions. Conservation principles are used to derive the full set of governing equations using the minimum number of assumptions and simplifications. Scaling arguments, based on estimates of physical parameters in marine sediments, show that both heat and mass transport are dominated by diffusive processes, so advection may be neglected in most formation environments. Analytical solutions to the leading-order set of equations are obtained for the case of a porous half-space cooled on its boundary. These solutions provide estimates of the growth-rate of a hydrate layer and the volume fraction of hydrate present. The model predicts that the layer grows on the thermal diffusion timescale with the phase-change interface moving at a rate which is proportional to the square root of time. The predicted hydrate volume fraction is determined by the rate at which compositional diffusion can supply gas to the moving interface. For the formation of a methane hydrate layer, the model generally predicts a hydrate volume fraction that is less than 1%. The modelled conditions are simulated in a reaction chamber constructed from a cast acrylic tube, 0.7 m in length, with an inner diameter of 0.14 m. To test the apparatus, n experiments were conducted in which the growth-rate of an ice layer in the sand-filled reaction chamber was monitored using RTD temperature probes. These experiments demonstrate that a simplified version of the hydrate formation model describes the for-mation of an ice layer in a porous medium. C02 was used as the hydrate former to test the predictions for the growth-rate of a hydrate layer and the layer's hydrate content. C02 has a high solubility in water, and the model predicts a much greater hydrate volume frac-tion for a CO2 hydrate layer than that for a low solubility gas such as methane. During hydrate formation experiments, the temperature probes were unsuccessful at detecting the position of the hydrate phase-change interface. At the end of some experiments, the pressure was reduced to dissociate the hydrate so that the presence of hydrate might be detected. However, the temperature probes failed to definitively identify the associated consumption of latent heat. This is despite the recovery of an iceThydrate mixture fol-lowing one such experiment. Additional instrumentation employing acoustic or electrical techniques may be necessary to quantitatively assess the model predictions. 111 Table of Contents Abstract ii List of Tables vii List of Figures viii Acknowledgement x 1 Introduction 1 1.1 Background • 1 1.2 Scope of Study 5 2 Theoretical Development 8 2.1 A Physical Model 8 2.2 Mathematical Formulation of Problem . . . 10 2.2.1 Conservation of Momentum '. 10 2.2.2 Conservation of Mass 11 2.2.3 Conservation of Energy 12 2.2.4 Conservation of Gas 15 2.2.5, Summary 18 2.3 Nondimensional Equations 19 3 An Analytical Model of Hydrate Formation 24 3.1 Problem Description 24 iv 3.2 Method of Solution 26 3.3 Theoretical Predictions . 33 4 Experimental Work 41 4.1 Experimental Apparatus 41 4.2 Formation of an Ice Layer 45 4.2.1 Experimental Method 46 4.2.2 Results . . 46 4.3 Formation of a Hydrate Layer 55 4.3.1 Hydrate Detection . 55 4.3.2 Experimental Method 58 4.3.3 Analysis of Results 59 5 Discussion 66 6 Conclusions 73 References 76 Appendices 81 A Nomenclature 81 B Interface Conditions 84 C Stability Analysis 90 C l Interface Stability 90 C.2 The Importance of Thermal Convection 98 D Effective Thermal Diffusivities 100 v D.l Fluid Saturated Sediments 100 D.2 Sediments Containing Ice 102 E Ice Saturation Level 106 F Data Acquisition 111 G Hydrate Layer Formation Above a Developing Icy Layer 114 vi List of Tables 3.1 Physical Properties 33 3.2 Dimensionless Model Parameters 34 4.1 Ice Layer Data 49 4.2 Properties used for Modelling Experiments 51 4.3 Ice Layer Predictions 52 E.l Data For Calculating Approximate Ice Volume Fraction 110 G.l Model Parameters for a Hydrate Layer Above a Developing Icy Layer . . 121 G.2 Sample Calculation of Growth Rates 122 vii List of Figures 2.1 Physical Model of Hydrate Formation 9 3.1 Modelled Formation Conditions 25 3.2 Dependence of h and A on c/, 36 3.3 Dependence of h and A on <b<x 37 3.4 Dependence of h and A on To 38 3.5 Dependence of h and A on STH • 39 3.6 Dependence of h and A on Teh 39 4.1 C02 Hydrate Phase Diagram 42 4.2 Experimental Apparatus . 44 4.3 Temperature Measurements During Ice Formation 47 4.4 A Comparison of Ice Interface Motion With Predictions 52 4.5 Predicted Temperature Records for Hydrate Formation Above 0 °C . . . 57 4.6 Temperature Records From a Hydrate Formation Experiment 60 4.7 Temperature Record From a Depressurization Experiment 61 C.l Wave Number for Marginal Stability as a Function of Time . 97 C.2 Wave Number as a Function of Growth Rate 98 C. 3 Wave Number for Marginal Stability as a Function of <j>a 99 D. l Rate of Temperature Change vs. Time in Fluid Saturated Sediments . . 102 D.2 Rate of Temperature Change vs. Time in Sediments Containing Ice . . . 104 D.3 Predicted Times of Maximum Rate of Temperature Change ....... 105 viii E.l Melting Behavior of an Ice Layer 108 E. 2 Melting Experiment to Determine Ice Saturation Level 110 F. l Circuit For Finding Changes, in Probe Resistances 112 G. l Physical Model of Hydrate Formation Above an Icy Layer . 115 G.2 Interface Positions of Ice and Hydrate Layers 122 IX Acknowledgement Many thanks go to Bruce for introducing me to the concept of freezing gassy liquids above 0 °C. I also owe a debt to several other people in the geophysics department for helpful discussions during the past two years. Dieter's efforts in constructing our apparatus were much appreciated. Finally I must acknowledge the support of many friends both inside and outside the department for reminding me of the interesting properties carbonated liquids that aren't frozen. Chapter 1 Introduction 1.1 Background Natural gas hydrates belong to the class of nonstoichiometric compounds known as clathrate hydrates (hereafter referred to simply as hydrates). These crystalline struc-tures consist of ice-like lattices of water molecules stabilized at high pressures and low temperatures by "guest" gas molecules occupying cavities between the water molecules. An important property of hydrates is their ability to store large amounts of gas in a relatively small volume. One cubic meter of methane hydrate, the most commonly oc-curring variety in the shallow subsurface, may contain the equivalent of up to 164 m3 of methane gas at standard temperature and pressure conditions [Kvenvolden, 1993]. Total worldwide hydrate reserves are estimated to contain 1019 g of carbon; about twice the carbon in all the conventional coal, oil, and gas reserves on land [Appenzeller, 1991]. Despite these impressive statistics, many questions remain unanswered as to the origin and distribution of naturally occurring hydrates. Upon discovery early in the last century, hydrates were regarded mainly as a labo-ratory curiosity, with initial studies concentrating on finding new hydrate formers and characterizing hydrate properties, lattice structures, and phase diagrams [see Jeffrey and McMullen, 1967, for a review]. Hammerschmidi, [1934] was the first to realize the po-tential impact of hydrates on the oil and gas industry when he showed that they could 1 Chapter 1. Introduction 2 form and disrupt pipeline flows. A further problem arose when exploratory drilling pene-trated methane hydrate deposits, causing them to dissociate and release large volumes of flammable gas. A field of research opened into the effects on hydrate phase equilibria of "inhibitors" which can help prevent disruption of gas flows through pipelines, and reduce the drilling hazard [Englezos and Bishnoi, 1988; Bishnoi and Dholabhai, 1993; Lederhos et al, 1993]. Recently, seismologists have identified vast reserves of hydrates located in the Arctic [Collett, 1993], and beneath the sea floor off continental margins [Collins and Watkins, 1985; Ginsburg et al, 1989; Hyndman and Spence, 1992]. Much of the current hydrate research is motivated both by the desire to extract the natural gas component for use as a clean-burning fuel [Cherskii and Bondarev, 1972; Collett, 1992], and to determine the effect that a possible release of the trapped greenhouse gas, methane, may have on the global climate [Nisbet, 1990; Paull et al, 1991; Loehle, 1993]. It has also been proposed that hydrates be used in such engineering applications as the storage of natural gas in subterranean reservoirs, the filtering of septic waste, and the desalination of drinking water [Sloan, 1990]. Before these potentials may be properly assessed, one must first understand the processes involved in hydrate formation. In order for formation to occur, there must be an adequate supply of gas and water molecules coexisting within the appropriate temperature and pressure regime. In the case of marine sediments, the sea floor temperature is relatively constant at a few degrees above zero. Since methane is the main guest constituent in naturally occurring hydrates, the methane hydrate phase diagram suggests that an excess of 300 m of hydrostatic head is necessary for hydrates to be stable at the sea floor. The pressure rises as the distance below the sea floor is increased, bringing about a corresponding increase in the hydrate equilibrium temperature. The base of the stability field is determined by the geothermal gradient; at some depth the temperature rises above the equilibrium temperature at that Chapter 1. Introduction 3 pressure, and hydrates may no longer exist. This depth is commonly associated with a discontinuity in acoustic impedance between partially hydrated sediment with a relatively high elastic velocity, and fluid saturated sediment, possibly including a small amount of free gas. The resulting seismic reflection, termed the BSR or bottom simulating reflector, is often used to map the extent of marine hydrate deposits [e.g. Hyndman and Spence, 1992]. Although most of the ocean floor is within the hydrate stability field, sufficient organic material for gas production is normally restricted to the continental shelf area. Natural gas may be produced either by the biogenic degradation of this organic material or by thermogenic production as a result of chemical reactions at greater depths [ Claypool and Kaplan, 1974]. One theory of hydrate formation holds that the gas contained in these deposits comes exclusively from in situ biogenic production, with negligible transport of gas into the stability field. Proponents of the in situ formation theory [e.g., Kvenvolden and Barnard, 1983] refer to geochemical data to justify their position. Isotopic analysis of methane from samples of hydrated sediments suggests that it is primarily biogenic in origin; furthermore, the concentrations of higher order hydrocarbons such as ethane and propane are lower than might be expected from natural gas of a thermogenic origin [e.g. Pflaum et al., 1985]. This suggests that gas migration from depths where thermogenic hydrocarbons are produced is insignificant. However, Hyndman and Davis, [1992], have argued that the relatively high hydrate volume fractions inferred from seismic obser-vations are incompatible with plausible levels of in situ gas production, implying that there must be some mechanism for concentrating the gas within the stability field. In their formation model, gas is carried by the pore fluid into the stability field where it is incorporated into the hydrate structure. While this theory appears to contradict the geochemical evidence, it is possible that some process such as isotope fractionation has biased the analysis of the recovered methane to make it appear to have a purely biogenic Chapter 1. Introduction 4 origin [Minshull et al, 1994]. In addition, the higher buoyancy of methane may allow it to be transported further than heavier hydrocarbons thus making the mixture of gases rich in methane. As yet no quantitative analysis of the implications of these competing theories on the properties of hydrate deposits has been reported. To be stable, a hydrate crystal must surpass some critical size or it will diffuse away. Once this size is reached, new hydrate may nucleate onto the existing structure much more easily. Some have argued that the low solubility of hydrocarbons in water requires that free gas be present for hydrate formation [e.g., Sloan, 1990]. However, while free gas should facilitate initial hydrate nucleation, its absence merely decreases the probability of nucleation. For systems operating on geological timescales this should not pose a problem. In addition, the porous media may play a role in collecting dissolved gas molecules. Laboratory research has revealed that once primary nucleation occurs, kinetic effects are not particularly important in subsequent hydrate growth [Sloan, 1990]. Indeed Mashirov et al., [1991] have reported success in simulating hydrate formation from dissolved gas. To date, aside from thermodynamic predictions of hydrate equilibrium conditions [e.g. Englezos, 1992], mathematical analysis in this field has been mainly restricted to pre-dicting the dissociation rates of hydrate deposits [Selim and Sloan, 1989; Tsypkin, 1991; 1992a; 1992b]. For the most part this research has been motivated by a desire to assess the viability of various methods of gas recovery in Arctic locations. The dissociation is forced to proceed rapidly so that large amounts of free gas are liberated and this gas phase becomes more important than the liquid phase to the dynamics of heat and mass transport. In contrast, formation probably occurs from small amounts of gas dissolved or suspended in liquid, so the system is expected to be dominated by heat and mass transport within the liquid. Major physical differences between the slow formation and the rapid dissociation of hydrate deposits necessitate the development of a new formu-lation to understand the generation and growth of a hydrate layer. In the dissociation Chapter I. Introduction 5 problem, the forcing conditions which promote the phase change are imposed and hence well known. The in situ hydrate formation conditions are not well known, but compar-isons between theoretical predictions and observations will help to resolve some of this ambiguity. 1.2 Scope of Study The goal of this work is to develop a quantitative model of hydrate formation and accu-mulation within a porous medium. Such models are needed to properly assess the amount and distribution of hydrate that may be expected to exist under varying physical condi-tions, and the time required for a layer to form. Seismic investigations are normally able to detect only the BSR, which marks the base of the stability field. Model predictions may be used to improve estimates of the total hydrate content above the BSR, which is largely unconstrained by seismic data, but vital in the assessment of the hydrate resource potential. The timescale required for layer formation has important implications for global change and ice-age dynamics. It has been argued that global warming might warm the sea floor, leading to a release of methane from hydrates, which would accelerate the greenhouse effect. On the other hand, global warming would melt polar ice caps, which would raise sea level and allow larger regions of the continental shelf to support hydrate growth, providing a sink for methane and a lessening of the greenhouse effect. During an ice age, the drop in sea level could destabilize hydrate deposits, and the corre-sponding release of greenhouse gasses could help to end the ice age. In another scenario, the increased pressure in polar regions beneath ice sheets might allow more hydrates to form, and less methane would escape from the tundra, resulting in lower concentrations Chapter 1. Introduction 6 of greenhouse gasses in the atmosphere, and a prolonged ice age. These competing hy-potheses motivate several questions. Which, if any, of these processes is important? Can a hydrate layer be reformed during the interglacial period? If the sea level drops by 100 m, what volume of hydrate becomes unstable? These are the types of problems that may only be properly addressed if the spatial distribution of hydrate in the sediment column and the time required for a hydrate layer to form are first determined. The approach taken towards investigating the problem of hydrate formation is to compare theoretical formulations based on conservation principles with data from labo-ratory simulations. Several difficulties are apparent from the outset. The mechanics of hydrate formation are not well understood. Is it an in situ formation process, or does the gas migrate from below? Is free gas required, or can hydrates form from gas dissolved in the liquid? Even if these processes do occur in nature, is it possible to simulate them in the laboratory and to accurately detect what is happening? Is the mathematical prob-lem even tractable? If these difficulties are overcome and the modelling is successful we will have developed quantitative tools to contribute to the study of more fundamental problems. The second chapter begins with a simple physical model describing a plausible hydrate formation environment. This is then translated into a set of mathematical conditions which must be satisfied in order to uphold the conservation laws. Simplifications are made based on scaling arguments, and dimensionless equations are presented; In Chapter 3 the solution to a reduced set of equations is presented for a scenario which may be easily replicated within the laboratory. A description of the laboratory apparatus and measurement techniques is found at the beginning of Chapter 4. Experimental methods and results are then produced for a related problem, the development of an ice layer in a porous medium. To complete the chapter, the results of laboratory efforts to simulate the modelled hydrate formation conditions are presented. Chapter 5 contains a discussion Chapter 1. Introduction 7 of the implications of the theoretical predictions and laboratory simulations. The work is summarized in Chapter 6, along with a few brief recommendations on the course of future research efforts. Chapter 2 Theoretical Development 2.1 A Physical Model A physical model of the hydrate formation process must be adopted prior to the derivation of the governing equations. Figure 2.1 shows a schematic diagram of a plausible hydrate formation environment based on conditions found in marine sediments. The pressure P increases with depth primarily due to the increased hydrostatic pressure, although there may also be a small non-hydrostatic component of the pressure gradient VP' which drives fluid flow. The temperature T increases with depth because of the geothermal gradient, and causes the temperature to surpass the equilibrium hydrate temperature at some level. A boundary results between sediments which may contain hydrate, above, and sediments which may not, below. Gas is transported through the sediment column by advection, along with or separate from the pore liquid. The gas will travel at a rate which differs from the liquid transport velocity u only if a significant portion of the pore space is occupied by free gas [Honarpour et al., 1986]. If It is dissolved, or if only a few percent of the pore volume contains free gas, it will travel at u. Since gas concentrations are normally thought to be quite low, the mixture of liquid and gas is modelled as an incompressible fluid, travelling at velocity u through the homogeneous, stationary porous medium. Gradients in gas concentration c also drive the transport of gas by diffusion and dispersion. Once the gas crosses the base of the stability field, hydrate nucleates and begins to fill 8 Chapter 2. Theoretical Development i Fluid and Sediment Mixture Moving Boundary Fluid, Sediment, and Hydrate Mixture Base of Stability Field * Fluid and Sediment Mixture Figure 2.1: Physical Model of Hydrate Formation The hydrate formation environment is modelled as a homogeneous sediment column. A hydrate layer begins to form at the base of the stability field where the pressure P and temperature T permit it to develop. Gas may be brought to the layer by the fluid moving with velocity u, and by chemical diffusion down the concentration gradient Vc. Chapter 2. Theoretical Development 10 a portion of the pore space. The hydrate saturation level h should affect both the bulk thermal properties, and the effective permeability of the fluid, sediment, and hydrate mixture. A hydrate layer develops as gas is transported up the sediment column. The position of the moving boundary and the hydrate pore fraction distribution are controlled both by the system's ability to discharge the latent heat of formation, and the speed with which gas can be supplied to the hydrate structure. 2.2 Mathematical Formulation of Problem The unknown fields in this problem, such as the fluid velocity u, temperature T, pressure P, and gas concentration c are governed by conservation equations. The conservation laws for momentum, mass, energy, and gas may be applied to the physical model to generate a mathematical representation of the system. An attempt has been made to keep the derivations as general as possible, leaving simplifications for section 2.3. A list of symbols is contained in Appendix A. 2.2.1 Conservation of Momentum The momentum conservation equation takes the form of Darcy's law [e.g., Phillips, 1991, p. 27], u = -g(h)-VP', (2.1) V which states that the fluid velocity u is proportional to the non-hydrostatic pressure gradient VP'. The constants of proportionality are the permeability of the sediment k and the dynamic viscosity of the fluid 77. The function g(h) is a relative permeability function which includes the effects of changes in the volume fraction of hydrate in the pore space h. No fluid may flow when the pore space is completely cemented by hydrate so #(1) must be 0. At the other extreme, when no hydrate is present the flow is restricted Chapter 2. Theoretical Development 11 by k, the intrinsic permeability of the sediment, so g(0) must be 1. In deriving the remaining three conservation equations, the balance relations are applied to a fixed arbitrary representative volume ft with surface element dS. The enclosed material is treated as a continuum of solid grains, pore fluid, and hydrate. 2.2.2 Conservation of Mass Any change in mass within the fixed volume ft must be accounted for by net mass flux across its boundaries. The total mass is simply the sum of the masses of the individual components contained within the volume, which is expressed as / [Pj<j>{i -h) + Ph<f>h + P,(i - $)] dn , - (2.2) where <f> is porosity, h is the fraction of the pore space occupied by hydrate, and pj, ph, and p, are the densities of the fluid, the hydrate, and the solid grains, respectively. Since the sediment matrix and the hydrate remain stationary, the only way for mass to escape the volume is through the fluid. The net mass transport across the bound-ary depends on the orientation of the interstitial fluid velocity vector v relative to the boundary. The surface element dS is a vector with unit outward normal, so the velocity with which the fluid exits ft is v • dS. Thus, the mass flux into ft is given by - Jspf(f>(l-h)v-dS = -Jspfu-dS , (2.3) where u is the macroscopic transport velocity appearing in Darcy's law, equation (2.1). This surface integral may be transformed into a volume integral using the divergence theorem to give - / pfu • dS = - f pfV • u dti . (2.4) J S vfl The mass balance becomes d_ dt f [Pf<f>(l -h) + Ph<f>h + pa(l -<f>)]dQ = - f PfV -udQ. (2.5) J o Chapter 2. Theoretical Development 12 The volume ft, is arbitrary, and the two integrands are continuous functions, so the integrands themselves can be equated, yielding jt[pfcb(l-h) + ph<bh + p,(l-<[>)] =-pfV-u. (2.6) The incompressibility condition means that the densities pf, ph, and pa are constant; thus the mass conservation equation reduces to <f>(pf - Ph)^ = PfV'U . (2.7) Equation (2.7) reveals that a density difference between the hydrate and fluid, will re-sult in a divergence of flow whenever hydrate is produced. However, when the density difference is small, and the rate of hydrate production dh/dt is small, the resulting fluid flow will be approximately solenoidal (e.g. V • u « 0). 2.2.3 Conservation of Energy The first law of thermodynamics [e.g. Denbigh, 1981, p. 19] dU = dQ-dW + fide , (2.8) states that changes in a system's internal energy dU must be balanced by heat flux into the system dQ, work done by the system dW, and the change in chemical energy p,dc associated with mixing system components. In the hydrate formation problem, the mixing term ndc is expected to be small since only small amounts of gas are available j and the change in gas concentration dc must therefore be small. If the work done against viscous forces is also negligible, only pressure work need be considered, and (2.8) may be rewritten as dU = dQ-PdV. (2.9) Chapter 2. Theoretical Development 13 Hydrate formation in nature occurs at roughly constant pressure; consequently it is useful to rewrite (2.9) in terms of the enthalpy of the system H which is defined as H = U + PV. (2.10) In differential form (2.10) becomes dH = dU + PdV + VdP (2.11) and, upon substitution into equation (2.9), it gives dH = dQ + VdP . (2.12) Thus, for constant pressure processes, conservation of energy may be written as dH = dQ , (2.13) which indicates that changes in enthalpy must be balanced by heat flux. This principle may be applied to the control volume ft where the change in total enthalpy is equivalent to the sum of the enthalpy changes for the individual components. With the specific enthalpies of the fluid, hydrate, and solid grains represented by Hf, Hh, and Ha, the total enthalpy change within ft is equal to ^ / [pfHfci>{l -h) + phHh<j>h + p,H.(l - <f>)] dCl . (2.14) at Jet Heat flux across the boundary dS may result from energy carried by the advecting fluid, or from conduction of energy through the bulk material. The amount of energy transported across the boundary by the advecting fluid is a function of the relative orientations of the velocity vector v and the surface element dS. The total advective transport into ft is - f pfHf(j){l - h)v • dS = - / V • (pfHfu) dn , (2.15) Js Jo Chapter 2. Theoretical Development 14 where the divergence theorem has been used once again to convert the surface integral into a volume integral. The thermal conductivity of the bulk mixture is a function of the volume fraction of hydrate in the pore space h. Denning the effective thermal conductivity as Ke(h) enables Fourier's law for the conductive energy flux into ft to be written as J Ke(h)VT • dS = J V • (Ke(h)VT) <fft . (2.16) Once again, the integrals in equations (2.14), (2.15), and (2.16), are calculated on the same arbitrary volume ft and the integrands are continuous, so equation (2.13) becomes p,4(l - h)^- + PHth?^ +Pt(l - 4 ) j j f + <f>[PH(Hh - H,) + H,(PK - PF)}^ +pfHfV • u = -pfVL • VHf + Ke(h)V2T + - n ^ V T • Vh . (2.17) da This version of the energy conservation equation, while complete, is not immediately useful because the enthalpy of a substance is not directly measurable. Instead it is desirable to rewrite (2.17) in terms of temperature. For a constant pressure process, the change in specific enthalpy is related to the change in temperature by dH = CpdT, , (2.18) where Cp is the specific heat capacity at constant pressure. The fluid, hydrate, and solid grains are characterized by different values of Cp. In addition, the heat capacity itself varies with temperature, however this thermal dependence may be ignored for small temperature changes. The difference between the enthalpy of the fluid and the enthalpy of the hydrate is equivalent to the latent heat of formation of the hydrate. If the latent heat per unit mass of hydrate is represented as Ln then Lh = Hf-Hh. (2.19) Chapter 2. Theoretical Development 15 Using equations (2.18) and (2.19), the energy balance becomes dT Oh dh [pfCf4>(l -h) + phCh<f>h + p.C.(l -<!>)]-£- ^PhLn-^ + <KPh - Ps)CiT—t -rPfCfTV-u = -pfCfU<VT + Ke(hW2T + ^VT-Vh, (2.20) dh where the specific heats of the fluid, the hydrate, and the solid grains at constant pressure are given by Cf, Ch, and C,. The notation may be simplified by defining the functions T(h) and Ke(h) as T { h ) _ p,c,fti--h) + Phch<iih + p.c.{i-4>) 2 2 1 ) and Ke(h) = ^ , (2.22) where T(h) is a normalized bulk heat capacity, and Ke(h) is the effective thermal diffusiv-ity. Furthermore, the divergence of velocity may be eliminated by using the expression for conservation of mass from equation (2.7). The final form of the energy conservation equation becomes r W f = « . ( 4 ) v ^ - . . v T + ^ » + . ^ v T . V » : (2,3) The term on the left is the change in heat content with time. The first two terms on the right account for diffusion and advection of heat, respectively. The third term is the heat production resulting from formation or dissociation of hydrate. The final term is present because any change in thermal properties due to variations in hydrate concentration will affect the manner in which heat is transported by conduction, much like changes in electrical resistivity affect the conduction of electrical energy. 2.2.4 Conservation of Gas The final governing equation is based on conservation of gas. In this problem, changes in the gas content within ft may occur due to advection, diffusion, and dispersion. If the Chapter 2. Theoretical Development 16 mass fraction of gas in the fluid is c, and the concentration of gas in the hydrate is c/», then the change in mass of gas within f2 is given by ^ Jn[pfC<b{\ -h) + phch<f>h] dtl . (2.24) In general, since hydrates are nonstoichiometric compounds, the ratio of gas to water molecules may vary within them. In practice, however, the ratio is nearly constant [Sloan, 1990] and thus Ch may be treated as a constant. The hydrate is stationary, but the pore fluid may carry gas. The term for advective gas flux is similar to equation (2.15) for advective energy flux, namely - / pfc<j>(l - hW • dS = - / V • (pfcu) dSl . (2.25) Gas may also be transported by molecular diffusion. When the gas concentration is small, the effects of the temperature gradient on the diffusive flux are small in comparison to those of the gradient of the chemical potential ft of the fluid mixture [Landau and Lifshitz, 1987, p. 232]. The diffusive flux i c may then be written as i c = -^(1 _ h)Vfi (2.26) where the gradient of the chemical potential is given by (2.27) The partial derivatives in (2.27) are replaced with the chemical diffusion coefficient Dc and the barodiffusion coefficient kp defined as and (2.28) (2.29) Chapter 2. Theoretical Development 17 Using these definitions, and neglecting the effects of temperature gradients gives \c = -Pf<t>{l-h)Dc{Vc+^VP) . (2.30) The term involving the pressure gradient is retained because of the large hydrostatic gradient. . The dispersive flux id depends on the dispersion coefficient Dd which is proportional to the fluid velocity, and to a correlation length which is comparable to the grain size of the sediment. This dispersive flux is expressed as id = -pf(f>(l - h)DdVc . (2.31) The combination of diffusive and dispersive flux of gas out of the control volume amounts to - f (ic + id) • dS = - / V • (ic + id) dtl Js Jn \ = / Pf<W • {(1 - h)[(De + Dd)Vc + -Dc-jjVP]} d£l: (2.32) The change in gas content, (2.24), is equal to the sum of the two transport terms from (2.25) and (2.32), so that /^alW1 + —chh] = -PfcV • u - pfu • Vc dt pf +Pf4>V • {(1 - h)[(Dc + Dd)Vc + £> cy VP]} . (2.33) Rearranging terms and using conservation of mass, (2.7), to eliminate the term involving the divergence of velocity gives (1 " h)^ = (1 - h)(Dc + 7Jd)V2c - i n • Vc+ ^ ( C - ch)fr (De + Dd)Vc-Vh + DekPV-{}-^VP). (2.34) Chapter 2. Theoretical Development 18 This equation is quite similar in form to the energy conservation equation. Changes in the fluid's gas content are caused by diffusion and dispersion, advection, transfer into the hydrate structure, and spatial variations in the hydrate volume fraction. One additional term appears, the final term on the right hand side which reflects the contribution of pressure driven gas diffusion. 2.2.5 Summary To summarize, the governing equations for the formation of gas hydrates in a homoge-neous porous medium are: u = -g(h)-VP' , (2.35) V dh PfV-u = 4>{P} - Ph)— , (2.36) r ( f t )f = K . w v* r - u • v r + + ^ f l v r • v , , ( 2 , r , and (1 - h)^ = (1 - h)(Dc + Dd)V2c - ±u • Vc + f{c - ch)^t -{Dc + Dd)Vc • Vh + DckPV • ( ^ ^ V P ) . (2.38) Analytical solutions to this set of four coupled nonlinear partial differential equations are not expected to exist in most circumstances, so in general numerical techniques would, be required to solve the system in this form. However, scaling arguments may be used to reduce the complexity of the system by identifying and eliminating insignificant terms. In the following section, simplifications based on estimated physical parameters are used to help generate the final set of leading order nondimensional equations used in subsequent analysis. Chapter 2. Theoretical Development 19 2.3 Nondimensional Equations In marine sediments, the temperature, gas concentration, pressure, and fluid velocity profiles are often only dependent on the depth below the sea floor. Restricting our attention to 1-D problems simplifies the mathematics considerably without neglecting any of the essential physics of the problem. The approximations which follow are used to further reduce the difficulty of the problem and produce a set of leading order equations to which analytical solutions may be found. Terms which are neglected are expected to have only a minor influence on the hydrate formation process, and their effects may be accounted for using perturbation techniques. As mentioned earlier, the mass conservation equation, (2.36), indicates that volume changes associated with hydrate production will produce fluid flow. However, the differ-ence in density between hydrate and water is less than 10%, so the volume change is not very large. Provided the rate of hydrate production dh/dt is not too high the flow will be solenoidal, and I can approximate (2.36) by V • u = 0 . (2.39) In one dimension this implies that the velocity does not depend on position (i.e. du/dz = 0). To further aid in the analysis, a number of characteristic, or representative scales are used to help assess the relative importance of the various terms in the energy and gas conservation equations (equations (2.37) and (2.38)). While the choice of scales is not unique, the physical processes involved in the problem render some scales more useful than others. Both energy and gas are transported by advective and diffusive mechanisms. The relative importance of advection and diffusion may be assessed by introducing a Chapter 2. Theoretical Development 20 diffusive time scale, «e(0) Ke and the corresponding velocity scale, _2 2 r - r (2.40) = ? = (2.41) Physically, ttP is the velocity required to traverse the length scale zr in the time scale tT\ if u is small compared with ur, then transportation by advection will be small. The associated nondimensional variables i, z, and it are defined as - t _ z , _ u t = — , z = — , and u = — . (2.42) <r : zr- uT The temperature and composition are made nondimensional using typical variations, AT and Ac, as representative scales. Thus the dimensionless temperature T and com-position c may be written as f = ~ ^ ' (2-43) and c = ^ , (2.44) where To and CQ are constant values which provide a convenient origin for these nondi-mensional quantities. With this choice of scales, the one-dimensional energy and gas conservation equations (refer to (2.37) and (2.38)) become df K.{h)&f • Of Oh 1 dKe(h) Of Oh m ^ = ^W-uM+SThTt+lTe^h-cn-bl> (2"45) and , / - . n ^ / - . 7 \ ^ 2 c ' -de . p i . . dh 1 - h)a-. = (1 - h)— - au- + * « - ( c - ch) dedh kPDc d 1-hdP dzdz (Dc-rDd)Acdr P 8z}'.. . 1 ' Chapter 2. Theoretical Development 21 where the nondimensional parameters Sxh, a, and cn are defined as: Ph<f>Lh T h ~ pfCfAT ' <b(Dc + Dd) ' and a = (2.47) (2.48) Ch - Co f Ch = -KT-- (2"49) The Stefan number Sxh measures the importance of the heat produced by hydrate pro-duction relative to the heat required to raise the system's temperature. If STK is large, for instance, the phase change should proceed slowly as the system struggles to trans-port away the latent heat. The Lewis number a indicates the relative ease with which heat and gas are transported by diffusion. The hydrate gas concentration parameter CH describes the amount of gas required to stabilize the hydrate structure. The importance of these three dimensionless quantities will be illustrated in Chapter 3. The coefficients of equations (2.45) and (2.46) reveal the relative importance of the various terms. In the temperature equation, for example, if u is much greater than ne(h)/Ke (e.g. the coefficient of the diffusive term), then advection dominates. Using the definition of the effective thermal diffusivity from (2.22), the coefficient of the diffusive term in (2.45) may be expressed as Ke(h) _ Ke(h) which is roughly order unity because the thermal conductivities of pure hydrate and pure water are quite similar [Sloan, 1990]. Thus, advection dominates the heat transport when u is large compared to unity; equivalently, the dimensional transport velocity u must be large compared with the reference velocity ur. Using the definitions of ur and tr from (2.41) and (2.40), advection dominates in the thermal equation when u » — . (2.50) Chapter 2. Theoretical Development 22 A similar argument for the gas conservation equation, (2.46), indicates that advection dominates if the coefficient, ocu, of the advective term is much greater than unity. Thus advection dominates if { > > 5 ( T b ) ' <2'51> or 4>(DC + Da) U zr(l - h) • ( 2 - 5 2 ) In the submarine environment, hydrate layers are of the order of 100 m thick. The typical changes in the temperature and gas concentration which are used to define T and c occur over this length scale. Typical thermal diffusivities are of the order of 10 - 6 m2/s. This implies that thermal diffusion dominates if u « — 10~8 m/s , and advection dominates if u » — « 10"8 m/s . Zr A typical porosity for unconsolidated marine sediments containing hydrates is 50% [ Yuan et al, 1994]. Experimental measurements of the chemical diffusivities of light hydrocar-bons in sedimentary rocks have been reported with values ranging from the order 10 - 6 to 10~12 m2/s [Kroos and Leythaeuser, 1988]. The diffusivity tends to be higher for smaller gas molecules and higher porosity sediments. Since methane is the smallest hydrocarbon molecule, and the porosity of marine sediments is high, a chemical diffusivity on the order of 10 - 8 m2/s is taken for this analysis; this suggests that provided the hydrate volume fraction is not close to one, chemical diffusion and dispersion are the dominant gas transport mechanisms if » « + D * ] * 5 x 10-"m/s , zr(l - h) Chapter 2. Theoretical Development 23 and advection dominates if m > > <f>(Dc + z y ^ 5 x 1 Q _ l l m / s z P ( l - ft) Values given in the literature for fluid transport velocities in marine sediments vary widely and hydrate layers are formed in different types of velocity regimes. Hyndman et al. [1993] report an upper bound of 4 x 10 - 1 1 m/s on the fluid velocity in the Cascadia accretionary margin where large deposits of hydrates are located. This suggests that diffusion certainly dominates in the temperature equation, and since this velocity is an upper bound, advection may be neglected in the gas conservation equation in most regions as well. The final term in (2.46) reflects the role of pressure driven diffusion of gas on the dynamics of the problem. Unfortunately the value of the barodiffusion coefficient kp is poorly constrained, and hence the relative importance of this process cannot be assessed. The following analysis proceeds on the assumption that compositional gradients dominate the diffusion of gas. In this case the nondimensional 1-D energy and gas conservation equations become df Ke(h)d*f dh I dKe(h)dfdh m ^ = ^w+SThWt+^e^k--dz-crz> ( 2 5 3 ) and ,/„ dc / H ,,92c . ph,_ _ . dh dcdh . ^ ( 1 _ A ) a _ , = ( 1 _ M _ + ^ ( c _ C k ) ( 2 . 5 4 ) In Chapter 3, solutions to these equations are presented for a problem which may be simulated in the laboratory to test the theory and the validity of the assumptions. Chapter 3 An Analytical Model of Hydrate Formation 3.1 Problem Description In this chapter, solutions to the governing equations are presented for the hydrate for-mation problem described below. This specific problem was chosen to satisfy a number of criteria. In order for the theoretical model to be quantitatively evaluated, it was im-portant that laboratory simulations of the problem were possible. Model conditions were chosen to be similar to those occurring in nature, so that the results would have practical merit. Finally, the boundary conditions and initial conditions were chosen to permit analytical solutions. These requirements are met by the problem depicted in Figure 3.1. I solve for the volume fraction of the pore space occupied by hydrate, and the rate of progression of the phase change boundary when a homogeneous porous half-space saturated with a gas-rich fluid is cooled on its boundary. Initially, there is sufficient gas present to form hydrate, but the temperature is too high to allow such structures to be stable. The initial conditions are uniform throughout the half-space with T = Too and c = C o o . At t — 0, the temperature at the bottom boundary is dropped below the equilibrium temperature and maintained at T = To. The developing hydrate-fluid interface at z = a(t) is in thermodynamic equilibrium, implying that its temperature Teh = Teh(P) is known as a function of pressure. Throughout the hydrate layer, the gas dissolved in the fluid is maintained at the equilibrium concentration ce = ce(P). If the concentration were to exceed this equilibrium value within the stability zone, additional 24 Chapter 3. An Analytical Model of Hydrate Formation 25 Homogeneous Halfspace Porosity High Initial Temperature T=T00 High Initial Gas Concentration c=c ft Moving Interface z = g(t) T = rehft Hydrate Layer Hydrate Pore Volume Fraction h c = ce Cold Boundary Figure 3.1: Modelled Formation Conditions The modelled hydrate formation environment is a homogeneous half-space with porosity <f>, initially at a constant temperature and high gas concentration CQQ. The boundary is cooled to temperature T 0 below the equilibrium temperature and a layer begins to develop with the moving interface at a(t) and a hydrate pore volume fraction of h. Chapter 3. An Analytical Model of Hydrate Formation 26 hydrate would be formed. Alternatively, if the gas concentration were to drop below ce, hydrate would dissociate until the system equilibrated. This problem is well-suited to the experimental apparatus described in Chapter 4. The uniform initial conditions and constant boundary temperature at z = 0 facilitate the derivation of analytical solutions. The conditions on the moving interface at z — a{t) are somewhat more complicated, but may be dealt with using the techniques described below. Pressure driven chemical diffusion is assumed to be unimportant, and the fluid is stationary within the porous medium so the problem is diffusively limited and equations (2.53) and (2.54) for conservation of energy and gas are valid. The implications of the predictions of this model on the natural formation of hydrates in marine sediments are discussed in Chapter 5. 3.2 Method of Solution The moving phase-change interface divides the half-space into two regions. The partially hydrated region between z = 0 and z = a{t) grows as the sediment column cools from below. In the region overlying z = a(t) hydrate is not permitted to form because the temperature exceeds the equmbriurn temperature Teh-Working with the nondimensional equations can simplify the algebra if the appropriate temperature and gas composition scales are chosen. The following definitions ensure that the nondimensional temperature T and gas concentration c introduced in (2.43) and (2.44), range between values of 0 and 1 throughout the half-space: (3.1) and c - ce (3.2) c = c, •oo Chapter 3. An Analytical Model of Hydrate Formation 27 The form of the nondimensional equations in (2.53) and (2.54) are unaffected by these definitions, but the expressions for the nondimensional parameters are dependent on (3.1) and (3.2). Explicitly, the Stefan number STH and the hydrate gas concentration parameter Ch become PfCfKloo - J-o) and eh = ^ — ^ . (3.4) C m Co The phase diagrams of common hydrates show that the equilibrium temperature is dependent on pressure, but for the pressure variations anticipated in the laboratory simulations the equilibrium temperature is nearly constant. The equilibrium gas concen-tration ce may also be a function of pressure but it will be treated as a constant. When ce is a constant, (2.54) reduces to | = 0 , 0 < i < * ( * ) . 4|> (3.5) so that for the hydrated region to be in chemical equilibrium, the hydrate volume fraction must be constant in time. Conservation of mass, (2.36), may be invoked to show that the flow velocity u must be solenoidal (i.e. du/dz = 0). Using this result and taking the spatial derivative of equation (2.35) for conservation of momentum reveals the following dependence of the hydrate volume fraction on position: dh_k(dg(h)yid*P> dz ~ r, { dh ) dz> • ^ Since the fluid is initially at rest throughout the half space, the non-hydrostatic pressure gradient must initially be uniformly zero. The only mechanism by which this gradient can change is through the volume change associated with hydrate production at the phase Chapter 3. An, Analytical Model of Hydrate Formation 28 boundary. However, when the volume change due to hydrate formation is sufficiently small to assume V • u = 0, the non-hydrostatic gradient will remain negligibly small for all time. Consequently, the hydrate volume fraction is constant throughout the stability zone. Applying this result to the conditions for conservation of energy and gas, (2.53) and (2.54), gives r 4 = a ? • » < * < « « • <3-7> dt d2f c = ce , 0 < z < a(i) (3.9) and , dc d2c .... ,„ *acTi = dT>' z > a { t ) (3-10) where the constants I\ and To are defined as and The initial conditions are: Ke(h) T0 = T(0) . (3.12) T(£,0) = 1, (3.13) c(5,.0) = 1 , (3.14) and 5(0) = 0 . (3.15) The boundary conditions at z = 0, a(t), and oo are: f(0,t) = 0 , (3.16) Chapter 3. An Analytical Model of Hydrate Formation 29 f(oo,i) = l , (3.17) c(oo,i) = 1 , (3.18) f(a,i) = f e h ^ ^ - ^ , (3.19) and c(a,i) = 0. (3.20) Two additional interface conditions are found from the application of the conservation equations in integral form across the moving phase-change boundary. The analysis, detailed in Appendix B, reveals that discontinuities in the temperature gradient and the gas concentration gradient must arise from the jump in the hydrate volume fraction from just below the interface at z = a" to just above it at z = a+. The volume fraction of hydrate in the layer h and the speed of the moving interface dd/di are related to the gradient discontinuities by df(a+,i) Ke(h)df(d-,i) _ da dz Ke dz T h di ' and S«£A = ^ a A § . (3.22) az pf dt Equations (3.7) through (3.10) along with the initial conditions, (3.13) through (3.15), boundary conditions, (3.16) through (3.20), and interface conditions, (3.21) and (3.22) form a complete mathematical description of the problem. General solutions for the temperature and gas concentration are used with the interface conditions to solve for the unknown hydrate volume fraction h and the growth rate of the hydrate layer da/dt. An examination of the equations and boundary conditions reveals that their is no inherent length scale in this problem. These are simply diffusion equations, and it should be possible to find similarity solutions in the variable £ defined by Chapter 3. An Analytical Model of Hydrate Formation 30 Furthermore, it is convenient to search for solutions where the interface motion is de-scribed by a(i) = 2\\fi, (3.24) such that the interface position in terms of the similarity variable is £ = .A. With this representation, the problem of determining a(t) is reduced to that of determining the unknown constant A. Substitution of the similarity variable reduces the energy conservation equations, (3.7) and (3.8), to second order ordinary differential equations, with associated boundary conditions r (oo) = l , (3.27) f(0) = 0, (3.28) and f(X) = feh, (3.29) and the interface condition mg) _ Kjvmg) = _ 3 af Ke d£ Integrating (3.25) and (3.26) twice and using the boundary conditions, (3.27), (3.28), and (3.29), to solve for the integration constants, yields the following temperature profile solution: f(0 = feheif^I , 0 < £ < A (3.31) erfTylfcA) Chapter 3. An Analytical Model of Hydrate Formation 31 and . f(0 = e r f ^ ) + r e " e ; f c e ( r J ^ A ) erfc(Vr~0 , <; > A (3.32) where the error function erf(a:), and complimentary error function erfc(a;) are defined by [Abramowitz and Stegun, 1972, p. 297] erffaj) = -4= [Xe-u2du , (3.33) y/ir Jo and 2 2 r°° a erfc(as) = -= / e"u <fo = 1 - erf(x) . (3.34) The interface condition, (3.30), relates the discontinuity in the temperature gradient to the hydrate volume fraction and the speed of the moving interface. Differentiating (3.31) and (3.32) allows (3.30) to be rewritten as Jj^-r,* 1 - _ J j T-r^ "eWTU = _ S T h h X . V 7r erfc(vr0A) V 7T K eerf(vl/iA) which is a transcendental equation for the two unknowns, A and h. The same similarity transformation is used to solve (3.9) and (3.10) for the gas con-centration. The equations become c = 0 , 0 < £ < A (3.36) and with boundary conditions and 2 ^ | + | f = 0 , £ > A (3.37) c(oo) = 1 , (3.38) c(A) = 0 , (3.39) Chapter 3. An Analytical Model of Hydrate Formation 32 and interface condition ^ l = 2 ^ c h f a h \ . (3.40) di Pi The gas concentration in the fluid can therefore be written as c(0 = 0 , 0 < i < A (3.41) and c(() = erft^aO - -i^jgL e r f c ^ ) , ( > A (3.42) The jump condition on the concentration gradient, equation (3.40), becomes [fa e - * * * 2 ph - , , , (n AO\ \ — — r / n-w = —ChfahX , (3.43) V 7r ertc(v <pct\) pf which is a second equation for the two unknowns, A and h. All that remains is to solve (3.35) and (3.43) for the hydrate volume fraction, and the parameter A, which is a measure of the interface velocity. Isolating h in (3.43) gives -<txx\2 pfe Phy^nfachX erfc(y^aA) Substituting (3.44) into (3.35) yields (1 - fe/OvTa Ke{h)fehy/Th PfS, (3.44) ero*2 erfc(Vr^A) ~ Keer^2erf(Vl\A) phy/fache^2 erfc(V^A) ~ ' 1 ' which is a transcendental equation involving A only. The value of A is found by plotting the left hand side of (3.45) and finding the zero crossing. Equation (3.44) gives the corresponding hydrate volume fraction, which is independent df position and time, as required. The temperature and gas concentration profiles described by equations (3.31), (3.32), (3.41), and (3.42), along with equations (3.44), and (3.45) constitute a complete solution to the problem. However, the complicated form of (3.45) makes it difficult to visualize Chapter 3. An Analytical Model of Hydrate Formation 33 Property Nominal Value Units Teh — T0 2.5 °G Too — To 10.0 °c Pi 1000 kg m"3 Ps 2650 kgm-3 aPH 930 kg m~3 hcf 4200 J kg-1 °c-j kg-1 °c-cc. 2200 dch 2080 j kg-1 °c-</> 0.5 dLh 430 x 103 J kg"1 10~6 m2 s"1 Dc + Dd 10"8 m2 s-1 ech 0.134 fce 0.0003 <-oo 0.001 Table 3.1: Physical Properties Sources of Information: a methane hydrate [Hyndman and Davis, 1992]; 6water at 5°C [Lide, 1990]; cmeasured in laboratory; ^ methane hydrate [Sloan, 1990, pp. 387, 389]; ccomposition CH^-h.lhH^O (8.4 moles methane/kg); ^take equilibrium composition 1/3 of solubility after Claypool and Kaplan, [1974]; 9solubility of methane at 880 PSI (6.06 MPa), 25 °C [Fogg and Gerrard, 1991, p. 116]. how the various parameters effect the dynamics of the process. In section 3.3, graphs are provided for the growth rate parameter A and the hydrate volume fraction h to test the sensitivity of the calculations to parameter values. 3.3 Theoretical Predictions Table 3.1 contains values for the physical properties which determine the nondimensional parameters appearing in equations (3.44) and (3.45). Table 3.2 includes both the nominal values for each of these model parameters, and a range of values which is used to test the sensitivity of the solution to parameter variations. Chapter 3. An Analytical Model of Hydrate Formation 34 Parameter Where Defined Teh (3.19) To (2.21), (3.12) Th (3.11) Ke(h) (2.22) Srh (3.3) <f>a (2.48) Ch (3.4) Definition Teh — Tp T — T PTCf<f> + p,C,(l-<f>) T(h)Ke Ke(h) Ke(h) Nominal Range Value 0.25 0-1 1.19 0.65-1.25 1.19 0.65-1.25 1 n.a. 4.6 1-20 100 10-1000 200 50-300 Ke(0) Ph<f>Lh PfCj(Too — T0) Dc + Dd Ch — Ce Coo Ce Table 3.2: Dimensionless Model Parameters Note that the values adopted for Th and T 0 in Table 3.2 are identical. Provided that the hydrate volume fraction is only a few percent, the phase change will have no noticeable effect on the thermal properties of the bulk matrix and Th should be approximately equal to To- By the same reasoning, the ratio of the effective thermal diffusivities Ke(h)/Ke which is defined to be equivalent to the ratio of the effective thermal conductivities, is taken to be 1. In each of the following five figures A and h are graphed as a function of one of the parameters in equations (3.44) and (3.45). The axis on the left side pertains to the growth rate parameter A, while the axis on the right pertains to the hydrate volume fraction h expressed as a percentage. In each case h is quite small; the calculation using the nominal parameter values from Table 3.2 yields h ?s 0.591% and A fa 0.205. Figure 3.2 shows the importance of Ch on the dynamics of the problem. The value of Ch is typically large because the concentration of dissolved gas in the fluid is much less r Chapter 3. An Analytical Model of Hydrate Formation 35 than the gas concentration in the hydrate. Lower values of Ch result when the initial gas concentration CQO is close to the saturation limit. This means that more gas is available for hydrate production and consequently h increases. The dependence of Ch on cn may be illustrated by considering the case where there is no gas transport and all of the in situ gas which exceeds the equilibrium concentration is converted to hydrate. In this limit, the gas balance condition is PfCoo = PhChh + PfCe{\ - h) , (3.46) so that the hydrate pore volume fraction becomes h= pfC°° ~ pfCe = , pf{cT - Ce) . (3.47) phCh — PfCe Ph\Ch — Ce) + (ph — Pf)Ce Since in practice the density difference between hydrate and water is small, and the hydrate gas concentration Ch is much greater than the equilibrium gas concentration ce, the second term in the denominator of (3.47) may be ignored, and h can be approximated by h « . (3.48) PhCh Equation (3.48) gives the minimum possible hydrate volume fraction for a given initial gas concentration. For the purposes of comparison, (3.48) yields a value for h of 0.538% using the nominal value of cV = 200, whereas the modelled value of h is 0.591%. The fact that (3.48) closely approximates the modelled value for h is a result of the low nominal value adopted for the Lewis number a. The growth-rate parameter A is also affected by the value of cv Since low values of Ch produce large values of h, the associated latent heat release retards the growth of the interface. Consequently low values of Ch cause low values of A as shown in Figure 3.2. The product <f>a is the ratio of the thermal and chemical diffusivities. Figure 3.3 demonstrates that as this ratio gets large, it becomes increasingly difficult to transport if Chapter 3. An Analytical Model of Hydrate Formation 36 0 . 2 0 6 0-205 -\ 0 2 C 4 O203 0202 O201 O 2 0 0 I - 25 Figure 3.2: Dependence of h and A on 4 The axis on the left is for the layer growth rate parameter A, where the interface position is a (i) = 2\Vt. The axis on the right is for the hydrate volume fraction h expressed as a percentage. Referring to equation (3.4), c/, is defined as Ch = (CH — c e )/(coo — c e ) . gas towards the stability field. The hydrate front moves quickly into the sediment column by converting the in situ gas into hydrate, with little transport from the overlying gassy region. The hydrate volume fraction h is limited by the in situ gas concentration and tends to the value given by (3.48) in the limit where <f>ct tends to infinity. In this extreme case, (3.45) reduces to 1 t-eh ° e - r 0 A 2 7T erfc(vTo~A) r f c _ Ke(h)ft eh PfSxhX (3.49) 7r Keerf(vTfcA) phCh which is simply the interface condition from energy conservation (3.35) with h replaced using (3.48). In (3.49), the layer growth rate parameter is independent of <j)ct, and hence the graph of A as a function of this product is flat on the right hand side of Figure (3.3). When the magnitude of the chemical diffusivity starts to approach that of the thermal diffusivity on the left hand side of the graph, gas is transported more efficiently towards the phase change boundary as it is converted to hydrate, and consequently h Chapter 3. An Analytical Model of Hydrate Formation 37 0.203 O.205 0.2C4 O-203 O202 0-2O1 a2or> 0.9 1000 Figure 3.3: Dependence of h and A on <j>ot The axis on the left is for the layer growth rate parameter A, where the interface position is a (i) = 2AV1 The axis on the right is for the hydrate volume fraction h expressed as a percentage. Recalling equation (2.48), 4>a = Ke/(DC + Da). must increase in order to maintain the layer's gas concentration at ce. The associated increase in latent heat is responsible for the modest drop in the growth rate parameter for low values of (j>a. To is proportional to the heat capacity of the bulk mixture. Figure 3.4 shows that higher values of T0 translate into a slower growth rate since a larger heat flux is required to reduce the temperature to the hydrate equilibrium value Teh- The slower growth rate results in a somewhat enhanced hydrate volume fraction because more gas can be transported by diffusion into the stability region. The Stefan number is a measure of the importance of the energy released due to the change in phase relative to the energy required to change the temperature. The larger this parameter becomes, the more difficult it is for the system to carry away the latent heat. Figure 3.5 shows how this translates into a slower rate of growth for the hydrate Chapter 3. An Analytical Model of Hydrate Formation 38 0.30-0.22-0.20-0 . 6 0 0.59 0.56 0.55 0.7 0.8 019 l.O 1 . 1 1 . 2 Figure 3.4: Dependence of h and A on To The axis on the left is for the layer growth rate parameter A, where the interface position is a (i) = 2\Vt. The axis on the right is for the hydrate volume fraction h expressed as a percentage. T0 is defined in (2.21) as T0 = {pfCf + p,C,(l — <f>)}/(PfCf)-layer and a slight increase in the hydrate volume fraction. Since the nominal hydrate volume fraction is so low, however, changes in Srh have little effect on either ft or A in this region of parameter space. The role that Ten plays on the dynamics of the system is shown in Figure 3.6. As the initial temperature is reduced towards the equilibrium temperature and Teh increases, the phase change is driven more rapidly because less heat must be extracted from the overlying porous medium. Consequently, the layer growth rate increases dramatically as Teh increases. If the initial temperature were equal to the equilibrium temperature (i.e. Teh — 1)) the entire region would be within the stability field and A would would only be limited by the need to remove latent heat. At low values of Teh the layer evolves slowly enough that despite the high nominal value of the Lewis number a, gas has time to diffuse towards the interface and bring about an increase in the hydrate volume fraction. Chapter 3. An Analytical Model of Hydrate Formation 39 0-2C6 - \_ n o m i n a l v a l u e \ v . A'^h = 4 .6 — O X 5 - — 02C4 " A 0.2C3 -— 02C2 -— O.201 J O2C0— 1 1 1 1 0 .595 0 .590 1 4 8 12 16 20 Figure 3.5: Dependence of h and A on Srh The axis on the left is for the layer growth rate parameter A, where the interface position is d(t) = 2\Vt. The axis on the right is for the hydrate volume fraction h expressed as a percentage. From (3.3), Srh = (ph<l>Lh)I'{PfC'/(Too — T0)}. O 0-2 04 _ 06 OS l O 2 i h Figure 3.6: Dependence of h and A on Teh The axis on the left is for the layer growth rate parameter A, where the interface position is d(i) = 2AV1 The axis on the right is for the hydrate volume fraction h expressed as a percentage. In (3.19), Teh is defined as Teh = (Teh - T0)/(Too — T 0). Chapter 3. An Analytical Model of Hydrate Formation 40 The problem has been solved with the assumption that the interface is planar and the geometry is completely one dimensional. In nature their are always small inhomogeneities which can cause perturbations in the planar interface. If these perturbations grow with time then a planar interface is unlikely to be realized in nature. To test how the solution would react to interface perturbations, a stability analysis was performed and is outlined in Appendix C. The results show that the solution is stable to short wavelength pertur-bations and only becomes unstable at long wavelengths. The fastest growing instabilities are those with the longest wavelengths; since the planar solution may be regarded as the infinite wavelength solution, the theoretical predictions presented here are expected to be realized in practice. To test whether this is indeed the case, an experimental apparatus was designed to simulate the conditions illustrated in Figure 3.1. This apparatus and the experimental method are described in the following chapter. Chapter 4 Experimental Work 4.1 Experimental Apparatus In order to test the hydrate formation model presented in Chapter 3, an apparatus was constructed to simulate the conditions considered in the theoretical analysis. Although this study is mainly motivated by an interest in naturally occurring methane hydrate, the analysis was carried out in a general manner, and should be relevent to the formation of any other type of clathrate hydrate. To avoid the logistical difficulties associated with the high pressures required to stabilize methane hydrate, a substitute hydrate former was sought for use in the experiments. Propane is favourable in this regard since propane hydrate is stable at approximately 30 PSI (0.21 MPa1) at temperatures slightly above 0°C. However, the low solubility of propane in water results in a modelled hydrate volume fraction that is well below 1%; thus to facilitate experimental detection of the hydrate phase, a guest molecule which dissolves more readily in water is preferred. Carbon dioxide was chosen because of its high solubility and relatively low hydrate stability pressure above the freezing point of water. C0 2 has the further advantage that it is not flammable and hence is less hazardous to work with than hydrocarbons. In designing the apparatus in which experiments are performed, a number of criteria had to be met. The CO2 hydrate phase diagram, shown in Figure 4.1 [from Bozzo et al, 1975], indicates an equihbrium pressure of 151 PSI (1.04 MPa) at —1.5 °C. This is the quadruple point with ice, so to keep ice from forming, the operating temperature is 114.7 PSI = 0.1013 MPa=l atm. 41 Chapter 4. Experimental Work 42 400 Pressure (PSI) 3 5 0 Pressure 2.5 (MPa) 1.0 Temperature (°C) Figure 4.1: CO2 Hydrate Phase Diagram This portion of the CO2 — H20 phase diagram is based on a best fit equation to experimental data in the temperature range —1.5 °C to +9 °C [Bozzo et al., 1975]. The reported average deviations from this line for experimental equilibrium pressures is 0.24%. The equilibrium line is based on CO2 concentration levels in excess of its solubility. The invariant point with ice is at -1.5 °C, and 151 PSI (1.04 MPa). constrained to be greater than —1.5 °C; thus the minimum operating pressure is 151 PSI (1.04 MPa). Since the theory is derived for a one dimensional half-space, the reaction chamber in which the experiments are performed has to be sufficiently large to keep this approximation valid; if the sediment column were too short, the temperature at "infinity" would change significantly over the duration of the experiment. The chamber also has to be wide enough that the effects of lateral heat flux are not appreciable. The temperature of the chilled boundary has to be maintained at a constant value for several hours despite variations in'the ambient temperature. In addition, although visual observations of the sediment column yield no quantitative information, they are helpful in some aspects of apparatus testing and may provide insight in future experimental work. With these concerns in mind, the reaction chamber was constructed from a 9.5 mm (3/8 inch) thick, clear acrylic tube, 0.7 m long with inner diameter 0.14 m, capped at each Chapter 4. Experimental Work 43 end by an aluminum lid in which coolant is circulated to maintain a boundary tempera-ture (see Figure 4.2). A pressure gauge situated on top of the container has shown that the vessel can be subjected to pressures of 300 PSI (2.07 MPa). The characteristic ther-mal diffusion time for 0.7 m of water saturated sand (from (2.40) with « e = 10 - 6 m2/s) is in excess of 5 days, so end effects should be unimportant for these experiments, which are typically completed within a 24 hour period. The acrylic tube is a poor thermal conductor and it is insulated with foam to further reduce lateral heat loss. Coolant is circulated with a Haake D8GH refrigeration unit capable of maintaining the temperature of its internal bath to within a few tenths of a degree Celsius. The hoses in which the coolant is transported, the coolant reservoirs, and the 6 aluminum bars which hold the chamber together are also insulated, to isolate the system from variations in the ambient temperature. To detect the position of the phase change interface, temperature is measured using RTDs (resistive thermal devices) which utilize the temperature dependence of the elec-trical resistance of platinum. The RTDs are held in place by aluminum sheaths which are glued to a plastic rod in the centre of the cylinder. The metal sheaths serve to in-crease the surface area over which ohmic heat is dissipated, thus minimizing the effects of temperature measurements on the in situ conditions. The leads, which are coated with a water resistant toluene-based compound, are ushered out through the base of the cylinder to some amphfying circuits and on to an A/D card connected to a 386 personal computer. The temperature at each RTD position is found by comparing the resistance across the leads once they are brought outside the chamber to that of a nominal precision resistor at ambient temperature. Temperature data is dumped to a file which is later accessed for analysis and plotting. A more complete description of the data acquisition system is found in Appendix F. Versatility was also a major consideration in the apparatus design. While the desired Chapter 4. Experimental Work 44 Reaction Chamber pressure gauge Mixing Chamber to gas cylinder agueous solution coolant reservoir aluminum bars aqueous solution cast acrylic tube variable speed motor coolant Lane Mtn. sand coolant reservoir 1 toA/D card Figure 4.2: Experimental Apparatus The reaction chamber is used in the CO2 hydrate experiments and the water freezing experiments. Pressure is maintained by a regulated gas cylinder which is normally con-nected to the flexible hose on the top of the chamber while an external valve is used to block the hose exiting its base. The mixing chamber is designed for experiments with low solubility gases and is not used in this study. Chapter 4. Experimental Work 45 conditions for this particular problem are relatively easy to simulate, in the future it would be convenient to use the same apparatus to observe hydrate formation characteristics in other experimental environments. For instance, it would be interesting to monitor the development of a hydrate layer when gas saturated fluid is forced to flow through the system at a constant rate, perhaps with a constant temperature imposed on both the top and bottom boundaries of the confined sediment column. Experiments could be performed with different hydrate formers which have low solubilities and need to be mixed into the water in a separate container. Although these types of studies are beyond the scope of the present work, the apparatus is flexible enough to be adapted to such uses. Depending on the valve configuration, fluid may enter or exit the reaction chamber through flexible hoses which attach to connectors that pass through both the top and bottom reservoirs and caps. If necessary, coolant from the refrigeration unit may be pumped to both ends of the reaction chamber and to the copper cooling coil of a separate mixing chamber. The mixing chamber is equipped with a small motor attached to a propeller that circulates the liquid to aid in the dissolution of low solubility gases. This mixing chamber was not used for the CO2 experiments. 4.2 Formation of an Ice Layer Before performing hydrate formation experiments, an investigation was made into a re-lated problem, the formation of an ice layer in a sediment column. An examination of this somewhat simpler process provides an opportunity to test the apparatus and re-fine the experimental method before introducing the additional complications of elevated pressures and the presence of gas. Moreover, these water freezing experiments enable an assessment of whether factors not considered in the theoretical analysis, such as pore scale phenomena, play a significant role in this related problem. The development of Chapter 4. Experimental Work 46 either an ice layer or a hydrate layer should obey the same conservation laws for momen-tum, mass, and energy. The freezing of water in sediments may be thought of as an end member of the hydrate formation problem with an ice pore volume fraction near 1. 4.2.1 ' Experimental Method To begin the water freezing experiments, water is added to the sand filled reaction cham-ber and allowed to sit for several hours to enable the bulk mixture to equilibrate at the ambient temperature Too. The refrigeration unit is then turned on to bring its internal coolant bath below the desired reaction chamber base temperature To. This ensures that once the coolant circulating pump is turned on, the time required for the base temper-ature to reach a constant value of To is minimized. A data acquisition program on the personal computer is used to collect temperature measurements at intervals as short as 3 s from the four RTDs which are located within the bottom 0.1 m of the container. Data is normally collected until all four probes are within the frozen layer, at which time the refrigeration unit is turned off and the reaction chamber is allowed to warm to ambient temperature once more. By repeating this process several times using different base temperatures, a database was acquired for comparison with theoretical predictions. 4.2.2 Results The data from a typical experiment is shown in Figure 4.3. Each line displays the temperature record at the distance above the cold base plate indicated in the bottom right-hand corner of the graph. When the temperature falls below the freezing point, ice formation and the associated release of latent heat result in an abrupt change in the slope of the temperature vs. time curve. The precise time at which freezing occurs is determined from a plot of the slope of the temperature record, which is found using the five point derivative filter described by Abramowitz and Stegun [1972, p. 883]. Chapter 4. Experimental Work 47 2 0 1 — 1 — 1 — 1 — i — 1 — 1 r i 1 1 1 r 0 200 400 600 Time (minutes) Figure 4.3: Temperature Measurements During Ice Formation Temperature measurements at the four RTD probes during an ice formation experiment. The initial temperature was uniform at 19 °C, and the base temperature T 0 was -7 °C. Probes were located at 20, 40, 60, and 80 mm above the base. Data from this experiment is labelled "run 4" in Tables 4.1 and 4.3 and in Figure 4.4. Chapter 4. Experimental Work 48 Table 4.1 summarizes the results of five experimental runs with different initial tem-peratures Too and base temperatures To. Problems with the cooling system were encoun-tered during runs 1 and 5, causing those experiments to be aborted before the frozen layer encompassed all four probes. The temperature probes were calibrated at 0 ° C , and slight differences in the electrical resistance of the circuitry connected to each probe are responsible for errors in the initial temperature Too. More confidence is given to the probes' ability to detect variations in the rate of temperature change than to the absolute accuracy of the temperature measurements themselves. The base temperature was cal-culated by assuming a linear gradient in the frozen layer several hours after the freezing front had passed all of the probes. Inaccuracies in the probe positions arise from the finite width of their metal casings, which tend to average the temperature measurements over their surface area. The interface crossing time becomes more difficult to pinpoint as the thermal gradient becomes shallower at later times. As mentioned previously, the dynamics of ice layer formation and hydrate layer forma-tion are governed by the same conservation equations for momentum, mass, and energy. A modified form of the interface condition from energy conservation, (3.35), with the hydrate pore volume fraction h replaced by the ice pore volume fraction i , enables pre-dictions of the interface velocity parameter A. The frozen volume fraction of the pore space i should be close to 1, so (3.35) becomes /To F n A 2 l-fei jTj r . A 2 Ke(i)fei __ VTT6 erfc(v^A) V*6 Ke rf (y^A) ~ T i ' 1 J where ne(i) is the thermal diffusivity of the icy sediments, and Tj has a definition similar to that of Th in equation (3.11), namely , r . _ PiCi<(> + P.C,(1 -<j>) Ke PfCf Ke(i) ' Chapter 4. Experimental Work 49 Run Initial Base Interface Freezing 104 a /(2y/i) Number Temperature Temperature Position Time -1/2) - ^ (°C) To (°C) a (mm) t (min.) (m s 1 19.0 ± 0.5 -4.3 ± 0.2 20 ± 2 70 ± 10 1.54 (1.30 - 1.83) 40 ± 2 280 ± 15 1.54 1.43 - 1.67) 2 13.8 ± 0.3 -4.4 ± 0.1 20 ± 2 95 ± 10 1.32 (1.13 - 1.54) 40 ± 2 270 ± 10 1.57 (1.47 - 1.68) 60 ± 2 550 ± 50 1.65 (1.53 - 1.86) 80 ± 2 1075 ± 50 1.58 (1.50 - 1.65) 3 18.5 ± 0.5 -5.6 ± 0.1 20 ± 2 67 ± 5 1.58 (1.37- 1.80) 40 ± 2 205 ± 10 1.80 (1.67- 1.94) 60 ± 2 460 ± 15 1.81 (1.72- 1.90) 80 ± 2 840 ± 15 1.78 (1.72 - 1.85) 4 19.0 ±0.5 -7.0 ± 0.2 20 ± 2 40 ± 5 2.04 (1.73 - 2.4) 40 ± 2 150 ± 10 2.11 (1.94- 2.29) 60 ± 2 330 ± 15 2.13 (2.02 - 2.25) 80 ± 2 545 ± 20 2.21 (2.12 - 2.31) 5 14.9 ± 0.2 -7.25 ± 0.1 20 ± 2 40 ± 5 2.04 (1.73 - 2.40) 35 ± 2 95 ± 10 2.32 (2.08 - 2.60) 50 ± 2 190 ± 15 2.34 (2.16 - 2.54) Table 4.1: Ice Layer Data The freezing times and probe positions for five experimental runs are presented here. The final column gives a/(2y/i) for comparison with the modelled value of Ay^. The brackets contain the range of values for a/(2y/i) using the errors on a and t. Chapter 4. Experimental Work 50 The Stefan number in (4.1) is defined by STi = PtCtfL\ To) ' (4-3) and the nondimensional ice equihbrium temperature is Tei = ^ ^ - . (4.4) -t oo — J-0 In Table 4.2 the physical properties of the materials used in the experiments are listed. The thermal diffusivities « e and Ke(i) were determined experimentally using the methods outlined in Appendix D. Table 4.3 summarizes the predictions of equation (4.1) for the growth-rate parameter A in each experimental run. The ice equihbrium temperature Tei was found to be 0.0 ± 0.1 °C, so the effects of pore scale phenomena on the water freezing point appear to be negligible. The theoretical results in Table 4.3 may be compared with the observational results in Table 4.1. The theoretical interface position is given by a = 2A\/M , (4.5) so that experimental data for a and t may be used to determine the factor v^ =^ - («) The advantage of this representation involving y/Kl is that the term on the right of (4.6) depends solely on the experimental observations and not on the material properties used in the model. Separate estimates of may be obtained from each of the four probes, which are located at different positions a. These experimental results are listed in the final column of Table 4.1 and may be compared directly with the theoretical results in the final column of Table 4.3. A chart displaying both sets of results is displayed in Figure 4.4. Chapter 4. Experimental Work 51 Property Nominal Value Units Pf 1000 kg m-3 Ps 2650 kg m"3 api 916 kg m-3 1110 kg m"3 <C} 4200 J kg"1 °C J kg"1 °C J kg"1 °C dca 2200 eCi 2100 fch 2080 J kg"1 °C d4> 0.4 au 333.6 x 103 Jkg" 1 448 x 10s. Jkg" 1 m2 s"1 (6.5 ± 0.6) x lO" 7 dKe(i) (9 ± 3) x 10-7 m2 s"1 To Ti 1.23 ± 0.05 (1.02 ± 0.05)/cE/Ke(i) 0.298 JCe 0.01 0.03 Table 4.2: Properties used for Modelling Experiments Sources of Information: aice at 0 °C, [Hobbs, 1974, p. 348]; 6 [Bozzo et al., 1975]; cwater at 5 °C [Lide, 1990] (the presence of dissolved gas affects the fluid heat capacity, but only by about 1%, hence this effect is neglected); ^measured in the laboratory; eat 0 °C [Hobbs, 1974, p. 361]; ^ assume same as methane hydrate [Sloan, 1990, p. 387]; flat 0 °C [Hobbs, 1974, p. 362]; hC02 hydrate at 3.9 °C [Takenouchi and Kennedy, 1965]; 'hydrate composition C0 2 • 5.75H20; a^ssume equilibrium composition 1/3 of solubility after Claypool and Kaplan, [1974]; Solubility of C02 at 250 PSI (1.72 MPa), 15 °C, [Stewart and Munjal, 1970]. Chapter 4. Experimental Work 52 Run Sri 104 X^Te (m s-1/2) 1 0.18 ± 0.01 1.25 ± 0.03 1.23 ± 0.40 2 0.24 ± 0.01 1.59 ± 0.02 1.48 ± 0.43 3 0.23 ± 0.01 1.21 ± 0.03 1.52 ± 0.46 4 0.27 ± 0.01 1.12 ± 0.03 1.76 ± 0.51 5 0.33 ± 0.01 1.31 ± 0.02 2.01 ± 0.53 Table 4.3: Ice Layer Predictions The nondimensional parameters, tei and Sri, in the second and third columns are for the data shown in Table 4.1. The values of Xy/Ke' are found by solving (4.1) using Tei, Sxi, and the physical properties listed in Table 4.2. 1 2 3 4 5 Run Number Figure 4.4: A Comparison of Ice Interface Motion With Predictions This chart displays model predictions to experimental data for the five ice formation runs summarized in Tables 4.1 and 4.3. The experimental data for each probe position is shown with light grey error bars while the model results are displayed with dark grey error bars. Chapter 4. Experimental Work 53 In each experimental run a/(2y/i) varies from probe to probe. Some of this variation can be explained by inaccuracies in the probe positions and in the timing of the interface arrival. The experimental data generally shows faster layer growth rates than the model predictions using the nominal parameter values, however in each case the model results match the observations within the error bars. Much of this offset could be explained by errors in the thermal diffusivity of the icy sediment «e(*)> though uncertainties in the other model parameters, such as Sri and Tei, may also play a role. Another possibility is that some of the liquid in the pore space might remain unfrozen so that the latent heat released is less than otherwise expected. An attempt to quantify the ice pore volume fraction using the melting rate is outlined in Appendix E. The results were inconclusive due to large error sources, however, since the Stefan numbers for these experiments are all relatively low (see Table 4.3), the presence of a small unfrozen volume fraction could only produce a modest effect on the results. The freezing time at the bottom probe's location is quite sensitive to the time required to bring the base temperature to To. In practice, the base temperature is not imposed as a step function in time (as is assumed in the model), but is applied over a finite length of time. The delay in achieving the final base temperature likely responsible for the comparatively lower growth rates inferred from its interface crossing times. Referring to Table 4.1, the initial temperatures for runs 2 and 5 (13.8 °C, and 14.9 °C respectively) were below room temperature. In both of these cases, the sediment column was not given time to reach a uniform temperature before the experiments were started. The warmer upper regions of the sediment column should have resulted in a slower layer growth rate than that which would have been realized if were uniform. The modelled growth-rates for runs 2 and 3 are almost identical (X^/K^ = (1.48 ± 0.43) x 10 - 4 m/s/s and (1.52 ± 0.46) x 10 - 4 m/y/s, respectively). The fact that the experimental data for run 2 shows the frozen layer to have grown significantly slower than in run 3 is easily Chapter 4. Experimental Work 54 explained by additional heat flux from above due to the nonuniform initial temperature in run 2. For run 5, the experimental data indicates a layer growth rate which is faster than the model prediction with the nominal parameter values to approximately the same degree as those runs with a near room temperature value of Too (runs 1, 3, and 4). The short duration of this experiment, only about 3 hours compared to about 18 hours in run 2, likely allowed the layer to grow with minimal influence from the effects of the nonuniformity of Too. This set of experiments has demonstrated that equation (4.1) adequately describes the formation of an ice layer in initially fluid saturated sediments. With the nominal parameter values, the model tends to underestimate the layer growth-rate, however, much of the difference can be ascribed to the value used for the thermal diffusivity of the ice layer Ke(i). This error source is not important for hydrate formation because the thermal diffusivity of the hydrate layer K e ( h ) is expected to be very close to the thermal diffusivity of the fluid saturated sediments ne since the volume fraction of hydrate is expected to be small, and the thermal conductivity of hydrate is close to the thermal conductivity of water (Ke(h) oc K(h)). The relatively low growth-rates derived from the data of the first probe highlight the fact that the base temperature is not dropped from to To instantaneously at the start of the experiment as the analysis assumes, but that a finite amount of time elapses during which the apparatus components cool down. This suggests that some aspects of the apparatus design could be improved upon, but if these limitations are kept in mind there is no reason not to proceed to the focal point of this study, the formation of hydrates. Chapter 4. Experimental Work .55 4.3 Formation of a Hydrate Layer 4.3.1 Hydrate Detection Before beginning the hydrate formation experiments, it is expedient to first perform some simple calculations to estimate the magnitude of the thermal signature that the RTD probes are required to detect. Equation (3.48) for the hydrate volume fraction h = -r%-PhCh approximates the hydrate volume fraction when there is negligible gas transport. As such, (3.48) yields the minimum hydrate volume fraction that should by present for a given value of cn, which was defined in (3.4) as Ch — ce -Ce ' where Ch is the gas concentration in the solid hydrate, ce is the equilibrium gas con-centration, and CQO is the initial gas concentration in the sediment column. Using the parameter values listed in Table 4.2, Ch is found to be approximately 14.4. Employing this result and using the densities of C0 2 hydrate and water given in Table 4.2, (3.48) returns h f« 0.06; thus the equipment must be able to detect hydrate occupying as little as 6% of the pore space. During experiments, the reaction chamber temperature is monitored so the obvious question is: "What effect will hydrate formation have on the thermal profile?". The model developed in Chapter 3 predicts that the temperature during hydrate formation will be governed by (3.31) and (3.32), namely ^ ' " Ket 2 T(z,t) =-(roo - T0)Teh ) ( + To ,0 < z < 2XVKet (4.7) erf (vTfcAJ Chapter 4. Experimental Work 56 and T(z,t) = (Too-T0) e ( fTo~z\ Te fc -er f (VF0\) I ffo~z | +T0 ,«' > 2X^/^el (4.8) where the similarity variable has been replaced using its definition, £ = z/(2y//cei), and the temperature is given in dimensional form. With an approximate hydrate volume fraction of 6%, A may be estimated from the interface condition from energy conservation, (3.35), which becomes The values of Ten and Srh depend on the inital temperature Too, the base temperature To, and the equihbrium hydrate temperature Teh. The initial temperature for most experiments is the ambient temperature of the laboratory, so T^ « 20 °C. To avoid forming ice, the base temperature is constrained by the freezing point of water, thus To ~ 0 °C. The cast acrylic tube is rated for a maximum working pressure of 350 PSI (2.4 MPa), so allowing a safety margin of 100 PSI (0.69 MPa) gives a gauge pressure of 250 PSI (1.72 MPa), for a total pressure (total pressure = gauge pressure + atmospheric pressure) of about 265 PSI (1.83 MPa). The phase diagram in Figure 4.1 indicates a corresponding equihbrium temperature of Teh ~ 3.4 °C. Therefore with the parameter values compiled in Table 4.2, equation (3.19) gives Ten ~ 0.17, and equation (2.47) gives STH ~ 1.91. Assuming that the small volume of hydrate has a negligible effect on the thermal properties of the bulk material so that Ke(h) 10"7, and ~ To ~ 1.23, results in an estimated growth rate parameter from (4.9) of A «' 0.13 (A-y/Kg" « 1.1 x 10 - 4 m/y/s). Figure 4.5 displays the corresponding modelled temperature records for probes at 20, 35, 50, and 65 mm above the base of the reaction chamber. Figure 4.5 demonstrates that for the conditions outlined above, the hydrate phase would produce only a very small thermal signal which would be difficult to detect using Chapter 4. Experimental Work 57 T _i ,—, , , •_ > 5 0 0 , 1 0 0 0 , 1500 2000 Time (mm.) Figure 4.5: Predicted Temperature Records for Hydrate Formation Above 0 °C This graph shows the predicted temperature records at 20, 35, 50, and 65 mm above the base of the reaction chamber during a hydrate formation experiment when To = 0 °C, Too = 20 °C, Teh — 3.4 °C, and Dc + D& « Ke. The equilibrium hydrate temperature is indicated with the horizontal bar on the right. The hydrate volume fraction is 6%. temperature measurements alone. The temperature record at each of the four positions shows no significant deviation from a diffusive cooling pattern as the hydrate interface passes. A protracted time period is required before the hydrate stability field encompasses all four locations; so the rate of temperature change is relatively low when the interface passes. During the ice formation experiments it was observed that the interface was more easily detected when the rate of temperature change was high. In Appendix G predictions are made for the hydrate formation rate when the base temperature is brought below the ice formation temperature. It was found that the hydrate phase should still produce only a small thermal signal. An alternative method of hydrate detection could be based on the pressure dependence of its phase equilibrium. If hydrate is initially present in the pore space, and the pressure is suddenly decreased to beneath the hydrate equilibrium pressure, the hydrate should Chapter 4. Experimental Work 58 dissociate. In order to accomplish this, it must consume latent heat, which should cause a reduction in the temperature of the surrounding material. The change in enthalpy caused by the dissociation of hydrate is AHdis = -h<f>PhLh , (4.10) and, the resulting change in temperature is given by \rpl- AHdis AH (All) A i d l S ~ PfCf<t> + p.C.(l -<j>)- ToPfCf ' l 4 ' i l J Substituting (4.10) into (4.11) gives < which for h « 0.06 gives ATdu « —1.9 °C. Such a large temperature change should be easily detected by the RTD probes. 4.3.2 Experimental Method Some care was exercised to obtain a high gas concentration in the fluid before conducting the hydrate experiments. To maximize the contact of the gas with the pore liquid, the pore space was filled with COz prior to the addition of water. This was accomplished by connecting the regulated C O2 cylinder to the top of the reaction chamber and applying a modest pressure to blow the fluid initially present out the bottom (refer to the diagram of the apparatus in Figure 4.2). Water was injected through the base of the reaction chamber so that some CO2 could dissolve as the chamber was filled; the remaining gas was either left trapped as bubbles in the pores or pushed out through the top. Once the liquid was added, the exhaust valve was closed and C02 gas was percolated up the sand column until the desired pressure was achieved. Gas-liquid contact was maximized by introducing the gas from below, to enable more of it to dissolve before it reached the Chapter 4. Experimental Work 59 top of the vessel. In addition, some of the gas should have been trapped in the pore space due to surface tension effects, and could thus provide an opportunity for the gas concentration in the fluid to increase as this trapped gas dissolved with the application of increased pressure and decreased temperature, and with the passage of time. Once gas bubbles were no longer seen rising in the fluid above the sediment column, the valve at the base of the chamber was closed, and the gas cylinder was reconnected to the top of the chamber to maintain the pressure at its elevated state while the CO2 dissolved. Because the solubility of gas increases at lower temperatures, coolant was cir-culated through the bottom reservoir at a few degrees above zero to aid in the dissolution of gas remaining in the pore space. After leaving the chamber in this state for several days, the refrigeration unit was shut off and the chamber allowed to warm to ambient temperature. With the temperature and gas composition uniform throughout the reaction chamber at Too and Coo, the remaining experimental procedure is much the same as for the ice freezing experiments. The temperature at the base was brought to a constant value of To, below the equilibrium hydrate temperature Te/,. Data was collected using the four temperature probes, located at 20, 35, 50, and 65 mm above the reaction chamber base. The thermal conditions were monitored until all four probes returned temperatures several degrees below the expected value of Te^ . 4.3.3 Analysis of Results Figure 4.6 displays the results of an experiment for which the base temperature T 0 ~ -7.6 °C. The gauge pressure was 225 PSI (1.55 MPa) for a total pressure of 240 PSI (1.65 MPa), which corresponds to an equilibrium hydrate temperature Teh of about 2.5 °C. By comparison, ice forms from the water-gas mixture at Tej « —1.5 °C, thus both ice and hydrate should have formed in this particular experiment. The ice interface crossing Chapter 4. Experimental Work 60 200 300 l i m e (min.) Figure 4.6: Temperature Records From a Hydrate Formation Experiment The temperature records from the bottom three probes are shown here for a hydrate formation experiment in which the base temperature To ~ —7.6 °C is less than the ice formation temperature Tei ~ —1.5 °C, which is in turn less than the hydrate equilibrium temperature Teh ?a 2.5 °C. An equipment problem corrupted the data from the fourth probe. times are clearly visible on the temperature record, at approximately 55, 130, and 255 min.. The hydrate interface crossing times, however, are not discernable. Repeated attempts to track the motion of the hydrate phase change interface failed, since the temperature measurements during these simple cooling experiments did not reveal the presence of hydrate. The initial experimental goal of testing the model pre-dictions of the layer growth rate was set aside and an effort was mounted to determine whether in fact any hydrate was indeed present in the pore space. To this end, ex-periments were performed in which the pressure in the reaction chamber was suddenly reduced to put the hydrate out of equilibrium and cause it to dissociate. The data in Figure 4.7 was collected when the pressure was dropped to the atmospheric level several hours after the data displayed in Figure 4.6 was collected. The particular record shown is from the second lowest temperature probe, the records from the other probes are similar Chapter 4. Experimental Work 61 T ( ° C ) T i m e (min . ) Figure 4.7: Temperature Record From a Depressurization Experiment The temperature data from the second lowest temperature probe is shown here for a depressurization experiment in which the gauge pressure was dropped from approximately 225 PSI (1.55 MPa) to 0 PSI (0 MPa) over a period of several tens of seconds. Any hydrate present should have dissociated at the reduced pressure. in form and amplitude. Since the temperature record for the depressurization experiment depicted in Figure 4.7 is for a probe that is located within the icy layer, the hydrate should dissociate into gas and ice rather than gas and aqueous solution. Thus the expected temperature change is given by a modified form of (4.12), namely (h + h*)4>Keph(Lh - Li) ATdis = (4.13) TiKe(i)pfCf where the latent heat required for dissociation is the difference between the heat of formation of hydrate from aqueous solution and gas, and the heat of formation of ice from aqueous solution, and the thermal diffusivity ratio Ke//ce(i) appears because of the definition of Ti in (4.2). The parameter h* w 0.028 is the additional hydrate volume fraction expected to form due to the exclusion of gas from the ice structure, as explained in Appendix G. Substituting the appropriate parameter values from Table 4.2 into (4.13) Chapter 4. Experimental Work 62 gives ATdu « —0.72 °C. The data in Figure 4.7 shows that the region surrounding the probe initially cools and reaches a minimum approximately eight minutes later of about 0.02 °C below its inital value, then warms to reach a peak of about 0.1 °C above its initial value some thirty minutes later, before beginning to cool again. Obviously this is not the behavior that was expected, so some of the essential physics of the problem must have been neglected. Aside from hydrate dissociation, other sources of temperature change that might be associated with a rapid pressure release in the reaction chamber include the exsolution of gas; the decompression of ice, aqueous solution, sand grains, or gas; and the advective mixing that accompanies the exsolution of gas. The temperature change associated with the adiabatic compression of a material is given by [Anderson, 1989, p. 90] ATcomp = a r ^ f P , (4.14) where ay is the coefficient of thermal expansion, AP is the pressure change, p is the density, and Cp is the specific heat. For ice that completely fills the pore space ATcompi becomes gTiTAP4>Ke A l c o m t n - T i K e { { ) p f C f • (4.15) For a pressure drop of 225 PSI (1.55 MPa), at a temperature of —5 °C (268 K), and with a thermal expansion coefficient of an w 60 x 10 - 6 K - 1 [Hobbs, 1974, p. 347], this gives ATcompi « —0.002 °C. If the pore fining fluid were water at 0 °C, ay/ ~ —68 x 10 - 6 K - 1 « —art [Weast, 1975], and the temperature would be expected to increase by a few thousandths of a degree. The enthalpy change associated with gas coming out of solution is calculated as _ AHex,Acpf<l> A r " - T0PfC, • ( 4 , 1 6 ) To get the largest possible effect, all the gas would have to come out of solution from Chapter 4. Experimental Work 63 a gas-rich pore fluid that completely filled the pore space. This is not a feasible sce-nario considering that the temperature record displayed in Figure 4.7 shows that the temperature is below freezing, but nevertheless it useful to get an idea of the magni-tude of temperature change that exsolution could cause. Using AHd%a « —441 kJ/kg (—19.43 kJ/mole) [Sloan, 1990, p. 76], and assuming that virtually all of the gas is ex-solved from a saturated solution, so that Ac « —0.03, then ATexs « 1 °C. Clearly, if there were even a small amount of fluid in the vicinity of the temperature probe that unfrozen, then exsolution of some of it's dissolved gas could raise the temperature by a significant amount. It is possible that a combination of mechanisms might be responsible for the tem-perature record shown in Figure 4.7. Thus the temperature signature due to hydrate dissociation may be present, but it is obscured by other effects. It is also possible, however, that there is no hydrate present to dissociate. In an attempt to clarify which processes are more important in these depressurization experiments and to improve vi-sual observations, the sand was removed from the reaction chamber for another set of hydrate detection experiments. Without the presence of sand, and with the gas cylinder connected to the top of the chamber to maintain the pressure, there is no possibility of free gas in the bottom of the chamber. An ice layer was formed by cooling the base temperature to below the ice freezing temperature Te{. In these experiments the thermal profile in the reaction chamber is not diffusive since thermal convection occurs. This was not the case when the sand was present because the Rayleigh number was lower since the fluid had to maneuver around the sand grains (see Appendix C). The hydrate crystals are expected to bind themselves to the ice at the base of the reaction chamber. Since the volume of hydrate which forms from dissolved gas is generally small, only a fraction of the ice surface will be covered with hydrate. Thus the hydrate should be intermixed with ice Chapter 4. Experimental Work 64 crystals as the icy layer grows from the base of the reaction chamber. Once the frozen layer was formed, the pressure was released to dissociate the hydrate. Bubbles of CO2 formed in the fluid as the pressure dropped and the gas solubility was lowered correspondingly. In many experiments, once the gauge pressure was reduced from 250 PSI (1.72 MPa) to about 100 PSI (0.69 MPA), the exsolution of bubbles was so vigorous that visual observation of the ice surface was no longer possible. The turbulent motion caused the fluid to circulate so that warmer water from the upper regions of the reaction chamber was brought down to the ice surface and facilitated the rapid melting of the icy layer. The temperature change near the melting boundary was quickly detected at the RTD probe locations. No definitive evidence of hydrate dissociation was acquired, but the large temperature change associated with the circulation of fluid may have hidden the effects of the dissociation enthalpy change. One final experiment was performed in which the chamber was left for one week with the gauge pressure held constant at about 290 PSI (2.00 MPa) to allow the gas to dissolve. After this time span, the base temperature was dropped from approximately 0 °C to —2 °C to allow the growth of a thin icy layer. Approximately 28 hours later, the acrylic tube failed causing a rapid pressure decrease to the atmospheric level, and a discharge of the reaction chamber's fluid contents. The frozen layer at the bottom of the chamber remained intact, and upon close inspection was found to be releasing gas. The observed white, icy structure was interpreted to be a mixture of hydrate and ice. Occasionally as the gas was released, pieces of the bulk structure were discharged several centimeters into the air. This was taken as evidence that as the gas component was freed from its water cages, it was sometimes trapped by the ice surrounding it and pressure was able to build up before it could be exhausted to the atmosphere. The substance continued to dissociate for at least ten minutes, after which observations were discontinued. In the following chapter the implications of the hydrate formation model and the Chapter 4. Experimental Work experimental results are discussed. Chapter 5 Discussion One of the questions posed in the introduction was whether or not free gas is required in order for hydrate to form. The experimental results presented here have shown that hydrate does form from an aqueous CO2 solution without the presence of free gas when the temperature is below the ice formation temperature. No clear evidence was found, however, for hydrate formation from solutions above the solution freezing temperature. This is not to say that the hydrate did not form, merely that if it did, it escaped detection. There are several plausible explanations that promote each of these possibilities. Several authors [eg. Sloan, 1990] have reported that kinetic effects can play a large role in the formation of hydrates from water and free gas. If the phase change is rate-limited by kinetic processes it is conceivable that once the equilibrium temperature was reached, an insufficient length of time was allowed for a significant volume of hydrate to crystallize. In addition, if phase change kinetics were to play a dominant role, the moving interface might be characterized by a gradational increase in the hydrate concentration, which would make it difficult to detect. Cha et al. [1988] observed that the formation rate of methane hydrate was greatly enhanced by the presence of a bentonite surface. They argue that the solid surface provides a nucleation sight to which the hydrate crystal can bind during its initial formation stages. In this study, Lane Mountain sand, which consists almost exclusively of quartz grains, was used as the sediment component, and it is possible that hydrate formation might have been facilitated by the presence of more clay. ) 66 Chapter 5. Discussion 67 The concentration of CO2 in equilibrium with its hydrate is not well constrained. Claypool and Kaplan [1974] report that for the methane-water system, the equihbrium concentration is about one third of the methane solubility limit immediately outside the hydrate stability field. The abruptness of the reported change in solution gas content once the methane hydrate stability field is reached suggests that it is reasonable to approximate ce as a constant. It was assumed in this study that the CO2 — H2O system has a similar behavior, however if in fact the equihbrium concentration were closer to the initial concentration in the experiments, the expected hydrate volume fraction would be correspondingly reduced. If future experiments are successful in determining the amount of hydrate that forms under different conditions, this data may be used to refine the estimates of ce. Before beginning each hydrate formation experiment, an attempt was made to ensure that a high concentration of gas was dissolved in the liquid. Despite these efforts, some of the experiments may have been performed with initial gas concentrations well below the solubility level. It is possible that in some cases the initial gas concentration was also below the equihbrium gas concentration, which would leave no opportunity for hydrate to form above the melting point of ice. To determine the initial gas concentration in future experiments, the pressure could initially be reduced slightly until gas bubbles are seen rising from the sediment column at which point the gas concentration may be taken from the solubility tables. The slow dissociation rate of the hydrate and ice sample which was removed from the reaction chamber may provide an explanation for the failure of the temperature probes to detect hydrate during earlier depressurization experiments. The observations suggest that the latent heat consumed during dissociation is removed from the hydrate's sur-roundings over a protracted time period. This would significantly reduce the magnitude of the associated thermal signature and allow it to be masked by the effects of other Chapter 5. Discussion 68 processes. Future depressurization experiments might be more successful at detecting hydrate if the pressure is gradually decreased to just under the stability pressure so that the mixing of pore fluids and the heat generated by the exsolution of gas would be minimized. An alternative method for hydrate detection using its dissociation characteristics could be based on the temperature dependence of its phase equilibria. In a series of ice melting experiments outlined in Appendix E, it was found that when the base tempera-ture was increased to above 0 °C, the temperature measured by each of the probes was constant at the ice melting temperature for an extended time period. In a similar set of experiments aimed at dissociating a hydrate layer, the probe measurements might remain constant at the hydrate equihbrium temperature for a significant length of time. In this type of experiment, the base-temperature and the equilibrium temperature may be used to estimate the heat flux that balances the latent heat; this could enable quantitative estimates of the hydrate volume fraction. The error bars on Figure 4.4, which compares the experimentally determined ice layer formation rate to the model results, highlight the need for more accuracy both in the values of the model parameters and in the experimental measurements. The most significant source of experimental error is the finite width of the aluminum sheaths which surround the temperature probes. These casings tend to average the measurements over their surface area and make precise interpretations of the interface crossing times ambiguous. Their purpose is to provide a sink for the approximately 1.4 mW of ohmic heat flux generated by the current that is driven through the RTD probes. They also serve to hold the probes in place and protect them when sand is dumped into the reaction chamber. A better method of holding the probes in place could be devised, and the probes themselves could be coated with a thermally conductive metallic paint so that only their width is an error source in interpreting the origin of the temperature data. The errors Chapter 5. Discussion 69 on the experimentally determined thermal diffusivities would also be greatly reduced by more accuracy in the probe locations. If hydrate is detected in future dissociation experiments, confirmation will be provided that hydrate formation from solution is feasible without the presence of free gas. In itself this is important since there are doubts as to whether there is sufficient organic material in the marine sediment column to produce concentrations of methane that exceed the solubility limit. However, for.the purposes of testing the model predictions, the most that can be hoped for from dissociation of experiments is a quantitative estimate of the volume fraction of hydrate present in the developing layer. Temperature measurements alone have thus far proved unsuccessful in tracking the progression of the moving phase change interface, thus further instrumentation of the reaction chamber's interior may be necessary to determine the layer growth-rate. Techniques which make use of changes in the bulk electrical properties or acoustic properties associated with the formation of hydrate are two promising alternatives. The success of such measurements would be greatly influenced by the spatial distribution of the hydrate crystals in the pore space. Acoustic velocities, for example, are controlled by the density and the elastic properties of the propagating medium. The density differences between hydrate and the pore fluid are minor, however the change in elastic properties could be significant if the hydrate formed at grain contacts and served to increase the rigidity of the bulk structure. Electrical properties are greatly influenced by the interaction between the surface charge on pore walls, and ions in the pore fluid. If hydrate were to coat a large portion of the pore surface area, its effect on the bulk resistivity could be substantial. Alternatively, if the pore throats were blocked by hydrate, the impediment to ion motion could produce a large electrical response. The model predictions might best be tested by combining the use of temperature measurements to determine the hydrate volume fraction, and electrical or acoustic techniques to monitor the interface position. Chapter 5. Discussion 70 The formation model presented here makes a number of simplifications to the com-plete set of equations. Scaling arguments were used to suggest that the thermal and gas concentration profiles are dominated by diffusive processes so that advection plays only a secondary role. Laboratory experiments revealed that the model is successful in describing the growth-rate of an ice layer. Contributions from neglected effects, such as those due to the volume change from water to ice, were found to be small. Similarly, the effect of the porous medium on the freezing temperature was found to be small. Even though hydrate was observed after one experiment, attempts to remotely detect the presence of hydrate using temperature probes failed. Further experimental efforts are therefore necessary before the validity of the model and its underlying assumptions may be quantitatively evaluated. At this point a review of the major features of the model is in order. The volume fraction of hydrate and the layer growth-rate are determined by the interplay between the need to remove the latent heat associated with the phase change and the efficiency with which gas can be supplied to the moving interface. Generally speaking, since heat transport is normally more efficient than compositional transport, the layer grows at a thermally controlled rate and converts the available gas into hydrate. The amount of gas available, and hence the hydrate volume fraction, is influenced by the ability of compositional diffusion to bring gas to the moving interface. Since hydrate production releases latent heat, which must be removed to allow the layer to grow, the hydrate volume fraction in turn affects the layer the growth-rate. Several nondimensional parameters that determine the relative importance of these processes have been identified. The Lewis number a is proportional to the ratio of the thermal and chemical diffusivities. Its high nominal value is reflected in the growth-rate of the hydrate layer, which occurs on the thermal diffusion timescale. The hydrate gas concentration parameter cn measures the amount of gas that hydrate crystals contain relative to the amount of gas in the pore Chapter 5. Discussion 71 fluid that is available to form hydrate. Since the solubility of gases in water is low in comparison to the gas concentration contained in hydrate, the hydrate gas concentration parameter must necessarily be much greater than one. With a high Lewis number, the layer develops more quickly than chemical diffusion can supply additional gas to its interface, thus the high value of cn imposes strict limits on the fraction of the pore volume that can be occupied by hydrate. The Stefan number Sxh compares the amount of latent heat that can be released by hydrate production to the amount of energy required to change the temperature. Despite the fact that the Stefan number can be quite high, the dynamics of layer growth are relatively insensitive to this parameter because of the low volume of hydrate that the gas solubility will allow. The specific heat ratio To is the heat capacity of the bulk mixture normalized by the heat capacity of water. Its value is related to the amount of energy required to change the temperature of the bulk mixture, so it acts as a scaling factor and the dependence of the layer velocity and hydrate content on To is nearly linear. The initial temperature, boundary temperature, and hydrate equilibrium.temperature together determine the value of Ten, the nondimensional equilibrium temperature. Its magnitude strongly influences the layer growth-rate, and at low values it can also have a large effect on the volume fraction of hydrate by retarding the motion of the interface and allowing chemical diffusion to supply additional gas for the phase change process. Equilibrium conditions are assumed to prevail on the interface of the developing hydrate layer, and the constant value of the equilibrium gas concentration implies a constant hydrate volume fraction throughout the layer. If the hydrate formation model presented here can be shown to accurately describe the growth of a hydrate layer under controlled laboratory conditions, an obvious question is: "What are the implications for hydrate formation in marine sediments?". Both the model environment and the marine sediment column are described by the same laws of conservation of momentum, mass, energy, and gas. Scaling arguments have shown Chapter 5. Discussion 72 that in both situations the fluid velocity is less important than thermal and chemical diffusion to the dynamics of heat and mass transport. Therefore the diffusion equation provides a leading order description of the thermal conditions and the gas concentration profile in both of these systems. Some of the physical conditions upon which the model is based are somewhat different from those which prevail beneath the sea floor. The initial state in the model is one with a constant temperature and gas concentration throughout. In continental shelf sediments, the thermal conditions are dominated by the geothermal gradient, with an essentially constant heat flux measured at the sea floor over extensive regions of hydrate-bearing sediments [Hyndman et al., 1993]. The distribution of gas in the sediments is often not well constrained, but it might be reasonable to suspect that it could be relatively uniform beneath the hydrate stability field. These differences are embodied within the initial conditions and boundary conditions for the two problems, however, and it must be emphasized that the physics involved in the two problems is the same. If the model presented here compares favourably with future laboratory simulations, the path will be cleared to apply the same equations to the more complicated problem of hydrate formation in marine sediments. Chapter 6 Conclusions The primary goal of this study was to develop a quantitative model of hydrate formation within a porous medium. This was accomplished for the case of a porous half-space which is cooled on its boundary. Prior to cooling, the temperature and gas concentration are assumed to be uniform throughout the half-space. As the boundary is cooled, a hydrate region grows into the porous half-space on the thermal diffusion timescale. The model predicts that the interface between the hydrated and fluid saturated sediments moves with the square root of time. The volume fraction of the pore space occupied by hydrate is constant throughout the developing layer, and is predominantly controlled by the initial gas concentration, though chemical diffusion may enhance the hydrate concentration by supplying additional gas to the hydrate region. The governing equations for this problem are derived from on conservation of mo-mentum, mass, energy, and gas. The governing equations contain transport terms which describe both advective and diffusive process. Scaling arguments show that advective transport is dominated by diffusive transport of heat and mass. Since advective trans-port may be neglected, there is no inherent length scale in the problem and similarity solutions may be obtained for the evolution of the temperature and gas concentration profiles as a function of time. The modelled values of the interface velocity and the hy-drate volume fraction are found by solving a transcendental equation which results from the conservation of energy and gas on the moving boundary. A laboratory apparatus was constructed to simulate the conditions considered in 73 Chapter 6. Conclusions 74 the model. The reaction chamber is made from a 0.7 m long, 0.14 m inner diameter cast acrylic tube. The ends are capped with aluminum lids through which coolant is circulated in order to control the boundary temperatures. The rate of growth of an ice layer was monitored using RTD temperature probes and shown to conform to the model predictions. This lends additional support to the validity of neglecting advection in the energy conservation equation. To test the model predictions for hydrate formation, carbon dioxide was chosen as the gas component because of its high solubility and relatively low hydrate stability pressure. Repeated attempts to detect the position of the hydrate interface using temperature measurements failed. This is interpreted to be a result of both the low hydrate volume fraction present during the experiments, and a combination of other thermal signals which were larger than the effects of the latent heat released during hydrate formation. Hydrate was removed from the reaction chamber after one experiment when the base temperature was low enough to promote ice formation, but the RTD probes failed to detect its presence. In order to test the formation model, improvements to the apparatus instrumentation will be needed to facilitate hydrate detection. Using the results of further laboratory tests, the validity of the model may be more properly assessed. If results are unfavourable, a review of the primary assumptions and simplifications inherent in the model may be in order. If, however, the model predictions agree with the experimental evidence, future analysis should focus on applying the governing equations to the task of determining the formation characteristics of hydrate in marine sediments. Developing a quantitative model of hydrate formation in marine sediments will be complicated by the geothermal gradient, which imposes a length scale on the problem so that similarity solutions of the governing equations will not be possible. In addition, hydrate layers in marine sediments occur over a vertical extent on the order of 100 m, and this may result in a significant Chapter 6. Conclusions 75 depth dependence for both the equihbrium temperature and the equihbrium gas concen-tration. These mathematical difficulties must be overcome to further contribute to our understanding of hydrate layer characteristics and their implications for more fundamen-tal problems. References [1] Abramowitz, M., and L A. Stegun, Ed., Handbook of mathematical functions, p. 297, Dover Publications Inc., New York, 1972. [2] Anderson, D. L., Theory of the Earth, pp. 79-102, Blackwell, Boston, 1989. [3] Appenzeller, T., Fire and ice under the deep-sea floor, Science, 252, 1790-1792, 1991. [4] Bishnoi, P.R., and P.D. Dholabhai, Experimental study on propane hydrate equi-hbrium conditions in aqueous electrolyte solutions, Fluid Phase Equil., 83, 455-462, 1993. [5] Bbzzo, A.T., Chen, H.-S., Kass, J.R., and A. J. Barduhn, The properties of hydrates of chlorine and carbon dioxide, Desalination, 16, 303-320, 1975. [6] Claypool, G.E., and I. R. Kaplan, The origin and distribution of methane in marine sediments, in Natural gases in marine sediments, edited by I. R. Kaplan, pp 99-139, Plenum, New York, 1974. [7] Cha, S. B., Ouar, H., Wildeman, T. R., and E. D. Sloan, A third-surface effect on hydrate formation, J. Phys. Chem., 92, 6492-6494, 1988. [8] Cherskii, N. V., and E. A. Bondarev, Thermal method of exploiting gas-hydrated strata, Sov. Phys. DoM., 17, 211-213, 1972. [9] Collett, T.S., Potential of gas hydrates outlined, Oil and Gas J., June 22, 84-87, 1992. 76 References 77: [10] Collett, T.S., Natural gas hydrates df the Prudhoe Bay and Kuparuk River area, North Slope, Alaska, Am. Assoc. Pet. Geol. Bull, 77, 793-812, 1993. [11] Collins, B.P., and J.S. Watkins, Analysis of a gas hydrate off southwest Mexico using seismic processing techniques and deep see drilling project leg 66 results, Geophysics, 50, 16-24, 1985.. . [12] Denbigh, K., The principles of chemical equilibrium, 4 ed., pp. 19-195, Cambridge Univ. Press, Cambridge, Great Britain, 1981. [13] Englezos, P., and P. R. Bishnpi, Prediction of gas hydrate formation conditions in aqueous electrolyte solutions, AIChE J., 34, 1718^ 1721,1988. [14] Englezos, P., Computation, of the incipient equilibrium carbon dioxide hydrate for-mation conditions in aqueous electrolyte solutions, Ind. Eng. Chem. Res., 31, 2232-2237, 1992. . [15] Fogg, P. G. T., and W. Gerrard, Solubility of gases in liquids, pp. 113-159, /Wiley; Chichester, Great Britain, 1991, , [16] Ginsburg, G. A., Kremlev, A. N., Grigor'ev, M. N., Pavlenkin, A. D., Larkin, G. V., and V. Ye. Tsar'kov, The discovery of solid gas hydrate in the rock at the foot of the continental slope of the Crimea, Trans. (Dokl.) USSR Acad. Sci. Earth Sci. Sec, 309, 83-85, 1989. , [17] Hammerschmidt, E. G-, Formation of gas hydrates in natural gas transmission lines, Ind. Eng. Chem., ^ 5, 851-857, 1934. [18] Hobbs, P. V., Ice physics, pp. 346-363, Clarendon Press, Oxford, 1974. References 78 [19] Honarpour, M., Koederitz, L., and A. H. Harvey, Relative permeability of petroleum reservoirs, 143 p., CRC Press, Boca Raton, Florida, 1986. [20] Hyridman, R. D., and G. D. Spence, A seismic study of methane hydrate marine bottom simulating reflectors, J. Geophys. Res., 97, 6683-6698, 1992. [21] Hyndman, R. D., and E. E. Davis, A mechanism for the formation of methane hy-drate and seafloor bottom-simulating reflectors by vertical fluid expulsion, J. Geo-phys. Res., 97, 7025-7041, 1992. [22] Hyndman, R. D., Wang, K., Yuan, T., and G. D. Spence, Tectonic sediment thicken-ing, fluid expulsion, and the thermal regime of subduction zone accretionary prisms: The Gascadia margin off Vancouver Island, J. Geophys. Res., 98, 21865-21876, 1993. [23] Jeffrey, G. A., and R. K. McMullen, The clathrate hydrates, Prog. Inorg. Chem., 8, 43-108, 1967. [24] Kroos, B. M., and D. Leythaeuser, Experimental measurements of the diffusion parameters of light hydrocarbons in water-saturated sedimentary rocks-II. Results arid geocheniical significance, Org. Geochem, 12, 91-108, 1988. [25] Kvenvolden, K. A., Gas hydrates - Geological perspective and global change, Rev. Geophys., 31, 173-187, 1993. [26] Landau, L. D., and E. M. Lifshitz, Fluid mechanics, 2nt* ed., pp. 227-237, Pergamon Press, Oxford, 1987. [27] Lederhos, J. P., Christiansen, R. L., and E. D. Sloan Jr., A first order method of hydrate equihbrium estimation and its use with new structures, Fluid Phase Equil., 83, 445-454, 1993. References 79 [28] Lide, D. R., Ed., CRC handbook of chemistry and physics, 71st ed., p. 6-8, CRC Press, Boca Raton, Florida, 1990. [29] Loehle, C, Geologic methane as a source for post-glacial, CO2 increases: The hy-drocarbon pump hypothesis, Geophys. Res. Lett., 20, 1415-1418, 1993. [30] Mashirov, Yu. G., Stupin, D. Yu., Ginsburg, G. D., and V. A. Solov'yev, Simulation of the generation of hydrates from dissolved gas, Trans. (Dokl.) USSR Acad. Sci. Earth Sci. Sec, 316, 176-179, 1991. [31] Minshull, T. A., Singh, S. O, and G. K. Westbrook, Seismic velocity structure of a gas hydrate reflector, offshore western Columbia, from full waveform inversion, /. Geophys. Res., 99, 4715-4734, 1994. [32] Nisbet, E. G., The end of the ice age, Can. J. Earth Sci., 27, 148-157, 1990. [33] Paul, C. D., Ussier III, W., and W. P. Dillon, Is the extent of glaciation limited by marine gas-hydrates?, Geophys. Res. Lett., 18, 432-434, 1991. [34] Pflaum, R. O, Brooks, J. M., Cox, B., Kenicutt, M. C. II, and D. D. Sheu, Molecular and isotopic analysis of core gases and gas hydrates, deep sea drilling project leg 96, pp. 781-784, Init. Repts. DSDP, 96, Ed. Bouma et al, U. S. Govt. Printing Office, Washington, 1985. [35] Phillips, O. M., Flow and reactions in permeable rocks, 285 p., Cambridge Univ. Press, New York, 1991. [36] Selim, M. S., and E. D. Sloan, Heat and mass transfer during the dissociation of hydrates in porous media, AIChE J., 35, 1049-1052, 1989. References / 80 [37] Sloan, E. D., Clathrate hydrates of natural gases, 641 p., Marcel Dekker, New York, 1990. • ' . •• ~.. . ' "' • [38] Stewart, P. B., and P. Munjal, Solubility of carbon dioxide in pure water, synthetic sea water, and synthetic sea water concentrates at —5 to 25 °C and 10-. to 45-Atm. pressure, /. Chem. Eng. Data, 15, 67-71, 1970. [39] Takenouchi, S., and G. C. Kennedy, Dissociation pressures of the phase C02 • 5 3/4 H20, J Geol., 73, 383-390, 1965. [40] Tsypkin, G. G.; Dissociation of gaseous hydrates in beds, J, Eng. Phys., 60, 556-561, 1991. . • • .; . . •; ' V' :y, . [41] Tsypkin, G. G., Appearance of two moving phase transition boundaries in the dis-sociation of gaseous hydrates in strata, Sov. Phys. Dokl. , 37, 126-128, 1992a. [42] Tsypkin, G. G., Effect of liquid phase mobility on gas hydrate dissociation in reser-voirs, Fluid Dynamics, 26, 564-572, 1992b. [43] Weast, R. C., Ed., CRC handbook of chemistry and physics, 56 t n ed., p. F-5, CRC Press, Cleveland, 1975. •[44} Yuan, T., Spence, G. D., and R. D. Hyndman, Seismic velocities arid inferred porosi-ties in the accretioriary wedge'sediments at the Cascadia Margin, J. Geophys. Res., 99, 4413-4427, 1994. U . : Appendix A Nomenclature The following is a list of symbols used throughout the text. a(t) interface position d(t) ;.' dimensionless interface position c \- fluid gas concentration Ce • • • • • equilibrium gas concentration CQO ••• ••• initial gas concentration Ch , hydrate gas concentration c dimensionless gas concentration Ch hydrate gas concentration parameter, see (2.49) C ; specific heat capacity at cbnstant pressure De chemical diffusion coefficient Dd chemical dispersion coefficient dS surface element outward normal unit vector g(h) relative permeability function, see (2.21) h hydrate pore volume fraction h* additional hydrate formed below ice melting temperature H . specific enthalpy AHexs enthalpy change on exsolution of gas ic • ' .V.. diffusive gas flux id • — dispersive gas flux 81 Appendix A. Nomenclature 82 k sediment permeability kp > barodiffusion coefficient Ke(K) , effective thermal conductivity L , latent heat of formation u macroscopic transport velocity u dimensionless transport velocity v interstitial fluid velocity P total pressure P' non-hydrostatic pressure component Q heat flux ST • ;. Stefan number, see (2.47) t time t ; dimensionless time T : temperature Te equihbrium temperature To base temperature Too • : initial temperature T dimensionless temperature Te dimensionless equihbrium temperature ATcornp temperature change from adiabatic compression ATdis temperature change from dissociation of hydrate ATex, temperature change from exsolution of gas U internal energy V volume W work z • • • • position Appendix A. Nomenclature 83 z .' dimensionless position a : Lewis number, see (2.48) OCT • • • • • • • • • coefficient of thermal expansion T(h) heat capacity ratio, see (2.21) To = r(0) heat capacity ratio for fluid saturated sediments Th — T T T - • heat capacity ratio for hydrated sediments Ke(h) Ti = ^ ^..e heat capacity ratio of icy sediments rj dynamic viscosity K e ( h ) effective thermal diffusivity, see (2.22) (0) ..' thermal diffusivity of fluid saturated sediments A ;.... growth-rate parameter, see (3.24) ft • • • • chemical potential £ : similarity variable, see (3.23) p -..i. density 4> porosity ft fixed representative volume subscripts / fluid h .hydrate s ; sediment grain i ice r representative quantity Appendix B Interface Conditions In Chapter 2, the conservation laws are applied over a continuum of solid grains, pore fluid, and hydrate to derive the governing equations. As a hydrate layer is formed, the interface at the top of the layer is defined by a discontinuity in the hydrate volume fraction h. This discontinuity in h will introduce discontinuities in the temperature gradient and the gas concentration gradient. The conservation laws may be applied in integral form across the interface to determine how the magnitudes of these discontinuities and the velocity of the moving phase change interface are related. In the following derivation, the interface is represented as a surface z = Zi(x,y,t) with unit outward normal ii. As shown in Chapter 2, the conservation laws may be represented in terms of surface integrals and volume integrals applied over some fixed region. To evaluate these inte-grals over a region where the integrands are discontinuous, consider a fixed box which is constructed to enclose part of a mobile interface in such a way that its top and bottom surfaces, cr+ and <r-, are parallel to the interface and equal in area. The other sides of the box have a much smaller surface area so their contribution to any surface integral over the box is negligible compared to the contributions from cr+ and cr~. The normal component of a vector field q may be integrated over the surface of the box to obtain ^ q - d s « f+ q-nd(T-rJ_<i-(-n)dcr . (B.l) If the box is made small enough that the component of q perpendicular to the interface 84 Appendix B. Interface Conditions 8 5 is uniform over each surface integral, then (B.l) may be simplified to j^q-dS w [ci-A]tJ d<r = [q-n]j> , (B.2) where the notation [q • n]t is used to denote the magnitude of the discontinuity in the normal component of q across the interface. The time derivative of the volume integral of a scalar function F over the box, may be rewritten as d_ 8t f F <Kl = 4- [if Fdw}d<r , (B.3) Jn at J<r Jw where the surface o~ is again parallel to the interface, and w is the width of the volume ft. The differentiation may be performed within the surface integral since its limits are fixed in time and the integrand is continuous along every surface parallel to the interface Zj. However, F may be discontinuous across the interface, so it is useful to split the integral over w into parts above and below z,- so that s/.(//^=/i(/.>i+r where w~ and w+ indicate the top and bottom of the box. The line integrals may be evaluated using Leibnitz' rule, which states that for a function $ which is continuous over the interval from ftj to b2, mi "•"•^ "•'"•St+X mdz- (B'5) The lower bound of the first integral and the upper bound of the second integral on the right hand side of (B.4) are fixed, only their common bound at Z{ moves. If the values of F on either side of the interface are given by F+ and F~, then Appendix B: Interface Conditions 86 The box is considered to be thin enough that F may be treated as a constant on each side of the moving interface, so (B.6) becomes 4- I ? dn = -[F}+^<r+[ ?f d£l . (B.7)' dtJo 1 J" dt Jet dt v ' As the width of the box tends to zero, the volume integral on the right hand side of (B.7) becomes much smaller than the first term since the volume shrinks more rapidly than the surface area, leaving Equations (B.2) and (B^ 8) may be applied to obtain interface conditions for the hydrate formation problem from the conservation laws. The integral form of the condition for conservation of mass may be stated as (see (2.4) and (2.5)) | - / [Pf<t>{l - h) + Ph<t>h + p,(l - <j>)} dSl = - [ pfu • dS . (B.9) ot Jn J s Evaluating these integrals over the box on the interface gives -[pf<t>(l -h) + ph4>h + ps{\ - tfHi^fr = -[Pfu • &]±<r ... (B.10) This equation may be simphfied to pf[n-n)t = 4>{PH-Pt)[h]t^, (B.ll) because the incompressibility conditions imply that the porosity and the densities of the three components cannot change. Notice that if the density difference between the hydrate and fluid is small, and the magnitude of the discontinuity in the hydrate volume fraction is small, the velocity field will be continuous across the moving interface. The same type of analysis can be used to find an interface condition from the energy conservation condition. The balance requirement in this case is (see (2.14), (2.15), and Appendix B. Interface Conditions 87 (2.16)) ^ [[pfHf4>(l-h) + phHh<f>h + P,Ha(l-<f>)} dfl = - f pfHfu-dS+ f Ke{h)VTdS . at Ja Js Js (B.12) The temperature must be continuous across the boundary since any discontinuity would require an infinite thermal gradient and, by extension, infinite heat flux. Thus the tem-perature dependence of enthalpy from (2.18) implies that the enthalpies of the individual components are continuous, so only the hydrate volume fraction, the temperature gra-dient, the fluid velocity, and the effective thermal conductivity may be discontinuous. Using (B.2) and (B.8), the equation can be re-written as <f>{PfHf - phHh)[h]±^o- = (-PfHf[u • n]+ + [Ke(h)VT • n]t)<r . (B.13) The definition of latent heat from equation (2.19), and the interface condition from conservation of mass (B.ll) may be invoked to show that P^[h]t^ = Mh)VT.n]t, (B.14) Pft'f at •where the effective thermal diffusivity Ke(h) is denned as in (2.22). Conservation of gas on the boundary requires that (see (2.24), (2.25), and (2.32)) / [pfC<f)(l — h) + phCh<f>h] dQ = — Pfcu • dS + fgPf<H(l-h)[(De + Di)Vc + De!^VP)]}'dS.. (B.15) The pressure and the concentration of gas in the fluid must be continuous across the boundary to avoid infinite gradients. Equation (B.15) may therefore be written as -<j>(phch ~ pfc)[h]t^ = -pfc[u • n]+ +p,<p[(l - h){(Dc + Dd)Vc + D^VP} • n]t . (B.16) 8_ 8t Appendix B. Interface Conditions 88 Substituting in the interface requirement from mass conservation, (B. 11), and dividing by porosity gives -Ph(ch - c)[h]t^ = P,[(l ~ h)(Dc + A , ) V c • n]+ -f pf[(l - h)Dc^VP' • f\]t. (B.17) In (B.17), the pressure gradient VP has been replaced by its non-hydrostatic com-ponent VP'. Using momentum conservation, (2.1), to replace VP' allows (B.17) to be rewritten as -Ph(ch - c)[h]+-§ = P,[(l - h)(Dc + Dd)Vc • A]t + ^ ^ [ ^ ^ « ' AH • (B.18) While the form of the relative permeability function g(h) depends on the location of the hydrate in the pore space, it might be reasonable to approximate it as a hnear function of h (i.e. g(h) = 1 — h). This function has the desired properties that it returns an effective permeability of 0 when the pore space is completely clogged with hydrate, and the intrinsic permeability k when there is no hydrate present. This choice for g{h) results in the final term in (B.18) being directly proportional to the jump in the normal component of the fluid velocity. Provided that either this velocity jump is small, or the coefficient of this final term is small, (B.18) may be approximated by -Ph(ch-c)[h]+^ = pf\(l-h)(De+^ (B.19) The jump conditions on the temperature gradient, (B.14), and the gas concentration gradient, (B.19), may now be modified for the geometry considered in Chapter 3 and transformed to dimensionless form using the variables introduced in (3.1) and (3.2). The interface is planar so the gradients of the temperature and gas concentration are normal to the moving phase change interface at Zi = a(t). The hydrate volume fraction h is zero above the moving boundary, and some finite value beneath it. The gas concentration within the hydrated region is at ce. With these substitutions the nondimensional interface Appendix B. Interface Conditions 89 conditions become dd' array"), Ke(h)df(d-,i) ~OThtl-r= = wz —r — , (B.20) dt dz «e(0) dz and —ch(j)ah—= = — — — , (B.21) Pf dt dz where the Stefan number Srh, the hydrate gas concentration parameter Ch, and the Lewis number a are denned by equations (3.3), (3.4), and (2.48), respectively. Appendix C Stability Analysis The hydrate layer formation model presented in Chapter 3 is based on a simplified view of nature. In reality a medium always contains small inhomogeneities which may cause in-stabilities that magnify with time. A stability analysis is necessary to determine whether perturbations to the phase change interface are likely to magnify or decay. Another effect that requires investigation is that of thermal convection driven by density contrasts in the interstitial fluid. If transport by convection were important, the model presented in Chapter 3 would require the addition of convective transport. C . l Interface Stability In some moving boundary problems, perturbations to a planar interface may lead to the growth of dendrites which invalidate one dimensional model formulations. To test the effect of perturbations on the model presented in Chapter 3, solutions to the governing equations are sought for a modified phase change interface geometry described by • z = A(x,t) = d(t) 4 : S d e i k i + f f i ( C . l ) The interface A is described by a planar component a and a component with a low amplitude periodic spatial dependence. If the growth parameter cr is positive for a given spatial wave number k, the amplitude of the perturbation will increase exponentially with time t suggesting that the planar interface model found in Chapter 3 is unlikely to be physically realized. If cr is negative, however, the disturbance will decay so the model 9) Appendix C. Stability Analysis 91 is stable to interface perturbations with that particular spatial wave number k. The nondimensional governing equations for this problem are r ^ - ^ + « P ••<*<**•'>' ( c - 2 ) . ' / ' % = W*W (c-3) c = 0 , 0 < i < A(x, i) (GA) and ± dc d2c d2c . , which are the same as (3.7) through (3.10) with extra terms to account for diffusion in the x direction. The boundary conditions and initial conditions become (see (3.13) through (3.20)) T(5,z,0) = l , c(x,z,0) = 1 , f(x,0,i) = 0 , (C.6) 2'(5,oo,£) = 1 , c(x,oo,i) = 1 , , (C.7) f(x,A,i) = feh , and c(x,A,i) = 0 . (C.8) The changing orientation of the outward normal fi causes the interface conditions, (B.14) and (B.19), to become - - - K. (h) - - - ft A VT(x,A+,i)-n-^^VT(x,A-,t)-n = -SThh-^ , , (C.9) l^e • at and Vc(x,A+,i)-h = ^ C h f a h ^ - , ( C I O ) Pt Ot where the symbol V is used to denote the nondimensional gradient. The strategy is to break the temperature and gas concentration fields into zeroth-and first-order components, and look for solutions for each component separately in the form • f(x,z,i) = f0(z,$ + 8f(z)eih*+*i , (C.ll) Appendix G. Stability Analysis 92 and c(x,z,i) = c0(z,i) + 8c(z)eik*+°i,^ (C.12) with the hydrate volume fraction given by h = ho + She****1 . (C.13) The perturbation quantities, 65,, ST, 8c,. and Sh, are assumed to be small enough that second-order contributions may be ignored. For boundary conditions on A(x,t) a Taylor series expansion about a(t) is used so that for a function /, f(A) = f(a) + Saeik^ai^^- + 0{Saf . (C.14) oz Substituting (C.ll) and (C.12) into (C.2) through (C.5) and collecting the zeroth-order terms gives r ^ = p^, 0<z<A(x,i) (C.15) am r)2T To~df = ~di ' (C16) co = 0, 0<z<A(x,t) (C.17) and dc0 d C o • ^a~dl = ~di' z>A(x,t) , (C.18) while the associated boundary conditions, (C.6) through (C.8), become ro(z,0) = l , 2,(5,0) = 1, f o(0,i) = 0, (C.19) fo(oo,*) = l , c0(oo,f) = 1 , (C.20) fo(a,i)-feh , co(x,a,t) = 0. (C.21) Equations (C.13) and (C.14) are used to produce the zeroth-order interface conditions from (C.9) and (CIO) giving , dfo(a+,i) Ke(h)df0(a-,i) da —T-^Z -OThho—- , (C .22J Oz Ke Oz dt Appendix C. Stability Analysis 93 and dc0(a+,i) Ph,.,da . — — = _ ^ a C h f e 0 _ . (C.23) az pf dt This problem is identical to that considered in Chapter 3, with the solution described by equations-(3.31), (3.32), (3.41), and (3.42), which may be written with the similarity variable £ replaced using its definition from (3.23) as - - - -t^tfo c0(z, i) = 0,0<z<A (C.26) and The interface conditions, (C.22) and (C.23), lead to ^° phy/ir<f>aChX eiic(y/fac\) ' ^ ^ and (1 - Te^yTo" _ Ke(h)Teh\/Th pfSrh Q 2 Q x eroA 2erfc(v^A) Keer^2erf(vT^A) + P f c y ^ e ^ e r f c t y ^ A ) ' . K ' J which determine h0 and A. The first-order terms from the substitution of (C.ll) and (C.12) into (C.2) through (C.5) describe the conditions on the perturbation quantities, namely Tk(T8f=<^Jf-k26f, 0<z<A(x,t) (C.30) dz2 T o a S f = y ± - k 2 6 f , z>A(x,t) (CM) d2_St dp 8c =0 , 0<z<A(x,t) (C.32) Appendix C. Stability Analysis 94 and d2 8c <j>cto-8c = - k28c , 2 > A(x,i) (CM) CLZ while the associated boundary conditions, (C.6) through (C.8), become ST'(co) = 0 , Sc(oo) = 0 , 8T(0) = 0 , (C.34) Sf(a+) = -Sa^p^- , 8c(a+) = - S a ^ ^ - , (C.35) az oz 8f(a-) = - S a ^ ^ , Sc(a~) = - S a ^ ^ - , (C.36) oz Oz and the interface conditions, (C.9) and (C. 10),.become i stji*) _ M M i ay-) = _ S T h ( h o „ s - a + S-A dz Ke dz dt OZ2 Ke OZ2 and d8c(a+) ph, cld°-\ c~d2c0(a+,i) . Tz—- = —ChpalhocrSa + Sh—f) - Sa ^ — - . (C.38) dz pf dt azi The solutions for ST and Sc. are found in terms of exponential functions and the coeffi-cients are evaluated using the boundary conditions, (C.34) through (C.36), to give and <9z eV^+iV- 5 - e-\A2+r„<r a = - ^ ( " M ^ - y W i V , 5 > i (C.40) dz 6c = 0 , 0 < z < A (C.41) Sc=-d°°(S'+'thae-T/ka++°" (*-«) . 5 > " i (C.42) dz All that remains is to evaluate the interface conditions, (C.37) and (C.38), to ob-tain two transcendental equations involving the unknown perturbation hydrate volume Appendix C. Stability Analysis 95 fraction Sh and the exponential growth factor a. Algebraic manipulations eliminate Sh to yield cr as a function of the wave number k. Evaluating (C.37) while making use of the zeroth-order interface condition from energy conservation, (C.22), and the location of a(i) from the zeroth-order solution (d(i) = 2\Vt), gives the following result: dT0(a-,t) K(h) A 2 + I > + ^ 2 + r V r coth(y/k* + rhcr~a) + (Th - T 0 ) f } oz KeSxhho dt /— ~—dd „ (dd\2 Sh dd . _ ,„x The interface condition from gas conservation, (C.38), may also be evaluated to give r~ ; dd , ( dd\2 Sh dd + = ^ - O T . (C.44) Equating the left hand sides of (C.43) and (C.44) produces dT0(d-,t) K(h) { v / j ^ ^ + > / f t r T ^ oz KeZ>Thho dt i- iJ~\2 +(yjk* + <f>a<r - y/k* + T0<r)-J - (<f>a - T 0 ) ^ = 0 , (C.45) which is a transcendental equation relating cr and k. Since the hydrate phase is not expected to occupy a significant portion of the pore space, the thermal properties of the hydrate layer and the overlying region should be relatively uniform. Making the approximation that r / l = r 0 = r (C.46) and Ke(h) = Ke (C.47) allows (C.45) to be rewritten as ^ " f ' ^ 1 , (Vk^T^+Vh^Y^cothWk^Yo'd) oz Sxhho v ' Appendix C. Stability Analysis 96 +(y/k* + <t>a<r- y/V + T<r)-J - (<f>a - T0) ^ JJ =0. (C.48) Equation (C.48) indicates a rather comphcated relationship between the exponential growth factor cr and the wave number k. Since the purpose of this analysis is only to determine whether the disturbance grows or decays for a given wave number, the next step is to determine the condition for marginal stability. This is found by setting a equal to zero in (C.48) giving d T ° % r t ) S™1 (1 + coth(\km\a)) - (*a - To) ( f = 0 . (C.49) The first term on the left hand side of (C.49) is always positive whereas the term involving the square of the zeroth-order velocity da/dt is always negative since <f>ct is much larger than T; thus there should be at most one root to this equation. However, using the nominal parameter values from Table 3.2 it was found that a plot of the left-hand side of (C.49) as a function of km was always positive for a given time t. This suggests that for the nominal parameter values, there is no interface perturbation wave number that is marginally stable. If the ratio of the thermal and chemical diffusion coefficients is increased so that <f>a takes a value of 5000, there is a solution to (C.49) which is plotted in Figure C.l as a function of time. The stable and unstable regions may be assigned by evaluating (C.48) using small positive and negative values of cr. Examining equation (C.48), the third term on the left-hand side is the only term which is always negative. The second term may be positive or negative, but it should be small for a small magnitude of cr. The first term involves both k and cr in groupings of the form \/k2 + Tcr, hence an increase in k above the marginal value km decreases the exponential growth parameter cr (e.g. makes cr negative) in order to balance the third term, which is not a function of k or cr. Figure C.2 shows a plot of the wave number A; as a function of cr at t — 2 using <j>a = 5000 which indicates that Appendix C. Stability Analysis 97 <JKX Unstable <j>0 o i 1 : I I T~ 0 0.5 1.0 1.5 2 0 Time't Figure C.l: Wave Number for Marginal Stability as a Function of Time The wave number for which cr = 0 is plotted as a function of nondimensional time i with the parameter <f>ct = 5000. The region above the line is stable (i.e. a < 0) and the region below is unstable (i.e. cr > 0). Decreasing the value of <f>a shifts the marginal stability line down to the left making smaller wave numbers stable at any give time. wavenumbers k < km correspond to unstable perturbations. When cr = 0, the solution to (C.48) defines km. For values of cr greater than 0, the wave number drops below km so the lower region in Figure C.l is unstable, while the upper region is stable. The increase in the growth rate cr with decreasing wave number means that the longest wavelength interface perturbations are the most unstable. Since the case where k = 0 corresponds to the infinite wavelength perturbation, and this mode has the fastest growth rate, the model of Chapter 3, which is in effect an infinite wavelength interface model should be valid at these high values of <f>a. The nominal value for 4>a from Table 3.2 was 100. The solution to (C.49) may be found as a function of <f>a to determine what happens to km as <f>a is lowered, or equivalently as the rate of gas transport by diffusion becomes closer to the rate of diffusive heat transport. Figure C.3 is a plot of km as a function of <f>a at i = 2. It shows that km decreases to 0 when <f>a ?a 2750. This means that the graph of km as a function of time in Figure C.l Appendix C. Stability Analysis 98 0.5 1.0 1.5 2.0 Growth Rate, a Figure C.2: Wave Number as a Function of Growth Rate The wave number is plotted as a function of growth fate for <f>ct = 5000 at t stability is achieved when cr = 0, at km 2.3. moves down and to the left as chemical diffusion of gas becomes more efficient, allowing increasingly longer wavelength perturbations to fall within the stability zone above the marginal stability line. Once <f>ct is reduced below 2750 at this particular time t = 2, perturbations at all wave numbers are stable. C.2 The Importance of Thermal Convection The criteria for the onset of thermal convection in porous media are well established [e.g. Phillips, 1991, pp. 140-148] and will not be rederived here. The Rayleigh number, defined by Sa=p,aTtgATklr VKe is a nondimensional parameter which represents the ratio of buoyancy forces caused by the temperature dependence of the pore fluid density, to dissipative effects which tend to retard the motion. In the laboratory experiments described in Chapter 4, cooling = 2. Marginal Appendix C. Stability Analysis 99 k .1 i i i i. 3000 3500 4000 4500 5000 <l>a Figure C.3: Wave Number for Marginal Stability as a Function of (f>a This graph shows how the marginal stability wave number decreases with decreasing <f>a at i — 2. If (f>a is less than about 2750, all perturbations are stable. the base of the reaction chamber below 4°C results in a buoyancy force because of the unusual coefficient of thermal expansion aj/ of water at these temperatures (e.g. ax/ « —68 X 10 - 6 K - 1 [Weast, 1975]. The temperature change AT across the region which can convect is —4 °C, and the depth of this region, denoted by lr is typically 10-1 m in the experiments. Using nominal values for the fluid viscosity rj, density pf, and thermal diffusivity Ke, and an upper bound on the sediment permeability k, the Rayleigh number in the experiments is . Ra < 5 . (C.51) The minimum Rayleigh number required for the onset of convection is [Phillips, 1991] Ra" = 4TT2 . (C.52) In this case Ra « Ra* so that the fluid is convectively stable, and the assumption of diffusive heat and gas transport should be valid. Appendix D Effective Thermal Diffusivities The effective thermal diffusivities of the various mixtures of sand, water, ice, and hydrate present in the experiments discussed in Chapter 4 can be approximated using mixing laws and published values for the thermal properties of each component. It is also possible to obtain the effective diffusivity directly from experiments. This appendix describes the experiments that were performed on sediments both above and below the freezing point of water in order to find the effective thermal diffusivity in each case. D. l Fluid Saturated Sediments To calculate the thermal diffusivity of fluid saturated sediments, consider a semi-infinite sediment column at constant initial temperature T(z,0) = Too. If the base temperature is changed to and maintained at some new value T(0,i) = T 0 the temperature profile will adjust with time in a predictable manner which is only dependent on the bulk thermal diffusivity «. Mathematically, this problem may be stated,as °— = K V 2 T , z > 0 (D.l) subject to T(0,t) = To , and T(oo,*) = T^ = T(z,0) . (D.2) The solution is found by integrating (D.l) twice and applying the initial condition and boundary conditions, (D.2), to give T(z,t) = TooeT{(^) + T0eIfc(-?=). (D.3) 100 Appendix D. Effective Thermal Diffusivities 101 Note that the bulk thermal diffusivity K used in (D.l) is somewhat different than the effective thermal diffusivity Ke used in the earlier analysis (i.e. refer to (2.37)). These two parameters are related through the heat capacity ratio T(h), so that in this case Ke = T(h)K = To* . (D.4) An obvious method for determining K is to compare experimental data to the pre-dictions of (D.3) and see which value gives the best agreement. It was found, however, that better results were obtained by observing the times at which probes at two different locations record their maximum rate of temperature change. This method has the ad-vantage that it is insensitive both to small calibration errors in the absolute temperatures of separate probes. It is also insensitive to the exact temperature difference between Too and T 0 . Setting the second temporal derivative of (D.3) equal to zero reveals that the maximum rate of temperature change will occur at . < d - 5 > Using data from two probes removes any time-zero offset from the experimental data giving z2 - z2 K = 6(tm(zl) - • ( D ' 6 ) A series of experiments was carried out to find the maximum values of the time derivative of temperature measured at fixed locations in the fluid saturated sediment column. An example of the time derivative of the data from one such experiment is shown in Figure D.l. The times f.m(zj) of the maximum rate of temperature change at each of the three probes. Using these times and the probe locations Z{ in equation (D.6), I obtain a value for a of (5.3 ±0.5) x 10~7 m2/s which best fits the data. The corresponding effective thermal diffusivity Ke of the fluid saturated sediments, usingT(0) = 1.23, is (6.5±0.6) x 10~7m2/s. Appendix D. Effective Thermal Diffusivities 102 dT/dt ( ° C / min) - 0 . 2 h -20 SO Time (minutes) W z v> Figure D.l: Rate of Temperature Change vs. Time in Fluid Saturated Sediments Derivative of temperature data during experiments to determine Ke. Equation (D.6) is used along with the times of the minimum rates of temperature change tm(zi), tm(zi), and tm(zs) at probes located at z\ = 20 mm, z 2 = 35 mm, and Z3 = 50 mm. D.2 Sediments Containing Ice In order to experimentally determine the effective thermal diffusivity of sediments contain-g ice, the base plate of the reservoir was first cooled to some fixed temperature To below freezing. Once the temperatures measured at the probes reached near constant sub-zero values, it was inferred that the frozen layer had achieved its maximum extent. The temperature profile in the layer at this time is approximated by T{z,0) = T o + T e h / " z , o (D.7) where Teh is the freezing point of the pore fluid, and a is the thickness of the frozen layer. The next step was to change the base temperature to some new value T 0 + AT and monitor the temperature response at each probe. The temperature within the layer is described by r)T ^- = /cV2T, 0 < z < a (D.8) dt . Appendix D. Effective Thermal Diffusivities 103 subject to T(0,t) = To + AT , T{a,t) = Teh, (D.9) and T(z,0) = T o + T e h T ° z . . (D.10) a To solve this problem, the steady state solution is first found and subtracted from the temperature T to produce homogeneous boundary conditions in the residual problem. Separation of variables is then used to express the temperature profile in terms of a sine series giving nVt Teh-T0-AT I ~ 2AT -« T(z, t) = T0 + AT + z + X) — e fl2 «zn(—-) . (D.ll) a ~; n7r a Of course in reahty, the boundary at a must move, but by comparing (D.ll) with the data from the lower probes at early times, the effects of this motion should not be appreciable. Once again it was found that comparing the times of thepeak rate of temperature change at two or more probes to the predictions using the time derivative of equation (D.ll) gave the best results. The maximum rate of temperature change should occur when the second temporal derivative of (D.ll) equals 0. This condition is satisfied when • , n2ic2t £ > V a 2 sin(^) = 0 . (D.12) A plot of the time derivative of the temperature data from the bottom two probes in one experiment is shown in Figure D.2. The time difference between the peak rate of change at the two probes is about 2.9 ± 1.2 minutes, and the ice layer thickness was inferred to be about 75 ± 3 mm. Figure D.3 is a plot of the predicted times of the maximum rate of temperature change found using (D.12) at the two probe positions, Z\ = 20 mm and z 2 = 35 mm, as a function of k. The value of « determined from a number of different experimental runs using temperature data from the bottom three Appendix D. Effective Thermal Diffusivities 104 - 0 . 0 2 1 ' J ' ' ' ' - 1 ' 1 ' ' . J 1 O 2 4 6 8 1 0 Time (minutes) Figure D.2: Rate of Temperature Change vs. Time in Sediments Containing Ice Data is from the bottom two temperature probes at Z\ = 20 mm and z 2 = 35 mm, with the interface at a «. 75 mm. The shaded areas denote the uncertainty in the times of the peak rate of temperature change dT/dt at the two locations. The elapsed time between these peaks is Ai = 2.9 ± 1 minute. probes is (9 ± 3) x 10 7 m2/s. Using T(i) = 1.02 this translates to a value for «e(i) of (9 ±3) x 10~7 m2/s. Appendix D. Effective Thermal Diffusivities 105 - 7 - 7 - 6 - 6 - 6 6 . 10 8 . 10 1 . 10 1 . 2 10 1 . 4 10 K (mA2/s) Figure D.3: Predicted Times of Maximum Rate of Temperature Change These predictions derive from equation (D.12) and are for probes at Z\ = 20 mm and z2 = 35 mm with the interface at a = 75 mm. The shaded region encompasses the range of « that corresponds to the observed time difference of At = 2.9 ± 1 minute as shown in Figure D.2. Appendix E Ice Saturation Level lit applying the mathematical model to the freezing of water in the pore space in Chapter 4, I assumed that ice would completely fill the pores. There is the possibility, however, that the effects of salts dissolved in the water or the curvature of the pore walls or some other pore-scale phenomenon could limit the fraction of the pore space occupied by ice. To test whether any of these processes might be important, experiments were conducted in which the base temperature of the reaction chamber was cooled below 0°C and the frozen layer was allowed to form. After the frozen layer had encompassed all 4 temperature probes, the base temperature was brought above freezing and the resulting temperature profiles were monitored. At each probe position the change in base temperature caused a rapid increase in temperature until the freezing point Tej was reached. The temperature then stayed relatively constant for several tens of minutes, before beginning to rise again. By modelling this data, it is possible to estimate the frozen volume fraction. The form of the experimental melting curves in Figure E.2 suggests the following series of events, which are shown in the diagrams of Figure E.l. Initially the temperature gradient in the icy layer is essentially linear. Warming the base causes the temperature to rise rapidly to the melting point and the layer begins to erode from both top and bottom. The temperature quickly adjusts to become nearly uniform throughout the layer so that virtually all of the heat transported to its upper and lower boundaries must be consumed in the melting of ice. The lower regions are melted by heat transport from below, and the upper regions are melted by heat transport from above. The low magnitude of the 106 Appendix E. Ice Saturation Level 107 thermal gradient within the layer implies negligible heat flux between the two interfaces, so the melting rates of the upper and lower boundaries are independent of each other. The heat flux available for melting the lower boundary may be calculated from the known base temperature and the measured melting temperature at the probe locations. Observations of the times at which two or more probes are freed from ice may then be compared to provide estimates of the frozen volume fraction. The temperature profile in the melted region will be approximately linear so that if the temperature of the base is Ti, the melting temperature is Tei, and the position of the layer bottom is s(t), then T(z,t) w Ti - T l ~T"z . 0<z<s (E.l) s '• • The heat flux Q is related to the temperature gradient via the thermal conductivity k of the fluid saturated sediments so that Q = -k^^kTl~Tei . (E.2) oz s This heat flux must balance the latent heat required to melt the layer, giving k^^i » Li4>pii%, (E.3) s at where i is the pore volume fraction of ice pi is the ice density, and Li is the latent heat of melting per unit mass. The position of the moving boundary may therefore be represented as V^ W1'""^ " : (E-4) where s(0) is the unknown position of the boundary when the layer reaches a nearly uniform temperature and its internal heat flux becomes small. Using data from two probes at positions s(ti) and 5(^2), the pore volume fraction of ice may be calculated as 2KePfCf(T1-Tei)(Vt~2-Vh)2 Li(f>Pi(s(t2) - sfa))* (E.5) Appendix E. Ice Saturation Level 108 Stage 1 ice layer / T • Stage 3 ice layer s(t) 1 ft i:Ti' • Stage 2 * ice layer / T • 1ei Stage 4 / ice layer s(t) ft * : 1 ft • \ i r Figure E.l: Melting Behavior of an Ice Layer The melting of an ice layer when the base temperature is warmed to T\ > Tei seems to proceed in a number of stages. To begin with, in stage 1, the thermal gradient in the ice layer is essentially linear, and the layer growth rate is small. The base temperature is warmed in stage 2, but has yet to achieve its constant value of T\. The layer begins to melt from below in stage 3, and the thermal gradient within the relaxes towards the melting temperature Te;. In stage 4, the gradient within the ice layer is negligible, and the temperature profile below the interface s(i) is essentially linear. Appendix E. Ice Saturation Level 109 where the thermal conductivity k has been replaced with the product of the effective thermal diffusivity Ke(i) and the heat capacity p/Cf. Defining a Stefan number as STl = rP(tLl rp x (E-6) pfCf{li - lei) allows (E.5) to be rewritten in the form: (E.7) ST1(s(t2) - sih))* • The.form of (E.7) indicates that larger values of STI will tend to separate the melting times, ti and t2. This is reasonable since STI is large if the base temperature Ti is only shghtly above Tei- The small temperature difference reduces the heat flux and slows the melting rate. In order to ensure that for at least two probe positions the ice has melted from below rather than above, there are some limitations on how close to the melting temperature To may be. Keeping these concerns in mind, the data shown in Figure E.2 was acquired. The fact that the third probe stayed at Tei longer than the second is inferred to mean that the bottom two probes did indeed melt from below. Using the parameter values summarized in Table E.l, in equation (E.7), the frozen volume fraction was found to be about 0,9. It was expected to be around 1, indicating a fully ice saturated sand, and this calculated value is certainly in the proper range considering that extreme parameter values using the errors listed in Table E.l give a range of 0.4 to 1.9 for i. The large error in the ice saturation level calculated in this manner is disappointing, however the long time during which the equilibrium temperature was maintained is suggestive of a possible basis for a means to detect small hydrate volumes. Appendix E. Ice Saturation Level 110 _g I * • • • I i i i—i—L_i i i i I i i—i—i—I—i i i i I—i—i—i—i—I—i—i—i—i—I—i—i—i—i—1—i—i—i—i 0 5 0 100 150 200 250 300 3 5 0 4 0 0 Figure E.2: Melting Experiment to Determine Ice Saturation Level The data shown here is from the bottom three probes. The ice surrounding the two lower probes melted before the ice around the third probe did, so the assumption that the ice at the bottom two probe positions was melted through heat transport from below is valid. The shaded areas indicated error bars on the melting times, ti and t2. Parameter Value Error Units T i - r « i 1.0 0.2 ° c STI • - r e t ) 29.5 « e 6.5 x 10"7 0.6 x lO" 7 m2/s 5 ( t 2 ) - 5 ( i i ) 0.02 0.002 m ti 46 8 min. t2 337 12 min. Table E.l: Data For Calculating Approximate Ice Volume Fraction Appendix F Data Acquisition The voltage drop across each pair of leads connected to an RTD probe must be converted into a temperature measurement and dumped to a data file. The first step in this: task is accomplished using the circuit shown in Figure F.l. The resistance of each probe is 100 ft at 0 °C, and it increases with temperature at the rate 0.39 ft/ °G. Each probe is in series with a 1.5 kft resistor and each resistor-RTD pair is connected in parallel to a battery with output voltage VB- A precision 100 ft ? . • • ' • ' . -resistor which is insensitive to temperature fluctuations is connected in the same way as the RTDs. Leads from each of nodes 1 through 4 on the circuit diagram are brought to one input of four separate AM 502 Tektronix differential amplifiers. A lead from node 0 is brought to the second input of each amplifier. The difference between the voltages measured at node 1 and node 0 with respect to the common ground is _ / i i i : 100 ft V 1 ~ ° ~ B W+ 1-5 kft ~. 100 ft + 1.5 kft) " ( } Since the RTD resistance is nominally 100 ft, and the temperature variations cause changes that are much smaller than the combined 1.6 kft series combination, (F.l) can be approximated as ' __ T r T r (Rx - 100 ft\ Equation (F.2) appfies to nodes 2 through 4 with R\ replaced by i?2, Rz, and R4 re-spectively. The differential amplifiers are used to boost these signals Gr> times, and the outputs are sent to channels 0 through 3 on an A/D card. '.' I l l ; '•. ' . Appendix F. Data Acquisition 112 1.5 kQ 1.5 Kli 1.5 kfl 1.5k£2 <1 . 5kQ % 1 , 5 k Q S 1 . 5 k Q I 3 R R R T. R4 T T "2 T D D D D 100 Q Figure F.l: Circuit For Finding Changes in Probe Resistances The four RTD probes are connected to this circuit via leads which pass through the base of the reaction chamber as shown in Figure (4.2). Nodes 1 through 4 are each connected to one input of four separate differential amplifiers, which have node 0 as their second input. Node 5 is connected directly to an A/D card to monitor the battery output voltage vB. ; To account for variations in the battery output voltage VB, the voltage at node 5 is measured relative to the common ground. The leads are attached directly to channel 4 on the A/D card, which has a maximum input voltage of 5 volts. A 6 volt gel-cell battery is used, hence the parallel 1.5 kfi resistors. A program on a 386 personal computer is used to collect data from the A/D card and convert the readings at channels 0 through 3 into temperature measurements. The leads which carry the current from the circuit board through the base of the reaction chamber and on to the RTD probes have some finite resistance rn (where n runs from 1 through 4). The resistance of each of the probes themselves behaves in a predictable way with temperature T so that for the first probe for example Rr = 100 Cl + 0.3m ft/ °C + rx . (F.3) If the signal at each of channels 0 through 4 is defined as Vcn (where n runs from 0 Appendix F. Data Acquisition 113 through 4), equation (F.3) maybe used along with (F.2) to show that V—TVUQo^ 1 6 x l 0 3 J , (F.4) which may be solved for Ti to give To find the lead resistance, the probe is put in ice water and (F.5) is solved for r\ using the observed values of VC0 and V^, and the known values of Ti and GD- Values for through r 4 were all less than 2 ft. The gain on the differential amplifiers GD was normally set at 100 times, so the temperatures at the four probes are and ^ o k t f r + ^ i V (F'9) The form of (F.6) through (F.9) reveals that the precise determination of the tem-perature at each probe is dependent not only on the precision of the RTD itself, but on the precision of the 1.5 kft and 100 ft resistors, the amplifier gain GD, and an accurate determination of the lead resistance rn. For measurements of the change in temperature with time, however, the lead resistance r n becomes unimportant, and only the circuit board resistors and the amplifier gain are relevent. In actual fact, the resistances of the supposedly 1.5 kft resistors varied between 1.494 kf2 and 1.498 kft, and the "100 ft" resis-tance was measured at 99.72 ft. These inaccuracies affect the magnitude of the inferred change in temperature slightly, but should not affect the determination of the rates of these changes. Appendix G Hydrate Layer Formation Above a Developing Icy Layer The analytical model of hydrate formation presented in Chapter 3 involved imposing a temperature To lower than the equilibrium hydrate temperature Teh on the boundary of a warm gas-rich porous half-space. If T0 is also below the freezing point of the aqueous pore solution however, this model is no longer valid since ice formation will occur, releasing latent heat and changing the thermal properties of the lower regions of the sediment column. A modified model which includes the development of an ice layer is outlined here. Figure G.l gives a physical description of the formation scenario. The conditions are much the same as those considered in Chapter 3. The initial conditions are uniform with T = Too and c == Coo throughout the half-space. At t — 0 the temperature at z = 0 is dropped to TQ. Here TQ is below the freezing point of the aqueous solution Tej, which is itself below the equilibrium temperature of the hydrate Te/,. The developing hydrate-fluid interface at z = a(t) is in thermodynamic equilibrium at T — Teh and the entire hydrate layer is in chemical equilibrium at c = ce. The interface at z = b(t) between the hydrate layer and the icy layer beneath it is also in thermodynamic equilibrium at the eutectic temperature T = Tei. The ice structure is assumed to reject the gas component so that all of the gas remaining in solution after the passage of the hydrate-fluid interface is incorporated into the hydrate structure to increase the hydrate volume fraction by some amount h*. This is equivalent to assuming that the ice solidus is vertical in the temperature composition phase diagram of the gas-water system. 114 Appendix G. Hydrate Layer Formation Above a Developing Icy Layer 115 High Initial Temperature T = High Initial Gas Concentration c = t z jft z = a ( t ) T = Teh ft Hydrate Layer Hydrate Pore Volume Fraction; h £ z~b( t ) T - T . , ft ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ Ice Pore VolurneFraqtion: ( l - W ) .A.>..;v;v..A.A.A..\,.v...l,.V..\.^ v^ '^ "vl^ '""V'V'V'V>^ ""-v,V'V'V%TA-'^ f. "V. C o l d B o u n d a r y z = 0 T = TQ Figure G.l: Physical Model of Hydrate Formation Above an Icy Layer At t = 0 the boundary is chilled to T0 < Tei < Teh < T^, a hydrate layer develops with its interface position at z'.= a(t) above an icy layer with its interface position at z = b(t). Appendix G. Hydrate Layer Formation Above a Developing Icy Layer 116 The method of solution is essentially identical to that for the model in Chapter 3. Once again the nondimensional temperature T and gas concentration c from (3.1) and (3.2), are used, namely T = — a n d c e Too To Coo ce This problem requires the use of both the hydrate and ice Stefan numbers from (3.3) and (4.3), which were ph<t>Lh PfCf(Too — T0) and • •- ' PfCf(Too - T 0 ) The hydrate gas concentration is defined as in equation (3.4) as Ch — ce Ch = —— , C<x> Ce and the nondimensional length and time variables are again I z z2 t — — , and z = —- , where tr — — . tr %r The nondimensional equation for conservation of energy in each of the three regions is the diffusion equation dt d 2 i 7 a F = d p ' ( a i ) where the constant 7 is defined as ^ = T i = m^.- o<z<i(i) . (G.2) Ke{l) 7 = r , = ^ ^ , b(i) < z < a(i) (G.3) and 7 = r 0 = r(o), z > a(i) (G.4) Appendix G. Hydrate Layer Formation Above a Developing Icy Layer 117 with r(i), T(h), and T(0) defined as the ratio of the heat capacity of the bulk mixture to that of water in the ice and hydrate, fluid and hydrate, and fluid containing regions respectively. The conditions for conservation of gas are as in (3.9) and (3.10), namely c — ce, b(t) < z < d(t) and , dc d2c _ .... , ^acTi = dP •>.'><&>.•• where the product of porosity with the Lewis number (f>a is again the ratio of the thermal and chemical diffusivities. Since it is assumed that no fluid is present beneath the ice interface and all the gas is contained in the hydrate, a condition on the gas concentration beneath z = b(t) would be meaningless. From the problem description come the following initial conditions and/boundary conditions: the initial temperature T(z, 0) = - T f(z,0) = 1;. (G.5] the initial gas concentration c(z, 0) = c(z,0) ='i; (G.6; the hydrate interface is initially at z = 0 o(0) = 0 ; (G.7; the ice interface is initially at z = 0 6(0) = 0; (G.8] the base temperature T(0,i) = T0 f(o,i) = 0; (G.9; the temperature at infinity T(oo,i) = Too r(oo,i )'=i-; ( G . I O ; the gas concentration at infinity c(oo,i). = CQO c(oo, t) = 1 ; (G.ir the hydrate interface temperature T(a, t) = = Teh => f(a,i) = Teh _ Teh ~ ~ T -"To -To ;(G.12; the ice interface temperature T(6, t) -— J-ex T(b,i) = fei T • -~ T --1 oo To To ;(G.i3] and the interface gas concentration c(a, t) = c(a, i) -= 0 . (G .14; Again, these are identical to those presented in (3.13) through (3.20), with (G:8) and Appendix G. Hydrate Layer Formation Above a Developing Icy Layer 118 (G.13) added to describe the temperature and initial position of the ice layer interface. There are interface conditions on the temperature gradient discontinuities which follow from (B.20) at both z — d and z = b, namely df(d+,i) Ke(h)df(d-,i) dd —dz ne dz— = - S T h h T i ' ( G - 1 5 ) — d T - - K e { h ) dz = - ( ^ + W ) ^ , (G.16) where i is the ice pore volume fraction, which is assumed to be (1 — h — h"). Note that the increased hydrate production on freezing necessitates the inclusion of a second term in the factor multiplying the ice layer growth rate in (G.16). The interface condition on the gas concentration gradient discontinuity at z = d is given by (3.22) which was, dc(d+,i) ph_ dd wz = —ch(f)ah—= . . oz pf dt The similarity variable £ from (3.23), is used to convert the partial differential equa-tions into ordinary differential equations, and solutions are sought for the case where the interface positions are given by 1 d(i) = 2\h\ft, (G.17) anc l(i) = 2\i\[i. (G.18) Under the similarity transformation (G.l) becomes ° ' ( a i 9 ) where (G.2) through (G.4) indicate the value 7 takes in each of the three regions. The temperature and gas concentration profiles are found in terms of error functions and complimentary error functions, with the coefficients determined from the initial con-ditions and boundary conditions given above in (G.5) through (G.14). Working through Appendix G. Hydrate Layer Formation Above a Developing Icy Layer 119 the algebra yields the following result for the temperature profile: and where f(t) = A erfCyT^O + B e r i c ^ O , K < (< A* and f r ierfc(yT7Ah-) - f e h erfc(vT^AO erfc(vT7A fc)erf(VTIAO - erfc(yTlA0erf(vT7A f c) ' rr ierf (yT7A.) - f e f cerf(vT^A0 er f (v^A h ) er fc (Vl \A, ) - erf(vTlAOerfc(vTlA f c) ' The gas concentration profile is given by (3.31) and (3.32), namely c(0 = o , o<t<\h A = B -and c(0 = erfCv^O - ^f^\h\^Hy/M) • t > Afc-(G.20) (G.21) (G.22) (G.23) (G.24) erfc(y/c5aA/l) Now substituting these expressions for the temperature and gas concentration pro-files into the interface conditions from energy conservation, (G.15) and (G.16), and the interface condition from gas conservation, (3.22), yields the following three transcenden-tal equations involving the unknown interface velocity parameters A/, and Aj, and the unknown hydrate volume fraction h: A'/fei ^ ~ \l^e-^lK^A ~ B ) =-SThhXh , (G.25) erfc(vToAh) • V , T " «e ' r ° : g - r o A i /I* A — R\ — i / r V - r . A ? Ke(l)Tei 7T V ^ ( ^ ( v T - A i ) -(5 , K(l-A-fe*) + 5'rfcfe*)A ij (G.26) Appendix G. Hydrate Layer Formation Above a Developing Icy Layer 120 and r r rr~\ \ — ch<f>ahXh . (G.27) 7r erfc(vci>Q:Ah) pf Equations (G.23) and (G.24) may be algebraically manipulated to show that eriCvT^AO-eriCvT^AO " K ' ' Equations (G.25), (G.26), and (G.27) are coupled through their dependence on h, Xh, and Aj. It was shown in Chapter 3 that if the chemical diffusivity is small in comparison to the thermal diffusivity (i.e. <f>a —>• co), the hydrate volume fraction may be calculated using (3.48), which was PhCh The amount of additional hydrate which must form at b(t) is found using conservation of gas, which requires that h * = ( i - h ) ^ . . (G.29) PhCh Adopting these values for h and h* leaves (G.25) and (G.26) dependent only on the two interface velocity parameters A^ and A^ . The resulting system is nonlinear, and must be solved iteratively. The strategy is to solve (G.25) with a starting value of Aj = A;i, calculated using the transcendental equation for ice layer formation from (4.1), namely V 7r erfc(v/ToA;i) V if /ce rf(vT7Aii) The root of (G.25) gives A^ = AM, which is used in (G.26) to solve for an updated value of A{ = Ai2, and the process is continued until further iterations no longer alter the values for Xh and Aj. Table G.l contains the parameters used for one such calculation of Aj and A/,. The intermediary values for the growth rates using the above procedure are tabulated in Table G.2. The resulting interface velocity parameters are Aj = 0.171 and A/, = 0.319; Appendix G. Hydrate Layer Formation Above a Developing Icy Layer 121 Parameter Value Units 6.5 x 10 Ke(h) 6.5 x 10 Ke(i) 9.0 x 10 IV 1.23 TH 1.23 T \ 0.74 Srh 1.73 Sn 1.08 Teh 0.385 0.204 h 0.063 h* 0.028 Table G.l: Model Parameters for a Hydrate Layer Above a Developing Icy Layer The physical properties of Table 4.2 were used to calculate the values tabulated above. The Stefan numbers and nondimensional equihbrium temperatures are for Too = 20 °C, Teh = 3.4.°C, Tei = -1.5 °C, and T0 = -7 °C. further iterations produced no changes past the 6th decimal place. Figure G.2 displays the interface positions for the ice and hydrate layers as a function of time. As time progresses, the two interfaces become increasingly separated so it should be relatively easy to distinguish between them. Appendix G. Hydrate Layer Formation Above a Developing Icy Layer 122 Iteration Aj 1 2 3 4 5 6 7 8 9 0.173897 0.171766 0.171043 0.170796 0.170684 0.170674 0.170670 0.170669 0.170669 0.321135 0.319339 0.318729 0.318522 0.318427 0.318419 0.318416 0.318415 Table G.2: Sample Calculation of Growth Rates Further iterations yielded no changes in Aj or Xh past the sixth decimal place. 1 7 5 Z 150 (mm) 1 2 5 h 100 75 50: 25 Hydrate Layer Interface Position . z = a(t) Fluid Saturated Sediments Ice Layer Interface Position z = b(t) Ice Layer 5 0 0 1 0 0 0 Time (min.) 1 5 0 0 2 0 0 0 Figure G.2: Interface Positions bf Ice and Hydrate Layers The interface positions for a developing hydrate layer above an evolving ice layer are shown as a function of time for the model parameters summarized in Table G.l.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Theoretical and experimental investigations into the...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Theoretical and experimental investigations into the formation and accumulation of gas hydrates Rempel, Alan 1994
pdf
Page Metadata
Item Metadata
Title | Theoretical and experimental investigations into the formation and accumulation of gas hydrates |
Creator |
Rempel, Alan |
Date Issued | 1994 |
Description | The substantial volumes of gas hydrates found in the Arctic and in marine sediments are both a possible source of global climate change, and a potential future energy resource. The rate at which a hydrate layer forms, and the spatial distribution of hydrate in the layer are controlled by the physical conditions of the formation environment. To better understand the physical conditions that affect hydrate layer characteristics, I present a quantitative model for the formation of hydrates in a porous medium. The theory is tested using the the results of laboratory simulations of the modelled conditions. Conservation principles are used to derive the full set of governing equations using the minimum number of assumptions and simplifications. Scaling arguments, based on estimates of physical parameters in marine sediments, show that both heat and mass transport are dominated by diffusive processes, so advection may be neglected in most formation environments. Analytical solutions to the leading-order set of equations are obtained for the case of a porous half-space cooled on its boundary. These solutions provide estimates of the growth-rate of a hydrate layer and the volume fraction of hydrate present. The model predicts that the layer grows on the thermal diffusion timescale with the phase-change interface moving at a rate which is proportional to the square root of time. The predicted hydrate volume fraction is determined by the rate at which compositional diffusion can supply gas to the moving interface. For the formation of a methane hydrate layer, the model generally predicts a hydrate volume fraction that is less than 1%. The modelled conditions are simulated in a reaction chamber constructed from a cast acrylic tube, 0.7 m in length, with an inner diameter of 0.14 m. To test the apparatus, experiments were conducted in which the growth-rate of an ice layer in the sand-filled reaction chamber was monitored using RTD temperature probes. These experiments demonstrate that a simplified version of the hydrate formation model describes the formation of an ice layer in a porous medium. CO₂ was used as the hydrate former to test the predictions for the growth-rate of a hydrate layer and the layer's hydrate content. CO₂ has a high solubility in water, and the model predicts a much greater hydrate volume fraction for a CO₂ hydrate layer than that for a low solubility gas such as methane. During hydrate formation experiments, the temperature probes were unsuccessful at detecting the position of the hydrate phase-change interface. At the end of some experiments, the pressure was reduced to dissociate the hydrate so that the presence of hydrate might be detected. However, the temperature probes failed to definitively identify the associated consumption of latent heat. This is despite the recovery of an ice Thydrate mixture following one such experiment. Additional instrumentation employing acoustic or electrical techniques may be necessary to quantitatively assess the model predictions. |
Extent | 6577658 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-01-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.0053009 |
URI | http://hdl.handle.net/2429/3501 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1995-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1995-0061.pdf [ 6.27MB ]
- Metadata
- JSON: 831-1.0053009.json
- JSON-LD: 831-1.0053009-ld.json
- RDF/XML (Pretty): 831-1.0053009-rdf.xml
- RDF/JSON: 831-1.0053009-rdf.json
- Turtle: 831-1.0053009-turtle.txt
- N-Triples: 831-1.0053009-rdf-ntriples.txt
- Original Record: 831-1.0053009-source.json
- Full Text
- 831-1.0053009-fulltext.txt
- Citation
- 831-1.0053009.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-0053009/manifest