FAR-INFRARED SPECTROSCOPY OF C02 CLATHRATE HYDRATE WITH APPLICATIONS TO THE MARTIAN NORTHERN POLAR REGION by Joseph Carl Landry A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in the University of Michigan 1995 Doctoral Committee: Professor Anthony W. England, Chair Professor Sushil K. Atreya Assistant Professor Kamal Sarabandi Professor Fawwaz T. Ulaby

TABLE OF CONTENTS LIST OF FIGURES iv LIST OF TABLES v CHAPTERS 1. INTRODUCTION 1 1.1 Overview 1 1.2 What is a Clathrate Hydrate? 4 1.3 Crystal Structure 7 1.3.1 Structure I Clathrate Hydrate 7 1.3.2 Structure II Clathrate Hydrate 9 1.4 Thermodynamic Properties 11 1.4.1 Phase Diagrams, Stability, and Reaction Rates 11 1.4.2 Thermal Conductivities 14 1.4.3 Latent Heat of Dissociation 16 1.5 Clathrate Hydrates in the Solar System 16 1.5.1 The Gas Giants 16 1.5.2 Titan 20 1.5.3 Mars 21 1.5.4 Earth 23 1.5.5 Comets 24 1.6 Summary 25 2. MECHANISMS FOR INFRARED ABSORPTION 26 2.1 Introduction 26 2.2 Lattice Absorptions 28 2.3 H20 Molecule Absorptions 31 2.3.1 Vibrational Absorptions 32 2.3.2 Librational Absorptions 37 2.4 C02 Guest Molecule Absorptions 38 2.4.1 Rotational Absorptions 38 2.4.2 Vibrational Absorptions 39

2.4.3 Translational Absorptions 42 2.5 Summary 49 3. A SIMPLE MODEL OF THE TRANSLATIONAL ENERGY MODES OF A GUEST MOLECULE IN A CLATHRATE HYDRATE CAGE 50 3.1 Introduction 50 3.2 Potential Model and Wave Equation 51 3.3 Solution of the Eigenvalue Equation and Boundary Conditions 53 3.4 Energy Spectra 57 3.5 Shortcomings of the Model 62 4. FAR-INFRARED MEASUREMENT OF CO2 CLATHRATE HYDRATE 65 4.1 Introduction 65 4.2 Far-infraredMeasurementTechnologies 66 4.3 Bomem DA8 FTS 69 4.4 Cryostat 74 4.5 CO2 Clathrate Hydrate Sample Cell Design 79 4.5.1 Sample Cell Design 79 4.5.2 FrostCell 81 4.5.3 Liquid/Solid Hydrate Cell 90 4.6 Measurement Results and Discussion 97 5. QUASI-CRYSTALLINE APPROXIMATION 102 5.1 Introduction 102 5.2 Dense Medium Scattering Equations 103 5.3 Average Field Approximations 106 5.3.1 Foldy's Approximation 106 5.3.2 Quasicrystalline Approximation 107 5.4 Physical Interpretation of Foldy's Approximation and QCA 110 5.5 Implementation of QCA 115 5.5.1 Incident Field Equations 115 5.5.2 Effective Propagation Constant and Average Fields 118 5.6 Pairwise Distribution Functions 124 5.7 Propagation Constants of Hydrate and Ice Particle Layers 125 * *

6. EMISSION MODEL FOR HYDRATE AND ICE LAYERS 128 6.1 Radiative Transfer and the Quasicrystalline Approximation 128 6.2 Rejection of the Distorted Born Approximation 131 6.3 Effective Single Scattering Albedo 133 6.4 Simplification of the Radiative Transfer Equation 136 6.5 Solution of the Radiative Transfer Equation 139 6.6 Scattering and Absorption of Atmospheric Radiation 147 6.7 Emissivity of Hydrate and Ice Particle Layers 149 7. ATMOSPHERIC TRANSMISSION MODEL 155 7.1 Particles 155 7.1.1 Scattering and Absorption by Small Particles 155 7.1.2 Water Ice Clouds, C02 Ice Clouds, and Dust 158 7.2 Gases 159 7.2.1 Molecules in the Model 160 7.2.2 Lineshape 161 7.2.3 Temperature, Pressure, and Number Density Profiles 163 7.3 Radiative Transfer Model 166 8. FEASIBILITY OF DETECTING MARTIAN C02 HYDRATE 172 8.1 Resolution 172 8.1.1 Earth-based System 172 8.1.2 Mars-based System 173 8.2 Sensitivity 174 9. SUMMARY, CRITIQUE, AND PROPOSALS FOR FUTURE WORK 176 APPENDIX 179 BIBLIOGRAPHY 184 iv

LIST OF FIGURES (1.1) Carbonated Water in a Water Ice Bath 5 (1.2) Cages of a Structure I Clathrate Hydrate 8 (1.3) Cages of a Structure II Clathrate Hydrate 10 (1.4) Phase Diagrams of CO2 Clathrate Hydrate 10 (1.5) Model of the Cloud Structure at Uranus 18 (1.6) Uranus and Neptune Predicted Methane Profile and Saturation Vapor Pressures and CH4 Hydrate Dissociation Pressure 19 (2.1) Mechanisms for Infrared Absorption. 27 (2.2) Lattice Mode Absorptions of Infrared Mulls 30 (2.3) Absorbance of Ethylene Oxide Hydrate and Deuterate at 100 K in the H20 Molecule Vibrational and Librational Absorption Bands 34 (2.4) Absorbance of Water Ice Based on Index of Refraction Data at 213 K in the H20 Vibrational and Librational Absorption Bands. 35 (2.5) H20 Asymmetric Stretching Mode Absorption in Hydrate and Water Ice 36 (2.6) C02 Asymmetric Stretching Mode Absorbance and Bending Mode Absorbance in C02 Hydrate 40 (2.7) C02 Asymmetric Stretching Mode Absorbance in C02 Ice at -80K 41 (2.8) CO2 Gas Line Strengths in the Asymmetric Stretching Absorption Region 43 v

(2.9) C02 Ice Spectra in the Bending Mode Absorption Region 44 (2.10) CO2 Bending Mode Line Strengths in CO2 Gas 45 (2.11) Far-infrared Spectra of CC14, CHC13, and CH2C12 Hydrate with H2S as a Help Gas 47 (3.1) Energy Spectrum of the Translational Modes of the C02 Guest in the Smaller Structure I Cage 58 (3.2) Energy Spectrum of the Translational Modes of the C02 Guest in the Larger Structure I Cage 59 (3.3) Energy Spectrum of the Translational Modes of CCI4, CHC13, 61 CH2C12 Guests in the Larger Structure II Cage (3.4) Sample Energy Eigenfunctions for Infinite and Finite Spherical 63 Square Wells (4.1) Setup for Ultrafast Laser Technique to Measure Far-infrared Dielectric Spectrum 67 (4.2) Bomem DA8 Fourier Transform Spectrometer 70 (4.3) Optical Path of Bomem DA8 Fourier Transform Spectrometer 72 (4.4) Cryostat 75 (4.5) Cryostat Cooling System 76 (4.6) Vapor Deposition Method for Clathrate Hydrate Formation 79 (4.7) Unassembled Frost Sample Cell 83 (4.8) Hydrate Frost Production System 84 (4.9) False C02 Hydrate Absorption Peak in Reflection Spectra 87 (4.10) The Effect of Substrate Wedging at High and Low Frequencies 89 (4.11) Experimental Setup for Observing Hydrate Formation at 0oC 91 (4.12) Exploded View of Far-infrared Cell 94

(4.13) Transmission Spectra of C02 Clathrate Hydrate with Water Ice for 98 Comparison (4.14) Imaginary Index of Refraction of C02 Clathrate Hydrate with Water Ice Index for Comparison 101 (5.1) Scattering and absorption in a dense medium 104 (5.2) Particle and field arrangement for calculation of unscattered field 111 (5.3) Particle and Field Arrangement for Calculation of Singly Scattered Field 112 (5.4) Particle and Field Arrangement for Calculation of Doubly Scattered Field using FA 113 (5.5) Particle and Field Arrangement for Calculation of Doubly Scattered Field using QCA 115 (5.6) Particle and Field Arrangement for Calculation of Triply Scattered Field using QCA 116 (5.7) Effective Complex Indices of Refraction of Dense Media Consisting of C02 Hydrate or Water Ice Particles 127 (6.1) Modeling of a Semi-infinite Layer 129 (6.2) Reciprocal of Emission Problem 132 (6.3) Calculation of Absorption Coefficient 134 (6.4) Emissivity Spectrum of Hydrate and Ice Particles Layers with Average Diameter 100 glm 151 (6.5) Emissivity Difference between 14 cm'1 and 24 cm-' Channels for Hydrate and Ice Layers with Varying Average Particle Diameters 152 (6.6) Emissivity Difference Between 14 cm'1 and 24 cm'1 Channels for Hydrate and Ice Layers with Varying Average Particle Diameters and Coverage Ratios 154 (7.1) Polar Temperature Profile Used in the Atmosphere Propagation and Emission Model 164 vii

(7.2) Pressure Profile Used in the Atmosphere Propagation and Emission Model 166 (7.3) Number Density Profile Used in the Atmosphere Propagation and Emission Model 167 (7.4) Far-infrared Transmission at Zenith for the Martian Northern Polar Region 169 (7.5) Hydrate and Ice Emissions Combined with Absorption and Emission in the Atmosphere 171 (A. 1) Optical and Electronic Configuration of a Fourier Transform Spectrometer 181 (A.2) Beams with Various Delays and Their Sums 182 (A.3) Electronic Processing of the Interferometer Signal 183 viii

LIST OF TABLES Table (1.1) Dissociation Pressures of Clathrates Hydrates at 0~C and Triple Points of Their Guest Species 13 (1.2) Thermal Conductivities of Water Ice and Clathrate Hydrates 15 (1.3) Heats of Dissociation and Sublimation of Hydrates and Ices 23 (3.1! izes, Masses, and Absorption Frequencies of Structure II Hydrate Guests 60 (4.1) Specifications and Components in Our Bomem DA8 System 69 (7.1) Model Mixing Ratios of Absorbing Consituents 161

CHAPTER 1 INTRODUCTION 1.1: Overview This thesis presents the results of research on the feasibility of detecting CO2 clathrate hydrate in the Martian northern polar region. Because the topic of clathrate hydrates is perhaps an esoteric one, we devote almost this entire introductory chapter to a discussion of their structure and thermodynamic properties and their probable and confirmed existence in the solar system. We specifically emphasize C02 hydrate and its importance to energy budget models of Mars. C02 hydrate is made up of a water molecule lattice that contains a regular pattern of voids that can each encage a single C02 guest molecule. This clathrate hydrate structure has an anomalously low thermal conductivity and could significantly affect the vertical flow of heat in the polar regions. CO2 hydrate also has a sizeable latent heat of dissociation that may result in a large thermal inertia that could buffer seasonal temperature variations. Given that the possible presence C02 hydrate on Mars is a worthwhile problem to study, we are motivated to try to detect it in a remote sensing experiment. To do this, we must find some region of its dielectric spectrum that is sufficiently interesting to result in a distinctive radiometric signature. Chapter 2 examines the mechanisms that give rise to features in the infrared dielectric spectra of clathrate hydrates. These mechanisms are compared to any analogous mechanisms that exist in water ice and CO2 ice and gas, which are important background constituents in the polar regions. If the spectra of CO2 hydrate and the background constituents are too similar, the emission signature of the hydrate 1

2 would probably not be sufficiently distinctive to unambiguously identify the hydrate in a remote sensing measurement Rattling mode absorption by the C02 guest molecule is the only hydrate absorption mechanism that does not have an analog in the background constituents. Chapter 2 argues that this is probably the only absorption that will result in an unambiguous signature. The distinctive far-infrared C02 hydrate absorption is due to transitions between translational modes of the guest C02 molecule inside its cage. In Chapter 3 we discuss a simple quantum mechanical model for the translational modes. The cage is modeled as a potential well with a uniform potential inside the well and an infinite potential outside the well, and the C02 guest molecule is modeled as a hard sphere. This crude model predicts that transitions between adjacent translational energy states occur at far-infrared frequencies in the range 10-20 cm-l. We see later in Chapter 4 that the models compares favorably with a measured absorption band at 17 cm-l. The model was also used to predict the frequencies of a few other clathrate hydrate rattling mode absorptions published in the literature. The model is in reasonable agreement with these data. Chapter 4 describes our measurement of the far-infrared dielectric spectrum of C02 hydrate. The measurement was performed with a Fourier transform infrared spectrometer configured with a mercury vapor lamp, a 100 pm Mylar beamsplitter, and a silicon bolometer operating at 4.2 K. The most difficult part of the measurement was designing a clathrate hydrate sample cell in which a well characterized sample could be grown. The cell had to have far-infrared access and be coolable to 150 K, which is the estimated wintertime temperature of the Martian northern polar region. In the final cell design, the hydrate was grown directly from the liquid water state, which also required the cell to be able to withstand several hundred psi of pressure. Once all the difficult design constraints were met, the far-infrared spectrum of CO2 hydrate was measured, and a clear absorption band centered at 17 cml was observed. Water ice and CO2 ice and gas do not exhibit a similar absorption.

3 Once the far-infrared properties of C02 hydrate are known, an emission model can be used to predict the radiobrightness of the hydrate on Mars. The hydrate is expected to be in the form of a layer of packed particles that have sizes comparable to the wavelength of observation. This requires a special model for the propagation of a wave through a dense medium. Classical radiative transfer models are not appropriate for dense media because they greatly overestimate the absorption due to scattering. This overestimation is due to the neglect of coherent scattering effects. The quasi-crystalline approximation (QCA) partially takes these effects into account by including the interdependence of the positions of pairs of particles when calculating the average field in a dense medium. Chapter 5 discusses the theory of QCA in a more intuitive approach than exists in the literature. A Monte Carlo simulation technique is presented to highlight the advantages of QCA over Foldy's approximation, which does not include pairwise interdependence of particle positions. The technique also reveals the conditions under which QCA breaks down. The chapter concludes with a description of the implementation of QCA to find the effective propagation constant of a dense medium and presents examples of hydrate particle and ice particle semi-infinite layers. A Rayleigh distribution is used to describe the particle sizes, and several different average particle size cases are evaluated. Given the propagation properties of the layer, we then develop a radiative transfer model in Chapter 6 to predict its emission. The differential equation of transfer is converted to a system of linear equations by substituting the scattering integral with a sum over a finite number of angles chosen by the Gaussian quadrature integration method. The model predicts a high emissivity for the hydrate layer at frequencies inside the rattling mode absorption band. At frequencies above the band, the emissivity can drop by 0.1 or more, depending on the assumed average size of the particles. The ice layer exhibits a much lower emissivity inside the hydrate absorption band, and for some cases its emissivity increases as the frequency is increased above the band. The hydrate therefore shows a large drop in emissivity, while the ice layer has a flat or slightly increasing

4 emissivity. Such a large difference in emissivity behavior should be easily differentiated by a sensitive detector. The emitted radiation passes through the atmosphere, where it can be attenuated and obscured by atmospheric emission. Chapter 7 discusses the mechanisms or scattering, absorption, and emission in the atmosphere and presents an atmospheric model that is coupled to the emission model. Rayleigh scattering and absorption by the small particles in water ice, C02 ice, and dust clouds is found to be insignificant and is neglected in the formal model. Of the atmospheric gases that absorb in the far-infrared, only CO was found to have an effect on the hydrate and ice layer emission spectra. It exhibits two strong absorptions in the 14-24 cm1l range, but ample windows exist for spectral observations away from these narrow lines. Chapter 7 concludes that the atmosphere will have little effect on the far-infrared surface emission spectrum. Chapter 8 briefly discusses the feasibility of building an instrument to remotely sense C02 hydrate on Mars. An Earth-based measurement is probably not practical due to the enormous telescope aperture that would be necessary to resolve portions of the Martian polar regions. A Mars based system would require an instrument of much more moderate size, and with a sufficiently sensitive detector, integration times of -1 sec should be possible for 1 K temperature resolution. Finally, chapter 9 sums up the research and makes recommendations for future work. The areas of the present research that can be improved and extended are emphasized. 1.2: What is a Clathrate Hydrate? Consider the experiment shown in Figure (1. 1). A bottle of water is highly pressurized with CO2 gas and placed in a water ice bath. The high background pressure reduces the freezing point of the water in the bottle, making it impossible for ice to form at

HAL ~~CO2 hydrate ~:,...~,... =.......... i... I.....t.......,j.:Figure:.....1) Cabnae Wae....Wtr cBt

6 the 00C equilibrium temperature. However, we do see an ice-like solid condense. Since ice is not stable at the temperature and pressure inside the bottle, the solid is not ice. It is instead a compound called a clathrate hydrate. When water freezes, it forms a lattice made up of water molecules. Within the lattice, voids can form which are not normally stable, collapsing due to attractive forces between oxygen molecules. The voids can be stabilized, however, by the introduction of a guest molecule. As long as the background pressure of the guest molecule gas is sufficiently high, this new hydrate structure is stable and will be preferentially formed over water ice. Furthermore, water ice sufficiently pressurized with CO2 gas will convert (although slowly) to C02 clathrate hydrate. If the partial pressure of C02 is lowered below the hydrate dissociation pressure, the clathrate hydrate structure will collapse, leaving liquid water or ice. Clathrate hydrates differ from other hydrates in that they are non-stoichiometric. For example, in ammonia hydrate, which is not a clathrate hydrate, the possible ratio of water molecules to ammonia molecules is 2:1, 1:1, or 1:2, depending on the crystal structure. In any of the three structures, the ratio is a definite number, and the ammonia hydrate is a stoichiometric compound. Ammonia does not form a clathrate hydrate. In a clathrate hydrate, a sufficient fraction of voids must be filled to stabilize the structure, and this fraction can be less than unity. One hundred percent occupancy has never been seen in nature. Instead, some range of smaller values is observed. Although the range may be narrow, no definite ratio of guest molecules to water molecules is required, and the clathrate hydrate is therefore non-stoichiometric. In this chapter we will discuss the general properties and the natural occurrence of clathrate hydrates. We will pay special attention to CO2 hydrate because it is the topic of this thesis, and we will focus on its properties that make it an interesting subject of a remote sensing experiment. However, since other clathrate hydrates are believed to occur in many places in the solar system and could also be interesting subjects of remote sensing

7 experiments, we will seek to maintain a broad outlook with the attitude that the techniques and models developed in this research can also be applied to other clathrate hydrate scenarios. 1.3: Crystal Structures The host lattice of a clathrate hydrate can assume different crystal structures depending on the size of the guest molecule and, to some extent, the method of formation. The'two naturally occurring forms are the structure I and the structure II hydrates. At least one other form, the tetragonal hydrate, has been formed in the laboratory with Br2 as the guest molecule. We will only discuss the naturally occurring forms here, since they are the only forms that will be observed in remote sensing experiments. 13.1: Structure I Clathrate Hydrate The lattice structure of structure I hydrates contains two cage types, which are shown in Figure (1.2). The smaller cage has a pentagonal dodecahedral shape and can accommodate molecules up to about 5.1 Angstroms in diameter, while the larger tetrakaidecahedral cage can fit molecules up to about 5.8 Angstroms in diameter. The general rule for forming clathrate hydrates is that as long as the guest molecules can fit in the cages, it will form a structure I hydrate. Otherwise, a structure II hydrate with larger cages sizes will be formed. The structure I lattice cell is cubic with a side dimension of 12 Angstroms. One cell contains 46 water molecules making up three small cages and two large cages. From the cell dimensions and the water molecule masses, we can calculate the mass density of an empty structure I hydrate, which turns out to be 0.79 g/cm3. Since water ice has a density of about 0.91 g/cm3, the empty hydrate is less dense than water ice, and ice will

ee~~ ~,, ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ eee~~~~~~~~~~~~~~~~~e~ ~ ~ f~..... ~ ee~~~~~~~~~~~~~~~e~ ~ ~ ~ e eee eeeee~~~~~~~~~~~~~~~~~~~~~ ee"~e ~.~. eee~~~~~~~~~~~~~~~~~~~~~~~~ -.::.~nm~~onal dod~~ahedron ~trakaidecahedron~~ Figure (1.2): Cages of a Smucmm I Cla~~~~~~ra~ Hydra~~

9 expand by about 15% when converting hydrate. This knowledge of expansion will be necessary when we calculate slab thicknesses in Chapter 4. Naturally abundant guest molecules that can form structure I hydrates include C02, CH4, C12, and H2S. Ar and K, which can satisfy the size restrictions for structure I hydrates, form structure I hydrates instead. The thermodynamic reasons for this are still not well understood. If all cages in a structure I hydrate were filled, the ratio of water molecules to guest molecules would be 5.75:1. However, the non-stoichiometric clathrate hydrates never achieve full occupancy, and the occupancy percentage depends on the guest molecule. For C02 hydrate, the approximate ratio is 6:1, while for CH4, the ratio is about 7.1:1. When mixed hydrates are formed, where different types of guest molecules occupy both sized cages, the occupancy issue becomes much more complicated. The host lattice structure, however, remains the same as in the case of single guest hydrates. 1.3.2: Structure II Clathrate Hydrate The structure II clathrate hydrate also has two cage types. The smaller of these has a pentagonal dodecahedral structure and is nearly identical to the smaller cage of the structure I hydrates. As in the structure I case, it can fit molecules up to about 4.7 Angstroms in diameter. The larger cage has a hexakaidecahedra structure that is shown in Figure (1.3). It can accommodate molecules with diameters up to 6.7 Angstroms. Because of the large disparity in the cage sizes of structure II hydrates, few single species structure II hydrates have been encountered. If a molecule is small enough to fit into the smaller structure II cage, it is small enough to fit in both structure I cages, and it is likely that it would form a structure I hydrate instead of a structure II hydrate. With a few exceptions, some of which were noted in the previous subsection, structure II hydrates are usually double hydrates. In a double hydrate, one molecule type occupies the smaller

10 t \~ pentagonal dodecahedron ~~~~~~~~. \~~~~~~~~~~~~~~~~~. i;~~~.~,~, i'... \. ~~~hcxakaidccahcdron ~~~~~~~~~~~~eeet ~e?~,rr' Figur (1.): Caes o a StU~t~e II lathate Hdrat...... /~ Fiur (13'Cgso tutue1 ltrt yrt

11 cages while another type occupies the larger cages. For example, Gerbaux et al. have formed a double structure II hydrate with the smaller cages occupied by the 4.1 Angstrom H2S molecule and the larger cage occupied by the 6.4 Angstrom CHC13 molecule. Besides Ar and K, which are known to anomalously form single species structure II hydrates, N2 and 02 are believed to form structure II hydrates as well. Their natural occurrence on earth, as well as the natural occurrence of other hydrates in the solar system, will be discussed in the next section. 1.4: Thermodynamic Properties In this section we discuss several important thermodynamic properties of hydrates. These properties make them significantly different from water ice or the guest molecule ice. The differences may have important effects on the energy balance models of their planetary environments. 1.4.1: Phase Diagrams, Stability, and Reaction Rates The long term stability of a clathrate hydrate can be expressed through its phase diagram. Two phase diagrams of C02 hydrate showing different regions of pressuretemperature space are given in Figure (1.4). The diagrams show that a sufficiently high pressure of CO2 gas can keep the hydrate more stable than water ice over the entire temperature range shown. More interestingly, at temperatures below 121 K, the extrapolation curves by Miller and Smythe (1970) show that CO2 ice becomes more stable than C02 hydrate. Therefore, a hydrate sample cooled to below 121 K would (slowly) convert to C02 ice. Long term experiments performed by Davidson et al. (1984) add support to this idea. If it were true, C02 would be the only known clathrate hydrate to have its guest molecule form a solid that is more stable than the hydrate itself. More

I/TI 0/ loOC\ to To ~ o 5t 1// 110 3 41 40 39 38 3I 36 35 34 33 I, l ao, 1 8;.1 1i...." - 1. P..... 80 CP 60F u40 C Hydrole D4 C02(1) 40 _ i-, JO I WOle, 4 CO2 (91 8 WeHydate C(02ll) Ice C02 19).1 14 1 i I,, [MIPRATUR[ (a): Low Temperature Phase Behavior (b): High Temperature Phase Behavior Figure (1.4): Phase Diagrams of CO2 Clathrate Hydrate (reproduced from Miller, 1985)

13 importantly, it indicates that CO2 hydrate cannot be formed below 121 K simply by pressurizing water ice with CO2 gas. We will use this fact when discussing experimental procedures for growing C02 hydrate in Chapter 4. Phase diagrams can tell us when hydrates will be stable, but they cannot directly tell us how fast a hydrate will form or dissociate. We can get an idea of how readily hydrates form on a relative scale based on their 00C dissociation pressures. In general, a hydrate will form more readily if it has a lower dissociation pressure. Triple Point Guest Molecule Dissociation pressure (atm) TemperaPn(K) Temperature (K) H2S 0.92 188 SO2 0.39 200 CHC13 0.065 210 CO2 12.5 195 CH4 26.0 90 N2 160 63 02 120 55 Ar 95.5 84 K 14.5 117 Table (1.1): Dissociation Pressures of Clathrate Hydrates at 0~C and Triple Points of Their Guest Species Table (1.1) shows the dissociation pressures for several structure I and structure II hydrates. We see that the polar molecules (H2S, S02, and CHCl3) tend to have lower dissociation pressures than the non-polar molecules (CO2, CH4, N2, 02, Ar, and K). Experimenters such as Fleyfel and Devlin (1988) and Kiefte et al. (1985) have found that polar molecules are much more easily enclathrated than non-polar molecules. We will not try to explain this phenomenon, but will simply say that since CO2 is a non-polar

14 molecule, it will probably be more difficult to enclathrate. This is confirmed in Chapter 3, and in fact posed the major experimental challenge of this research. Table (1.1) also gives the triple point temperatures for the guest species of each clathrate hydrate. At these low temperatures, ice formation competes with hydrate formation, though hydrates are generally more stable than the ices if water is available in the background. When we discuss possible scenarios for clathrate hydrate formation in the solar system, it will be important to consider the existence of ices in the same environments. 1.4.2: Thermal Conductivities While most of the thermodynamic properties of clathrate hydrates are reasonably well understood, the most striking anomaly is the strange behavior of the thermal conductivity. The hydrogen-bonded lattice of hydrates is similar to that of water ice, and once the difference in the densities is taken into account, the hydrate and ice lattices have many similar properties. For example, they have similar heat capacities, sound velocities, and infrared absorption. However, their thermal conductivities are radically different. Table (1.2) shows the thermal conductivities of a few hydrates above 200 K along with that of hexagonal ice. We see that the thermal conductivities of the hydrates are about a factor of 5 less than that of ice. This disparity increases with decreasing temperature. The thermal conductivity of liquid water increases slowly with decreasing temperature, while the conductivity of hydrates decreases. At about 40 K, the thermal conductivity of hydrates is only about 2% of that of ice. Thermal conductivity depends on the properties of lattice wave propagation, with scattering causing the conductivity to decrease. Since hydrates contain guest molecules and ice does not, it is natural to suspect that guest molecules are involved in the

15 conductivity decrease. Coupling between lattice waves and the motion of the guests has been suggested as a mechanism, for increasing scattering. thermal conductivity compound (structure) temperature (K) t ml K1) water ice (Ih) 260 2.35 xenon hydrate (I) 245 0.36 methane hydrate (I) 213 0.45 ethylene oxide hydrate (I) 263 0.49 cyclobutane hydrate (II) 260 0.47 tetrahydrofuran hydrate (II) 260 0.51 1,3-dioxolane hydrate (II) 260 0.51 Table (1.2): Thermal Conductivities of Water Ice and Clathrate Hydrates The thermal conductivity of clathrate hydrates has been studied by many investigators. Ross and Anderson (1982) have measured the thermal conductivity of several hydrates with polyatomic guests in the temperature range 100-260K. They attributed the decrease in thermal conductivity to the coupling of lattice waves to rotational modes of the guest. Later, however, Handa and Cook (1987) measured the thermal conductivity of Xe hydrate, which also showed the decrease found for the polyatomic case. Since Xe guest atoms have no rotational modes, rotational coupling cannot account for the decrease. Handa and Cook have suggested that phonon coupling to translational modes of the guest may increase phonon scattering. Tse et al. (1983, 1984) have performed lattice dynamics simulations of methane hydrate, which indicate that translational modes have a small effect on the phonon density of states, but they propose that this may be sufficient to radically change the phonon scattering mechanism and cause a substantial change in thermal conductivity.

16 Whatever the cause of the decrease in thermal conductivity in clathrate hydrates, the decrease could have a strong effect on the energy balance of a planetary environment if the hydrate is sufficiently abundant. We will return to this concept when we discuss the Martian polar regions in section 1.5.3. 1.4.3: Latent Heat of Dissociation Forming a clathrate hydrate in a water + guest gas system requires the removal of energy from the system, and dissociating the hydrate necessitates the addition of energy. In the case of C02 hydrate, 5.86 kilocalories are required to dissociate one mole of hydrate. We will see that this may be important in some planetary environments when we look at the Martian problem in section 1.5.3. 1.5: Clathrate Hydrates in the Solar System Clathrate hydrates can form wherever the guest molecule partial pressure is high enough and the temperature is low enough for the hydrate to be stable. Many places in the solar system have conditions favorable for their formation. Miller, who has been the major proponent of the study in planetary systems, has given a few reviews of their possible locations (Miller, 1961, 1985). We will briefly discuss the more likely scenarios. 1.5.1: The Gas Giants Uranus and Neptune are colder and more dense than the other gas giants in the solar system. For example, at an altitude corresponding to 1 bar of pressure, Jupiter and Saturn are believed to have ambient temperatures above 200 K, while Uranus and Neptune have temperatures below 150 K. If we were to assume equivalent mixing ratios for non

17 saturated constituents, then hydrates would be more likely to form on Uranus and Neptune because the guest species would have a higher partial pressure than on Jupiter or Saturn. The mixing ratio of the constituents is not the same for all of the gas giants. In particular, CH4 on Uranus and Neptune is believed to be enhanced to several times the Jovian and Saturnian ratios. Atreya and Romani (1985) developed models for the clouds of Uranus and Neptune, and the results of their Uranus model are shown in Figure (1.5). The model predicts a CH4 ice cloud layer in the vicinity of the 1.4 bar level at a temperature of about 90 K, as was also inferred from a radio occultation experiment on Voyager 2 (Tyler et al., 1986). (90 K is near the triple point of CH4, as given in Table (1.1).) The exact location and density of this cloud depends on the assumed CH4 mixing ratio. The saturation vapor pressure at the bottom of the cloud is -0.1 bars and equals the CH4 partial pressure at that level. By extrapolating hydrate phase data above 100 K down to the temperature of the clouds, Miller concludes that the CH4 dissociation pressure at these low temperatures is a few orders of magnitude lower than the saturation vapor pressure of CH4 ice. Thus, CH4 hydrate would be more stable than the CH4 ice clouds. He further states that CH4 at lower altitudes (and higher temperatures and partial pressures) must form CH4 hydrate. This argument may have a problem because as one goes deeper into the atmosphere of Uranus or Neptune, the temperature rises quickly enough so that the dissociation pressure soon exceeds the CH4 partial pressure. Furthermore, according to the cloud model, the atmosphere is largely free of water vapor until we get down to altitudes corresponding to temperatures above 200 K. At higher altitudes the amount of water in the atmosphere is exceedingly low because much of the water has frozen out in lower altitude water ice clouds. An interesting question is whether or not the temperature of the water ice clouds is low enough for them to be more stable as CH4 hydrate clouds. Figure (1.6) shows the partial pressure of CH4 as a function of temperature in the Uranus or Neptune atmosphere

19 4 1 I' I'' i'.______ Uranus and Neptune CH4 partial pressure. —---- CH4 hydrate dissociation pressure 3 2 C.. 100 5 150 75 200 2.. I 100 125 150 175 200 225 250 temperature (K) Figure (1.6): Uranus and Neptune Predicted Methane Pressure Profile (Atreya, 1986) and CH4 Hydrate Dissociation Pressure. The methane pressure exceeds the hydrate dissociation pressure at temperatures below about 205 K, thereby allowing a stable hydrate. However, most of the water will have frozen out at lower altitudes.

20 using the cloud model. A generous CH4 mixing ratio of 0.05 is assumed, which is 60 times the solar ratio. (Neptune is expected to have a CH4 mixing ratio of 30-75 times solar, while the preferred value for Uranus is 25-30 times solar (Atreya, 1986).) The figure also shows the dissociation pressure of CH4 hydrate as a function of temperature. According to this model, the hydrate should be stable down to an altitude corresponding to about 205 K. The water ice cloud layer shown in Figure (1.5) extends to the 205 K level, but it is very tenuous at that point because much of the water has frozen out at lower altitudes. If this is the case, then very little water will be in the CH4 hydrate form on Uranus and Neptune. If a more conservative estimate of the CH4 mixing ratio is used, then the amount of hydrate can be much lower. While uncertainties in the cloud models allow for hydrates to be present, the presence of large clouds of hydrates on Uranus and Neptune is by no means a certainty. Because Jupiter and Saturn have a smaller methane and have higher temperatures for a given pressure level, CH4 hydrate clouds do not form on these gas giants. 1.5.2: Titan Titan, Saturn's largest moon, has an atmosphere primarily of N2, with CH4 and probably Ar as additional major constituents. Miller (1985) proposes a mixed hydrate of N2, Ar and CH4, and Owen (1982) uses such a model to explain the lack of Ne in the Titan atmosphere. The temperature of Titan's surface (-80 K) is sufficiently low to stabilize a N2 clathrate hydrate given the currently accepted value of the N2 partial pressure (-1 bar). However, the rate of hydrate formation at this temperature is exceedingly slow. Fishermen, if ice on Titan is in the form of thick layers, hydrate formation is greatly impeded. If large deposits of clathrate hydrate exist on Titan, they may have formed in an earlier age of the moon's history when thermodynamic conditions were markedly

21 different, or during the creation of the solar system. Atreya (1986) has discussed the merits and the problems of the Titan clathrate hydrate scheme. 1.5.3: Mars The Martian northern polar region undergoes profound seasonal changes (James et al., 1992). In the winter, the temperature drops to about 150 K, lowering the saturation vapor pressure of CO2 to below the summertime partial pressure. As a result, large amounts of CO2 condense out of the atmosphere, and the polar region is blanketed with a layer of CO2 ice. In the spring and summer, the CO2 ice layer retreats to expose a residual cap of water ice. A similar cycle occurs in the southern hemisphere. The south pole is covered by a layer of CO2 ice in the winter, but at least in some years, the seasonal cap does not completely retreat, and the residual cap is C02 ice instead of water ice. The south pole is colder than the north pole. A water ice cap may lie beneath the CO2 cap, but this cannot be determined directly from a remote sensing measurement. From the phase diagrams in Figure (1.4), we see that the saturation vapor pressure of CO2 is above the hydrate dissociation pressure over the entire range of the measurements. Therefore, when CO2 ice is forming in the Martian northern polar regions, the water ice in the residual cap should be converting to C02 hydrate because the hydrate is more stable that the CO2 ice. It is thus very likely that CO2 hydrate exists on Mars. The low temperature phase diagram of Figure (1.4) was presented in a paper by Miller and Smythe (1970). In the same paper they also presented reaction rates for the conversion of finely divided ice samples into C02 hydrate. They found that at Martian northern polar winter temperatures, the ice samples became completely converted to hydrate within a few hours. The water ice on Mars is also expected to be finely divided, as

22 will be discussed in Chapter 7. A substantial amount of it is therefore likely to become hydrate in the winter. Since Miller and Smythe proposed the existence of Martian clathrate hydrates in 1970, there has been some speculation on the effect their presence would have on the atmosphere and the CO2 cycle. At the time, temperature ranges in the polar regions were not well known. Dobrovolskis and Ingersoll (1975) proposed that the north polar cap could contain a permanent deposit of hydrate that serves as a reservoir of C02, stabilizing the ambient pressure against large changes. Since then it was determined that the summertime temperatures in the northern polar are higher than previously suspected and are much too high for C02 hydrate to remain stable. It was also discovered that the south pole probably has a permanent C02 cap, which has a much greater buffering capability than a hydrate deposit. As a result, the hydrate reservoir concept has fallen out of favor. The presence of C02 hydrate still has importance for Mars. If a large deposit of hydrate temporarily exists in the northern polar region, it could have a profound influence on the energy budget. Clathrate hydrates have an anomalously low thermal conductivity, as discussed in section 1.2.3. A hydrate deposit could therefore be a more effective insulator than a similar deposit of water ice. This might result in a much slower transmission of heat to and from the lower strata and significantly change polar heat flow models. Clathrate hydrate in the polar region can also affect the polar energy budget through the energy it emits as it forms and absorbs as it dissociates. The latent heat of sublimation is shown in table (1.3) along with the latent heats of sublimation of water ice and CO2 ice. We see that a mole of dissociating C02 hydrate consumes about the same amount of heat that a sublimating mole of CO2 ice consumes and about half the heat that a sublimating mole of water ice consumes. Since the sublimation of water and CO2 ices are important sources and sinks of energy in polar energy balance models, a large deposit of CO2 hydrate might also be important.

23 molar heat of sublimation or dissociation(kcal/mole) CO2 hydrate 5.86 water ice 12.20 C02 ice 6.3 Table 1.3: Heats of Dissociation and Sublimation of Hydrates and Ices The current Martian climate could be undergoing episodic changes in its geological history. The present surface and atmospheric thermal structure of the planet are by no means what the planet has had since its last known vulcanism. Since the mid-1970's the average temperature on the planet has dropped 20 K, largely due to the removal of dust particles from the atmosphere (Lee et al., 1995). Studies of the energy budget of the planet and its polar regions are especially interesting in light of rapidly changing Martian climate. 1.5.4: Earth Temperatures at the surface of the earth and in the atmosphere are much too high to stabilize clathrate hydrates given the ambient partial pressure of possible guest molecule gases. Beneath the surface, however, the pressures are high enough to compress gas to form clathrate hydrates. A mixed hydrate of N2 and 02, the so-called air hydrate, has been found in antarctic ice cores (Gow et al., 1968; Miller, 1969), where air bubbles were highly compressed over time. CH4 hydrates have been found in permafrost in Siberia and later in North America, as well as in ocean sediments (Miller, 1985).

24 1.5.5: Comets Comets are believed to be icy conglomerates of mainly water ice with several important minor constituents, including CO, HCN, NH3, and CH4 (Whipple, 1950). Spectra from the comae and tails of comets have revealed the presence of some of these minor constituents or of radicals that are probably their by-products. Since comets are thought to have formed in the primordial gas cloud at the same time the solar system did, their composition may hold important information on the early solar system's makeup and creation. In their highly elliptical orbits around the sun, comets undergo large temperature changes. At aphelion, a comet's temperature is on the order of 10 K, while at perihelion it can rise to above 200 K. Some of the highly volatile constituents of the comet, such as CH4, are thought to be too unstable to survive a close pass by the sun, yet they or their byproducts are observed in cometary gas emissions. The existence of CH4 and other constituents as clathrate hydrates in comets was thus suggested by Delsemme and Swings(1952) to explain their stability. Furthermore, the relative amounts of gas or gas by-by-products of each volatile should be related to the stability of the volatile. Since the relative abundance of the constituent gases and by-products is fairly constant with radial distance from the sun, Delsemme and Swings also proposed that this equalization is the result of mixed clathrate hydrates. Since 1952, the necessity of clathrate hydrates in comets has been largely downplayed. Whipple has pointed out that cometary emission phenomena can be explained by a comet nucleus with an outer surface of water ice and minor constituents trapped inside. The sublimation is then controlled by the water ice, and minor constituents are freed as layers of water ice are removed. While the observed gas emissions are consistent with this ice model, the model does not support or defy the existence of internal clathrate hydrates.

25 One problem with a clathrate hydrate composition in comets is that clathrate hydrates may not be able to form in the cold regions where comets are thought to originate. This problem has been studied by Blake et al. (1991), who grew H20 and CH30H mixtures at very low temperatures (-4 K) and observed the phase changes as the samples were heated. Using transmission electron microscopy and selected area electron diffraction, they verified that upon warming, their initially amorphous ices could rearrange themselves to form clathrate hydrate. Therefore, a comet composed primarily of water ice could partially convert to hydrate as it approaches the sun. 1.6: Summary Clathrate hydrates have several interesting thermodynamic properties. They are more stable than water ice or the guest molecule ice if the background pressure is sufficiently high. As a result, they can occur in many places in the solar system where water and the guest constituent are present. Their large latent heat of sublimation and low thermal conductivity make them important in energy budget models, particularly in situations where the temperature or pressure changes dramatically enough for them to change phase. Outside the earth, their presence is especially likely in the Martian northern polar region in winter and spring. CO2 hydrate has never been spectroscopically confirmned there, but instead has been suggested through indirect thermodynamic arguments. In order to directly verify its presence and follow its seasonal evolution through a remote sensing experiment, its dielectric spectra must be known. In the next chapter, we examine the dielectric properties of CO2 hydrate and determine which frequency range is most suitable for a remote measurement

CHAPTER 2 MECHANISMS FOR INFRARED ABSORPTION 2.1: Introduction To determine the feasibility of remotely detecting C02 clathrate hydrate on Mars, we need a model for its emission. Such an emission model would require knowledge of the dielectric properties of CO2 hydrate. The dielectric spectrum of the hydrate must be significantly different from the spectra of other environmental constituents to result in a distinctive radiometric signature. The constituents that are likely to have similar spectra to CO2 hydrate are water ice and C02 ice and gas. They each share some structural property with the hydrate that may result in nearly the same energy states and transition frequencies. Furthermore, they are each abundant in the Martian polar regions and are likely to have a strong effect on its emission spectrum. In this chapter we will examine the mechanisms that give rise to infrared absorption in C02 hydrate. We will compare the resulting spectra to the analogous spectra in water ice and CO2 ice and gas to determine if the differences might be sufficiently large to produce distinctive hydrate signatures. The mechanisms are shown pictorially in Figure (2.1) and include lattice vibrations, individual H20 molecule motions, and guest molecule motions. The different mechanisms are not completely uncoupled, but the coupling is small enough to assign a dominant mechanism to each absorption. 26

27 (a): Host Lattice Vibrations vibration libration D..../.... /.00. (b): Individual H20 molecule motion vibration translation rotation (c): CO2 Guest Molecule Motion Figure (2.1): Mechanisms for infrared absorption. (Simplified 2-D lattice structures and cages are not accurate and are for illustration only.)

28 2.2: Lattice Absorptions Ices belong to a family of solids called molecular crystals (Califano, 1981). The lattice structure of these solids is characterized by strong intramolecular bonds and relatively weak intermolecular bonds. The molecules in such a solid maintain their individual identities while still being part of a large lattice. A water lattice is held together by weak O-H bonds between adjacent H20 molecules, as illustrated by the dotted lines in Figure (2. lb). The hydrogen atom involved in the weak bond is more tightly bound to the oxygen atom in its own molecule than to the oxygen in the neighboring molecule. The stronger bonds are illustrated by the solid lines in Figure (2. ib). The wavefunction of the hydrogen atom is localized near the oxygen to which it is more strongly bound, making the strong and weak O-H bonds asymmetric. Because of the structure of the ice crystal, the asymmetric bonds need not form a regular pattern in the lattice. The ice crystal is thus a disordered crystal. The orientations of water molecules can be made to change through the application of an electric field, but the process is very slow. The Debye resonance frequency resulting from these reorientations in water ice is in the kilohertz range. If a hydrate has a polar guest species, then additional Debye relaxation features are present in the dielectric spectrum. Acetone hydrate, for example, has guest molecule Debye absorption in the range 37-140 GHz (Davidson, 1973, p. 209), which is an ideal range for remote sensing radio technologies. Unfortunately, C02 is not a polar molecule, and if it follows the pattern of other non-polar guest clathrate hydrates, its hydrate has a dielectric relaxation spectrum similar to that of water ice. A shift in its low frequency spectrum might be detectable in a laboratory, but the kilohertz range is not suitable for a remote sensing measurement because of the exceedingly large antennas that are required to attain a reasonably fine resolution.

29 A lattice can also absorb electromagnetic waves through coupling to higher energy lattice phonons. The disordered nature of ice crystals leads to a relaxing or elimination of the selection rules that normally apply to ordered crystals. This results in rather broad ice lattice absorption spectra, absorption bands over 100 cm'1 wide. Clathrate hydrates also have a disordered structure, and we therefore expect their lattice transition spectra to also be broad. Since the structure of hydrates is different than that of hexagonal ice, we might expect that the lattice spectrum would occur over a different range and have a different shape. However, because of the large-scale disorder, lattice mode transition frequencies are largely governed by small scale order, and the lattice absorption band has the same range and general shape as that of hexagonal or cubic water ice. The lattice absorption spectrum of water ice and ethylene oxide (C2H40) structure I clathrate hydrate was measured by Bertie and Othen (1973) at 100 K, and ethylene oxide hydrate and other hydrates were later measured by Bertie and Jacobs (1977) down to 4.3 K. We will consider the higher temperature measurements here because they are closer to the 150 K temperature estimated for the Martian northern polar region. The samples were prepared by grinding ice or hydrate into a powder, verifying its structure through X-ray diffraction, and mulling the powder with liquid propane or Freon-13. The results of the higher temperature measurement are shown in Figure (2.2). The broader features of the hydrate spectrum were attributed to the irregular hydrate structure 0-0 nearest neighbor separations of 2.74, 2.75, 2.79, and 2.92 Angstroms, while ice only has a fixed separation of 2.75 Angstroms. The two spectra are said to be similar because the majority of the hydrate 0-0 separations have lengths of 2.75 Angstroms. The breadth of the hydrate spectrum might also be due to the quality of the mull. Bertie and Jacobs (1977) also prepared samples from mulls, though with a different mulling agent. They found that less finely ground samples resulted in broader spectra in at

30 2.0 A..... - 40 30 o 0 v/ / cm1 1.5 - - 3.0 0350 300 250 200 15 0 100 50 Figure (2.2): Lattice Mode Absorptions of Infrared Mulls (reproduced from Bertie and Othen, 1972). The solid lines are hydrate (top figure) and deuterate (bottom figure) of ethylene oxide. The dotted line is undeuterated water ice. Letters identify different mulls, with absorption scales indicated with arrows.

31 least the cyclopropane hydrate case. It is difficult to generalize this to the Bertie and Othen (1973) data, since most of the 1977 data were reported at 4.3 K. The sharp peak in the ice and hydrate spectra occur at 229.2 cm-l in both cases. Its position is therefore not a parameter suitable for distinguishing between ice and hydrate. The broadening of the peak could be used, but its effect might be too subtle to avoid confusion with ice emission, and the hydrate peak is likely to be broader at the higher Martian temperatures. The dip in the ice absorption at 180 cm'l and the shoulder in the hydrate absorption at 200 cm-l might also be considered, but they are relatively weak and become less pronounced at the higher 150 K temperature expected in the Martian polar region. They are not ideal because they occur in the same absorption band. A much better feature would be a band that occurs in one case and not in the other. The lattice absorptions of structure I clathrate hydrates apparently will not provide us with such a feature. The hydrate spectra exhibit an absorption plateau below 100 cml, while the water ice absorption continues to decrease with frequency. This plateau is not due to lattice absorptions, but is instead a result of guest molecule motions. The plateau does not exist in C02 hydrate, and the reason relates to the C02 molecule being nonpolar. The ethylene oxide absorption will be discussed when we consider guest molecule absorption in section 2.4. 2.3: H20 Molecule Absorptions Since the intramolecular forces in an ice or hydrate lattice are much stronger than the intermolecular forces, the intramolecular motion of the H20 molecule gives rise to its own absorptions that are distinct from the lattice absorptions. The mechanisms for H20 absorption are vibrations and librations.l Intermolecular forces are so weak that the vibrational absorptions of the H20 molecule in the lattice are similar to those of an H20

32 molecule in a gas. The librational motion is somewhat analogous to the rotation of a free molecule, but the absorption is different because the molecule in a lattice is subject to a restoring force due to intermolecular bonds. The free molecule is under no such influence. 2.3.1: Vibrational Absorptions If the modes of a particular water molecule vibration (for example, the symmetric bending mode v2) can be modeled with a purely harmonic potential, then the energy levels of the modes will be equally spaced, and the only purely vibrational dipole transitions that are allowed are between adjacent energy levels. (We assume here a mode where the dipole moment of the molecule changes while it vibrates. Otherwise, no transitions at all would be allowed.) The purely vibrational spectrum is thus a single line at the characteristic frequency of the oscillator. If the molecule is in a gas, it is also free to rotate, and transitions simultaneously involving vibrations and rotations are allowed. The resulting spectrum then has many spectral lines in addition to a purely vibrational transition line. The purely vibrational line is known as the Q branch, while the lower frequency and higher frequency vibrational rotational lines are known as the P and R branches, respectively. In an ice or hydrate, the H20 molecule is not free to rotate. However, it will not exhibit a sharp spectral line because it can couple to its neighbor molecules and to the lattice as a whole. The position of an ice or hydrate vibrational band is mainly determined by the frequency of the intramolecular vibration, and its width is largely determined by coupling to the lattice directly or through librations. Because of the weakness of the intermolecular bonds in both hydrates and ice, their H20 vibrational bands occur in essentially the same places and have similar shapes. This is shown in the hydrate data 1. A libration is a constrained rotation. The entire molecule twists in place while its intramolecular bond lengths remain fixed. The motion is illustrated in Figure (2.lb).

33 measured by Bertie and Othen (1973) at 100 K and the water ice data compiled by Warren (1984) at generally higher temperatures (> 200 K). Their data are shown in Figures (2.3) and (2.4), respectively. The sharp features in the hydrate date are either due to the ethylene oxide guest molecule or the Freon-13 mulling agent used to prepare the samples. Bertie and Othen state that the hydrate water molecule features are broader and less finely resolved than water ice features. They attributed this to the multiple 0-0 separations present in the hydrate structure. The broadening might be used in a radiometric signature of hydrate, but it is likely to be too subtle to convincingly confirm its existence. In an effort to find distinctive spectra of solar system hydrates, Smythe (1975) measured the infrared reflection spectra of CH4 and C02 hydrate frosts in the 1-6 gm range and compared them with similar spectra of water ice frost. The only change he found in the H20 molecular absorption was perhaps a slight broadening of its 3 pm feature due to stretching modes. This, and the lack of strong guest molecule features, led Smythe to make disparaging comments on the possibility of detecting clathrate hydrates spectroscopically. Richardson et al. (1985) found large increases in the bandwidths of vibrational features of isolated HDO and D20 molecules in clathrate hydrate lattices when compared to similar spectra of the same molecules isolated in cubic ice. However, since the water on Mars is not expected to be highly enriched with deuterium, this difference is not important to this work. Apparently they did not find similar broadening of H20 molecular vibrational features because none were reported. They did, however, report a change in the shape of the V2 absorption feature of H20 as their sample went from a mixture of amorphous H20 ice and ethylene oxide to ethylene oxide clathrate hydrate. The figure from their paper is reproduced in Figure (2.5a). Unfortunately, the change in the shape of the absorption is not large. Furthermore, it turns out to be very close to the v2 absorption band of hexagonal ice at -7 C (Warren, 1984), which is given in Figure (2.4b). The high

34 0 60 / 40 20 2 Z Z 60 a 3 2020- 3600 2800 2000 1600 1200 800 400,/cm'l Figure (2.3): Absorbance of Ethylene Oxide Hydrate (Top) and Deuterate (Bottom) at 100 K in the H20 Molecule Vibrational and Librational Absorption Bands (reproduced from Bertie and Othen, 1973). Samples were prepared in an infrared mull.

0.() (OX 4(XX) 36M(K 32Mn 2 8M 24MK 2MnK) AMn 18Mn 16MW 14MK 12M HXX(~) 9M NX6() 4M, c wavenumberbr s (cm-I' waveiiumlbers (cml') Figure (2.4): Absorbance of Water Ice Based on Index of Refraction Data at 213 K (Warren, 1984) in the 1i20 Vibrational adLba tional Absorption Bands.

1300 1270 1240 C I _ _ _ _ _ _ _ _ _ _ _ |________ _ _1800 1600 14(X) 12() 1650 1450 1250.1 I cm wavenumbers (cm- ) a) Ethylene Oxide Hydrate Absorbance (Reproduced (b) Water Ice Absorbancc at 213 K (Warfrom Richardson, 1988). Top curve is a mixture of ren, 1984) amorphous water ice and ethylene oxide. Bottom curve is ethylene oxide clathrate hydrate. Figure (2.5): H20 Asymmetric Stretching Mode Absorption in Hydrate and Water Ice

37 temperature hexagonal ice data has a broad peak a few hundred wavenumbers higher than that of the 10 K C02 hydrate. At lower temperatures, Tsujimoto et al. (1982) found that the peak position shifts to lower frequency. While we could find no ice data at 10 K to compare with the Richardson et al. data, it is likely that the peaks' positions are close. In any case, the broad shapes and ranges are similar enough for the v2 band to be unlikely to give rise to a sufficiently distinct emission spectra. 2.3.2: Librational Absorption Librational absorption data for clathrate hydrates are much scarcer in the literature than vibrational data. This is perhaps because the librational motion of a molecule does not have a direct analog in the free gas state and, hence, is not directly comparable. Bertie and Othen's (1973) spectra of ethylene oxide hydrate given in Figure (2.3a) includes the librational modes in the 800 cm1l region, and it is virtually identical to the librational data compiled by Warren (1984) shown in Figure (2.4b). The similarity is somewhat surprising because the restoring forces that govern the librational modes should depend on the 0-0 distances, which, as previously stated, are fixed at 2.75 Angstroms in water ice but can take on four possible values in a structure I hydrate. Harvey et al. (1964) have published spectral data of SO2 hydrate including the librational region. They claimed a large frequency shift of the librational peak, which they attributed to the large dipole moment of the S02 guest. If the shift were actually due to coupling to the dipole moment of the guest, then a similar shift would not be present in the CO2 hydrate because CO2 is a non-polar molecule. In any case, the group retracted their data saying that their samples were not actually clathrate hydrate (Hardin and Harvey, 1971).

38 More careful study of the librational modes seems worthwhile. The present data suggest that the librational absorptions are unlikely to produce a distinctive emission signature. 2.4: CO2 Guest Molecule Absorptions The guest molecule in a clathrate hydrate cage is not bound to the host lattice and has some freedom to rotate, vibrate, and translate within its cage. We will discuss these three mechanisms separately to determine their likelihood of producing distinctive hydrate spectra. 2.4.1: Rotational Absorptions Bertie and Jacobs (1977) measured the far infrared spectra of the structure I clathrate hydrates of ethylene oxide, cyclopropane, and trimethylene oxide in the spectral range 330 - 15 cm'l. The polar molecules ethylene oxide and trimethylene oxide each have a major axis about which the rigid molecule can be rotated without reorienting the dipole moment. Transitions corresponding to rotations about these axes are forbidden when the molecules are in the free gas state. Not surprisingly, corresponding absorptions by the polar guests were not seen in the hydrate spectra. However, each polar guest did exhibit two rotational absorption bands corresponding to the two axes of rotation that reorient the dipole moment. At high temperatures (-100 K), the absorption features smear together to form a low frequency plateau below 100 cm-l, which is clearly seen in Bertie and Othen's (1973) spectra in Figure (2.2). Gerbaux et al. (1975) observed rotational absorption bands in mixed hydrates that contained the polar H2S molecule as a guest. Because CO2 is not a polar molecule, CO2 gas does not exhibit pure rotational absorption bands. We therefore do not expect it to exhibit pure rotational bands in its

39 hydrate spectrum. This is an advantage for remotely sensing C02 hydrate. C02 has a high partial pressure on Mars (-7 mbars at the surface). If it had a permanent dipole moment, it would no doubt have many rotational absorption lines in the far-infrared. These lines would obscure the rotational emission spectrum of the hydrate, probably making it impossible to detect. More importantly, the rotational lines would obscure other parts of the far-infrared that might be used in remote sensing. We will see the importance of this when we consider the absorption spectrum of the Martian atmosphere. 2.4.2: Vibrational Absorption Fleyfel and Devlin (1988) measured the vibrational spectra of several guest species including C02. C02 was the only non-polar molecule they examined, and they found that their vapor deposition technique could only enclathrate it if they mixed it with a polar help gas. They used ethylene oxide as a help gas, and its presence affected the C02 vibrational spectra. Nevertheless, it is still possible to determine the spectral range of the absorption and whether other Martian environmental constituents would interfere with it. The asymmetric stretching mode absorption of the C02 guest is shown in the main part of Figure (2.6). In curves (a), (b), and (c), the H20/C02/C2H40 ratios used for vapor deposition are 5/0.66/0.33, 5/0.5/0.5, and 5.2/0.25/0.75. Two peaks are plainly visible. Fleyfel and Devlin claim that each of the peaks corresponds to absorption by a C02 molecule in a different cage type. The two different sized structure I cages result in different shifts in the vibrational absorption. The reduction in the lower frequency hydrate peak is the result of competition between the ethylene oxide and CO2 molecules for occupation of the smaller structure I cages. The double peak result is unique to the hydrate. In C02 ice, for example, only one peak is present, as shown in Figure (2.7). It would seem that the hydrate vibrational spectrum would give rise to a distinctive emission. Unfortunately, it lies well within the

40 I I 1 Figure (2.6): CO2 Asymmetric Stretching Mode Absorbance (Main Figure) and Bending Mode Absorbance (Inset) in C02 Hydrate (Reproduced from Fleyfel and Devlin, 1988). The letters identify samples prepared with different ratios of CO2 to the help gas ethylene oxide.

41 2.80 2400 2350 2300 2250 wavenumbers (cm-') Figure (2.7): CO2 Asymmetric Stretching Mode Absorbance in CO2 Ice at -80K (Warren, 1986)

42 range of the analogous C02 gas absorption band, whose line strengths are shown in Figure (2.8). A vibrational emission signature due to the hydrate absorption band would be obscured by Mars' atmospheric CO2 absorption. There is a similar problem with the doubly degenerate C02 guest bending mode. The data of Fleyfel and Delvin (1988) are shown in the inset of Figure (2.5), and the C02 ice spectra (Warren, 1986) and the gas live strengths (Rothman et al., 1987) are shown in Figures (2.9) and (2.10), respectively. In this case, the CO2 ice has two lobes and the hydrate has a single lobe. However, we cannot see through the Martian atmosphere at these frequencies to observe the hydrate. We need an absorption mechanism that occurs in the hydrate that has no analog in any of the other surface or atmospheric constituents. We will find just such a mechanism in the translational modes. 2.3.3: Translational Absorptions In addition to vibrating and rotating, the guest molecule can rattle inside its cage. Due to the confinement of the cage, this translational motion is quantized into energy eigenstates. A simple model (chapter 3) for a caged molecule indicates that these absorptions occur in the far-infrared. Due to coupling to other types of guest motion and to the lattice, guest molecule translational mode absorption features are expected to be broad and weak. Unlike the other absorption mechanisms, the translational modes do not have analogs in the pure solid, liquid or gaseous phases of water or the guest. As a result, we expect the translational absorptions to be unique to the clathrate hydrates. If translational absorption features are sufficiently strong and are not largely overlapped by other absorption features, they may be useful mechanisms for emission in a remote sensing measurement.

43 -20.0 -2 1.0 -22.0. -23.0 -24.0 M 2200 2250 2300 2350 2400 wavenumbers (cm'l) Figure (2.8): C02 Gas Line Strengths in the Asymmetric Stretching Absorption Region (Rothman, 1987)

44 7.1. 680 670 660 650 640 wavenumbers (cml') Figure (2.9): CO, Ice Spectra in the Bending Mode Absorption Region (Warren, 1986)

W-mle V >Oa~ll4~ ~45 -21.0 * -23.0 ~, -24.0 640 6_50 660 670 wavenumbers (cm-') Figure (2. 10): CO,) Bending Mode Line Strengths in CO2 Gas (Rothman, 1987)

46 Gerbaux et al. (1975) measured the far-infrared spectra of CC14, CHC13, and CH2C12 structure II clathrate hydrates with H2S added to each sample to aid in the clathration process and to fill up the smaller cages. Their spectra is given in Figure (2.11). The solid lines are the absorption coefficients of their hydrate samples, and the dotted line is the result of their model of the rotational absorption of the H2S molecule. We see that the far-infrared spectra are dominated by an absorption at about 39 cm1l, with weaker absorptions in the cases of CC14 and CHCl3 hydrates and no additional absorption in the CH2C12 case. All three measured spectra also include an increase in absorption at higher frequencies, which is due to the low frequency tail of the lattice absorption of the disordered hydrate crystal. The absorption at 39 cm-l is common to all three spectra and matches the rotational absorption model closely. It is therefore probably due to H2S rotations. Gerbaux et al. say that the weak peaks at 52 cm-1 and 57 cm-1 can be explained by translational mode transitions. We shall see in chapter 3 that the frequencies of the peaks are consistent with a simple model for the translational energy states of a guest molecule. We shall see that the transition energies for the CH2C12 guest are lower than those of the CC14 and CHC13 guests, implying that a CH2C12 translational mode absorption would occur at a lower frequency and could be obscured by the strong H2S absorption. This would explain why no CH2C12 peak has been observed. The model in chapter 4 predicts that translational mode absorptions will occur in the far-infrared, but this spectral region is also the realm of rotational absorptions. We should consider the possibility that the data of Gerbaux et al. is due to rotational modes. CHC14 is a polar molecule with allowed rotational transitions in the free gas state. It may have guest species rotational transitions. However, CC14 has a symmetric tetrahedral shape with no permanent dipole moment. Hence, the weaker peak in its spectrum is not due to a CC14 rotational transition. It could be due to H2S rotations about a different axis

47.- 1720 C H2s,M) 45K -— H2 S Th. ) C Help Gas (reproduced from Gerbaux et al., 1975) 200_ _.1200 200 - 20 40 60 so 100 Help Gas (reproduced from Gerbaux et al., 1975)

48 than the axis corresponding to the 39 cm- l absorption, but then the same peak should occur in the spectra of the other hydrates, which it does not. Other than Gerbaux's data, we could not find any other spectra in the literature showing (probable) translational absorptions. The spectra of Bertie and Jacobs (1977) (also discussed in section 2.3.2) show the far-infrared spectra of structure I hydrates of ethylene oxide, trimethylene oxide, and cyclopropane. Both ethylene oxide and trimethylene oxide exhibit rotational absorptions, which may mask any translational absorption. Cyclopropane, on the other hand, is non-polar, and one might expect that its translational mode spectrum would be clearly visible. However, the cyclopropane molecule is so large that it forms a structure I hydrate with some difficulty. Its size is about 5.5 Angstroms, which is in the transition region between structure I and structure II hydrates (Davidson, 1973). It must be contorted in some way to fit in the smaller 5.0 Angstrom structure I cages and probably will not rattle appreciably. It is about the same size as the larger structure I cages, and if it rattles, it will do so at higher frequencies. Such a rattling mode would probably be obscured by the strong lattice mode absorptions. We will discuss translational mode frequencies of cyclopropane as well as ethylene oxide, trimethylene oxide, and other guest species using a simple model in Chapter 3. CO2 is not a polar molecule, and we therefore do not expect the far-infrared spectrum of its hydrate to have rotational absorptions that might mask translational absorptions. If its rattling mode absorption occurs at a low enough frequency, it will not overlap the strong absorption in the tail of the lattice mode absorption. The 4.7 Angstrom C02 molecule is much smaller than cyclopropane molecule, and it should have a lower translational absorption frequency. From the model in the next chapter, we will see that possibility of a distinctive C02 guest translational absorption peak is quite promising.

2.4: Summary In this chapter we have considered the mechanisms for C02 hydrate absorption. For all mechanisms except translational vibrations, there is an analogous mechanism in water ice or in C02 gas or ice that causes similar absorption. As a result, emission of C02 hydrate would be difficult to distinguish from the emission of one of these other background constituents. If C02 gas has a similar absorption, it is not even possible to see through the atmosphere to the surface where the hydrate might form. Translational motion of the C02 guest molecule is the only mechanism that does not have an analog in a background constituent. The cage that causes the confinement of the guest leads to translational energy states, and no cage for C02 molecules exist in its ice or gas phases. We therefore expect the translational absorptions to be unique to the hydrate.

CHAPTER 3 A SIMPLE MODEL OF THE TRANSLATIONAL ENERGY MODES OF A GUEST MOLECULE IN A CLATHRATE HYDRATE CAGE 3.1: Introduction It is well known that confinement causes quantization in a wave system. In a rectangular metal resonator, for example, only certain types of waves that have well defined spatial frequencies can exist. These allowed spatial frequencies are quantized and correspond to particular modes of the reservoir. Mathematically, it is the application of the boundary conditions to the solution of the wave equation that gives us discrete values. Clathrate hydrate cages are similar to resonators in that they confine the wave function of the guest particle. If the walls of the cage are rigid, we expect physical quantities (momentum, energy, etc.) associated with particular modes of the guest to be quantized. In particular, we expect that the energy levels of the particle associated with the movement of its center of mass inside the cage to be quantized.This is strictly true only if the walls of the cage are perfectly rigid and if the rattling of the guest within the cage does not couple to its intramolecular motion. Because the translational motion of the guest can couple to the host lattice motion and to its own internal vibrations and rotations, the rattling mode energies of the guest are smeared, and we expect rattling mode absorptions to be broad. In this chapter we present a simplified model of the translational motion of the guest. The size of the cage and guest are based on their actual physical sizes, but their shapes are chosen to simplify the calculations. Furthermore, we assume that the cage walls 50

51 are rigid, and that the motion of the guest's center of mass does not couple to its rotation or intramolecular vibration. The discrete energy separations predicted are therefore only rough estimates. They are in reasonable agreement with observed spectra. 3.2: Potential Model and Wave Equation CO2 clathrate hydrate is a structure I clathrate hydrate and has two roughly spherical cage types, which were shown in Figure (1.2). The corners of the cages are oxygen molecules that are bound to each other and to the rest of the lattice through hydrogen bonds. X-ray crystallography measurements show that the smaller of the two cages has a pentagonal dodecahedron structure that can accommodate guest molecules up to about 5.1 Angstroms in diameter. The larger cage has a tetrakaidecahedral structure and can accommodate molecules with a size of about 5.8 Angstroms. The shapes of these cages are close to spherical, and we will assume that they are spherical in our model. We will also assume a zero potential inside the cage and a large constant potential Vooutside the cage. Furthermore, we will treat the guest molecule as a hard sphere that is incapable of rotating or vibrating. The potential U (r) of a molecule with radius rm in a cage with radius rc is then given by O, r < rc- rm U(r) = (3.1) Vo, r > rc rm This potential only accounts for motion of the center of mass and does not consider intramolecular energy.

52 We now consider the energy eigenvalue equation for the guest molecule in the cage. fe' (r) = ETf (r) (3.2) Here Hi is the Hamiltonian operator (energy operator) and T (r) is the wavefunction of the guest. v (r) is an eigenfunction of H and has energy eigenvalue E. The Hamiltonian has the form H - 2m P+ U(r) (3.3) where p is the momentum vector operator, m is the mass of the guest molecule, and 2 2m is the kinetic energy operator. The momentum operator in position coordinates is p = -ihV (3.4) p ~ is thus given by op = -h2V2 (3.5) where V2 is the Laplacian operator. Using (3.3) and (3.5), the energy eigenvalue equation (3.2) can be written as 2m - 2m V2V (T) + U (F) = (F ) EV (3.6) For the regions r > rc - rm and r < rc - rm, the potential is constant. Therefore, using a generic constant potential UO, we get the following general form of our energy eigenvalue equation:

53 2m -2m V2~ (2 ) = (E- UO ) N(37 This is the equation to be solved with appropriate boundary conditions to get the allowed translational energies of the guest molecule. 3.3: Solution of the Eigenvalue Equation and Boundary Conditions We can rewrite (3.7) slightly to get,V2i() 2m (E- U) (3.8) V2~ (2) =- () (3.8) h2 This is the eigenvalue equation for the Laplacian operator and has been solved in many books for different coordinate systems. In order to facilitate the application of the cage boundary conditions, we consider (3.8) in spherical coordinates. d2 2d 1 d!A2 cot~ 1 d21 2) 2m(E-Uo) - 4 4- -- + - Coto + - 2 —f (r) 2 - g(r) (3.9) r2 r 2d02 r2 cotO r2sin20 d ) h2 If E > Uo, we get the general solution (2m (E U=A( 2 ) r,)Y7( ) +By(m 2E- U0) r lI' (0,, ) (3.10) Ain h 2 l n h2 where j, and y, are spherical Bessel functions of the first and second kind, respectively, r' is a spherical harmonic1, and A and B are arbitrary constants that will be determined by boundary conditions and normalization. If E < Uo, we get solutions of the form 1. We could express the angular dependence in a more general way than a spherical harmonic. However, since the spherical harmonics are a complete set and the angular dependence does not explicitly affect the energy eigenvalues, we will not bother using a more general form.

54 y (r) - Ci 2m(Uo A Ed ~)Y2 ( +=Dk( O( - E)2 rmE) (0), (3.11) where in and kn are modified spherical Bessel functions of the first and third kind, respectively, and C and D are arbitrary constants. We now apply boundary conditions to quantize the energies. We start with the solution Nfl in the region r > rc - rm and the boundary at infinity. If the particle is bound to the cage, the integral of the probability density v1 (r) 2 over the region must be finite. J d~rlv l1(r)I2 (3.12) r> r~- rm The asymptotic behaviors at infinity of the unmodified and modified spherical Bessel functions under consideration are jn (z) isinz, -cosz Z Z 1. 1 Yn (z) sinz, - cos (3.13) Z Z Z --- ~ (3.13) in (z) - eZ kn (z) -eZ The integral in (3.12) will only be finite when kn is used in the solution. Therefore, for bound solutions, the wavefunction in the region r > rc - rm has the form (r) = Dk, ( 2m (V - E)rm (3.14)

55 We see that E < Vo for a bound state. The physical explanation for this is quite intuitive. If E > VO, the particle can cross the barrier of the cage and retain enough kinetic energy to move away from the cage indefinitely. We now consider the inner region r < rc - r,. We require the energy E to be greater than or equal to the potential energy inside the cage. If this were not true, then the average kinetic energy of the system would be negative, which is non-physical. Hence, inside the well, E> UO and we need only consider solution involving unmodified spherical Bessel functions. Furthermore, we constrain the wavefunction to be finite inside the cage. Therefore, we use only spherical Bessel functions of the first kind in the inner region, where the solution Wo is to (i) = Aj (2lE- U) )r, ) (3.15) We will now consider the boundary conditions at r = - rm. The energy eigenvalue equation for the problem is given by 2m 2mV2v (i) + U (F) V (F) = EV (F) (3.16) Since U (F) has a discontinuity at r = rc - rm, V2N (F) is not continuous there either, and hence d~r- (F) is not as well. However, if we integrate (3.16) with respect to r, we get an equation that has no discontinuities. Therefore, dd — (F) is continuous at r = rc - rm, as is Wr (F). These boundary conditions can be expressed as 2mE2m ( LE- (rUO) r))/(,) =Dkn(t 2m( (rc -rm)y) "' (0,) (3.17)

56 A2m (E- U) 2m (E- U) A h2 (rh- r. ) Z, =D'| h2 k' ( (up) (3.18) dr We must have I = l' and m = m'" to be able to match W (F) and dr- (T) at the boundary. The boundary condition now reduce to Ai( 2m (E- U0) 8(rc-r) Dk( 2m (rc - rm) (3.19) 2m(E- Uo) ( 2m (E- UO) A h2 h Y (rc - rm)) =D 2m (V- E) k( h2 (rc-r ) (3.20) h 2m(V2-E) The unknowns in these two equations are A, D, and the energy. A and D can be reduced to one unknown through normalization, leaving us with two unknowns to be determined by the two equations. This leads to a discretization of the energy eigenvalues. If the cage barrier energy V0 becomes large, the guest molecule wavefunction will be highly localized in the inner region. As we raise the barrier, the magnitude of the wavefunction decreases in the outer area. As V0 -e o, we get a new boundary condition at

57 Ain( 2m(E- (rc-rm)) = (3.21) We will use (3.21) to predict the energy spectrum of the translational modes of the guest. If zn, i is the ith zero of jn, then the spectrum is given by h2 En 2m(rrm) i (3.22) 3.4: Energy Spectra. Figures (3.1) and (3.2) show the lower levels in the energy spectrum for the guest molecule in each of the two possible structure I cage types. The frequencies of several transitions are also shown in the figures. The larger cage has lower energies because the guest has further to move in the cage and thus will rattle at a lower frequency. Typical transition frequencies between states are 10-20 cm1l. We will see that the low frequency transitions for the larger cage case are quite comparable to the frequencies of the absorption band centered at 17 cml1 shown later in Figure (4.12). It is therefore likely that the band is due to absorption by the guest in the larger cage. Given its crudeness, the model has done exceedingly well in predicting the region of absorption. It is interesting to see how well the model can predict rattling mode absorptions in other clathrate hydrates. Gerbaux et al. (1975) used H2S as a help gas to form structure II mixed hydrates of CC14, CHC13, and CH2C12. They measured their farinfrared spectra and found absorption peaks in two cases that they attributed to rattling modes. We will examine how well the model predicts the absorptions they measured. Below is a table summarizing the sizes and masses of their guest molecules and the frequencies of their measured absorption peaks. A strong, broad absorption due to the H2S

500 cm'1 69 cm'1 89 cm'1 250cm-1 124 cm'1 98 cm-1 0 cm'l Figure (3.1): Energy Spectrum of the Translational Modes of the C02 Guest in the Smaller Structure I Cage

59 100 cm4' 10 cm'1 23 cm'1 12 cm-l 50 cm' _,, 24 cm'1 16 cm-l 13 cm'1 0 cm' Figure (3.2): Energy Spectrum of the Translational Modes of the C02 Guest in the Larger Structure I Cage

60 help molecule was present in all measurements and tended to obscure the rattling mode absorptions. In the case of CH2C12, where no rattling absorption peak is seen, the H2S absorption may have completely hidden the rattling mode absorption. absorption diameter mass absorption frequency (angstroms) (amu) (cm1) stucture II 6.7 large cage CC14 6.45 152 57 CHC13 6.4 118 52 CH2C12 6.2 84 obscured? Table (3.1): Sizes, Masses, and Absorption Frequencies of Structure II Hydrate Guests The predicted energy spectra of the three molecules are shown in Figure (3.3). We see that the transition frequencies predicted for CC14 and CHC13 compare favorably to the observed absorptions of 57 cm-1 and 52 cm-1 given in the table. The model correctly predicts that the transition frequencies of CHC13 are lower than those of CC14. The physically intuitive reason for the lower frequency of the smaller CHC13 molecule is that it requires extra time to rattle back and forth through a cycle and thus rattles at a lower frequency. The lower rattling frequency of the guest will couple to a lower frequency photon. We also see that the model predicts the smallest molecule CH2C12 to have the smallest transition frequencies, which is consistent with its rattling absorption being hidden by the low frequency H2S absorption.

400 cm'1 300 cm66 cm'1 300 cm-i 59 cm'1 200cm'1 92 cm'1 83 cm'1 65 cm' iO0cm-i' 73 cm'1 00~~CM~~ I _.65 cm'1 " - 51 cm'1 CHC14 CHC13 CH2C12 0 cm'1 Figure (3.3): Energy Spectrum of the Translational Modes of CC14, CHC13, CH2C12 Guests in the Larger Structure II Cage

62 3.5: Shortcomings of the Model The model presented in this chapter is simplistic, yet it appears to do a good job of predicting the range of frequencies of rattling mode absorptions. It should not, however, be trusted to predict exact absorption frequencies, and its neglect of coupling makes it useless in predicting linewidths. In this section we discuss these and other shortcomings of the model. The energy eigenfunctions in any potential well problem are a strong function of the boundary conditions used at the well edges. In our model, we assumed an infinitely sharp and infinitely high barrier, whichled to an infinite number of bound states, all of which must go to zero at the cage wall. Making the barrier of finite height no longer forces the wavefunction to go to zero at the cage boundary, and this generally results in lower energy eigenvalues for bound states. The effect is illustrated in Figure (3.4). The energy of a state is related to the number of oscillations the wavefunction has inside the well, with a higher energy state having more oscillations. In the figure, the wavefunction inside the finite well can spread out into the classically forbidden region, allowing lower frequency oscillations in the well region, and, hence, lower energies. A finite barrier also results in a finite number of bound states. More importantly, if the barrier has a finite width, there are no bound states because a guest with any finite kinetic energy can eventually tunnel through the barrier and escape. Each quasi-bound state would therefore have some mean lifetime. Our infinite potential model does not provide for quasi-bound states, and therefore cannot predict lifetimes. When the guest molecule is very large and barely fits in its cage, its translational energies will be very high. More importantly, a small change in the guest or cage size in the model can lead to a large change in the predicted energy spectrum. For example, in the case of the structure II hydrate of CCL4, we used the value of 6.45 Angstroms as the diameter of the guest, which resulted in a lowest energy of 64 cm1l. Had we used a value

63 T a oU3 = / o CUC radius 4't va 0 U o CU radius -e Figure (3.4): Sample Energy Eigenfunctions for Infinite and Finite Spherical Square Wells

64 of 6.5 Angstroms, the lowest energy would have been over 100 cm-l. Therefore, our model can be exceedingly sensitive to guest and cage sizes when the two are close to each other. Perhaps the most serious limitation of our model is its lack of coupling to the lattice walls and to the intramolecular vibrations and rotations. Our approximation results in discrete energies and sharp absorption lines, which are not observed spectroscopically. More importantly, the lack of coupling prevents us from saying anything about which transitions are allowed and how strong such transitions would be. To adequately answer this question, we have to calculate matrix elements between the various states for the type of excitation. This is exceeding difficult because it requires knowledge of the change in polarization of the guest molecule as it moves in its confining potential. Our hard sphere model of the guest is not capable of predicting changes in polarization, and hence is not useful in predicting selection rules. We will not address the issue of selection rules in this study, but will instead assume that enough asymmetries exist in the system to allow for far-infrared transitions.

CHAPTER 4 FAR-INFRARED MEASUREMENT OF CO2 CLATHRATE HYDRATE 4.1: Introduction Arguments in Chapter 2 led to the conclusion that guest molecule translational modes are probably the most suitable absorption mechanism for producing a distinctive emission signature of C02 hydrate. Measurements of other hydrate translational absorptions (Gerbaux et al., 1975) and the simple model in chapter 3 indicate that such absorptions will occur in the far-infrared. The dielectric spectrum resulting from these absorptions is necessary to model the emission from a layer of hydrate particles. To determine the spectrum, we need a device that can measure far-infrared dielectric properties of materials and a method for preparing the sample for measurement. This chapter begins with a survey of the available far-infrared measurement technologies and concludes that a Fourier transform spectrometer (FTS) is best suited for our purposes. It then describes the FTS and the cryostat used in our measurements, followed by a long discussion of the clathrate hydrate sample cell design process. This includes a discussion of a failed vapor deposition technique and more successful frost and slab sample measurements. Commentary is given on the pitfalls of doing far-infrared spectroscopy and working at low temperature. The chapter ends with a discussion of the results of the measurement. 65

66 4.2: Far-infrared Measurement Technologies Since it is difficult to determine a priori the frequency ranges of the C02 translational absorption bands, it is preferable to use a wideband far-infrared instrument to make hydrate dielectric measurements. This will allow for determination of absorption peak positions and widths and will also show how competing absorptions might obscure translational spectra. A reasonable requirement for the frequency range of the instrument is 5-100 cm-l. Above 100 cm-l, a translational peak would probably be overwhelmed by the strong tail of the hydrate lattice absorption. Below 5 cm1l, the guest would have to be very small to rattle at such a low frequency. An ideal instrument would make a continuous measurement across this range so as not to miss any important spectral features. Unfortunately, there are no tunable coherent sources that can continuously span this range or even a sizeable fraction of it. Schottkey diodes together with multipliers can produce moderate power in this range, but only at discrete points or narrow ranges. Lasers can produce high power far-infrared radiation, but suffer from the same nontunability problem. As a result, we must rely on broadband, low power sources. A relatively new technique for measuring far-infrared properties of materials uses ultrafast laser technology. We visited the University of Michigan Ultrafast Sciences Laboratory (USL) to determine if the technique was suitable for our purposes. A simplified geometry of the technique is pictured in Figure (4.1). Antennas A and B consist of two metal pieces connected by photoconductive material. Antenna A is biased as shown, and an ultrafast laser pulse is directed at the photoconductor, shorting the poles of the antenna and causing a brief current surge. The accelerating charge gives rise to a broadband electromagnetic field that propagates through the test material and impinges on antenna B. The field biases antenna B, and this varying bias field can be sampled by shorting the antenna's photoconductive gap with a delayed version of the pulse that excited antenna A. The experiment is run over and over with the delay between the two

fixed mirrors ultrafast laser pulse generator movable mirrors biasing Antenna A Antenna B test sample Figure (4.1): Setup for Ultrafast Laser Technique to Measure Far-infrared Dielectric Spectrum

68 pulses stepped so that a uniformly sampled sequence of the antenna B biasing field is obtained. Fourier transformation of the sequence gives the spectrum of the field. By collecting spectra with and without the test sample, a transmission spectrum can be constructed by ratioing the two raw spectra, and the dielectric spectrum can be extracted. The ultrafast technique has been used successfully to measure spectra up to a few terahertz and offers a signal to noise ratio superior to FTS for frequencies below 1.5 THz. Unfortunately, a system is not commercially available in a single package and requires significant expertise to build and maintain. While it is available to us at USL, its access is restricted, and its geometry is not convenient for integrating with a clathrate hydrate sample cell. We found it inconvenient for our purposes. FTS has been used by many investigators to measure the infrared transmission and reflection spectra of samples of all phases and over wide temperature and pressure ranges. The nominal range of commercially available systems is from 4 cmml up to the ultraviolet, with resolutions as fine 0.0015 cm-l. It is a mature technology that is standard equipment in many university and industry laboratories. A simple explanation of FTS theory is given in Appendix A. FTS would seem at first to be ideal for our purposes. Unfortunately, its performance in the very far-infrared is limited by the available radiation sources. The typical broadband FTS source for the far-infrared is a hot black body (-1500 K). The power of such a source is quite low at the low frequency end of the far-infrared spectrum, especially below 50 cm-l. To help compensate, a much hotter mercury vapor lamp is typically used. The quartz envelope of the lamp is opaque in much of the far-infrared, but below about 100 cm-l, it becomes more transparent, and the hot mercury vapor (-10,000 K) emission power surpasses that of the standard globar ceramic source. The power is still quite low, however, and to make low noise measurements requires a very sensitive detector. Liquid helium cooled silicon bolometers operating at 4.2 K or 1.6 K are available that can measure below 10 cml when appropriate filters are

69 used to remove the higher frequency power. These bolometers require several hours of setup time and have limited helium hold times (-1 day at 4.2 K, -6 hours at 1.6 K). They can greatly limit the flexibility of an experiment, especially one that requires the complex procedures necessary for growing a clathrate hydrate. Despite its drawbacks, FTS is the most reliable technology available for measuring far-infrared dielectric properties. The far-infrared measurements of other hydrates discussed in chapter 2 were all made with FTS instruments. In the next section we describe the FTS procured for this experiment. 4.3: Bomem DA8 FTS Our FTS is a Bomem DA8 high resolution spectrometer and is shown in Figure (4.2). Some of the nominal system specifications are given in table (4.1). Fine? Resolution 0.0026 cm'l Coarsest Resolution 32 cm-1 Frequency Range 4 - 8500 cm'1 Sources Mercury Vapor Lamp Globar Ceramic Source Quartz Halogen Bulb Beamsnlitters Mylar Films KBr CaF Detectors Silicon Bolometer Deuterated Triglycine Sulfate Manganese Cadmium Telluride Indium Antimonide Table 4. 1: Specifications and Components in Our Bomem DA8 System

Scan motor Scan tubs Upper section Si Source cover plate Scan control indicators.-Bamrspllttsr cover plate sapl COR4110oe CV9 Emission port -jj-ItJlI emalplon beam 105to 10mm aperture r'ZirBltv — ~I% m dismeterl8esm 32w f I Middle if~~ LJ r)~ eu~ eance, to 32Mi J dsection eletor modules *Ci Front panel 1I4 been laI*mosed of senple COuPhels06U Lower section Figure (4.2): Bomem DA8 Fourier Transform Spectrometer (From Bomem DA8 Training Manual)

71 The system is used by other investigators at the University of Michigan to study mid- and near-infrared solar emission spectra, which explains why the system can measure up to 8500 cm'l and has a maximum resolution of 0.0026 cm-l. For our experiment, we only require a resolution of -1 cm'1. The optical path of the instrument is shown in Figure (4.3). It is physically divided into three regions: the source compartment, the interferometer and transfer optics compartment, and the sample compartment Each of these compartments can be independently accessed through plates bolted to the sides of the main housing. The source compartment contains three sources to cover the infrared, visible, and part of the ultraviolet spectrum. (The high frequency limit of our instrument is dictated by our detectors.) A rotating parabolic mirror under software control selects the appropriate source. For our far-infrared work, we use exclusively the mercury vapor lamp. Radiation from the lamp is focused down through a rotatable filter wheel and a variable aperture to the interferometer and transfer optics compartment. Transfer optics redirect the beam to a collimating mirror, which sends a parallel beam to the interferometer section. The interferometer consists of a beamsplitter, a fixed retromirror, and a moving retromirror. Our far-infrared beamsplitteres are stretched Mylar films of varying thicknesses. The beamsplitter divides the collimated beam, sending a subbeam to each of the retromirrors. The two sub-beams are reflected back and recombine at the beamsplitter to form a new beam with a modulation that depends on the path difference of the two subbeams. This modulated beam is directed by a series of transfer optics to the sample compartment. The sample compartment can accommodate samples of moderate size. The modulated beam can be focused through samples on either side of the compartment and onto a detector. We used a liquid helium cooled silicon bolometer operating at 4.2 K for our sensitive measurements, and a room temperature deuterated triglycine sulfate (DTGS) detector for crude measurements. The signal from the detector is sampled by an A/D

72 SOURCE 3 ( SOURCE 2 SOURCE SELECTION MIRROR SOURCE I t t, RAED WROR FIAT TPANSFER MOUATED,BEA~ME O',UTPUT U,' "" ROR L~~~~~~~E~~~~~~R SAMPLE DETECTOR?.EECTON FOCU DEDETECTO WIRRYMIRR OR Figure (4.3): Optical Path of Bomem DA8 Fourier Transform Spectrometer MIRROROF

73 board, Fourier transformed by the FTS vector processor, and transferred via an ethernet cable to our VAX 3100 control computer. Because water vapor exhibits high absorption in the far-infrared, the entire optical path is evacuated during measurements. The source and transfer optics can be isolated from the sample compartment by a vacuum flap that rotates up to the sample compartment portal and seals against it with an O-ring. This allows a sample to be changed or removed without disturbing the interferometer sections, thus insuring that the sample and reference spectra are acquired under the same optical conditions. The interferometer has a dynamic alignment control system that insures that the two sub-beams of the divided beam are parallel when they recombine. As the moving retromirror scans, it tends to rattle and tilt slightly. To correct for this, a HeNe laser beam is introduced into the interferometer, and the modulated laser beam output is sent to three detectors that measure phase differences across the beam cross section. If the detectors show too large a phase variation, a pair of motors tilt the fixed retromirror to realign the two beams. Our system suffers from a 15 Hz mechanical vibration that introduces spurious spectral features. The interferometer modulates the far-infrared radiation in this frequency range, and as a result, the mechanical vibration causes an artificial spectral feature in our far-infrared spectrum. Because the vibration is not in phase with the motion of the interferometer mirror, the feature that it causes will diminish as more interferometer scans are averaged. Unfortunately, the rate at which it diminishes is very slow. In theory, the amplitude of the spurious feature is proportional to 1/1R, where N is the number of scans in the average. Thus, reducing the spurious feature by a factor of 10 can require 100 times more scans, which can increase the time span of an experiment by an unacceptable amount. Instead of trying to get rid of the feature through averaging, the apparent farinfrared frequency of the feature can be shifted by changing the mirror speed. This is in

74 fact how such a feature is identified. If a feature's frequency changes when the mirror speed is changed, it is not a real feature and must be due to some external source, such as a mechanical vibration or 60 Hz line noise. By changing the mirror speed, the feature can be moved into some unimportant part of the spectrum, and the interesting part of the spectrum is left unaltered. 4.4: Cryostat In order to simulate the conditions in the Martian northern polar regions, the hydrate sample must be cooled to about 150 K. This is accomplished with an Oxford CF1 104 continuous flow cryostat. A simplified drawing of the cryostat is given in Figure (4.4). We will briefly discuss its operation. In the cryostat's most basic configuration, a sample is thermally contacted to a coldfinger, as shown in the figure. (We had to modify this configuration several time during the hydrate cell design process. This will be discussed in the next section.) The cold finger is typically made of deoxygenated copper or some other metal with a high thermal conductivity. The detector either measures the power reflected or transmitted by the sample. For reflection measurements, special transfer optics are required to collect the reflected beam and route it to the detector. For transmission measurements, the cold finger has a hole to allow the beam to pass through. The cold finger is attached to a copper block, which is directly cooled by cryogen. The cryostat cooling system is shown in Figure (4.5). While the system was originally intended for use with liquid helium, we used liquid nitrogen and found that it performed adequately. The nitrogen was pumped up from our MVE Lab 30 dewar, through the center pipe of the horizontal span of the gas flow shielded transfer tube, and into the cryostat. The exhaust from the cryostat was pumped back up into the transfer tube, through the outermost tube of the horizontal span, and out through a needle valve to the transfer tube's

75 cryogen inlet cryogen vent heater ____ copper block > cold finger Wsample spectrometer beam Figure (4.4): Cryostat

vacuum space needle valve cryogen vented cryogen gas flow shielded transfer tube flow meter and needle valve vacuum jacket dewar optical access ports pump Figure (4.5): Cryostat Cooling System

77 exhaust port. This port was connected through a flow meter and a manual needle valve to a diaphragm pump, which provided the pumping action for the system. The temperature was measured with a temperature sensor mounted either on the sample or the cold finger. We used a platinum resistance temperature detector (RTD), which is a resistor with a resistance that is nearly linear with temperature over a wide range. Our Oxford ITC4 temperature controller provided a constant 1 mA current across the RTD and independently measured the resulting voltage drop at the RTD contacts. This four wire measurement thus eliminated errors that a two wire measurement would incur due to voltage drops in the long lead wires. The RTD was calibrated at the beginning of each day of measurements by immersion in water ice and liquid nitrogen baths. The calibration was checked on a few occasions at the intermediate temperature of dry ice. In each case it was found to be accurate to better than 1 K. The temperature could be controlled by either heating the sample or restricting the flow of cryogen. Above the cryostat's copper plate were two ohmic heaters whose voltages could be varied by the ITC4 controller. The needle valve at the transfer tube exhaust contained a small opening and closing stepping motor that was also under ITC4 control. The manual needle valve in the flow meter box could be used for temperature control. We found it most convenient to run with the motorized needle value all the way open at all times. If it were less than half way open, it tended to freeze in place and clog, which requires about an hour of cooling and purging to remedy. We thus set the ITCr4 for manual control of this valve and left it open. The temperature control method that seemed to work the best was as follows: 1) Set the ITC4 thermostat for the desired temperature and turn on the heater control.

78 2) Open the manual needle value all the way. After a few minutes, the weight in the flow meter should rise and oscillate about the mark for 1.5 liters of helium per hour. (This of course corresponds to a different rate for liquid nitrogen.) 3) When the sample temperature is about 50 K above the desired temperature, close the manual valve until the flow reading drops to about 0.1 liters per hour. In about 20 minutes, the sample should stabilize to within 1 K of the desired temperature. The 20 minute time to stability was accurate when a simple sample was used (e.g., a 0.5 mm wafer of silicon). When a more massive sample was used, such as the clathrate hydrate cell discussed in the next section, the sample temperature could take 40 minutes or more to stabilize. 4.5: CO2 Clathrate Hydrate Sample Cell Design In this section we will discuss several attempts at designing a working hydrate sample cell and how failures or inadequacies led to changes and overhauls in the cell design. Designing the cell was by far the most challenging experimental part of this research and represents over a year of work. Because the cell design was intimately related to the measurement process, the two are discussed concurrently in this section. 4.5.1: Vapor Deposition The most straightforward and common way of preparing a sample is by vapor deposition. This technique requires no sample cell at all. The basic geometry of a vapor deposition experiment is shown in Figure (4.6a). An ejection nozzle deposits a cooled gas onto a cooled substrate to form a sample film. To make a clathrate hydrate, a mixture of water vapor and the host species with the correct ratio is used.

IRA t3~111 FINE El gas ejection nozzle mass of CO2 ice (a): "Classic" Vapor Deposition Geometry (b): Results of Pressurizing the Cryostat with C02 Gas Figure (4.6): Vapor Deposition Method for Clathrate Hydrate Formation

80 Fleyfel and Devlin (1988) attempted to make form a pure CO2 clathrate hydrate using the vapor deposition technique. They were growing their samples at low temperature (<100 K), where reaction rates are low. They found it impossible to form pure C02 hydrate in this way, attributing this difficulty to the nonpolarized nature of the C02 guest. They instead used ethylene oxide as a help gas and grew a mixed hydrate. Measurements of ethylene oxide in the far-infrared (Bertie and Jacobs, 1977) show two absorption bands due to rotation. (See discussion in chapter 2.) These absorptions, which are likely to occur in the far-infrared spectrum of any hydrate with polar guest molecules, could easily obscure the C02 guest absorption. Removing the effects of polar guest absorption from far-infrared spectra would be difficult. We therefore attempted to grow pure C02 hydrate through vapor deposition, and, like Fleyfel and Devlin, were unsuccessful. Miller and Smythe (1970) grew C02 hydrate from water ice frosts by pressuring them with CO2 gas. At 150 K, the required pressure was 4.5 mbar, and the hydrates took a few hours to form. In general, the formation time depends on the temperature and size of the ice particles, with thinly divided samples converting to hydrate much more quickly than thick slab samples. We decided to try a similar technique with a vapor deposited water ice film. We grew a film of water ice on a cooled polyethylene substrate using the classic vapor deposition technique and checked to make sure its far-infrared spectrum was correct. Through the pump out port of the cryostat we then introduced a few millibars of C02 gas and checked the far-infrared spectrum every few minutes. The distinctive lattice mode absorptions of C02 ice began to form at 68 cml and 114 cm-l. The transmission quickly dropped to nearly zero and the raw spectrum was at the noise level. A look into one of the unused optical access ports of the cryostat revealed that a large mass of dry ice had formed around the cold finger and copper block, as shown in Figure (4.6b).

81 While the cold finger may have been above the sublimation temperature for CO2 ice, the copper plate and the metal above it apparently were not. Even when the temperature was raised well above the sublimation temperature of C02 for the given pressure, the mass of C02 remained. The temperature of the metal above the copper plate is apparently always close to the cryogen temperature. As a result, the C02 mass formed and obscured the cold finger even when the cold finger was quite warm. Thus our modification of a simple deposition experiment seemed doomed to failure. The problem resulted from different parts of the cryostat being at vastly different temperatures. The logical way around this seemed to be the use of a cell, attached to the cryostat, that was at a single temperature. This would allow the cell to be pressurized with C02 gas without causing a buildup of C02 ice. We set out to design and build such a cell. 4.5.2: Frost Cell Growing a C02 hydrate sample from the liquid water state requires a few hundred psi of pressure and a robust sample cell. A much simpler cell can be constructed if the sample is grown by pressurizing water ice. The problem with preparing a sample in this way is the long waiting time required for the ice to completely convert to hydrate. As we shall see in the next section, it may be impossible to convert a thick slab of water ice (thickness - 1 mm) in a reasonable amount of time. Apparently C02 molecules have great difficulty diffusing through the surface of a solid hydrate sample. Following the lead of Miller and Smythe (1970), we decided to grow samples from a frost. While they may be relatively easy to prepare, it is important to recognize from the outset that hydrate frost samples will not give quantitative dielectric data. They can give the general shape of the spectrum, including peak positions, but they cannot give the quantitative values of the dielectric constant because they are difficult to characterize and because scattering losses cannot be separated from dielectric losses. We only attempted

82 this method to determine if any absorption features are present and at what frequencies they occur. Quantitative determination of the dielectric spectrum had to await measurements of samples that are easier to characterize. After several design iterations, we finally arrived at the cell configuration shown in Figure (4.7), which is designed for a reflection measurement. The cell is made of polyethylene because it provides low loss transmission in the far-infrared, and because it is easy to machine. The sample is grown on copper foil glued to a 20 mil thick piece of polyethylene sheet. One end of the foil is wrapped around a platinum resistance temperature detector (RTD) that measures the temperature of the sample. The foil stuck well enough to the sheet to form a good thermal contact, but glues do not stick well enough to polyethylene to form vacuum seals. Fortunately, polyethylene melts at a low temperature and is easily welded with a soldering iron to another piece of polyethylene. We were able to make vacuum tight welds that were sufficiently strong to withstand the moderate pressures experienced by the cell. The sheet joined to a shallow tray that was machined from polyethylene stock. The tray was much thicker than the sheet and could easily withstand the moderate pressures of this experiment. 1/8 inch holes were drilled in three corners of the tray, and polyethylene tubes were inserted and welded in place. Two tubes were used to run water vapor through the cell to deposit frost and to pressurize the cell with CO2 gas. The other tube served as a conduit for the wires connected to the RTD. The sheet and the tray were welded together and attached to the cold finger with a metal clamp, as shown in Figure (4.8). The polyethylene sheet was thin enough to allow the copper foil and the sample to be effectively cooled. Because the sheet was in contact with the much stronger aluminum cold finger, it did not need to be able to withstand high pressures. The polyethylene tube containing the wires for the RTD was connected to a 1/8 inch copper tube through a stainless steel Swagelok union. The ferrule clamping down on

83 polyethylene tray polyethylene sheet C02 gas and water copper foil vapor input; "'l.^.t'/gC~L~LL*.>~ *c` —A- | < -bt —t. ___i _S * * w...-i a -s...........;.'~.i.. ". ~~~~~~~~~~~~~~~~~aw.o.~:J *- w....-J- ^... ^ee- s_..............'... i'. _, - >, ^ i,:.~~~ s....:'.'.....'~'_. e...:i';> _. " 7i; $ *S.,..,....;.,~~~~~~~~~~~~~~~~~~~~.~~wate ~~,r vapor output'-s~~ Y~.- / |.';~ I -.. I... - % -% lo~. / / I Platinum RTD RTD) wires Figure (4.7)' Unassembles Frost Sample Cell

84 RTD wires copper tubing polyethylene tubing Swagelok union -optical access port feedthrough C02 belljar pump spectrometer sample Swagelok unions compartment vacuum plate Figure (4.8): Hydrate Frost Production System

85 the polyethylene tube was nylon to provide a better seal. Epoxy was pushed into the far end of the copper tube, and the wires were pushed through the tube before the epoxy had cured. When it did cure, a vacuum tight seal was formed. The protruding wires were then run through connectors to the cryostat's built-in vacuum feedthrough, which connected to the temperature controller. The tubes carrying the water vapor and C02 were also coupled to copper tubes through Swagelok unions. The two copper tubes fed through a brass disk in one of the spare optical access ports of the cryostat and entered the spectrometer sample compartment. The tubes were soldered to the disk to form a vacuum seal. Once out of the cryostat, the two copper tubes were expanded to 1/4 inch polyethylene tubes through another pair of Swagelok unions. The larger, flexible tubing connected to two more vacuum feedthroughs and exited the sample compartment. The outside tubes were then connected as shown in Figure (4.8), with one going to a bell jar containing a source of water vapor (a beaker full of deionized, degassed water) and to a CO2 gas tank, and the other going to a vacuum pump. The water ice frost samples were prepared by cooling the cell to its lowest possible temperature (- 90 K) and opening the valve of the bell jar and vacuum pump to allow water vapor to flow through the cell and condense on the foil. The reflectance spectra of the frost was checked with a the DTGS detector in the region of the water ice lattice mode absorption to insure that the frost was present. After a few minutes of deposition, the cell was sealed off from the bell jar and pump, and the cell was warmed to about 210 K. The valve from the CO2 gas tank was then opened, and the pressure was adjusted to 45 psi, which is well above the 210 K dissociation pressure of C02 hydrate, but well below the saturation vapor pressure of CO2. The lattice mode absorptions of CO2 ice were not present in the reflection spectrum, confirming that CO2 ice was not forming. The cell was left pressurized over night, and the following morning the far-infrared spectrum of the sample was obtained with the liquid helium cooled bolometer operating at 4.2 K.

86 To get a reference spectrum, the sample had to be removed. The CO2 gas valve was closed, the cell was heated to above the freezing point of water ice, and the vacuum pump removed the resulting liquid water. The temperature of the cell was then lowered to 210 K and a reference spectrum was obtained. The transmission spectrum was obtained by ratioing the hydrate and reference raw spectra. The only clearly interesting feature was an "absorption peak" at about 75 cm1l, which is shown in Figure (4.9). The absorption is somewhat higher than what is expected from the model in chapter 3. Furthermore, it was much stronger than expected, dwarfing the tail of the hydrate lattice absorption. Unfortunately, a careful analysis of transmission spectra of polyethylene revealed that it has an absorption in almost exactly the same place and that its absorption increases rapidly as its temperature decreases. What looked like a C02 hydrate absorption peak was apparently only an imperfect cancellation of the polyethylene absorption peak in the sample and reference spectra. What may have happened is that the polyethylene tray was cooled more effectively in the hydrate case due to the C02 gas. No other features were present in the spectrum, but the data below about 30 cm-1 was too noisy to be reliable. Too much power was lost in the complicated path that the spectrometer beam had to take in order to get in and out of the sample for the very low frequency data to be measured. As a result, there could still exist a spectral feature at low frequency that would be useful. To see it would require a transmission measurement. The frost cell can be easily modified for a transmission measurement. If we remove the copper foil from the polyethylene sheet, the spectrometer beam can pass through the entire cell, but it cannot penetrate the metal cold finger. This problem is remedied by drilling a honeycomb pattern of holes in the cold finger. The holes in our honeycomb coldfinger were 1/8 inch in diameter, and the coldfinger gave at least 40% transmission down to a least 20 cm-l. (We only made measurements of it with the DTGS

87 1.10 1.00 0.900 0 E 0.800 0.700 0.600 0.500...... 65.00 70.00 75.00 80.00 85.00 Frequency [I/cm] Figure (4.9): False CO, Hydrate Absorption Peak in Reflection Spectra. The absorption is actually due to polyethylene

88 detector, which is not capable of going much lower in frequency.) It also provided more than adequate cooling of the frost. The major problem with this new configuration is that the resulting spectrum exhibits large oscillations due to constructive and destructive interference in the polyethylene sheet. Though these oscillations perhaps should cancel each other when the sample and reference spectra are ratioed, in practice they never do. This may be due to the slightly different boundary conditions between the frost and sheet in the sample measurement and the vacuum and sheet in the reference measurement. The problem of this "channel spectra" is usually handled by wedging the material causing the problem, as shown in Figure (4.10). At high frequencies, a small angle of wedging will scramble the phase of a wave in the transverse direction, thereby smoothing out interference effects. At longer wavelengths, the angle must be larger, and the wedge works by either scrambling the phases or by bending successively bounced waves at different angles. Wedging the polyethylene sheet at a sufficiently large angle turned out to introduce more problems. The large wedge required made the sheet thicker and reduced the rate at which some parts of the cell could be cooled. This introduced a large temperature gradient in the cell, undoubtedly resulting in only part of the frost being sufficiently cold to convert to hydrate. We decided to end the frost measurements at this point. Because the frost reflection measurement worked down to 30 cm-l, and because we eventually had to make a measurement of a well characterized sample anyway, the effort required to make the frost transmission experiment work did not seem like an efficient use of energy. We instead devoted our time to developing a new cell that could produce a uniform slab of hydrate whose transmission properties could be well understood.

(a): High frequency phase scrambling in a wedged substrate. The thickness of the slab must be greater than or comparable to a wavelength for sufficient scrambling. Because the wavelength is small, the wedge angle is small, and the multiply scatter waves in the substrate are only slightly bent. (a): Low frequency phase scrambling in a wedged substrate. The longer wavelength requires a larger wedge angle. The phase of the incident wave is scrambled, and multiply scattered waves are deflected at different angles. Figure (4.10): The Effect of Substrate Wedging at High and Low Frequencies

90 4.5.3: Liquid/Solid Hydrate Cell One of the most important lessons derived from the frost measurements is that it is important to know whether or not the sample is actually hydrate. In the frost experiment we had no direct evidence of hydrate formation. If verification of formation could be incorporated into the new cell, the results of measurements could be presented with more confidence. We set out to build a cell that would allow for verification. There are a few ways to verify the hydrate structure. Bertie and Jacobs (1977) ground their samples into powder and certified that they were hydrates through X-ray diffraction. A powdered sample, much like a frost, is difficult to characterize and is not suitable for our experiment. Fleyfel and Devlin (1988) found distinctive mid-infrared behavior of the C02 hydrate, particularly in the guest molecule asymmetric stretching mode absorption. Building a cell that will allow far-infrared and near-infrared transmission is difficult given the dissimilar window materials needed for the two spectral ranges. Furthermore, their samples were quite thin (-10 gm), while our samples would probably need to be much thicker (- 1 mm). A sample this thick would probably be too opaque in the C02 asymmetric stretching region to reveal the distinctive features. Perhaps the simplest way to confirm that a sample is a hydrate is to observe its phase behavior. At the beginning of chapter 2, we described an experiment where a bottle of water pressurized with C02 was cooled, and C02 hydrate was observed forming. We did almost exactly this experiment to confirm that we could make a hydrate from liquid water. Figure (4.11) shows the setup of this simple experiment. A cylindrical cavity was machined out of a copper block, and a hole was drilled along the axis of the cavity to form a viewing port. This formed the top piece of the sample cell. A circular plexiglass window was inserted into this piece and sealed with five minute epoxy. A hole was drilled into the

91 C02 gas input port plexiglass window water in aluminum foil cup water ice bath epoxy seals Figure (4.11): Experimental Setup for Observing Hydrate Formation at 0~C.

92 side of the piece, and a copper tube was soldered there to provide a port for C02 pressurization. The bottom piece was a thinner piece of copper on which an aluminum foil cup of water was placed. The top and bottom pieces were sealed together with more five minute epoxy, and the two pieces were tightly bound together with four screws fed through holes drilled through the corners of the two pieces. The entire cell was immersed in a water ice bath. To form a hydrate, the gas port was connected to a CO2 tank, and the cell was pressurized to 415 psi. The water solidified into a clear material which could not have been water because of the high pressure and temperature. We concluded that it must be C02 hydrate. This was confirmed when we slowly lowered the pressure. Atjust above 180 psi, the sample began to melt and bubble as C02 gas was released. This pressure is in very good agreement with the 00C dissociation pressure of 184 psi given by Miller (1985). We therefore have a cell in which we could verify the formation of C02 hydrate. However, the cell is not suitable for far-infrared spectroscopy. It does not have a sample substrate or windows that are transparent in the far-infrared, and its epoxy seals are likely to crack at low temperature. The cell must clearly be modified before a measurement can be made. We will first consider the sample substrate and the window material. Polyethylene, polymethylpetane (trade name TPX), quartz, sapphire, and silicon are commonly used substrates for far-infrared work. Of these materials, polyethylene and TPX are easy to machine and adapt to a cell as its design changes to meet problems that are inevitably encountered. This turned out to be very important because the cell went through several design iterations, and if a new window material had to be ordered after each iteration, the projected would have been delayed by several months. The choice of a plastic material as a substrate has disadvantages. The thermal conductivity of plastics is quite low. Fortunately, the temperature of the Martian polar

93 region is only 150 K, which is not difficult to reach using 77 K liquid nitrogen as a cryogen. Because we made hydrate from a liquid, the cell had to seal at a few hundred psi, as well as at a lower pressure at the 150 K measurement temperature. This provided a significant challenge. We experimented with several high quality epoxies, none of would stick to polyethylene or TPX at low temperatures. One epoxy stuck to silicon and the copper cell at low temperatures, but the long curing time (-1 day) and the even longer time required to dissolve the epoxy in trichloroethylene (-3 days) made the use of such epoxies impractical. Fortunately, if polyethylene windows are used, the sealing problem solves itself. When a polyethylene piece is sandwiched between two pieces of copper and sufficiently squeezed, we found that the seal did not leak at liquid nitrogen temperatures. If the polyethylene piece is sufficiently thick, the same gasket could be used several times without failing. We therefore decided to use polyethylene as a window, a sample substrate, and a gasket. We tried using TPX as a gasket as well. If it had worked, it would have been preferred over polyethylene because it is transparent, thus allowing visual inspection of the sample. Unfortunately, the TPX tended to crack at low temperatures. Our TPX windows cracked in about half the runs we made, and a window that survived one run might crack in the next run even though the setup was exactly the same and the temperature was lowered more slowly. The final cell design is shown in Figure (4.12). It consists of copper clamps holding a pair of polyethylene disks over the ends of a short copper pipe. The disks serve both as windows and as gaskets. The lower, thicker window is also a substrate on which the hydrate sample is grown. We milled 0.5 mm from the center of this window to form a shallow cup to hold the water. It was necessary to keep the water away from the walls of

94 thin polyethylene windo short copper pipe short copper pipe copper clamps gas inlet iL cryostat coldfinger thick polyethylene window/substrate water/hydrate sample Figure (4.12): Exploded View of Far-infrared Cell

95 the pipe. In a previous design iteration, the water was in contact with the copper, and it flowed up the copper as it cooled and froze. The screw threads holding the cell together and to the cold finger had to be strong because the cell had to be clamped very tightly to prevent leaks. Furthermore, the screws had to have a higher coefficient of expansion than the cell so that the squeezing would remain tight as the temperature was lowered. Stainless steel screws were therefore used because of their strength and their favorable thermal expansion properties. The largest uncertainty in this measurement is due to uncertainty in the sample thickness. We determined the thickness from the amount of water in the cell. the geometry of the cell, and the expansion properties of the water to hydrate phase change. At 0IC, the C02 hydrate has a cubic unit cell containing 46 water molecules with a lattice parameter of 1.207 nm. Using the molecular weights of the water molecules, we calculate an empty hydrate density of 0.80 g/cm3. Liquid water thus expands by 25% when forming a structure I hydrate. Using a micropipet, we placed 40 p1 of distilled, deionized water on the substrate. Based on the dimensions of the cell, this corresponds to a slab of water 0.45 mm thick. When converting to a hydrate, the sample should expand to 0.56 mm. We somewhat arbitrarily place a generous 0.05 mm uncertainty on the thickness. To test the evaporation rate of the water, we closed the cell, left the C02 port valve open, and let the system sit for 16 hours. We observed no measurable loss of water. We therefore assumed that no water was lost to evaporation in the relatively short time of about one hour during which the hydrate was prepared. We performed extensive experiments to determine a feasible method for preparing a high quality sample. Initial tests were made with a plexiglass window replacing the upper polyethylene window so that we could see the state of the sample. The cell was loaded with water, attached to the coldfinger, and pressurized with C02 gas to 415 psi. The pressurized cell was slowly cooled toward 00C. At about 30C the water froze, which is in

96 reasonable agreement with published hydrate phase behavior (Berecz and Balla-Achs, 1983). The existence of the hydrate was further confirmnned by reducing the pressure and observing the sample bubble profusely and melt. The solid, pressurized sample at 3 C was clear, but, as the temperature was lowered a few degrees below 00C, the sample became cloudy. This behavior suggests that only a thin layer at the top of the sample had converted to hydrate, and that the liquid below it was freezing into a cloudy sample of water ice at the lower temperature. To confirm this hypothesis, low resolution, far-infrared transmission measurements of the pressurized sample were made over the range -100C to 00C. At temperatures above about -40C the sample was quite opaque, indicating the presence of liquid water. At lower temperatures, the sample showed a spectrum with absorption increasing with frequency, which is characteristic of water ice (Mishima, et al., 1983) and clathrate hydrate (Gerbaux et al., 1975) far-infrared spectra. The hydrate layer forming at the top of the sample slowed migration of the C02 gas to the underlying water. To effectively cool the bottom first, thin strips of copper foil were run across the substrate and thermally joined to the short copper pipe. Samples prepared in this manner did not become cloudy and were not opaque to far-infrared radiation in the range -100C to 00C. Clathrate samples were cooled to -5 0C, and the pressure was reduced to 300 psi to prevent liquid C02 from condensing. The cell was then sealed and cooled to 150 K, the temperature of the Martian pole in winter, and another spectrum was measured. The spectrometer was configured with the mercury lamp, a 100 glm mylar beamsplitter, and the silicon bolometer operating at 4.2 K. The spectrum was ratioed with the spectrum of the empty cell at the same temperature to form a transmission spectrum. Significant amounts of C02 ice did not accumulate in the spectrometer beam path because there was no evidence of the 68 cml C02 ice lattice absorption. We will discuss the results of the measurement in the next section.

97 4.6: Measurement Results and Discussion The transmission spectrum of the clathrate hydrate at 150 K is given by the solid curve in Figure (4.13). The decrease in transmission at higher frequencies is characteristic of the low frequency tail of the broad water ice and hydrate lattice absorption peak with a maximum above 200 cm-l. There is a clear low frequency absorption from about 22 cm'l to below the minimum frequency our system could reliably measure. The position and the weak two lobe structure are reproducible between measurements. The absorption is not due to rotational transitions of the guest molecule because C02 has no permanent dipole moment. We attribute the absorption to transitions between guest translational modes. The frequency band of the absorption compares favorably to the transition frequencies predicted by the model in Chapter 3 for a C02 molecule in the larger structure I clathrate hydrate cage. The apparent feature at 37 cml moves with mirror speed and is an artifact caused by mechanical vibration, as discussed at the end of section 4.2. The data between 25 cm-l and 33 cm-l is unreliable because destructive interference in the 100 pm mylar beamsplitter causes the spectrometer to have a power minimum there. Crude measurements made with a thinner beamsplitter showed no features in this gap. Figure (4.13) also shows simulated transmission data of 0.56 mm water ice samples for comparison. The higher frequency curves are based on measurements made of 100 K and 200 K ice samples by Whalley and Labbe (1969). The transmission spectrum of an ice sample of 150 K ice presumably would lie between the 100 K and 200 K curves and would be similar to the hydrate sample. Since our hydrate transmission curve lies between the two ice curves, the assumed thickness of the sample is plausible. The lower frequency ice curve is based on measurements and theory by Mishima et al. (1983) of 150 K water ice. The ice transmission clearly does not exhibit absorption features in this region of the far-infrared. C02 gas and CO2 ice do not have strong

98 wavelength [mm] 1.0 0.4 0.2 1.0' ~- hydrate 150 K N E 0.5 ----- ice 200K cn I --- ice 100 K. -.-.- ice 150 K 0.0 I 10 20 30 40 50 wavenumber [cm ] Figure (4.13): Transmission Spectra of C02 Clathrate Hydrate with Water Ice for Comparison

99 absorption here either. (We know of no far-infrared data for C02 ice away from its lattice absorptions. This is probably because the attenuation is very weak and thick samples are difficult to make. See Warren's review (1986).) To model the emission from a Martian hydrate layer, we must know the complex hydrate index of refraction or, equivalently, the dielectric constant. The dielectric constants of structures I and II hydrates have been measured in the water molecule and guest molecule relaxation regions (Davidson, 1973, p. 183). For non-polar guests, the dielectric relaxation spectra asymptotically approach 2.9 with increasing frequency. Since there should be no strong absorptions in these hydrates at higher frequencies until the lattice starts absorbing, the dielectric value of 2.9 should extend into the far-infrared. We therefore assume the real part of the dielectric is 2.9, which corresponds to a real index of 1.7. We used the measured transmission spectrum to estimate the imaginary index. For this, we considered reflections at the hydrate-substrate and hydrate-gas interfaces. We determined the index of refraction of the substrate to be 1.6. Since this is so close to the index of the hydrate, we neglect reflections at the clathrate-substrate interface. Furthermore, the reflection at the clathrate-gas interface in the sample measurement is similar to the reflection at the substrate-gas interface in the reference measurement. We thus assume that these losses cancel out when the raw sample and reference spectra are ratioed to form the transmission spectrum. These simplifications lead to a transmission given by — 4o nid T= e (4.1) where a is the free space wavenumber, n is the imaginary index, and d is the sample thickness. The imaginary index is obtained by In (T) ni = - (4.2) n1 4i:o'ad

100 Figure (4.14) gives the spectrum of the imaginary index for the hydrate at 150 K as well as that of ice, as reported by Mishima et al. (1983). The dotted lines above and below the hydrate curve are error bounds resulting form the 0.05 mm uncertainty in the hydrate sample thickness. Whether or not the observed absorption is sufficiently strong to provide an unambiguous radiometric signature depends on the physical characteristics of the hydrate layer being measured and the atmosphere. In the next few chapters we present layer emission and atmospheric absorption models in an attempt to predict the emission signature of hydrate as well as water ice.

101 wavelength [mm] 1.0 0.8 0.6 0.4 0.02, hydrate C 0.01- 1, E 0.00 10 15 20 25 30 wavenumber [cm ] Figure (4.14): Imaginary Index of Refraction of C02 Clathrate Hydrate with Water Ice Index for Comparison

CHAPTER 5 QUASI-CRYSTALLINE APPROXIMATION 5.1: Introduction This chapter discusses the quasi-crystalline approximation (QCA) and how it is applied to the problem of predicting the effective propagation constant of a dense medium consisting of hydrate or ice particles. The use of QCA in geophysical problems was popularized by Tsang, Kong, and Shin (1985), and an exhaustive derivation of the method is given in their book. The purpose of this chapter is not to rehash their complicated equations and methodologies, but instead to present the theory of QCA in a form that someone with a basic background of electromagnetics can (I hope) grasp intuitively. This more physical approach should allow for greater insight into the strengths and weaknesses of QCA. Section 5.1 gives a derivation of the average field in a dense medium. Section 5.2 applies Foldy's approximation (FA) and QCA to the average field to give simplified average field expressions. Section 5.3 gives a physical discussion of FA and QCA based on a Monte Carlo simulation. Section 5.4 outlines how QCA is numerically implemented to give the effective indices and average fields. Section 5.5 discusses the particle pairwise distribution functions commonly used in conjunction with QCA. Section 5.6 gives some results of QCA applied to hydrate and ice particles. Some of the details of QCA, such as application of the translation theorem of vector spherical waves, are properly left to the book by Tsang et al (1985). 102

103 5.2: Dense Medium Scattering Equations Consider the problem of scattering in a dense medium consisting of N particles immersed in a background medium, as shown in Figure (5.1). The field EaPP is applied to the medium and scatters to form the total field E!o. N -tot — pp E;+ x ca (5.1) i= 1 In the above equation, kca is the scattered field from the ith particle.l Eca can be determined from the field Fk' incident on the ith particle by applying the appropriate boundary conditions. We define the transition operator Ti that operates on Einc to produce the scattered field from the ith particle.2 ~.c a -Y-nc tca = Ti (5.2) In general, Ti depends explicitly on the position, shape, size, orientation, and composition of the ith particle. It does not explicitly depend on the other particles. In this discussion we will only consider spherical particles with identical indices of refraction. Therefore, Tj only depends on position and radius. Now consider the incident field Enc. It is equal to the sum of the applied field and the fields scattered by all the particles except the ith particle. 1. Tang et al. use Green's functions instead of fields in their derivation of QCA. Fields are used here because of their more intuitive feel. 2. In Tsang et al.'s derivation of QCA, a transition operator T is introduced that operates on the inci dent field to produce an equivalent source to produce the scattered field. This is later replaced by a T-matrix, which operates on the incident field to give the scattered field. Tsang at al. may use this approach for historical reasons. Here instead a transition operator operating on the incident field to give the scattered field is used from the start.

104 rapp X T~~~~~~~.o T4 00 Figure (5.1): Scattering and Absorption in a Dense Medium.

105 N inc = iapp+ X TCa (5.3) j=1 j.* Substituting (5.2) into (5.3) gives a system of equations in the incident fields only. N z. nc =app + j nc5 j = I j=1 j*i This is a systems of N equations in N unknowns and can in principle be solved. By substituting (5.2) into (5.1), we get an expression for the total field in terms of the incident fields. N ot = app + T_-nc (5-5) i= 1 By recursively substituting (5.4) into (5.5), we get an infinite expansion for the total field. N N N Et0 =Epp+ S Tipp + Ta, 7T)Ea +.. (5.6) i= 1 i= lj= 1 It is informative to examine the physics of (5.6). The first term on the right hand side is the unscattered field. The second term is the field that has been scattered by exactly one particle. The third term is the doubly scattered field, etc. It is possible that under some conditions, particularly when the particles have no losses, the above series expansion will not converge. We will not address this question, but will instead assume that the series converges. In practice, the series has appeared to converge in our simulations when the particles had some finite (though at times very

106 small) loss. In cases where the particles were lossless, convergence was generally slower. However, no detailed study of the convergence properties of the series was undertaken. Many dense media inr nature, such as snow, sand, or molecules in a fluid, have random particle positions. Such complicated media usually require statistical treatment, and the average field is of interest. Applying the expectation operator E to (5.6) gives the average field. N N N E[E'~'] = aPP + E Ti] E PP + E [TTj] EaPP +..(5.7) i= 1 i= lj= 1 jri In the next section we will look at ways of simplifying the average field in (5.7) through approximations. 5.3: Average Field Approximations 5.3.1: Foldy's Approximation Consider the third tenn on the right hand side of (5.7). N N A E[ TiITjJ jPP (5.8) i= lj= 1 jei This term corresponds to the average doubly scattered field. If the positions of the ith and jth particles are not independent, then we cannot separate the expected value of the product of their respective transition operators into the product of the expected values of the operators. E [ TiTj] E [ Ti] E [ T] (5.9)

107 In some cases, particularly when the medium has a sparse distribution of particles, the particles positions are approximately independent, and we can apply Foldy's approximation (FA). FA can be written concisely by replacing the "not equal to" sign in (5.9) with an "approximately equal to" sign. We will introduce a bit more formality and m.- nery in getting to the approximation. While this does not develop much new insigh: the approximation, it does motivate the presentation of the quasi-crystalline approximation (QCA) in the next section. Since QCA is the method we will use in dense media modeling, the extra effort we expend here is worthwhile. We can rewrite the expected value in the doubly scattered term given by (5.9) as E[TiTj]l E[TjE[Tjiri] (5.10) Here E ri] means "the expected value of Tj given ri." (The expansion of the expectation of the product of two variables as in (5.10) is often referred to as "iterating the expectation.") If the positions of particles i and j are approximately independent, then E[Tj ril] E[Tj] (5.11) This is a formal way of expressing FA. It is interesting and illuminating to see what the total average field looks like when FA is used. By iterating the expectations in (5.7), we get the following expression for the average total field.

108 E[F~t] = p+ E E[T EPP+ E E[TE[TJ ri,]]aPP i= 1 i= lj= 1 Nj i (5.12) N N N + E E[TjE[TjE[T ri, rj]r]]PP + i= lj= lk= 1 j i ktj Applying (5.11) in (5.12) gives the following approximation for the total field. N N N E[r~t] aPp+ Z E[ T PP + E[Ti] E[Tj]'PP i= 1 i= lj= 1 jNi (5.13) + C ] E[T]i E[Tj] El[T ZIaPP +... i= lj= 1 joi This is the average total field expansion under FA. We see that an infinite number of scatter terms are present. However, since the approximation neglects the effects of particle position interdependence, the multiple scatter terms do not match the exact theory in (5.7) and (5.12). Only the no scatter and single scatter terms match the exact theory. For dense media, particle positions become strongly interdependent. In particular, particles cannot overlap, and this is not included in FA. Clearly we must include some correlation between particles in a good dense medium scattering model. The QCA model introduced in the next section includes pairwise dependence of particle positions. 5.3.2: Quasi-crystalline Approximation (QCA) We consider again the average total field expression given by (5.12). In particular, consider the expectation in the triple scatter term.

109 E[TE[TE[Tl rig rj] rj] (5.14) If the positions of the first and the third particles (i.e., the kth and ith particles) in the triple scatter are approximately independent, we have E[ T rim rj] = E [rji] (5.15) This is the quasi-crystalline approximation. Comparing (5.15) to (5.11), we see that QCA is a one order higher approximation than FA. To see the effect of applying QCA on the average total field, we apply (5.15) to (5.12) to get N N N E[E'~'] = E' + E[Tifl Ef + E E[TjE[T ri]]El ] i= 1 i= lj= 1 joi N N N + E E[T)TjE[Tkrl rj]] i j= kl = 1 joi k*j N N N N < AE[TTjE[ IE[Tlr ]-PP+ k r(5.16) i = lj= = I = 1 j~i k j Ilk As in the case of FA, QCA contains terms corresponding to all orders of scattering. Since QCA makes an extra order of correlation approximation, the double scatter term as well as the single and no scatter terms in (5.16) match the exact theory in (5.12).

110 5.4: Physical Interpretation of Foldy's Approximation and QCA In this section we describe a method for detennining the average field given by FA and QCA through Monte Carlo simulations. This method is probably too inefficient to be used in practice, but it does give more insight into the assumptions used by the two approximations. Figure (5.2) shows the initial setup of our thought experiment. Consider a volume V of space in which we will place particles, which for now we keep in a bag outside the region of interest. For simplicity, we let the background medium of the volume be free space. We let a field EPP be incident on the volume. The field in the volume is then quite trivially just ErPP. This is the unscattered field, as given by the first tennrm in (5.12). To simulate the singly scattered field, we draw a random particle Pi from the bag, place it at a random position in the volume V, as shown in Figure (5.3), and calculate the total field El. The singly scattered field E1 is just the total field minus the applied field. 1 = otl. faPP (5.17) Both FA and QCA calculate this field correctly. To simulate the doubly scattered field, we allow only the singly scattered field El to excite the medium. Under FA, we put particle Pi back into the bag and draw a particle out. We put this second particle Pj (which could be the same as particle Pi) at a random place in V, independent of the position that Pi had occupied, and calculate the total field E' resulting from the new applied field. The doubly scattered field ~2 is then given by E2 = 2~ - E1 (5.18) Figure (5.4) shows a possible position of Pj along with the position of P,.. Since the positions of the two particles are completely independent, it is possible for the two

111 rpp bag of particles volume V Figure (5.2): Particle and field arrangement for calculation of unscattered field

112 kapp 0gofprilsvlmV bag of particles volume V Figure (5.3): Particle and Field Arrangement for Calculation of Singly Scattered Field

113 particle Pi was put back in the bag before particle Pj was drawn Q -aPP = t E0 =Ei12+E1 bag of particles volume V Figure (5.4): Particle and Field Arrangement for Calculation of Doubly Scattered Field using FA

114 particles to interpenetrate, which is not physical. Furthermore, since Pi and Pj can be the same particle, it is possible for one particle to scatter from different places. QCA does not suffer from any of the double scattering inadequacies that are inherent to Foldy's approach. To simulate QCA double scattering, we draw a new particle Pj from the bag and place it at some position r1 that depends on the position ri of P and is chosen according to the conditional probability p (rJ ri). This is illustrated in Figure (5.5). Given a reasonable p (rj ri), the two particles do not interpenetrate. We then put Pi back in the bag, calculate the total field E2' excited by El, and get the doubly scattered field E2 via (5.18). QCA calculates the doubly scattered field correctly. To see where QCA begins to fail, we consider the particle arrangement for the triply scattered field, as shown in Figure (5.6). The doubly scattered field E2 is now the exciting field, and we draw another particle Pk (which could be the same as Pi) from the bag and place it in V at some position rk that depends on rj, but that does not explicitly depend on ri. As a result, Pi and Pk can interpenetrate, and if they are the same particle, they can scatter from different places at the same time. To simulate higher order terms in the FA and QCA total field expansions, we would continue to draw and replace more particles. We would repeat this experiment over and over again to get an average field. 5.5: Implementation of QCA 5.5.1: Incident Field Equations When implementing QCA, it is easier to work with the fields incident on the particles, rather than the total field. If we apply the expectation operator to our expression (5.5) for the total field in terms of the incident fields, and if we iterate the expectation, we get

115 particle Pi was put back in the bag after particle Pj was drawn pp = E2+El ag of particl2 es volume V bag of particles volume V Figure (5.5): Particle and Field Arrangement for Calculation of Doubly Scattered Field using QCA

116 Particle Pj was put back in the bag after particle Pk was drawn. Particle Pi was put back in the bag after particle Pj was drawn but before P. was drawn. 0 \\ bag of particles volume V Figure (5.6): Particle and Field Arrangement for Calculation of Triply Scattered Field using QCA.

117 E[E"'] = PP + E [TiE[EiC ri]] (5.19) i= 1 The conditional expectation E[i cI ri] can be expanded by applying the conditional expectation to (5.4). N E[ECi ri] = +PP+ E[TjEIj ri] (5.20) j=1 j.i By iterating the expectation, we get N E [njri = EaPP+ E[TE[C r,,rj r,] (5.21) j*i Now consider one of the conditional expectations in (5.21). E [ Nlri, rj] (5.22) QCA says that if a wave scatters off particles k, j and i in succession, the positions of particles k and i are independent. In (5.22) we have a wave E incident on particle j that is the scattered field from particle k and other and other such particles (and also includes the applied field, which is independent of all particles.) The wave scatters off particle j and then off the third particle i. Under QCA, E' is approximately independent of the position of the third particle (i.e, particle i)1. Therefore, 1. The argument that (5.22) is equivalent to the QCA approximation given in (5.15) is not trivial and is perhaps a bit nebulous. To verify that it is the same, one can recursively plug (5.24) into (5.5) to reproduce the QCA approximation for the total field given by (5.16). This is left as an exercise.

118 E[ ri r,,rj]= E[lI "r-] (5.23) Plugging (5.22) into (5.21), we get N E[Enc r] = + E[TE[ irjl r] (5.24) j=1 jsi We will use this equation when implementing QCA numerically. 5.5.2: Effective Propagation Constant and Average Fields This section describes how to use QCA to calculate the effective propagation constant kiff of a dense medium. The method requires some restrictions on the joint probability distribution of particle positions. Let p (ri, rj) be the joint probability distribution for the positions of the ith and jth particles. Using Bayes rule, we can expand the distribution as p (r,, rj) = p (rjf ri)p (ri) (5.25) We assume that the marginal distribution of the position of the ith particle is uniform over the volume of the dense medium. p(ri) = V (5.26) Furthermore, we assume that the conditional probability density p (ri ri) is only a function of the distance between the two particles. p (rIri) = p (r3- ri) (5.27)

119 The method of finding keff that follows is essentially the same as that given by Tsang et al., but with many of the details omitted for clarity. For example, in their derivation, fields and operators are expanded in terms of vector spherical waves and the translation theorem for these waves is used to related incident fields at different particles. While their explanation is complete, it is difficult to follow because very complicated expressions and formulations are introduced. The point of this section is to give a simplified explanation that allows the reader to see the spirit of the derivation. The details are in sections 6.3.3 and 6.6.1 in their book. We consider the problem of an applied plane wave with propagation vector ko and polarization p obliquely incident on a semi-infinite half space filled with dielectric spherical scatterers. The applied field is then given by aPP = Ee ir (5.28) Using QCA, we will outline the calculation of the average field propagating down into the medium and briefly discuss the average field reflected from the medium. We will find that the medium can be replaced by an effective lossy dielectric that will contain an average evanescent transmitted wave and an average specularly reflected wave. Consider equation (5.24) from the previous section, which we repeat here for convenience. N E[Enc rL= + E[-TE[E rjl r,] (5.29) jri Define w (rd) as the expected value of the incident field EiC evaluated at the position ri given that the ith particle is at ri. w (ri) = E [EZTC (r = i | r]i (5.30)

120 If the dense medium is behaving like an effective lossy dielectric, we would expect that w (ri) will be an attenuated plane wave with some effective propagation vector. We shall see that this is indeed the case. Using (5.30) and (5.28), equation (5.24) can then be written as w (r) = EOpei T+ E [ T (rj)w (rj) ri] (5.31) 1=1 jri Expressing the expectation in (5.31) as an integral, we get w(ri) = Eoek~ + E fdrTj (rj)(rj)p(rjri) (5.32) j= 1V joi Using our restriction on the conditional probability density of the particle positions given by (5.27), equation (5.32) becomes w(ri, = Eopeik'+ JdTj(rj) w (rjp (rj - ri) (5.33) j= 1V jti In view of the random nature of the medium, we postulate that the field w (ri) is also a plane wave with some effective propagation constant keff. w(ri) = ae (5.34) Plugging (3.31) and (5.34) into (5.33), we get Iikr = EP iko r+ A ldrjT (irj) ae i) (5.35) j= = jt i

121 Defining the pair distribution function g by p(rj-ri) = (5.36) we can insert (5.36) into (5.35) to get aei = eikozi + ldr(r) aeiffZg (r- ri) (5.37) i=1 v joi Assuming a large number of particles and that all particles are the identical, this simplifies to - A = krOz, r.. - - ikfflz aeik e +no dr I T (r - ri) (5.38) V where no is the particle number density. In general the particles can be of different types. We assume there are m different particle types with number densities nj, n2,..., nm. Equation (5.37) then becomes -ikrtlt,, ~ikozi - _ ik./z -3_ ae = e + njnkdrTj (r)ae gjk (r - ri) ( j=lk= V where gjk is the pair distribution function relating a particle of type j to a particle of type k. Consider the operator Aj given by Aj = njldrTj (r) e ik, (r - r,) (5.40) Tsang et al. show that this operator can be expressed in the form Aj(keff) = Bj(kff)e + C (kff) ei=klZ (5.41)

122 Substituting (5.40) and (5.41) into (5.39) gives ae' = ye'ik + nj (ff) e + Cj (k) ei (5.42) ae ye~rl = 9e'~"Z' + n. B (kffl)e + C (k~ff) e a'~]h(5.42) j= 1 The solution ae for a wave propagating in the medium due to a applied wave yei ~z0 corresponds to the medium extinguishing the applied wave and creating a new wave. All coefficients of e ikz on the right hand side of (5.42) must sum to zero to satisfy the equation. This expresses the extinction of the applied wave. m 9+ njBj (kff) a = 0 (5.43) 1=1 The above equation is known as the Ewald-Oseen extinction theorem. By equating coefficients of eikf' on each side of (5.42), we get the generalized Lorentz-Lorenz law. ( nj Ci(kff) -I)a = 0 (5.44) m a has a non-trivial solution if the operator A njCj (keff) - I is singular. It is singular if j=1 and only if its determinant is zero. det njCj(kff) -I = 0 (5.45) J=1 The above non-linear equation is solved to find keff. Solving (5.45) is not necessarily an easy task, particularly when particle dimensions are close the wavelength. In this case, twenty or more terms must be included in the spherical wave expansion of the fields and operators. The computer code developed for this research can handle an arbitrary number of expansion terms, the speed of the code drops precipitously as the number of terms increases. The code uses Muller's method to

123 solve (5.45) iteratively, as suggested by Tsang et al. Muller's method requires three initial guesses, and the quality of these guess can determine whether or not convergence occurs in a reasonable time. Solving for kff via (5.45) gives us an effective propagation for the dense medium, but leaves the Lorentz-Lorenz law in (5.44) as a system of equations with one extra degree of freedom. This extra degree of freedom is constrained by the Ewald-Oseen extinction theorem in (5.43), and the fields can be computed. They obey equations similar to those of an ordinary dielectric. The average propagation vector obeys Snells law (keff)x = (ko) X (5.46) The average field scattered up is only in the specular direction and has an amplitude proportional to the incident field through the effective reflection coefficient r, which depends on the polarization. ko tan O - kefftan Oeff rh= kOtanOO + kefftanOeff (5.47) ko tan o+ keff tan0e r = kotan0ff- ksfftanO~ (5.48) k ktan fff+ kfftan El We give the caveat that the angle Eeff is slightly complex. We have now specified the field everywhere. The only remaining loose end to tie up is the issue of the pair distribution function used by QCA to get kff. We will discuss pair distribution functions in the next section.

124 5.6: Pairwise Distribution Functions The pairwise distribution function chosen in a scanttering model can have a strong effect on the effective propagation constant calculated. For sparse media, a uniform distribution may be suitable when the forces between particles are short range, as in the case when the particles are simple hard spheres. An improvement over this distribution is the hole-correction function, which is zero when the particles overlap and unity otherwise. These two distributions are suitable for circumstances where the particles are sparsely distributed, as are molecules in a gas. However, when the volume density of the particles becomes significant, the particle positions become correlated, and these two simple functions are no longer suitable for use in scattering models. The Percus-Yevick pair distribution function (PY) was designed to handle the case of dense particle distributions. It was originally applied to the problem of molecules in a fluid, which are mainly under the influence of interparticle forces. PY not only takes into consideration the direct force on one particle by the other particle in the pair, but also takes into account the indirect forces of all the other particles. As a result, it does a significantly better job than the simpler pair distribution functions in predicting scattering. The use of PY in QCA has been used by Tsang et al in geophysical applications. They use a hard sphere potential to describe the forces between the particles, which produces a pair distribution function that is sufficiently simple to calculate that it can be used in efficient scattering computer programs. They report reasonable results in their comparison of QCA PY model to actual snow measurements. To help verify the suitability of PY and QCA in dense medium scattering problems, Tsang et al have performed Monte Carlo computer simulations. Spheres were randomly placed in a box until the desired particle volume density was achieved. The effective propagation constant of the medium was then computed numerically, and the

125 problem was repeated many times to get an average propagation constant The results showed reasonable agreement with the propagation constant predicted by QCA using PY Because of the method used in assembling each particle realization, it is not surprising that PY worked in this case. The only forces acting on the particles were the hard sphere forces. In geophysical media, gravity can have an important effect on the pair distribution function. Since gravity is normally unidirectional, geophysical media will typically have a pair distribution function that is not spherically symmetric. Despite its drawbacks, PY is still a useful approximation to the spatial statistics of a dense medium. It satisfies the physical constraint of impenetrability, and it takes into consideration the indirect effects of particles outside the particle pair. When the pair distribution function of a dense medium is not known or is too difficult to measure, PY supplied a reasonable guess that can be used in simulations. Since it is impractical for us to measure the pair distribution function for Martian snows, we will use PY in our simulations. 5.7: Propagation Constants of Hydrate and Ice Particle Layers Using the QCA model described in this chapter and the dielectric data from chapter 4, we estimated the effective propagation constants of C02 clathrate hydrate and water ice particle layers. In general, the effective propagation constant keff depends on the particle size relative to a wavelength. In the resonant region, we found that the kff is an oscillating function of the size. To reduce this effect, we assumed a Rayleigh size distribution of sizes with 11 discrete sizes. Kieffer (1990) modeled the water ice particles in the Martian northern polar region and estimated them to be 100 gm or more in diameter. We used average particle diameters ranging from 50 to 400 plm in our Rayleigh distribution. The real and imaginary parts of

126 the effective index of refraction n, + ini are shown in Figure (5.7) for the case of a 100 pgm average particle diameter. The effect of the hydrate absorption feature is not apparent in the spectrum of ni. This is because the loss depends on both absorption and scattering. As frequency increases above the range of the hydrate feature, absorption decreases, but scattering increases. The net effect is a increase in the total loss. What is more important in radiative transfer models is the amount of loss due to absorption and the amount of loss due to scattering. We will see in the next chapter that the effective single scattering albedo is lower in the region of hydrate absorption and higher away from the absorption. This will have a clear effect on the modeled emission. The high frequency ends of the n, spectra show small oscillations. These are probably due to resonant scattering effects that were not completely smoothed by our particle size distribution. The oscillations are much larger when a single particle size is used. To reduce the oscillations, we could have included more sizes in our discrete Rayleigh distribution. However, this rapidly increases the computer time and memory required to calculate the effective indexes. Our choice of 11 particle sizes is a reasonable compromise between realistic simulation and computational expediency. We have not included dust in our model. Currently there is some discussion over how dirty the northern polar ices are (Kieffer, 1990). In any case, the dust particles are believed to be much smaller than the ice or hydrate particles and should scatter less. There composition is mainly oxides of iron, which are unlikely to have significanat absorption in the far-infrared. Dust may have a significant effect on the far-infrared spectrum below 25 cm1 if it is abundant enough to change the background dielectric in which the ice and hydrate particles reside, or if it coats these particles. It is a weakness of this thesis that dust could not be well-handled.

127 2.5 2.0 1.51- W 3 1.0 d L o. - -:: —- hydrate nr ice nr 0.5 - hydratenixl _ - e- - --- ice n. x 10 0.0 14 16 18 20 22 24 wavenumbers [cm'] Figure (5.7): Effective Complex Indices of Refraction of Dense Media Consisting of CO2 Hydrate or Water Ice Particles. Average particle diameter is 100 Wtm, and volume fraction is 0.3.

CHAPTER 6 EMISSION MODEL FOR HYDRATE AND ICE LAYERS Chapter 4 described the measurement of the far-infrared C02 clathrate hydrate refractive index, and chapter 5 presented a model for estimating an effective propagation constant in a hydrate particle layer given the hydrate index. This chapter describes how the index and effective propagation constant data can be incorporated into a radiative transfer model to predict emission from a layer of hydrate particles. The same model is used to predict emission from an ice particle layer. Figure (6.1) illustrates the problem of modeling a semi-infinite particle layer. Emission, absorption, and scattering must be included in the radiative transfer model, and boundary conditions must be considered. The results of the model show a large spectral gradient in the hydrate emission spectrum that is not present in the ice emission spectrum. In the next chapter the results of the emission model are coupled to an atmospheric propagation model to predict the spectrum that would be observed in a radiometric measurement. 6.1 Radiative Transfer and the Quasicrystalline Approximation Classical radiative transfer theory is based on an equation of the form m = v, h d-(,). __Kex t (r) + Kext(.r )l d, P n m(sip'1;kss ) ( 28 128

129 incident waves z=0 boundaryc O 0 0 0O Z c~ o emission scattering and absorption Figure (6.1): Modeling of a Semi-infinite Layer

130 In (6.1), 1, (r, 3) is the specific intensity with polarization n at position r in the direction 3, Kt,, is the extinction coefficient, Co is the single scattering albedo, Pnm (s, 3 ) is the scattering phase function from direction I and polarization m to direction 3 and polarization n for a particle with propagation constant ks and size rs, and B (r) is the emission source term. Radiative transfer has been used successfully to model propagation and emission in media with a sparse distribution of particles, such as an ice cloud. When the particles become densely packed, however, classical radiative transfer can greatly overpredict extinction due to scattering. This can result in much lower modeled emission than is experimentally observed. One way to circumvent the overestimation of scattering is to use an effective particle volume fraction that is significantly smaller than the true physical volume fraction. By choosing an appropriate effective value, radiative transfer models have successfully matched measured snow data. However, this non-physical method is not useful in predicting emission because the effective volume fraction is not determined until after the emission data is acquired. A more reasonable approach has been advanced by Tsang and Kong. Using the quasicrystalline approximation (QCA), which is discussed in the previous chapter, they obtain better predictions of scattering attenuation in dense media. As the volume of the particles in a dense medium increases, the particle positions become strongly interdependent. This interdependence leads to coherent multiple scattering effects that result in effectively less extinction of an incident wave. The QCA method takes into consideration the interdependence of particle positions, while classical radiative transfer does not. We use QCA to obtain an effective propagation constant keff and hence an effective extinction coefficient Ke. Chapter 5 gives a physically based description a QCA and an outline of how these parameters are calculated. More detailed discussions of

131 the calculations can be found in Tsang et al. The parameters are incorporated in the radiative transfer equation (6.1) to model emission. I (i = - cf() (I(, 3)) +KeAt (4) Iag | PnM (s, s;ks, rs) Im (r, 3 m = v,h + K/e (r) (1 - m) B (f) (6.2) 6.2 Rejection of the Distorted Born Approximation Before proceeding with the radiative transfer model, it is worthwhile to discuss an alternate model that has been used by many researchers in radar problems, but that we found to be unsuitable for our purposes. This alternate model is the distorted Born approximation (DBA), which is a variation of the Born approximation. We have found that it can give emissivities that are much lower than those typically found in nature. In its simplest form, the Born approximation consists of approximating the field incident on the particles in a medium as the incident field. This is modified in DBA by including extinction of the field due to absorption and scattering by the particles. The scattered field from a particle can then be computed directly from the incident field and the properties of the particle. The total scattered field is then the sum of the scattered field from all particles, and how these fields are added together depends on how particle positions are correlated. QCA predicts the average field inside the medium and thus supplies DBA with the input it needs. If we assume that the average field is the only field incident on the particles, then the problem of predicting the emissivity of a half space of particles can be addressed by considering a reciprocal problem, as shown in Figure (6.2). Instead of considering the intensity emitted from the layer of particles, we consider a plane wave incident on the slab

132 incident power scattered POwer 0~.. 009ii) O of O i (0 0 r'~ 0''..... o c0-C~ r O0 0..,I.............. ~~o oOo0o ~~0 0 ~o ~0 o o 0 Figure (6.2): Reciprocal i

133 from above, and for simplicity, we let the wave be normally incident. The emissivity e (OinC) of the isothermal slab is equal to the ratio of the power absorbed in the slab to the power incident on the slab. By power conservation, we could also calculate the total power scattered up and divide by the incident power to get the reflectivity r ((inc), which is related to the emissivity by r (inc) = 1-e (8ic). (6.3) If we were to solve the reciprocal problem rigorously, the two methods would produce the same value for the emissivity. DBA does not conserve energy, however, and will not give the same results. Furthermore, for dense media we have found that DBA can also give non-physical results. For example, using QCA to calculate the internal incident field and then using DBA to calculate the power scattered up can lead to reflectivities greater than unity. Similarly, using QCA and DBA to calculate the power absorbed in the medium has led to emissivities less that 0.01, which are seldom observed in nature for dielectric media. Because of its unrealistic or physically impossible predictions, we have chosen not to use DBA in our modeling of emission from hydrate and ice layers. We instead cast the problem in the form of classical radiative transfer, which gives us physically realistic results. 6.3 Effective Single Scattering Albedo Unfortunately, QCA does not give us any information on how much of the extinction is due to scattering and how much is due to absorption. This information is necessary for specifying the effective single scattering albedo Co in (6.2). We will estimate the contribution of absorption to extinction by the following simple method. Consider a wave incident on a thin layer containing N particles, as shown

134 Intensity I _u, 0 0.0 oOQ..:...........................-.................. Particle with Cross Sectional Area Ai Area A FiguSlab Are (6.3): Calculation of A Figure (6.3): Calculation of Absorption Coefficient

135 in Figure (6.3). The wave has intensity I and the total power incident on the slab is P. The slab has area A and thickness Az. The power AP absorbed by the slab is approximated by AP = IA Qbs (6.4) where Ai is the cross sectional area of the ith particle, and Qab is its absorption efficiency. Since the intensity is equal to the total power incident on the slab divided by the slab area, (6.4) can be written as AP = XA.i (6.5) If the slab contains M different particle types with number densities n, n2,..., nM, then AP can be rewritten as M P abs AP = P A njVAjeQby' (6.6) j= 1 where AJ and Qjabs are now areas and efficiencies of the jth type of particle and V is the volume of the slab. Since V = AAz (6.7) (6.6) simplifies to M AP = PAz A njAj Qbs (6.8) j= 1 The absorption coefficient is immediately recognized as M Kabs = - n-A-Qb (6.9) j=1

136 Then effective single scattering albedo is then given by Kabs o ff 1- (6.10) Kext We could have taken a similar approach using scattering efficiencies Qsca to find the scattering coefficient and hence a new scattering albedo. However, this would lead to an overestimate of the scattering because it would not take into account the coherent scattering effects that have been found to reduce extinction. 6.4 Simplification of the Radiative Transfer Equation Equation (6.1) still requires specification of a scattering phase function and an emission source term. An averaged Mie scattering phase function is used based on the particle size distribution and frequency. Let n (r) dr be the number density of particles with radii in the small range from r to r + dr, and let Pnm,, (s, s';k,, r) be the Mie scattering phase function for a particle with radius r and propagation constant ks. The average Mie scattering phase function is given by pave (,) j'dr n (r)pnm (';k, r) (6.1) Pnm s, ) (6.11) Jdr n(r) The emission source term B is in general a function of depth because temperature and particle characteristics are in general a function of depth. In the present case, however, the extinction turns out to be high enough so that far-infrared radiation can only penetrate to about 2 cm in the case of minimum attenuation. We assume that over this short distance, the temperature is uniform. We also assume that the particle dielectric properties and statistics are uniform with depth. If the layer is sufficiently porous, the particles a few centimeters deep should readily convert to COC hydrate when the tempertature drops low enough to stabilize

137 hydrate. Miller and Smythe (1970) grew CO2 hydrate from water ice frost at temperatures and pressures similar to those in the Martian northern polar winter. The time required for their thinly divided frost samples to convert to hydrate was on the order of a few hours. If the particles are thinly divided, as high frequency emissivity data suggests (Kieffer, 1990), then the long Martian winter should provide sufficient time for the water ice to convert to hydrate. Hence, kef, eff', Oeff, Pave' and B are uniform with depth. The physical temperature of the layer is taken as 150K in accordance with Martian observations and models. Inserting the effective parameters derived from QCA and the averaged scattering phase function into (6. 1), we get a new radiative transfer equation. TdSIZ' m = t xt 4 v, hn (6.12) +K ef ( 1 ( Beff)B Equation (6.12) can be further simplified. Clearly the intensity depends only on the depth z and polar angle 0. We can make the substitution in (r, s) JoIn (Z9 A) (6.13) where Cp = cosO (6.14) The differential ds can be related to dz by ds = d (6.15)

138 With these simplifications, the radiative transfer equation in (6.12) becomes 1eff 1 eff ) ff v, Z' d (Z-.) = —K + -K d, l) I (, ) /z, St) = Iext ext, = v,hx1 leff (6.16) +-K (1 - eff)B ex exf We now get rid of all azimuthal dependence in the equation. If we integrate the entire equation over 4 and divide by 2i, we get d I =1 K— I(z, t) = K dz' ext eff C I pae ( + CK 2 Pai (' dodo I;', ~') )I (z, ) eXt 2rl = v, hn 1 eff +-K (1 -o ff)B (6.17) ext The double azimuthal integral in parentheses can be evaluated to give a function q that depends only on polar angles. = 4 m *,Ith",ave (,, 04/;,, 4)) (6.18) m = v, h Equation (6.17) then becomes the simpler equation d (z' t) = -!Keft(z, ) +dl qnm ( )I(z, ) d ~ext ext 2 m=vh (6.19) +-Kez ( 1 2-,Oe ff ) B

139 To perform a final simplification, we introduce the optical depth parameter t defined by (6.20) Making the substitution In (z, g) - - In (~, g) (6.21) we get our final simplified radiative transfer equation J g =-In(, n + tt |d qnm (v, h ) In ( 6.22) m = v, h (6.22) + (1 -coeff) B 6.5 Solution of the Radiative Transfer Equation Following the approach of Chandresekhar (1951), we use Gaussian quadrature to approximate the scattering integral in (6.12). The integral is replaced by a weighted sum of integrand values at the quadrature angles. If the number of angles is even, then the angles occur in positive-negative pairs. =-i = -gi (6.23) Taking advantage of this symmetry, we can write the integral approximation as ldg' T qnm,(g,'). l,(z,') m = v, h N = aj( qnm(g, j)1I(,, j) + Z qnm(g'-~j)I(" — >)) (6.24) j= m=v,h m=v,h

140 where the weighting coefficient aj is known as a Cristoffel number. The method of Gaussian quadrature to approximate an integral is exact if the function being integrated can be expressed as an Mth order polynomial, where M is less than or equal to twice the number quadrature points used in the approximation. We assume that the specific intensity and phase function vary slowly enough so that the approximation is valid. No significant change in the results was present when more than about 20 angles were used. If we use the quadrature approximation above in the equation of transfer in (6.22), and if we only consider solutions at a quadrature angle gis then we get a discretized equation of transfer given by dzIn -(Z' i) In (Z, 1i) N + aj( E qnm (ll.j) I ('.9 j) +; qnm (l, -9j) I (x, -.j)) j = 1 m = v,h m= v,h + - (1 - O) B (6.25) (6.25) is a linear, inhomogeneous system of differential equations in 4N intensity unknowns. To solve the equation, we first do some work to get it in a more compact form. We define the following vectors and matrices.

a-l -1141 = ~21 m = (6.27) [m 0 M = m (6.28) [Pnm (iU41,jil).... Pnm (igil,jLN) ~ pi,.. (6.29) nm nt, m h, v Pnm[ ('RN'9 ill I)... Pnm (i5N, iJLN) =1, 1 =1, I1 =1,-1 =1,-1 Pvh Pvv P Pvh =-, 1 =-1, 1 =-1, -1 =-1, -1 LPIv hh Phv Phh B (r) = [B,..., B] (6.31) The system of differential equations can then be rewritten more compactly as - (X) = (-M+ MP) I () + (1 - co)MB (6.32)

142 Defining further the matrix A and the vector b, A = -M+ - MP (6.33) b= (1-o) MB (6.34) we have a very compact form of (6.25). I (t) = AI (X) + b (6.35) dt We will solve the above radiative transfer equation using standard linear algebra techniques. Consider the homogeneous equation of (6.35). 4-ro (') -AIo (I) = 0 (6.36) dt It has solutions of the form Io(t) = aves (6.37) where v is an eigenvector of A, X is an eigenvalue, and a is an arbitrary constant. The general solution to the homogeneous equation is 4N IO() = Eaivie i (6.38) i=1 Because of the physics of the problem, half of the eigenvalues have positive real parts and half have negative real parts. We rewrite (6.38) separating the sum into two sums, one of which includes terms with positive real part eigenvalues XP? and the other with negative real part eigenvalues X,.

143 2N 2N Io (X) = aive + iwe (6.39) i=1 i=1 The surface of the semi-infinite layer is chosen to be at z = 0, and the layer extends to negative infinity. We require the intensity to be finite as the depth goes to negative infinity. Therefore, all terms in (6.39) containing negative eigenvalue must vanish. Hence, [3i = 0, (1 < i < 2N) (6.40) and the homogeneous solution is reduced to 2N Io(~r)= a vie (6.41) i= 1 We now consider a particular solution I, (X) to the discretized radiative transfer equation (6.35). If Ip (X) is independent of t, then (6.35) reduces to Alp (r) = -b (6.42) A particular solution to (6.35) is then given by Ip,() = -A b (6.43) The total solution is the sum of the homogeneous and particular solutions given by (6.41) and (6.43). 2N I(T) = I, aivie A-A b (6.44) i= 1

144 Deep inside the layer, we expect the intensity to equal the source term B. As z - a, the homogeneous part of the solution vanishes leaving only the particular part. It turns out that the particular part is equal to B, and we can write (6.44) as 2N 1](X) =,aivie i +B (6.45) i=l 1 We now apply the boundary condition at the top of the layer. We assume a specular surface and compute Fresnel reflection coefficients based on the effective propagation constant kcff derived from QCA, the propagation constant ko of free space, and the angle of incidence. The reflection coefficients for vertical and horizontal polarizations are given by 2 2 2 RV (gi) = kogi -keff kerk(> 1 )2) (6.46) k|Iik k k- (1 _2) 0 I e 0 eff 2'22 Rh (d ~1) eff- k2k(1 i2) (6.47) keff Li + gO f - kf(1) - ) For convenience, we break the vector of intensities I (t) into an upgoing part iUP (t) and a downgoing part idown (t). ip (X) = [Iv (x, l1),..., IV (, gN) Ih (, g)... Ih ( ] T (N)] (6.48) and we ) [Iv ( -g IV (, -N) Ih ( -I )' Ih ( -define) ] reflection matrix(6.49) and we define a reflection matrix by

145 -RV (1) RV (EN) R = (6.50) Rh (gl) Rh (hN) The specular boundary condition can then be compactly expressed by Jdown (O) = RI"P (0) (6.51) We have assumed for the moment that there is no radiation incident on the layer form above. Such incident radiation is more conveniently treated separately, and we address it in a later section. We also split the eigenvectors into upgoing and downgoing parts Vi = [ (Vi) I, (Vi) 2,-", (Vi) 2N (6.52) (6.52) -down = [ (vi) 2N+ 1' (vi)2N+ 2' ". (vi) 4N] (6.53) and put their coefficients in a vector. a = [al, a2.., a'2N] T (6.54) = up =down We put the upgoing and downgoing eigenvectors into matrices V and V V = [vvup up, v (6.55) = [vI, v2,..., vN(6.55)

146 = down -down -down -down T V = [vl,...,V2N ] (6.56) Using (6.45), we can express 1P (0) and Idown (0) as?'P (O) = V p+B (6.57) = down Jdown (0) =V a + B (6.58) Plugging (6.57) and (6.58) into the specular boundary condition equation given by (6.51), we get = down =up V U+B = R(V a+B) (6.59) (6.59) can then be rearranged into a more convenient form =down ==up _ - _ (V -RV ) = (I-R)B (6.60) where I is the identity matrix. The linear equation (6.60) can be solved for ai. Substituting this solution into (6.45), we have a final solution for the intensity the layer. We strive to find the intensity emit emitted from the layer. To get this, we multiply the upgoing intensity IPP (0) by a Fresnel transmission matrix. ~emit = TI (0) (6.61) The transmission matrix is defined in terms of the reflection matrix by T = I- R (6.62)

147 For layer angles larger than the critical angle, the transmission matrix correctly contains a zero diagonal element. Cm"i is defined as ri (0) = [I='v (0, j),..., m,' (0, )N), hr (0 l)..., rehm (,O 1N)]T (6.63) where 1,mit (0, gi;) is the emission intensity at angle 1ti, and i. is related to the layer angle gli through Snell's law.'2 (6.64) For angles layer angles larger than the critical angle, the corresponding emission angles are not defined. In this model, however, the emission intensities at these angles are zero because the transmission coefficient is zero. This completes the derivation of the emission model for the layer. The results of the modeling are discussed in detail in Chapter 6. The scattering and absorption of atmospheric radiation by the layer will be discussed in the next section. 6.6 Scattering and Absorption of Atmospheric Radiation by the Layer The hydrate and ice layers affect observed radiobrightness not only due to their emission, but also due to their scattering and absorption of atmospheric radiation. Fortunately, the methods used to calculate emission can be used to calculate scattering and absorption of atmospheric radiation with only a few modifications to the boundary conditions. A strong caveat regarding the angular distribution of the incident atmospheric radiation is appropriate before proceeding. Since the modeling presented here relies on the accuracy of the Gaussian quadrature integration technique discussed in the previous section, it is important that the incident radiation vary slowly enough to be well

148 approximated by a reasonably low order polynomial. This is not the case for a unidirectional source of radiation, such as direct sunlight. A delta function cannot be well represented by a low order polynomial, and our method breaks down for the unidirectional case. Furthermore, our modeling assumes an azimuthally symmetric intensity distribution, which is also not satisfied in the unidirectional case. With these warnings in mind, we proceed with the model. The boundary condition given by (6.51) expresses the condition that the intensity coming down from the layer-atmosphere interface is equal to the reflection of the intensity impinging on the interface form below. When atmospheric intensity latm is incident from above, the downgoing radiation just below the interface also includes a term due to the transmitted part or "atm.?d""o () = RiP (0) + Tpfm (0) (6.65) The vector ]atm (0) is given by a~m (01) = [m (0N-).tm (h N)tm (, -_).tm (0, _. N) ] T (4.66) where Pntm (0, -L i) is the atmospheric intensity with polarization n incident on the layer at angle -4li. Angle -gli is related to the angle -.i inside the layer through Snell's law, as was discussed in the emission model of the previous section.For angles -pil corresponding to layer angles beyond the critical angle, we simply set atm (0, -L;) equal to zero. Since we are only concerned with the effect on incident atmospheric radiation, we do not include emission in this model. The equations for I"P (X) and Idown (X) given by (6.57) and (6.58) are changed to up n = v_ (6.67)u.

149 = down down (O) = V a (6.68) Plugging these into the boundary condition equation (6.65), we get =down UP. V a = RV +Tr' + (0) (6.69) Rearranging, we get a linear equation convenient for solving for a. =down ==up (V -RV )T = Tratm(O) (6.70) With the vector a, we can calculate Sup (0) using (6.57). The scattered intensity rca (O) can then be calculated using an equation analogous to the emission equation (6.61). rca= TI (0) (6.71) We have now completely characterized the layer through its emission and its response to atmospheric radiation. We can now calculate its emissivity over the region of the hydrate absorption feature. Its coupling to the atmosphere will be explored in the next chapter. 6.7 Emissivity of Hydrate and Ice Particle Layers Using QCA and the values of the index of refraction derived from the measurements, we can compute the emissivity of hydrate and ice particle layers. In general, the propagation constants depend on the particle size distribution. As stated in the last chapter, we have assumed a Rayleigh size distribution. We have considered distributions with average particle diameteres spanning the range 50-300 ptm.

150 The emissivity spectrum for the 200 4m average diameter case is shown in Figure (6.4). We see that the hydrate absorption feature causes a large spectral gradient in the hydrate emissivity spectrum. The hydrate emissivity has a large decrease with increasing frequency, while the ice emissivity increases somewhat. A reasonable method for differentiating hydrate from water ice would be to observe the difference in brightness at 14 cm-l and 24 cml1. In Figure (6.5) we have plotted the emissivity differencel between the two channels for several average particle diameters between 50 and 300 pm. We see that above about 60 plm, the difference between the hydrate and ice curves is more than 0.1 and that above about 150 glm, the ice curve becomes negative. Sensitive far-infrared radiometers should easily be able to discern such large differences in spectral gradient. A large emissivity difference is indicative of hydrate, while a small or negative difference is indicative of ice. Below we suggest a reasonable set of criteria for differentiation. 1) If the emissivity difference is greater than 0.07, then the hydrate presence is confirmed. 2) If the emissivity difference is 0.007 or less, then the ice presence is confirmed. 3) If the emissivity difference is between 0.007 and 0.07, then the data is inconclusive. With the above criteria, we can differentiate hydrate from ice if the average particle diameter is greater than about 150 pLm. It is possible that a telescope viewing the Martian polar surface would simultaneously view hydrate patches and ice patches, particularly as the temperatures 1. At 20 cmnl and 150 K, the Rayleigh-Jeans Law is not accurate. Brightness is not a linear function of the physical temperature, and we thus work with the emissivity instead of the brightness temperature. Brightness temperature can still be defined, but it less intuitive. In the 150 K far-infrared cases we consider here, a change in emissivity of 0.07 corresponds to a change in brightness temperature of about 10 K.

151 1.0 0.8 *> 0.6 0.4 0.2 0.0, I 14 16 18 20 22 24 wavenumbers [cm'l] Figure (6.4): Emissivity Spectrum of Hydrate and Ice Particles Layers with Average Diameter 100 pLm.

152. --- 0.1 hydrate _. —------ ice'- 0.0 -0.1 I I 50 100 150 200 250 300 average particle diameter [urm] Figure (6.5): Emissivity Difference between 14 cm-1 and 24 cm-l Channels for Hydrate and Ice Layers with Varying Average Particle Diameters

153 increase in the springtime and hydrate particles dissociate into water ice and C02 gas. We consider the emissivity differences for several cases in Figure (6.5). Large changes in the emissivity difference are still possible for sufficiently large particles when coverage goes from mostly hydrate to mostly ice, and the seasonal conversion may still be observable. Using our simple criteria above, we can differentiate 60% hydrate and 40% ice coverage from 100% ice coverage for average particle diameters above 150 gLm. A 20% hydrate and 80% ice ratio still leads to a negative emissivity difference for sufficiently large particles, but when the hydrate coverage reaches 40%, the difference becomes significantly positive, and the signature is ambiguous. The simulated emission of a single patch containing a mixture of hydrate and ice particles might be useful. We would expect the emissivity to be somewhere between the hydrate and ice emissivities. Unfortunately, applying the quasi-crystalline approximation to particles having different dielectric constants is much more complicated than it is for the single dielectric case. We have not addressed this mixing problem.

154 0.2 *,1 0.0 - 100% hydrate * -0.1 -------- 80% hydrate, 20% ice ------ 60% hydrate, 40% ice X- | — 40% hydrate, 60% ice -0.2 ~ - -.. 20% hydrate, 80% ice 100% ice -0.3. I 50 100 150 200 250 300 average particle diameter Figure (6.6): Emissivity Difference Between 14 cm'1 and 24 cml1 Channels for Hydrate and Ice Layers with Varying Average Particle Diameters and Coverage Ratios

CHAPTER 7 ATMOSPHERIC ABSORPTION, EMISSION, AND SCATTERING Radiation emitted by hydrate and ice in the Martian polar regions must pass trough the Martian atmosphere before being received by spaceborne or earth-based detection systems. We must consider the effect of absorption and scattering in the atmosphere on the surface emission, and we should also consider emission by the atmosphere itself. In this chapter we shall consider the atmospheric particles and gases that might affect the surface emission signature. We shall see that there are ample far-infrared spectral windows available for remote sensing. 7.1: Particles 7.1.1: Scattering and Absorption by Small Particles Gas molecules are quite small and are effectively incapable of scattering electromagnetic radiation with sizeable wavelength. Since dust and ice clouds are made up of macroscopic particles, they are capable of scattering infrared radiation, and this effect should be considered in a model of the atmosphere. We shall see, however, that the ice cloud and dust particles are too small and sparsely distributed to significantly affect the surface emission signature. We will also dispose of the absorption by ices and dust as an unimportant mechanism. In considering the scattering and absorption by particles, we shall show that even with overestimates of particle sizes and total volume, the effect on the surface emission is small. 155

156 In all cases, the particles are much smaller than a wavelength. It is therefore appropriate to use Rayleigh particle scattering and absorption models. We will assume spherical dielectric particles, which leads to the following expressions for the scattering and absorption cross sections (Ishimaru, 1978): 128t'ra 64 -r1 2 ___ 1287~5a6 E l, +2-1 i(7.1) 128a-2 CTa = rxi Er | i3 V (7.2) Here ac is the scattering cross section, ca is the absorption cross section, a is the radius, Er = Err + i~ri is the particle's complex dielectric constant relative to the non-absorbing background, V is the particle volume, and X is the wavelength. We see in (7.1) that the scattering cross section is proportional to the sixth power of the radius, or, equivalently, the square of the volume. Hence, if we fixed the total volume of the particles and make them larger (thereby decreasing the number density), we get more scattering. Therefore, if we know the volume of material is in a cloud, and if we overestimate the particle size, then we will overestimate the amount of scattering. We will overestimate both the total volume and the particle size and thus further overestimate the scattering. We see in (7.2) that the absorption cross section is proportional to the volume. For a fixed amount of material, the absorption is independent of the particle size, as long as the Rayleigh size criterion is satisfied. An overestimate of the total volume, however, will lead to an overestimate of the absorption. If the particle is too large to satisfy the Rayleigh

157 criterion, a Rayleigh scattering model will overestimate both the scattering and the absorption. The clouds on Mars are not a dense medium and should not require a dense medium propagation model such as the one discussed in chapter 5. We will assume the positions of the particles are independent, and the multiply scattered waves add incoherently. With this simple model we get the following expression for the scattering and absorption constants: Ks = asn (7.3) Ka = Can (7.4) The total attenuation factor fart is then given by fat: = e(+a) (7.5) where d is the length of the path through the uniform cloud. If the path is normnnal through the cloud, then the quantity nd is just the number of particles per unit area. We can write nd in the more convenient form nd = - (7.6) vp where I is the thickness precipitated material in the cloud, and vp is the volume of a particle. The attenuation factor is then fatt = e VP (7.7) If the path through the horizontal cloud is at some angle 0, then the attenuation factor becomes

158 -( + oo) 1-ce fartt = e (7.8) We will examine this attenuation coefficient for the different clouds and dust in the Martian polar atmosphere. 7.1.2: Water Ice Clouds, CO2 Ice Clouds, and Dust Viking observations indicate that the Martian polar atmosphere contains water ice clouds that accompany the retreat of the seasonal C02 cap in the springtime. Christensen and Zurek (1984) estimate that the ice clouds amount to 1-2 gm of precipitated water and particles sizes of about 0.5 gm. We will use overestimates of 10 glm of precipitated water and 10 glm particles. We consider absorption over the region of the hydrate attenuation feature. The high frequency limit is about 24 cm1l, where the complex dielectric constant is 3.17 + iO.0 16. Because this is the highest frequency under consideration, the highest scattering will be exhibited here. The imaginary part of the index of refraction is maximum here, and the cloud should exhibit the most absorption here as well. Plugging the overestimates in the previous paragraph into (7.8), we get an transmission coefficient of 0.997. If a 450 path through the cloud is considered, the transmission coefficient is still greater than 0.99. The multiple overestimates used have led to small effect by the springtime water ice cloud. We therefore will neglect its effect in our model. Water ice clouds are also present in the polar hood. For convenience, we model them exactly as we modeled the springtime clouds, though the amount of water ice in the polar hood is believed to be less. Since the springtime clouds showed a small effect, the polar hood water ice clouds also show a small effect, and we neglect it.

159 The polar hood is also believed to contain low level C02 ice clouds with precipitated CO2 of about 10 gm and particle radii up to 25 gm. We will use 20 gm of precipitated CO2 and a fixed 25 grm radius. At 24 cm-l the real part of the dielectric constant is about 2.0, but the imaginary part has proven to low for investigators to measure. We assume a generous value of 0.004. Using these values, we get a transmission coefficient of 0.998 at normal incidence and 0.997 at 45 0. These number are so close to unity that we neglect the effect of the C02 clouds. Even in the worst dust storms, the total dust in a column of atmosphere will only precipitate to a few microns. Models and data analysis of Toon (1977) and Pollack(1979) show an average particle radius of about 2 glm. We shall assume 10 glm of precipitable particles and 10 glm radii. Since the composition of the dust is not known, we will assume a very generous value of 0.1 for eri and a small value of 1.0 for Err. (Smaller values of Err result in more absorption, and absorption dominates over scattering in this regime. We are therefore overestimating extinction by picking a small value of Err.) This leads to a normal transmission coefficient of 0.992 and a 45 transmission coefficient of 0.99. We will therefore neglect the effect of dust particles on atmospheric propagation and emission. 7.2: Gases Table (7.1) lists the most abundant gases identified in the Martian atmosphere. The concentrations of each species can depend on temperature and season. In this section we consider the effects of C02, CO, H20, 02, and 03 on the far-infrared transmission and emission spectrum of the atmosphere. Kr, Xe, and Ar have no mechanisms for infrared absorption and hence do not play a role. We use a simple model of the pressure and temperature profiles to model the thermodynamical properties of the atmosphere.

160 7.2.1: Molecules in the Model When C02 is usually made of one carbon-12 atom and two oxygen-16 atoms. In this case it is a linear molecule symmetric about the carbon atom and thus has no permanent dipole moment. Electromagnetic waves can thus not apply a direct torque to the equilibrium molecule, and C02 exhibits no rotational lines. Carbon is also naturally abundant as carbon-13 (1.1% on earth), and oxygen is naturally abundant as oxygen-17 and oxygen-18 (0.038% and 0.2% on earth). If we replace the carbon-12 atom in the common C02 molecule with a carbon-13 atom, the molecule is still symmetric and has no dipole moment. However, if we replace one of the oxygen-16 atoms with an oxygen-17 or oxygen-18 atom, the symmetry is destroyed. The heavier oxygen atom will change the C=O bond length slightly, resulting in a small dipole moment. The C02 molecule should then exhibit rotational lines. Apparently the dipole moment is quite small. We could find no data for the low pressure rotational lines in the literature. For example, the Hitran database contains spectral lines for the different C02 isotopic mixtures in the regions of their asymmetric stretching modes and bending modes, but it does not show any purely rotational lines. Rotational lines have been reported at high pressure (-100 atm), but these are collision induced lines and would not be present in the much lower pressures in the Martian atmosphere. CO has a permanent dipole moment and exhibits rotational lines in the farinfrared. It is assumed to have a constant mixing ratio with altitude, with the value given in Table (7.1). We shall see that it is the strongest far-infrared absorber. The mixing ratio for H20 in the polar region is significantly lower than that given in Table (7.1) because the temperature is lower at the poles. We assume that its partial pressure is equal to its saturation vapor pressure. The H20 saturation vapor pressure is calculated using the Clausius-Clapeyron equation.

161 molecule mixing ratio CO 7 x 10 H20 saturation 02 1.3 x 10-3 03 2x 10-' Table (7.1): Model Mixing Ratios of Absorbing Consituents 02 has no permanent electric dipole moment, but it does have a weak magnetic dipole moment due to the addition of spins of unpaired outer shell electrons. The absorption due to magnetic moments is typically much lower than that of electric moments, and we shall see that the 02 spectrum is not visible in the Martian atmosphere. We assume that it has a constant mixing ratio, as given by Table (7.1). The abundance of 03 varies greatly on Mars with location and season. We highly overestimate its abundance by assuming a constant mixing ratio of 10 ppm. We shall see that under this assumption, it is barely visible in the far-infrared Martian atmospheric spectrum. The line positions and strengths for all molecules considered are taken from the Hitran database. We next consider the lineshapes. 7.2.2: Lineshape In the lower atmosphere, the pressure is high enough so that the dominant mechanism for absorption line broadening is collisions. As one goes to higher altitudes, the pressure drops off, and Doppler shifts due to the distributed velocities of the molecules

162 becomes the dominant mechanism. Doppler broadening results in a Gaussian line shape with a half width at half maximum Yd given by Yd =3.581 x 10-7 v (7.9) Here T is the temperature in Kelvin, m is the molecular weight in grams, and v is the frequency. Collisional broadening is more complicated, and many investigators have proposed lineshapes based on both theory and data (Townes and Schawlow, 1955). The shape depends on the types of molecules involved in the collision. Following Johnson (1993), we will assume a Lorentzian lineshape, with a halfwidth given by T -x yc = aOP(3 (7.10) Unfortunately, very little data exists showing the broadening of lines by collision with CO2 molecules, and we are forced to make a reasonable guess. As in Johnson's model, we will assume a temperature parameter x = 0.75 and a broadening parameter o = 3 MHz/mbar. In regions where Doppler and collisional broadening are comparable, we use the convolution of the Gaussian and Lorentzian lineshapes. The new lineshape is called a Voight lineshape, and a fast method of computing it is given by Drayson (1976). The absorption coefficient at a particular frequency is calculated by summing the contributions from all the lines within 5 cml1. In order to determine the absorption linewidths, we need a model for the pressure and temperature. We also need to know the number density of each species to determine the strength of the absorption each will exhibit. All of this is discussed in the next section.

163 7.2.3: Temperature, Pressure, and Number Density Profiles Lindal et al. (1979) have reported temperature profile data for the Martian atmosphere at latitudes between 740S and 73~N. Their data is based on S-band (2.3 GHz) and X-band (8.4 GHz) Viking radio occultation measurements. The measurements at 73'N in the springtime show surface temperatures of about 150 K increasing with altitude up to some point, and then decreasing up to the maximum altitude measured. We will use a temperature model consistent with these measurements for our simulations of atmospheric propagation and emission. The adopted temperature profile is shown in Figure (7. 1). The temperature is 150 K at the surface and rises linearly with altitude up to 160 K at 10 kmn. It then decreases linearly to 130 K at 40 km and is fixed at 130 K for higher altitudes. The same profile was used by Johnson (1993) in his atmospheric modeling. To get the pressure profile, we use the ideal gas law and the hydrostatic law. The ideal gas law is given by PV = nRT (7.11) where P is the pressure, V is the volume, n is the number of moles of molecules, R is the universal gas constant, and T is the temperature. Equation (7.9) can be rewritten as p pRT (7.12) m where m is the average molecular weight of the gas and p is the mass density. The hydrostatic law is given by dP d-' Pg (7.13)

40 30 20 10 130 I 150 160 130 140 150 160 temperature (K) Figure (7.1): Polar Temperature Profile Used in the Atmosphere Propagation and Emission Model

165 where z is the altitude and g is the acceleration due to gravity. Solving for p we get 1dP p = gdz (7.14) g dz Substituting (7.12) into (7.9), we get a differential equation for the pressure: dP mg = -z p (7.15) dz RT(z) Because of the temperature dependence on altitude, equation (7.13) must in general be solved numerically. However, for our simple piecewise linear model of the temperature profile, (7.13) can be solved exactly, given a surface pressure of 6.1 mbar. Once the pressure is known, the number density p, can be determined from the ideal gas equation (7.9) by noting that Pn VNA (7.16) where NA is Avogadro's number. Figures (7.2) and (7.3) give the resulting pressure and density profiles, respectively. We have now characterized the absorption properties of the important gas species in the Martian polar atmosphere. These parameters can now be incorporated into a radiative transfer model. 7.3: Radiative Transfer Model Since we are neglecting scattering in the atmosphere, the radiative transfer equation simplifies greatly. Using a similar notation to that used in the previous chapter, we have

166 50 40 30'- 20 10 0 0 2 4 6 8 10 pressure [mb] Figure (7.2): Pressure Profile Used in the Atmosphere Propagation and Emission Model

167 50 40: 30..m 20 10 0 0 1 2 3 molecule number density [cm'3 x 1017] Figure (7.3): Number Density Profile Used in the Atmosphere Propagation and Emission Model

168 4-l(Z lo) a= _-KaI(bsZ p) +KabsB (T(z) (7.17) dz' p p It is convenient to divide the intensity into upgoing and downgoing parts IU and Id, respectively. Their solutions are given by z Iu(z, g) = Iu(O, g)e + - Jdz'Kab(Z')B (T(z') ) e (7.18) 0 Id(ZP.) = 1 = Jdz'Kabs(Z')B (T(z')) e- (Z Z)/ (7.19) z where z (z, z') = dz" KabS(Z") (7.20) The transmission coefficient e- ("0 o) of the atmosphere at zenith was first calculated to determine whether any atmospheric windows might exist. The transmission values are plotted in Figure (7.4) for the range 10-30 cml1. We see that there are ample windows over the range for performing remote sensing observations. A spaceborne or earth-based probe would observe the upgoing intensity at infinity, which can be calculated via the following method: (1): Solve for the downgoing intensity to get the intensity Id (0, p.) incident on the polar surface. (2): Use Id (0, p) in the semi-infinite layer model in the previous section to get the reflected and emitted intensity Iu (0, p). (3): Use I. (0, p) in equation (7.18) to solve for the upgoing flux at infinity.

169 1.0'' - -"' 0.8 *. 0.6 0.4 0.2 0.0 0.0..!, I. I. I., I. I. I. I, 10 12 14 16 18 20 22 24 26 28 30 wavenumbers [cml'] Figure (7.4): Far-infrared Transmission at Zenith for the Martian Northern Polar Region

170 The results of these calculations are shown in Figure (7.5) for the water ice and hydrate cases for average particle diameters of 200 gLm. We see that the hydrate signature is hardly obscured by atmospheric absorption. As long as a measurement is made away from the strong CO lines, the hydrate and ice emission signatures are practically identical to those shown in Figure (6.4) in the last chapter. The atmosphere does not present a problem, and we conclude that the hydrate should be observable in a remote sensing experiment.

CHAPTER 8 FEASIBILITY OF DETECTING MARTIAN CO2 HYDRATE In this short chapter we discuss the feasibility of detecting Martian C02 hydrate from platforms orbiting Earth or Mars. The key parameters to examine are resolution and sensitivity. We will see that the requirements for a Mars orbiter are not exceptionally stringent, while building an Earth-based instrument presents considerable challenges. 8.1: Resolution 8.1.1: Earth-based System As seen from the Earth, the Martian disk at opposition subtends an angle ranging from about 14 arcseconds in 1980 to about 25 arcseconds in 1971. The difference is due t( the ellipticity of the Martian and Earth orbits. The northern polar region subtends a significantly smaller angle (- 1 arcsecond). Building a single antenna to resolve 1 arcsecond at far-infrared frequencies is difficult because of the size that is required. If we assume a diffraction limited, circular aperture, the half power beamwidth 01/2 is given by p1/2 ~= D(8.1) If we are operating at 15 cm-1 with P1/2 = 1 arcsecond, then the diameter of the antenna is 13.8 m. Building such a large antenna and getting it above the highly absorbing trace gases in the Earth's atmosphere is a challenging task. 172

173 NASA has proposed the construction of the Large Deployable Reflector (LDR) telescope (Swanson et al., 1983). The reflector is specified to have a diameter of 20 m and to be diffraction limited at 200 cm'l. It will contain heterodyne receivers with superconducting junctions, Schottkey diodes, or bulk detectors to cover the range 14-50 cm"l. If the system is ever launched, it could be useful in detecting CO2 hydrate on Mars. The European Space Agency (ESA) has proposed the Far-infrared and Submillimeter Space Telescope (FIRST) (Lamarre, 1994). Its designed calls for a 3-4 m aperture and an operating range of 10-100 cm-l. Its photoconductor and bolometer detectors will be cooled to liquid helium temperatures. This smaller antenna may not be large enough to resolve the polar regions on Mars. Unfortunately, a resolution of 1 arcsecond may not be sufficient to detect CO2 hydrate. In the northern Martian spring, temperatures rise and the season CO2 cap retreats to expose the permanent cap, which is likely to be mostly CO2 hydrate. As the temperature continues to rise, the hydrate converts to water ice, probably beginning at lower latitudes and proceeding north into the summer. Some lower latitude hydrate regions might convert to ice before higher latitude hydrate regions are uncovered. Resolving small hydrate regions would probably be impractical from earth. 8.1.2: Mars-based System Resolution restrictions are much less strict for a Mars orbiting system. For example, if a platform is orbiting 500 km above the surface and requires a resolution of 10 km (1.15 ) at 15 cm1l, a circular aperture with this half power beamwidth has a minimum diameter of 3.3 cm. Such a telescope would be much cheaper and easier to build than a large Earth-based telescope. Deploying at Mars, however, is much more expensive and much less routine.

171 oj 1.0 0.8 --------------- - r - - 1 0.6 C 0.4 hydrate ice 0.2 5 0.0 I 14 16 18 20 22 24 wavenumbers [cm-1] Figure (7.5): Hydrate and Ice Emissions Combined with Absorption and Emission in the Atmosphere. Average particle diameter is 200 km.

174 There are several missions planned to Mars in the next ten years (Pirard, 1994). In November of 1996, the Mars-96 spacecraft is scheduled to be launched carrying the payload of Mars-94. This mission includes an orbiter and a landing system capable of taking core samples. In the same month, Mars Global Surveyor is scheduled for launch, and it will place in Martian orbit an array of cameras and sensors and a data relay station for further missions. A month later the Mars Pathfinder mission should be launched, which will land several instruments that will characterize the Martian surface layer. More missions are planned that include ballonborne instruments positioned within 1 km of the surface and landrovers designed to cover up to 100 km of terrain. The resolution requirements of such instruments would be less strict than that of a high flying instrument. The possibility of detecting hydrates directly from core samples will be discussed with other future work in the next chapter. 8.2: Sensitivity The power received by a telescope from a source that exactly fills the half power beamwidth is given by P= IAvAaQ 1/2 (8.2) where I is the specific intensity of the source, Av is the bandwidth of the receiver, Aa is the area of the aperture, and Q1/2 is the half power solid angle. For a circular aperture, Q1/2 is given by X 2X2 p1/2 = (4) A (8.3) Equation (8.2) then simplifies to

175 P= () IvA2 (8.4) and is independent of the source and telescope sizes. As long as the source fills the beamwidth, we can use equation (8.4) for both Mars-based and Earth-based systems. Based on the emission and atmospheric propagation models of chapters 6 and 7, to differentiate C02 hydrate from water ice through its spectra gradient, we need a resolution of thermal signals of - 1%. At 150 K and 15 cm-1, the specific intensity of a grey body with 0.8 emissivity is - 10-15 W mi2 sr Hz-l. Assuming a modest bandwidth of 1 cm-1, equation (8.4) gives a received power of 8 x 10-12 W. High sensitivity far-IR measurements on spacecraft can be performed with liquid helium cooled bolometers, which typically have a noise equivalent power of -10-18 W/HzL/2. With a 1 cm"l bandwidth, this corresponds to an noise power - 5 x 10-14 W. The signal to noise ratio for a 1 second integration time is then -150. We see that reasonably short integration times are necessary to achieve the specified 1 K sensitivity. Less sensitive detectors can be used with longer integration times.

CHAPTER 9 SUMMARY, CRITIQUE, AND PROPOSALS OF FUTURE WORK This thesis presents an experimental and theoretical study of the far-infrared dielectric properties of C02 clathrate hydrate and its emission signature in the Martian northern polar regions. The goal of the study was to determine the feasibility of detecting the presence of the hydrate in a remote sensing experiment. The models indicate that it may be possible to differentiate a hydrate layer from a water ice layer with an instrument of moderate size and sensitivity orbiting Mars. Detection from an earth-based satellite is limited by the resolving power of large far-infrared telescopes and may not be feasible. Detection of C02 hydrate may also be possible from core samples taken by Martian landers. Davidson (1973, p. 183) reports a -30 C structure I clathrate hydrate DC permittivity of 61 and a similar permittivity of hexagonal ice of 107. This large dielectric difference could be exploited by a low frequency device to differentiate C02 hydrate form water ice. If the experimental methods developed in this thesis are to be used in far-infrared studies of other clathrate hydrates, a few improvements might be considered before proceeding, particularly in the sample cell design. The measurements required the construction of a high pressure, low temperature sample cell with far-infrared access. The final cell design was satisfactory, but suffered from a few flaws. The cell had to be carefully balanced at all times to avoid tipping the sample cup and losing water. The cell also had no optical access ports for viewing the quality of the hydrate as it formed. The quality of the sample had to be inferred from spectroscopic measurements and from earlier mock runs made with clear cell windows. Perhaps the biggest problem with the cell is that 176

177 it does not allow for a direct measurement of the sample thickness. Ideally, the sample would be grown from liquid water filling the space between two parallel window with know separation. Unfortunately, this would probably leave too little of the water in direct contact with CO2 gas. If no better cell design can be conceived, then the current design might be used with minor modifications to measure the far-infrared dielectric properties of other nonpolar clathrate hydrates requiring high pressures to form from liquid water. Methane hydrate has a 00C dissociation pressure of 26 atm, which is about twice the 00C dissociation pressure of C02 hydrate. If it forms in the gas giants, it would be at temperatures -90 K, considerably lower than the 150 K temperature of the Martian northern polar region. The cell would therefore have to be able to withstand higher pressures at 0C and remain sealed at lower temperatures in order to simulate the conditions of the outer solar system. The models in this study should also be scrutinized before they are applied to other remote sensing problems. In particular, the QCA dense medium propagation model should be better tested to evaluate its performance. We used it in conjunction with the PercusYevick pair distribution function, which was originally developed for molecules in a fluid. A more realistic correlation function might produce considerably different results. Since we do not know the packing properties of hydrate particles on Mars, we used the PercusYevick function as a first guess. The code developed for this study can easily be modified to handle an arbitrary, spherically symmetric pair distribution function. (Non-spherically symmetric functions significantly complicate the formulation.) Coupling QCA with the distorted Born approximation led to unrealistic predictions of the emissivity, while our simple radiative transfer model gave much more reasonable results. We favored the simpler model, incorporating into it the propagation constants predicted by QCA. We therefore used "physically based" QCA parameters in the heuristically based radiative transfer theory. A more rigorously derived theory of

178 emission might have been considered. The dense medium radiative transfer theory was developed to put dense medium emission on a more physical foundation. Unfortunately, its development in the resonant region becomes intractable. Atmospheric propagation is much better understood than dense media propagation, and our atmospheric models are probably more accurate. We used overestimated the concentrations of several gases, yet the effect on far-infrared atmospheric transmission was minor. The Martian polar atmosphere in winter and early spring apparently has large far-infrared windows that can be used for surface remote sensing. C02 hydrate as well as C02 ice, which has lattice mode absorptions at about 68 cml and 114 cm'1, may be detectable with far-infrared instruments. The far-infrared is a largely untapped spectral region in astronomy and planetology. With the launching of orbiting platforms holding much higher resolution farinfrared telescopes, the need for measuring and modeling the far-infrared properties of materials becomes increasingly necessary.

APPENDIX SIMPLIFIED THEORY OF A FOURIER TRANSFORM SPECTROMETER Consider the optical configuration in figure (A. 1). The source emits a beam of electromagnetic radiation toward an ideal 50/50 beamsplitter. Half of the beam power is reflected up to a movable retromirror, and the other half is transmitted through the beamsplitter to a fixed retromirror. The two beams are reflected back to the beamsplitter and recombine. One beam is a delayed version of the other, and the delay depends on how much further the movable mirror is than the fixed mirror from the beamsplitter. The difference in the two differences is represented in the figure by z. The extra distance that one beam has to traverse in going to and from the movable mirror is 2z. Hence the time delay between the two beams is 2z/c, where c is the speed of the radiation. After recombining, the beams are absorbed by a square law detector, the output of which is high pass filtered and Fourier transformed to give the spectrum of the source. We will eiamine why the output is the source spectrum by considering very simple sources. Let the source be a simple sine wave, and consider the input to the square law detector for different values of z. Figure (A.2a) shows the case for z=0. The two beams recombine in phase, and the output of the high pass filter is a maximum. Figure (A.2b) shows the case for z=X/4, where x is the wavelength of the source. The beam going to the movable mirror must travel half a wavelength further than the other beam. The two beams thus add out of phase, and the output of the high pass filter is a maximum. When z=./2, as 179

180 shown in figure (kA2c), the path difference is one wavelength, the beams are back in phase, and the output of the square law detector is a maximum again. The output of the detector is shown in figure (A.3). It is a raised sine wave whose frequency depends on the wavelength depends on the source wavelength. If x is small, the movable mirror does not have to travel very far to produce an oscillation in the detector output. If X is large, the mirror must travel a long way to produce an oscillation. As shown in figure (A.3), the DC component of the detector output is removed by the high pass filter, resulting in a pure sine wave. The Fourier transformer converts this to a spike in the frequency domain, with a position that depends on the frequency of the input, which depends on the frequency of the source. The spike position is linearly related to the frequency of the source, and we thus have a device that finds the frequency of a monochromatic source. Using superposition arguments and the orthogonality of sine wave of different frequency, we can show that a two frequency source will result in a two spike output. Similarly, an N frequency source results in an N frequency output Generalizing, if we have a source with a continuous spectrum, we can see that the output will be this spectrum.

181 s(t) source s(t) + s(t-2z/c) 2 detector high pass filter Fourier transformer F/T source spectrum Figure (A.1): Optical and Electronic Configuration of a Fourier Transform Spectrometer

s(t) s(t-2z/c)) s(t)+s(t-2z/c)) (a): z = 0 (b): z X= /4 (c): z = /2 Figure (A.2): Beams with Various Delays and Their Sums. In (a), the two beams have the same delay and up in phase to give a maximum sum amplitude. In (b), one beam has a total delay of half a period, resulting in deconstructive interference and a minimum sum amplitude. In (c), the total delay is a full period, and the beams again add in phase. 9 9~~~~~~~~~~~~~~~~~~ 9 9~~~~~~~~~~~~~~~~~~ 9 9~~~~~~~~~~~~~~~~~~~~o (a): z = 0 (b): z = 3,/4 (c)'z = 3,/2 9 9~~~~~~~~~~~ 9~~~~~~~~~~~~~~~~n a9 n

183 2 detector filter Fourier transformer F T I spectrum Figure (A.3): Electronic Processing of the Interferometer Signal

BIBLIOGRAPHY Atreya and Romani, (1985), "Photochemistry and Clouds of Jupiter, Saturn, and Uranus," in Planetary Meteorology, G. E. Hunt, Editor, Cambridge University Press, London, pp. 17-68. Atreya, S. K., (1986), Atmospheres and Ionospheres of the Outer Planets and Their Satellites, Springer-Verlag, New York Berecz, E. and M. Balla-Achs, (1983), Gas Hydrates, Elsevier, New York Bertie, J. E., and E. Whalley, (1964), "Infrared Spectra of Ices Ih and Ic in the Range 4000 to 350 cm1l," Journal of Chemical Physics, 40, 1637-1645. Bertie, J. E., and E. Whalley, (1967), "Optical Spectra of Orientationally Disordered Crystals. II. Infrared Spectrum of Ice Ih and Ic from 360 cml to 50 cm1l, Journal of Chemical Physics, 46, 1271-1284. Bertie, J. E. and D. A. Othen, (1972),'"The Infrared Spectrum of Ethylene Oxide Clathrate Hydrate between 360 and 20 cm-1, at 100 K," Canadian Journal of Chemistry, 50, 3443-3449. Bertie, J. E. and D. A. Othen, (1973), "'The Infrared Spectrum of Ethylene Oxide Clathrate Hydrate at 100 K between 4000 and 360 cm-1," Canadian Journal of Chemistry, 51, 1159-1168. Bertie, J. E. and S. M. Jacobs, (1977), "Far-infrared Absorption and Rotational Vibrations of the Guest Molecules in Structure I Clathrate Hydrates Between 4.3 and 100 K," Canadian Journal of Chemistry, 55, 1777-1785. Bertie, J. E. and J. P. Devlin, (1983), "Infrared Spectroscopic Proof of the Formation of the Structure I Hydrate of Oxirane from Annealed Low-temperature Condensate," Journal of Chemical Physics, 78, 6340-6341. Blake, D., L. Allamandola, S. Sandford, D. Hudgins, and F. Freund, (1991), "Clathrate Hydrate Formation in Amorphous Cometary Ice in Vacuo," Science, 254, 548-551. Califano, S., (1981), Lattice Dynamics of Molecular Crystals, Springer-Verlag, New York. Chandrasekhar, S., (1960), Radiative Transfer, Dover, New York. 184

185 Christensen, E. and R. Zurek, (1984), "Martian Northern Polar Hazes and Surface Ice: Results from the Viking Survey/Completion Mission," Journal of Geophysical Research, 89, 4587-4596. Consai, K. and G. C. Pimentel, (1987), "Infrared Spectra of the Clathrate Hydrates of Acetylene and of Acetylene/Acetone," Journal of Physical Chemistry, 91, 289293. Davidson, D. W., (1973), "Clathrate Hydrates," in Water: A Comprehensive Treatise, volume 2, F. Franks, Editor, pp. 115-234, Plenum Press, New York. Davidson, D. W., S. K. Garg, S. R. Gough, Y. P. Handa, C. I. Ratcliffe, J. S. Tse, and J. A. Ripmeester (1984), "Some Structural and Thermodynamic Studies of Clathrate Hydrates," Journal of Inclusion Phenomena, 2, 231-238. Delsamme, A. H. and P. Swings, (1952), "Hydrtaes de Gaz dans les Noyaux Cometaires et les Grains Interstellaires," Annales d'Astrophysique, 15, lff. Delsemme, A. H. and A. Wegner, (1970), "Experimental Study of Snows in a Cometary Environment," Planetary and Space Science, 18, 709-715. Dobrovolskis, A. and A. P. Ingersoll, (1975), "Carbon Dioxide-Water Clathrate as a Reservoir of C02 on Mars," Icarus, 26, 353-357. Drayson, R., (1976), "Rapid Computation of the Voight Profile," Journal of Quantitative Spectroscopy and Radiative Transfer, 16, 611-614. Fink, U. and G. T. Sill, (1982),'"The Infrared Spectral Properties of Frozen Volatiles," in Comets, L.. Wilkening, Editor, University of Arizona Press, Tucson, Arizona, pp. 164-20_. Fleyfel, F. and J. P. Devlin, (1988), "FT-IR Spectra of 90 K Films of Simple, Mixed, and Double Clathrate Hydrtaes of Trimethylene Oxide, Methyl Chloride, Carbon Dioxide, Tetrahydrofuran, and Ethylene Oxide Containing Decoupled D20," Journal of Physical Chemistry, 92, 631-635. Fleyfel, E. and J. P. Devlin, (1991), "Carbon Dioxide Clathrate Hydrate Epitaxial Growth: Spectroscopic Evidence for Formation of the Simple Type-II C02 Hydrate," Journal of Physical Chemistry, 95, 3811-3815. Gerbaux, M. M. X., C. Barthel, and A. Hadni, (1975), "Absorption de Clathrates de Glace dans L-infrarouge Lointain," Spectrochimica Acta, 31A, 1901-1903. Gow, A. J., H. T. Ueda, and D. E. Garfield, (1968), Science, 164, 1011-1013.

186 Handa, Y. P. and J. G. Cook, (1987), Thermal Conductivity of Xenon Hydrate," Journal of Physical Chemistry, 91, 6327-6328. Hardin, kA H. and K. B. Harvey, (1971), "Comments on the Infrared Spectra of SO2 Clathrate-hydrate Formed by Vapor Condensation," Canadian Journal of Chemistry, 49, 4114ff. Harvey, K. B., F. R. McCourt and H. F. Shurvell, (1964), "Infrared Absorption of the S02 Clathrate-HydrateMotion of the S02 Molecule," Canadian Journal of Chemistry, 42, 960-963. Ishimaru, A., (1978), Wave Propagation and Scattering in Random Media, Volume 1, Academic Press, New York. James, P. B., H. H. Kieffer, and D. A. Paige, (1992), The Seasonal Cycle of Carbon Dioxide on Mars," in Mars, Kieffer et al., Editors, pp. 934-968, University of Arizona Press, Tucson. Johnson, B., (1993), Laboratory and Theoretical Study of the Far Infrared Spectra of Martian Ices, Ph.D. Thesis, University of Michigan, Ann Arbor Kieffer, H. H., (1990), "H20 Grain Size and the Amount of Dust in Mars' Residual North Polar Cap," Journal of Geophysical Research, 95, 1481-1493. Kiefte, H., M. J. Clouter, and R. E. Gagnon, (1985), "Determination of Acoustic Velocities of Clathrate Hydrtaes by Brillouin Spectroscopy," Journal of Physical Chemistry, 89, 3103-3108. Lamarre, J., (1994), "FIRST: Far-infrared and Submillimeter Space Telescope, a major Scientific Project of the European Space Agency," Optical Engineering, 33, 785790. Lee, S., M. Wolff, J. Bell, P. James, T. Clancy, and L. Martin, (1995), "Mars Mosiac," EOS Transactions, 76, 129. Lewis, J. S., (1968),'"The Clouds of Jupiter and the NH3-H20 and NH3-H2S Systems," Icarus, 10, 365-378. Lindal, G. F., H. B. Hotz, D. N. Sweetnam, Z. Shippony, J. P. Brenkle, G. V. Hartsell, and R. T. Spear, (1979), "Viking Radio Occulatation Measurements of the Atmosphere and Topography of Mars: Data Aquired During 1 Martian Year of Tracking," Journal of Geophysical Research, 84, 8443-8456. Lunine, J. and D. J. Stevenson, (1985),'Thermodynamics of Clathrate Hydrate at Low and High Pressures with Application to the Outer Solar System," Astrophysical Journal Supplement Series, 58, 493-531.

187 Miller, S. L., (1961),'"The Occurrence of the Gas Hydrates in the Solar System," Proceedings of the National Academy of Science, 47, 1798-1808. Miller, S. L., (1969), "Clathrate Hydrates of Air in Antarctic Ice," Science, 165,489-490. Miller, S. L. and W. D. Smythe, (1970),'Carbon Dioxide Clathrate Hydrate in the Martian Polar Cap," Science, 1970, 531-533. Miller, S. L., (1985), "Clathrate Hydrates in the Solar System," in Ices in the Solar System, J. Klinger et al., Editors, pp. 59-78, D. Reidel, Boston. Mishima, O., D. D. Klug, and E. Whalley, (1983),'The Far-infrared Spectrum of Ice Ih in the Range 8-25 cml1: Sound Waves and Difference Bands, with Applications to Saturn's Rings, Journal of Chemical Physics, 78, 6399-6404. Owen, T., (1982), The Composition and Origin of Titan's Atmosphere," Planetary and Space Science, 30, 833-838. Pirard, T., (1994), "Mars Armada: 1996-2003," Spaceflight, 36, 315. Pollack, J., D. Colburn, D. Flasar, R. Kahn, C. Carlston, and D. Pidek, (1979), "Properties and Effects of Dust Particles Suspended in the Martian Atmosphere," Journal of Geophysical Research, 84, 2929-2945. Pouchert, C. J., Editor, (1989), The Aldrich Library of FT-IR Spectra, edition 1, volume 3: Vapor Phase, Aldrich Chemical Company, Milwaukee. Richardson, H. H., P. J. Woolridge, and J. P. Devlin, (1985), "FT-IR Spectra of Vacuum Deposited Clathrate Hydrates of Oxirane, H2S, THF, and Ethane," Journal of Chemical Physics, 83, 4387-4394. Ross, R. G. and P. Anderson, (1982), "Clathrate and Other Phases in the Tetrahydrofuranwater System: Thermal Conductivity and Heat Capacity Under Pressure," Canadian Journal of Chemistry, 60, 881-892. Rothman, L. S., R. R., A. Goldman, L. R. Brown, R. A. Toth, H. M. Pickett, R. L. Poynter, J.-M. Flaud, C. Camy-Peyret, A. Barbe, N. Husson, C. P. Linsland, and M. A. H. Smith, (1987),'"The Hitran Database: 1986 Edition," Applied Optics, 26, 40584097. Smythe, W. D., (1975), "Spectra of Hydrate Frosts: Their Application to the Outer Solar System", Icarus, 24, 421-427. Swanson, P. N., S. Gulkis, and T. H. Kulper, (1983), "Large Deployable Reflector (LDR): A Concept for an Orbiting Submillimeter-infrared Telescope for the 1990s," Optical Engineering, 22, 725-731.

188 Toon, O., J. Pollack, and C. Sagan, (1977), Physical Properties of the Particles Composing the Martian Dust Storm of 1971-1972, Icarus, 8, 227-258. Townes, C. and A. Schawlow, (1955), Microwave Spectroscopy, McGraw-Hill, New York. Tsang, L., J. A. Kong, and R T. Shin, (1985), Theory of Microwave Remote Sensing, John Wiley and Sons, New York. Tse, J. S., M. L. Klein, and I. R McDonald, (1983), "Molecular Dynamics Studies of Ice Ic and the Structure I Clathrate Hydrate of Methane," Journal of Physical Chemistry, 87, 4198-4203. Tse, J. S. and M. L. Klein, (1984), "Computer Simulation Studies of the Structure I Clathrate Hydrates of Methane, Teraflouromethane, Cyclopropane, and Ethylene Oxide," Journal of Chemical Physics, 81, 6146-6153. Tsujimoto, S. A. Konishi, and T. Kunitomo, (1982), "Optical Constants and Thermal Radiative Properties of H20 Cryodeposits," Cryogenics, 22, 603-606. Tyler, G., D. Sweetnam, J. Anderson, J. Campbell, V. Eshleman, D. Hinson, G. Levy, G. Lindal, E. Marouf, and R. Simpson, (1986), "Radio Science Observations of the Uranian System with Voyager 2: Properties of the Atmosphere, Rings and Satellites," Science, 233, 79-84. Warren, S. G., (1984), "Optical Constants of Ice form the Ultraviolet to the Microwave," Applied Optics, 23, 1206-1225. Warren, S. G., (1986), "Optical Constants of Carbon Dioxide Ice," Applied Optics, 25, 2650-2674. Weidenschilling, S. J. and J. S. Lewis, (1973), "Atmospheric and Cloud Structures of the Jovian Planets," Icarus, 20, 465-476. Whalley, E. and H. Labbe, "Optical Spectra of Orientationally Disordered Crystals. III. Infrared Spectra of Sound Waves," Journal of Chemical Physics, 51, 3120-3127. Whipple, F. L., (1950), Astrophysical Journal III, 375ff.

UNIVERSITY OF MICHIGAN 390111111111111115 03466 3990 3 9015 03466 3990