UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Molecular mechanisms of methane hydrate dissociation and inhibition Bagherzadeh Hosseini, Seyyed Alireza 2015

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2015_may_bagherzadehhosseini_seyyedalireza.pdf [ 21.41MB ]
Metadata
JSON: 24-1.0166109.json
JSON-LD: 24-1.0166109-ld.json
RDF/XML (Pretty): 24-1.0166109-rdf.xml
RDF/JSON: 24-1.0166109-rdf.json
Turtle: 24-1.0166109-turtle.txt
N-Triples: 24-1.0166109-rdf-ntriples.txt
Original Record: 24-1.0166109-source.json
Full Text
24-1.0166109-fulltext.txt
Citation
24-1.0166109.ris

Full Text

MOLECULAR MECHANISMS OF METHANE HYDRATE DISSOCIATION AND INHIBITION  by Seyyed Alireza Bagherzadeh Hosseini  M.A.Sc., The University of British Columbia, Canada, 2010 B. Sc., Sharif University of Technology, Iran, 2008  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  DOCTOR OF PHILOSOPHY in  THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Chemical and Biological Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  March 2015  © Seyyed Alireza Bagherzadeh Hosseini, 2015   ii Abstract Gas hydrates are crystalline compounds with cage-like structures formed by hydrogen-bonded water molecules hosting guest molecules such as light hydrocarbons and CO2. They are known to: • represent a potential reserve of natural gas embedded in seabed and permafrost sediments • pose a flow assurance challenge to the oil and gas industry Molecular dynamics simulations are employed to study the processes of gas hydrate decomposition and inhibition. To mimic the porous environment of the real gas hydrate reservoirs, hydroxylated silica surfaces are included in the simulations and placed in contact with hydrate and water.  Water molecules wet the silica surfaces and form a meniscus, confirming the hydrophilic properties of the hydroxylated silica surface. It is found that the silica surface alters the characteristics of the confined water up to ~6 Å away from the surface. The decomposition of methane hydrate in the presence of silica surfaces, 34 to 40 Å apart, follows a concerted behavior where layers of hydrate cages at the curved dissociation front collapse almost simultaneously. The rate of hydrate dissociation in contact with a silica surface is faster compared to that of a hydrate phase just in contact with bulk water. Additionally, the decomposition leads to the formation of methane-rich regions (nano-bubbles) in the liquid water phase. In more realistic simulations, gas reservoirs are added to the simulations to determine whether the formation of nano-bubbles is a general feature of the hydrate decomposition process. It is found that the nano-bubbles can form under simulation conditions where the dissociation rate is faster than the diffusion rate, thus generating dissolved methane mole fractions of greater than 0.044 that would lead to bubble nucleation. Finally, the binding mechanism of the alpha-helical 37 amino acid residue winter flounder antifreeze protein, which is a candidate as a kinetic hydrate inhibitor to methane hydrate, is determined to be the result of cooperative anchoring of the pendant methyl groups of the threonine and two alanine residues, four and seven places further down in the protein sequence, to the empty half cages at the hydrate surface.   iii Preface The research undertaken in this dissertation including the identification of the research questions, design and performing of the simulations/experiments, the analysis of the research data as well as the writing of the first draft of the dissertation and the six manuscripts arising from the work presented here are solely done by the candidate. The research supervisor and other authors of the papers contributed by providing comments on the first draft and thoroughly examining the content. A list of published/submitted manuscripts for the work presented in this dissertation is given below: • A version of Chapter 3 has been published. S. A. Bagherzadeh, P. Englezos, S. Alavi, and J. A. Ripmeester, Influence of Hydrated Silica Surfaces on Interfacial Water in the Presence of Clathrate Hydrate Forming Gases, Journal of Physical Chemistry C 2012, 116(47): 24907-24915 • A version of Chapter 4 has been published. S. A. Bagherzadeh, P. Englezos, S. Alavi, and J. A. Ripmeester, Molecular Modeling of the Dissociation of Methane Hydrate in Contact with a Silica Surface, Journal of Physical Chemistry B 2012, 116(10): 3188-97 S. A. Bagherzadeh, P. Englezos, S. Alavi, and J. A. Ripmeester, Molecular Simulation of Non-equilibrium Hydrate Decomposition Process, The Journal of Chemical Thermodynamics. 2012, 44(1): 13-19 • A version of Chapter 5 has been submitted/published. S. A. Bagherzadeh, S. Alavi, J. A. Ripmeester, and P. Englezos, Formation of Methane Nano-bubbles during Hydrate Decomposition and Their Effect on Hydrate Growth, submitted (under review). S. A. Bagherzadeh, S. Alavi, J. A. Ripmeester, and P. Englezos, Evolution of Methane During Gas Hydrate Dissociation, Fluid Phase Equilibria. 2014, 358: 114-120 • A version of Chapter 6 has been published. S. A. Bagherzadeh, S. Alavi, J. A. Ripmeester, and P. Englezos, Why Ice-Binding Type I Antifreeze Protein Acts as a Gas Hydrate Crystal Inhibitor, 2015, DOI: 10.1039/C4CP05003G. V. K. Walker, H. Zeng, H. Ohno, N. Daraboina, H. Sharifi, S. A. Bagherzadeh, S. Alavi and, P. Englezos, Antifreeze Proteins as Gas Hydrate Inhibitors, accepted for publication in the Canadian journal of Chemistry, 2015, dx.doi.org/10.1139/cjc-2014-0538    iv The following list summarizes the contributions arising from this work to various conferences: • S. A. Bagherzadeh, S. Alavi, J. A. Ripmeester, and P. Englezos, Mechanism of Binding of Type I Antifreeze Protein as an Inhibitor to the Hydrate Surface, in Proceedings of the Industrial Chemistry Conference, Edmonton AB, Canada, November 12-14, 2014 • S. A. Bagherzadeh, P. Englezos, S. Alavi, and J. A. Ripmeester, The Evolution of Methane Gas During Hydrate Dissociation, in Proceedings of the 13th International Conference on Properties and Phase Equilibria for Products and Process design (PPEPPD2013), Iguazu, Argentina, May 26-29, 2013 • S. A. Bagherzadeh, P. Englezos, S. Alavi, and J. A. Ripmeester, The State and Fate of Water in a Slit-like Silica Pore in the Presence of Clathrate Hydrate Forming Gases, in Proceedings of the 62nd Canadian Chemical Engineering Conference, Vancouver, Canada, October 14-17, 2012 • S. Alavi, S. A. Bagherzadeh, R. Kumar, P. Englezos, and J. A. Ripmeester, Molecular Dynamics Simulations of CH4 Clathrate Hydrate Dissociation Adjacent to Hydrated Silica Surfaces, in Proceedings of the 7th International Conference on Gas Hydrates, Edinburgh, Scotland, July 17-21, 2011 • J. Ripmeester, S. Hosseini, P. Englezos, and S. Alavi, Fundamentals of Methane Hydrate Decomposition, in Proceedings of Canadian Unconventional Resources and International Petroleum Conference; Calgary, Alberta, Oct. 19-21, 2010.     v Table of Contents Abstract ..................................................................................................................................... ii Preface ..................................................................................................................................... iii Table of Contents ...................................................................................................................... v List of Tables ............................................................................................................................ ix List of Figures ............................................................................................................................ x List of Symbols ........................................................................................................................ xx List of Abbreviations .............................................................................................................. xxi Acknowledgements ............................................................................................................... xxii Dedication ........................................................................................................................... xxiii Chapter 1: Introduction ............................................................................................................. 1 1.1 What Are Gas Hydrates? ................................................................................................................ 1 1.2 Significance of Gas Hydrates ......................................................................................................... 3 1.2.1 A Potential Source of Natural Gas ........................................................................................... 3 1.2.2 A Flow Assurance Problem ....................................................................................................... 3 1.2.3 Other Applications of Gas Hydrates........................................................................................ 4 1.3 Hydrate Research ............................................................................................................................. 4 1.4 Basics of the Molecular Dynamics Simulation ............................................................................ 6 1.4.1 Calculation of Forces .................................................................................................................. 8 1.4.2 Calculation of Physical Properties ............................................................................................ 9 1.4.2.1 Temperature ........................................................................................................................ 9 1.4.2.2 Pressure.............................................................................................................................. 10 vi 1.4.3 Treatment of the Long-Range Electrostatic Interactions ................................................... 11 1.4.4 Ensembles Used ........................................................................................................................ 12 1.4.5 The General Algorithm for MD ............................................................................................. 12 Chapter 2: Literature Review ................................................................................................... 14 2.1 Molecular Simulations of Hydrate Formation and Growth .................................................... 14 2.2 Molecular Simulations of Hydrate Dissociation ....................................................................... 15 2.3 Methane Water Solubility Issue ................................................................................................... 18 2.4 Molecular Simulations of Hydrate Inhibition ............................................................................ 20 2.5 Thesis Objectives ........................................................................................................................... 23 2.6 Thesis Organization ...................................................................................................................... 25 Chapter 3: State of Water in the Presence of Silica ................................................................. 26 3.1 Introduction .................................................................................................................................... 26 3.2 Computational Methods ............................................................................................................... 27 3.3 Results and Discussion ................................................................................................................. 29 3.3.1 Water and Gas Distribution .................................................................................................... 30 3.3.2 Contact Angle ............................................................................................................................ 34 3.3.3 Water Orientation ..................................................................................................................... 37 3.3.4 Residence Correlation Function.............................................................................................. 38 3.3.5 Hydrogen Bond Analysis ......................................................................................................... 39 3.3.6 Effect of the Separation between the Silica Surfaces ........................................................... 40 3.3.7 Discussion .................................................................................................................................. 43 3.4 Conclusions .................................................................................................................................... 43 Chapter 4: Methane Hydrate Dissociation in the Presence of Silica ...................................... 45 vii 4.1 Introduction .................................................................................................................................... 45 4.2 Computational Methods ............................................................................................................... 46 4.3 Results and Discussion ................................................................................................................. 49 4.4 Conclusions .................................................................................................................................... 59 Chapter 5: Nano-bubble Formation during Methane Hydrate Dissociation ......................... 61 5.1 Introduction .................................................................................................................................... 61 5.2 Simulation System and Methods ................................................................................................. 61 5.3 Results and Discussion ................................................................................................................. 63 5.3.1 Dynamics of Hydrate Dissociation ........................................................................................ 64 5.3.2 Nano-bubble Size ...................................................................................................................... 66 5.3.3 Critical Concentration of Dissolved Methane Required for Bubble Formation ............. 69 5.3.4 Effect of Nano-Bubbles on the Growth Rate of the Hydrate Phase ................................ 73 5.4 Conclusions .................................................................................................................................... 75 Chapter 6: Mechanism of Binding of wf-AFP to the Hydrate Surface ................................... 77 6.1 Introduction .................................................................................................................................... 77 6.2 Computational Methods ............................................................................................................... 77 6.3 Results and Discussion ................................................................................................................. 79 6.3.1 AFP/Water Simulation ............................................................................................................ 79 6.3.2 AFP/Water/Hydrate Simulation (case I) ................................................................................ 81 6.3.3 Protein Deformation during the case I Simulation ................................................................ 82 6.3.4 AFP/Water/Hydrate Simulation (case II) .............................................................................. 83 6.3.5 The F3 Order Parameter ........................................................................................................... 84 6.4 Hydrogen Bonding Analysis ........................................................................................................ 90 viii 6.5 Conclusions .................................................................................................................................... 92 Chapter 7: Conclusions and Recommendations ..................................................................... 94 7.1 Contributions to the Hydrate Knowledge ................................................................................. 94 7.2 Recommendations for Future Work ........................................................................................... 95 Bibliography ............................................................................................................................ 97 Appendices ............................................................................................................................ 108 Appendix A ............................................................................................................................................... 108 Appendix B ............................................................................................................................................... 110 Appendix C ............................................................................................................................................... 111 Appendix D ............................................................................................................................................... 112 Appendix E ............................................................................................................................................... 114 Appendix F ................................................................................................................................................ 116  ix List of Tables Table  1-1. Different ensembles used in this thesis ..................................................................................... 12 Table  2-1. Guest concentration (mole fraction) enhancement ................................................................. 19 Table  3-1. Parameters of the force field applied in this chapter ............................................................... 27 Table  3-2. The temperature, number of molecules in the water and gas phases and the calculated contact angles (degrees). The simulations all have a total run time of 2 ns. ........................................... 28 Table  3-3. The density of CH4 and CO2 molecules at the interface and gas phase of the simulation cell for different simulation set-ups. All values are reported with three significant digits. The standard deviation of the mean for gas density is approximately in the order of 10-7 mol.cm-3 for all simulations. A sample value is given for simulation 3 in brackets. ........................................................... 33 Table  5-1. Number of water molecules, nW, dissolved methane, nCH4 (aq), and gaseous methane, nCH4(g), for different simulations performed to determine the CDMM. ............................................... 71 Table  6-1. Number of molecules present in each simulation set .............................................................. 78 Table  6-2. Average numbers of hydrogen bonds between individual THR residues and neighboring water molecules. ............................................................................................................................................... 92 Table  A-1. Unit cell parameters for the original and the new silica unit cell ....................................... 108 Table  F-1. Heavy atoms of the side chain of all residues of wf-AFP around which the F3 is calculated. The wf-AFP molecular structure is given in Figure  2-5. .......................................................................... 116  x List of Figures Figure  1-1. (A) A burning methane-propane hydrate. (B) A schematic of the small cage of sI hydrate accommodating a methane molecule. Red lines represent hydrogen bonding between the water molecules. (C) Structure I hydrate. The large (51262) and small (512) cages are colored as red and green, respectively. Methane molecules are represented by spheres. Source: http://chemwiki.ucdavis.edu ........................................................................................................................... 1 Figure  1-2. The phase boundary of methane hydrate in pure water (blue). Adding methanol (MeOH), as a thermodynamic hydrate inhibitor, pushes the hydrate formation conditions to higher pressures and lower temperatures (red). ......................................................................................................... 2 Figure  1-3. The length and time scales over which hydrate research is being carried out. ..................... 5 Figure  1-4. The general algorithm of molecular dynamics simulations ................................................... 13 Figure  2-1. The common configuration for hydrate decomposition simulations. The hydrate phase, in the middle, is in contact with bulk water. ................................................................................................ 17 Figure  2-2. The configuration for hydrate decomposition simulations in the presence of porous media. ................................................................................................................................................................. 18 Figure  2-3. The configuration for hydrate decomposition simulations including separate gas reservoirs to allow the released guest reach the gas phase. ........................................................................ 19 Figure  2-4. Sample AFPs. (A) Type I, Winter flounder (Pseudopleuronectes americanus): wf-AFP, (B) Type III, Ocean pout (Zoarces americanus): op-AFP, (C) Snow flea (Hypogastrura nivicola): sf-AFP (D) Yellow mealworm beetle (Tenebrio molitor): mb-AFP. Color scheme: H, white; C, cyan; N, blue; O, red; S, yellow. The ice-binding surfaces of these proteins are relatively flat and rich in threonine. .... 20 Figure  2-5. The primary and three-dimensional structures of wf-AFP. Alanine (ALA, A): white bars. Threonine (THR, T) residues are circled. The (i+4) and (i+7)ALA residues with respect to the (i)THR residue are pointed by arrows. Other amino acid residues are shown with ball-and-stick models. Color scheme: carbon, cyan; oxygen, red; hydrogen, white; nitrogen, blue. ............................ 21 xi Figure  2-6. Chemical structures of (A) the THR residue and (B) the ALA residue. The side chains are emphasized by transparent volumes around them. THR has one hydroxyl and one methyl side chain whereas ALA only has one methyl side chain. Color scheme: carbon, cyan; oxygen, red; hydrogen, white; nitrogen: blue. .................................................................................................................... 21 Figure  3-1. (a) Initial set-up and (b) final configuration for simulation (2) at 275 K with 1584 water molecules between the silica slabs at the top and bottom of the simulation cell. The formation of the meniscus illustrates the hydrophilic character of the silica surfaces used in this work. Three distinct regions designated as G (in gas phase), I (at the water-gas interface) and C (in water phase) are used to calculate average gas concentration in the simulation box. ................................................................... 29 Figure  3-2. The final (t = 2 ns) snapshot of the simulations at T = 275 (simulations 3 and 5) and 300 K  (simulations 9 and 12) in the presence of  CH4 and CO2 gaseous phases. CH4 molecules are illustrated as cyan spheres and CO2 as linear cyan-red molecules. The enhanced concentration of the gas near the interface is evident. The numbers in parentheses correspond to simulation numbers indicated in Table  3-2. The asterisk marks a potential nucleation site for hydrate crystals. ................. 30 Figure  3-3. Water distribution along the distance between the two silica surfaces. Two major peaks of enhanced water density are present near the surface. For clarity, only half of the symmetric profile is shown. ............................................................................................................................................................ 31 Figure  3-4. Concentration surfaces for (a) methane and (b) water density projected on the yz-plane at 275 K (simulation #3). The enhanced concentrations of methane at the silica surface-gas and water-gas interfaces can be observed. Water density is also enhanced at the silica surface. ............................ 32 Figure  3-5. Distribution of hydrate forming gases adjacent to a flat water interface along the length of simulation box (no silica surface) for (a) the methane – water system and (b) the carbon dioxide – water system. The behavior of the methane-water system is similar to that of Figure  3-2. ................. 34 Figure  3-6. A typical distribution of water molecules between the silica surfaces obtained by binning of the areas in the yz-plane of the simulation. The higher local concentration of water molecules near the silica surfaces is apparent. The outer boundary of the water meniscus is used to estimate the xii contact angle. Four tangent lines were drawn on four points of contact of water to the silica surfaces and the contact angles are reported in Table  3-2. A schematic of the method for drawing tangent lines on the meniscus to determine the contact angle is shown in the bottom panel. An arc of approximate length of 2.5 Å is used to draw the tangent to the meniscus near the silica surface. ...... 35 Figure  3-7. (a) Schematic of the orientation of water molecules in the layer adjacent to the silica surface. The blue dotted arrow shows the direction of dipole moment vector M and the vector normal to the silica surface, n. (b) the average angle φ between the polar moment vector of water molecules and normal vector to the silica surface. The black lines are the L1 layer values, the red lines are the L2 layer values and the blue lines are the C layer values. ............................................................... 37 Figure  3-8. The residence correlation parameter for the water layers in different simulations. The numbers in parentheses correspond to the simulation numbers in Table  3-2. The correlations for the L1 water layer adjacent to silica surface are shown in black, the correlations for L2 water layer are shown in red, and the correlations for the water layer in the center of the meniscus C are shown in blue. Top row is at 275 K and bottom row at 300 K. The gas phase space is empty for left panel, filled with 306 methane molecules in the middle panel and 306 carbon dioxide molecules in the right panel. .................................................................................................................................................................. 39 Figure  3-9. The hydrogen bond correlation functions for water molecules in the L1, L2, and C layers (each 3 Å thick). (a) hydrogen bonding dynamics in layers near the silica surface (L1 and L2) and the central layer (C) for simulation 4. The hydrogen bonding in L1 is separated into two contributions: hydrogen bonding between water molecules (WW) and hydrogen bonding between water molecules and OH of the silica surface (SW). (b) hydrogen bonding dynamics in layers L1 and L2 for simulations at 275 K (5) and 300 K (11). ..................................................................................................... 40 Figure  3-10. Snapshot from Simulations 3, 3a, and 3b of Table  3-2 after 2 ns showing the meniscus and methane gas phase for systems with silica layer separations of y = 62 Å (top), 50 Å (middle), and 38 Å (bottom). The contact angle remains relatively the same in these three cases. ............................. 42 xiii Figure  4-1. The initial configurations for simulation of methane hydrate dissociation; (a) HW scenario; (b) SHW scenario and (c) SWHW scenario. The 10 slabs perpendicular to the z-direction and 6 and 8 slabs perpendicular to the y-direction used to characterize the system behavior during methane hydrate dissociation are also shown. For clarity, sampling layers are only shown along the z-direction in HW scenario and along y-direction in SWHW scenario. The length of the simulation box along the z- and x-direction are Lz = 145 Å and Lx = 36 Å, respectively. ...................................... 47 Figure  4-2. The configurations at the end of the 2 ns runs with Tinit = 283, 293, and 303 K in the HW scenario. Methane molecules, liquid water and hydrate are represented by cyan spheres, red-white bonds, and, red dotted hydrogen bonds, respectively. .................................................................... 49 Figure  4-3. The F3,Li order parameter for the water and hydrate phases as a function of time for different layers perpendicular to the z-direction (see Figure  4-1) for simulations with Tinit = 283 (top row), 293 (middle row), and 303 K (bottom row). In the HW scenario at 283 K, no large scale dissociation of the methane hydrate phase occurs and the incomplete cages at the interface (slabs 3 and 8), show F3,Li values of 0.04 which are between the hydrate-like and liquid-like states of water. . 50 Figure  4-4. Temperature profiles for the simulations at 283 K (top), 293 K (middle) and 303 K (bottom) in the NVE ensemble for three different scenarios introduced in the text, HW (left), SHW (middle) and SWHW (right). The drop in temperature is an indication of hydrate dissociation. ........ 52 Figure  4-5. Pressure profiles for the decomposition of methane hydrate. The first, the second and the third rows are at 283, 293 and 303 K. The first, second and third columns correspond to the HW, SHW and SWHW scenarios, respectively. ......................................................................................... 53 Figure  4-6. Snapshots of the NVE simulation of the HW scenario with methane hydrate (water molecules colored with red and white) in contact with liquid water (water molecules colored as blue) at 283 K. Some distortion in the interface of the methane hydrate phase and exchange of water molecules between the liquid phase and intact methane hydrate cages is observed. ............................. 54 Figure  4-7. Simulation configurations at the end of the 2 ns NVE runs for decomposition of methane hydrate in the SWHW scenario at initial temperatures of 283, 293, and 303 K. ................... 55 xiv Figure  4-8. Contours of the local F3 order parameter showing the dissociation trend under SWHW (left) and HW scenario (right) at Tinit = 303 K. ............................................................................................ 56 Figure  4-9. Bubble formation during the decomposition of methane hydrate at initial temperature of 303 K after 2 ns simulation times. Top: HW scenario, middle: SHW scenario, bottom: SWHW scenario. Methane molecules are shown by cyan spheres. The water molecules initially in the hydrate and liquid water phases are represented by hydrogen bonds in red and blue, respectively. ................. 57 Figure  4-10. Profile of the number of methane molecules along the z-axis of the simulation box. The simulations are all at Tinit = 303 K. ................................................................................................................ 59 Figure  5-1. Projection of the initial configuration of G-W-H-W-G three-phase simulation system on the yz-plane. A 3×3×6 super cell of hydrate phase (H) is placed between two bodies of water (W) and is sandwiched by gas reservoirs (G). For analysis, liquid water and hydrate phase are sliced into 20 layers, each 6 Å thick along the z-direction such that layers 5 and 16 lie at the water/hydrate interface. Methane molecules are shown as cyan spheres. Oxygen and hydrogen of liquid water are shown in red and white. Hydrate water is shown as hydrogen bonds in red evident the middle section of the simulation box. ........................................................................................................................ 62 Figure  5-2. Potential energy of the system during hydrate decomposition. The increase in potential energy indicates the dissociation of hydrate. As expected, rate of hydrate decomposition is faster at higher temperatures. It takes approximately 1.9, 4.3, 14 and 41 ns at 360, 350, 340 and 330 K for the hydrate block to melt completely. ................................................................................................................. 64 Figure  5-3. The local order parameter, F3, as a function of time for different layers in z-direction during the simulation at 330 K. For layer assignment refer to Figure  5-1. Hydrate decomposition progresses in a concerted manner where rows of cages parallel to the hydrate/water interface collapse simultaneously within a 15 ns time period. This is associated with the increase of F3 value from ~ 0.01 to ~ 0.09. .................................................................................................................................... 65 Figure  5-4. A sample snapshot of the decomposing hydrate-water-gas system at 350 K shown in projections on the xz- and yz-planes. Methane molecules initially in the hydrate and gas phases are xv shown in pink and cyan spheres, respectively. Water molecules are shown by hydrogen bonds in red. Methane molecules from the decomposing hydrate enter into the gas reservoirs and remain in the liquid phase as spherical and cylindrical nano-bubbles. The projections on the xz and yz planes show the cross sections of the nano-bubbles on these planes. As seen, the right-hand nano-bubble has a roughly spherical shape whereas the left-hand bubble is roughly cylindrical. A small amount of methane remains dissolved in water. ............................................................................................................. 66 Figure  5-5. Water density distribution for simulation at 350 K at 14 ns, after the dissociation of the hydrate phase. Blue area represents the gas phase (water density of zero) and green area corresponds to liquid water with average density of ~ 1 g·cm-3. ..................................................................................... 67 Figure  5-6. (A) yz-snapshot of the mass density distribution of methane molecules for the simulation at 350 K at 5 ns where the hydrate phase has decomposed entirely. Nano-bubble 1 (nb1), 2 (nb2) and 3 (nb3) are shown by circles. (B) The bubble radius profiles for simulations at 340 K and 350 K. At 350 K, initially three bubbles formed. The smaller bubbles, nb2 and nb3, gradually shrank by dissolving in water. The larger nano-bubble remained stable for about 20 ns and it also eventually began to dissolve in water and collapsed. ..................................................................................................... 68 Figure  5-7. Number of methane molecules in different phases during simulation at 350 K. Overall, the number of molecules in the gas phase increases, indicating release of methane to the gas reservoir as a result of hydrate decomposition. ........................................................................................... 69 Figure  5-8. (A) yz-snapshot of the initial state of the simulation with methane mole fraction of 0.074 in a water-methane mixture. (B) yz-snapshot of the simulation after 1 ns of NPT simulation at 350 K and 400 bar with the initial configuration shown in (A). The formation of bubbles and the subsequent depletion of methane in the liquid phase are evident. (C) yz-snapshot of the simulation with dissolved methane mole fraction of 0.042. (B) yz-snapshot of the simulation after 5 ns of NPT simulation at 350 K and 400 bar with the initial configuration shown in (C). No bubbles formed in this simulation and the methane molecules diffused into the gas reservoir. ........................................... 72 xvi Figure  5-9. (A) yz-snapshot at 4 ns of the partially decomposed hydrate at 350 K. Water molecules initially in the hydrate and liquid phase are shown by hydrogen-bonds in red and blue, respectively. For clarity, only methane molecules in the nano-bubble are shown. One of the methane molecules in the hydrate phase is also shown as a reference of position inside the hydrate phase. (B) yz-snapshot of the hydrate growth simulation at 280 K after 50 ns. As seen, hydrate cages have grown on both sides of the initial hydrate template. However, more hydrate cages are formed on the right hand side where initially a bubble was present. ........................................................................................... 74 Figure  5-10. F3 profile of the hydrate growth simulation at 280 K and 400 bar. L4 to L10 and L11 to L17 refer to the slices on the non-bubble side and bubble side, respectively. As can be seen, three layers of hydrate cages have grown on the bubble side whereas there are only two layers of hydrate growth on the other side of the initial hydrate template. ........................................................................... 75 Figure  6-1. yz-projection of the initial configurations of the case I (A) and case II (B) simulations. The AFP is initially aligned with the long axis along the x-direction. wf-AFP, methane in the gas phase and water in the liquid phase are represented by a cyan ribbon, cyan spheres and blue lines, respectively. Hydrate water and methane are shown as red hydrogen bonds and cyan spheres in the middle of the box. Simulation box sizes along the x, y and z-directions are 8.421, 8.421 and 15.639 nm, respectively. (C) Relative position of wf-AFP with respect to the hydrate block in case I simulation. The protein is roughly 0.1 nm away from the hydrate surface. Same relative distance is used for the case II simulation. The cyan helix represents the wf-AFP. For clarity only side chains of the threonine residues are shown. Liquid water molecules and gas phase methane are not shown. Left pane: xy view; top right: xz view; bottom right: yz view. ................................................................... 78 Figure  6-2. Root mean square deviation of the backbone of the protein during the wf-AFP/Water simulation. RMSD remains below 0.3 nm suggesting no significant change in the α-helical structure of wf-AFP under simulation conditions of 275 K and 800 bar. ................................................................ 80 Figure  6-3. (A) Initial (0ns) and (B) final (5 ns) configurations of the wf-AFP/Water simulation. α–helical structure of wf-AFP is evident in the middle of the simulation box. The size of the cubic box xvii is 8.421 nm and blue lines represent water molecules. x, y and z axes are colored as red, green and blue, respectively. (C) Radial distribution functions for water oxygen with the THR hydroxyl oxygen (OW-OG1, black), the THR terminal methyl carbon (OW-CG2, red) and other water oxygens (OW-OW, blue). The first peak of the OW-OG1 RDF occurs at 0.282 nm, similar to the OW-OW first peak, while the first peak of OW-CG2 occurs at 0.37 nm due to hydrophobic characteristics of THR methyl group. .................................................................................................................................................... 81 Figure  6-4. (A) yz- and (B) xy-projection and (C) zoomed in yz-view of the case I simulation at 200 ns. The binding residues, 2THR, 6ALA and 9ALA, are shown in yellow in panel (B) and as thick bars in panel (C). Liquid water is not shown in panel (B) for clarity. The entrapment of binding residues in the empty half-cages is evident in panel (C). Hydrogen bonds are shown as red dashed lines. ........... 82 Figure  6-5. (A) The structure of wf-AFP at different times during the case I simulation. Cyan, red, black, orange and transparent yellow correspond to 5, 80, 120, 140 and 170 ns, respectively. The blue line on the left shows the approximate position of the hydrate/water interface while the one on the right shows that of the water/gas interface. Other molecules are not shown for clarity. (B) RMSD of the wf-AFP during the case I simulation. There is a significant change in the backbone structure of the protein after 60 ns. This substantial structural change can be seen as unfolding of the section of the protein which is not adjacent to the water/hydrate interface. ................................... 83 Figure  6-6. (A) yz- and (B) xy-projection and (C) zoomed in yz-view of the protein 1 in the case II simulation at 100 ns. The binding residues, 24THR, 28ALA and 31ALA, are shown in yellow in panel (B) and as thick bars in panel (C). Liquid water is not shown in panel (B) for clarity. The entrapment of binding residues in the empty half-cages is evident in panel (C). Hydrogen bonds are shown as red dashed lines. .............................................................................................................................. 85 Figure  6-7. (A) F3 values for the side chains of all residues of wf-AFP for the case I simulation at 172.5 ns. 2THR, 6ALA and 9 ALA have F3 values below 0.05 suggesting they are trapped within the half-cages at the hydrate/water interface. Refer to  Appendix F for the naming scheme of the atoms. (B) F3 vs. time profiles of the methyl side chain of the four THR residues and (C) ALA 3, 6,7 and 9 xviii of the wf-AFP in case I simulation. The adsorption of wf-AFP to the hydrate surface is a collaborative action between THR and ALA residues. In this case: 2THR, 6ALA and 9 ALA. (D) A snapshot of the case I simulation at 195 ns which shows that the methyl group of 2THR is trapped inside a half-cage at the hydrate surface and its hydroxyl side chain forms two hydrogen bonds, shown as green lines (emphasized by arrows), with local water molecules. ........................................................................ 86 Figure  6-8. F3 values of all side chains of wf-AFP residues for case II simulation at 67.5 ns. As seen 24 THR, 28 ALA and 31 ALA (designated by red circles) as well as 35 THR for protein 1 have F3 values well below 0.0525 (red broken line). This suggests these residues have an ordered water neighborhood and are trapped within the half cages at the hydrate surface. Protein 2 is not effectively anchored to the hydrate surface at this time. ............................................................................ 88 Figure  6-9. F3 vs. time profiles of the methyl side chain of (A) all of the THR residues and (B) ALA 28, 31, 32 and 36 of protein 1 in case II simulation. The adsorption of wf-AFP to the hydrate surface is a collaborative action between THR residues and ALA residues, in this case 24THR, 28ALA and 31 ALA. F3(t) profiles of the methyl side chain of (C) all of the THR residues and (D) ALA 3, 7, 14 and 32 of protein 2. During the 135 ns of the simulation, this protein is not effectively bound to the hydrate surface. The protein 2 THR residues do not illustrate F3 values which stably remain below the interfacial value of ~0.05. ......................................................................................................................... 89 Figure  6-10. The yz-projection of three snapshots of the case I simulation at (A) 195.418 ns, (B) 195.608 and (C) 195.800 ns. The 2THR methyl group is trapped in the center of the half-cages at the water/hydrate interface while there are 1, 0 and 2 hydrogen bonds (shown in green lines) between the hydroxyl group and nearby water molecules. ........................................................................................ 91 Figure  A-1. The original hexagonal unit cell (parallelogram) and the new selected orthorhombic simulation cell (yellow rectangle). Si (blue), O (red). ............................................................................... 108 Figure  A-2. The hydroxylated cubic unit cell of silica. Si (blue), O (red), H (white). The (010) basal plane is hydroxylated. ................................................................................................................................... 108 xix Figure  B-1. The initial setup configurations for the simulation 3 of Table  3-2 (water and methane gas) and 5 (water and carbon dioxide gas). ............................................................................................... 110 Figure  B-2. Snapshots of CH4 and CO2 gas distribution adjacent to a flat water surface. The gas concentration enhancement near the interface is evident. ...................................................................... 110 Figure  C-1. The F3 order parameter along the layers in the y-direction, see Figure  4-1. .................... 111 Figure  D-1. zz-component of the pressure tensor during methane hydrate decomposition simulation (G-W-H-W-G) at 350 K and 400 bar. Red line represents the 20-points moving average trend line. ......................................................................................................................................................................... 112 Figure  D-2. zz-component of the pressure tensor during the simulation with dissolved methane mole fraction of 0.042 at 350 K and 400 bar. Red line represents the 20-points moving average trend line. .................................................................................................................................................................. 113 Figure  E-1. Initial (0 ns) and final (50 ns) configurations of the system. The pressure of the simulation is 800 bar and hydrate is grown at 305 K. Water molecules are shown by red (oxygen) and white (hydrogen) bonds. The hydrogen bonding between water molecules in the hydrate phase is also shown by red dashed lines and methane molecules are shown by cyan spheres. The size of simulation box at t=0 ns along the x, y and z directions are 0.24, 0.24 and 0.665 nm. ...................... 114 Figure  E-2. Potential energy profile for hydrate growth simulation at 305 K and 800 bar with different time steps. The energy trends show similar behavior and the end state at 50 ns is essentially the same. ......................................................................................................................................................... 115 xx  List of Symbols ࢘௜ Position vector of the particle i in space ݎ௜௝ Distance between particles i and j ࢜௜ Velocity vector of particle i ࢇ௜ Acceleration vector of particle i ࡲ௜ Pair-wise sum of forces acting on particle i arising from all possible pairs containing the particle i  ݉௜ Mass of particle i N Number of particles ݐߜ Time step ܷ Potential of the system ݇௕ Force constant for the harmonic oscillation of bonds ݎ௘௤ Equilibrium bond distance ݇ఏ Force constant for the harmonic oscillation of bond angles ߠ௘௤ Equilibrium bond angle ߶ Dihedral angle ߝ L-J potential depth ߪ Distance at which the inter-particle L-J potential is zero ݍ௜ Partial charge of particle i ݇௘ Coulomb’s constant ௗܰ௙ Number of degrees of freedom ௖ܰ Number of constraints ܸ Volume ݇஻ Boltzmann constant ܲ Pressure ܶ Temperature ܨଷ Order parameter ߠ௝௜௞ The angle formed by a triplet of water oxygen atoms where i is in the center Å Angstrom (10-10 m) xg Mole fraction of dissolved methane in liquid water xgc Critical mole fraction of dissolved methane in liquid water  xxi List of Abbreviations sI Structure I hydrate sII Structure II hydrate sH Structure H hydrate STP Standard temperature and pressure conditions vdW Van der Waals L-J Lennard-Jones MD Molecular dynamics NVT Constant volume constant temperature molecular dynamic ensemble NPT Constant pressure constant temperature molecular dynamic ensemble NVE Constant volume constant energy molecular dynamic ensemble ns Nano (10-9) second ps Pico (10-12) second fs Femto (10-15) second nm Nano meter HBCF Hydrogen bonding correlation function HW Hydrate/Water scenario SHW Silica/Hydrate/Water scenario SWHW Silica/Water/Hydrate/Water scenario G-W-H-W-G Gas-Water-Hydrate-Water-Gas simulation set up CDMM Critical dissolved methane mole fraction KI Kinetic inhibitor AFP Antifreeze protein wf-AFP Winter flounder antifreeze protein ALA or A Alanine residue THR or T Threonine residue VAL Valine residue SER Serine residue IBS Ice binding surface RDF Radial distribution function RMSD Root mean square deviation  xxii Acknowledgements I wish to express my sincere appreciation to the following people for their professional support as well as friendly assistance, without which it would have been very hard to accomplish this thesis: My supervisors, my family, my research colleagues, my friends and the CHBE staff I would like to acknowledge my supervisor, Prof. Peter Englezos, for encouraging and supporting me to grow as an independent researcher and for giving me the freedom to be involved in the extracurricular professional activities and advancing my soft skills as well. I am also very thankful to Dr. Saman Alavi and Dr. John Ripmeester for their co-supervision and in particular for challenging me with many stimulating feedbacks. I owe a special thanks to Dr. Alavi. You have been indeed very kind to meet with me through skype even after your working hours in Ottawa. I learned a lot from you both scientifically and professionally. I owe a great deal to my research colleagues, Hassan Sharif, Duo Sun, Negar Mirvakili and Nagu Daraboina. You all have been very helpful and friendly to me. I have had a very enjoyable time with you and meaningful discussions in different respects. I also wish to thank my friends, Mehdi Bazri, Amir Dehkhoda, Mohammad Masnadi and Hafiz Rahman for their true friendship and often chats which have created many pleasant memories from the graduate school for me. My debt to my family is substantial, particularly my wife and my mother who have always been supportive and encouraging in every way. To all of you, I wish to say I truly appreciate your unconditional love and support. Finally, I would like to thank the University of British Columbia for the Four Year Doctoral Fellowship award (4YF) and Compute/Calcul Canada for providing the computational facility through the Westgrid cluster. xxiii Dedication  To my family My dear wife, Motahareh, and my lovely son, Seyyed Mohtadi Seeing you every day is delightful and energizing for me and will always be…ChapteThis disseprocessestheir sign1.1 WhGas hydrmoleculeand near undersea such as CencapsulaSloan & Kinterestin Figure  1-1accommomoleculegreen, reshttp://chHydrate cguest mostructure r 1: Introdrtation focu by conductiificance and at Are Gaates are crysts known as tfreezing temsediments, aH4, C2H6, C3ting the gueoh 2008; Ggly have sim. (A) A burndating a mes. (C) Structupectively. Memwiki.ucdrystals have lecule. ThreeI (sI), cubic uction ses on the stng moleculathe basics ofs Hydratesalline solid she guests. Unperature), usnd permafroH8, and COst molecules ale & Steed ilar physico-ing methanethane molecre I hydrateethane moleavis.edu ordered stru common crstructure II udy of the fur dynamic si the molecu? olutions comder favorabually met in st deposits, 2, freeze to fo(Davidson 12012). Gas hchemical pro-propane hyule. Red line. The large (cules are repctures consiystalline stru(sII) and the ndamentalsmulations. Tlar dynamicsposed of wle thermodyhydrocarbonwater molecurm a solid c973; Englezydrates resemperties excepdrate. (B) A s represent h51262) and smresented by sting of diffectures have hexagonal st of gas hydrahis chapter i (MD) technater and appnamic condi transport ales in the prage-like netwos 1993; Ripble the phyt for the th schematic oydrogen boall (512) cagespheres. Sourent types obeen found ructure H (ste dissociatintroduces thique. ropriate hydtions (relativnd processinesence of guork, via hydmeester 200sical appearermal conduf the small cnding betwes are coloredrce: f cages depenfor hydrates:H). These ston and inhibe gas hydratrate former ely high presg facilities, est moleculerogen bond0; Sloan 200ance of ice activity. age of sI hyden the water as red and ding on the the cubic ructures diff1 ition es, sure s ing, 3; nd rate   er in the size aexample, moleculeThe smalpentagonstoichiomIn other wgas solubFigure  1-the pressuformers. that evenimmediattakes for Figure  1-2(MeOH),pressuresnd the geommethane fors containing l and large cas plus 2 hexaetry of the sords, the mility in liquid2 shows the re-temperatHydrates are when the P/ely. The onsthis to happe. The phase as a thermo and lower teetry of the cms the sI clatwo small cages of the sIgons (51262),I methane hole fraction  water (This gas-liquid waure plane. T stable only T conditionet of hydraten is known  boundary ofdynamic hydmperatures ages and thethrate hydrages (12-sided clathrate hy respectivelyydrate with iof methane wis discussed ter-hydrate he general tron and abovs favor the h crystal nuclas the induc methane hyrate inhibito(red).  number of cte and its cu) and six lardrate are for, with water ts cages fullyould be 0.1in greater dethree phase send is similae the equilibydrate formeation is a sttion or nucledrate in purr, pushes thages forminbic unit cell ige (14-sidedmed by 12 poxygen at th occupied w48 which is tail in sectiotability bounr for other lirium curve. Hation, crystaochastic pheation time.  e water (bluee hydrate forg the unit ces made up o) cages (see Fentagons (5e vertices. Tould be 8:46much highern  2.3 of Chadary for meght hydrocarowever, it ls of hydratenomenon an). Adding mmation condll. As an f 46 water igure  1-1).12), and 12 herefore, the (methane:w than its regpter 2). thane hydrabon hydrateshould be no may not nucd the time thethanol itions to hig2  ater). ular te in  ted leate at her 3 Interestingly, formation of hydrate from a solution that has previously experienced a hydrate formation or freezing is generally associated with a shorter induction time compared to the formation of hydrate from the fresh solution. This observation is referred to as the “memory effect” in gas hydrate context.  1.2 Significance of Gas Hydrates In practice, gas hydrates represent an opportunity as well as a challenge. 1.2.1 A Potential Source of Natural Gas Gas hydrates were first synthesized in 1800’s (Davy 1811). After the discovery of naturally occurring methane hydrates (Makogon et al. 1971), it was realized that this unconventional source of methane and natural gas can potentially serve as an energy source for the future (Holder et al. 1984; The Council of Canadian Academies 2008; Committee on Assessment of the Department of Energy, Methane Hydrate Research and Development Program: Evaluating Methane Hydrate as a Future Energy Resource 2010). Worldwide estimates of methane stored in the form of hydrate suggest that the hydrate deposits have the potential to produce natural gas (Hester & Brewer 2009) for at least 300 years based on the current global consumption rate. However, the challenge is to find safe, controllable, and economical strategies to exploit the gas from these vast reserves. Very recently, Japan has successfully performed the first extraction of natural gas hydrates from beneath the seabed (Tabuchi 2013). So far, depressurization, thermal stimulation, inhibitor injection and guest exchange (Ohgaki et al. 1996; Park et al. 2006; Baldwin et al. 2009) have been proposed as extraction techniques. 1.2.2 A Flow Assurance Problem It was reported in 1930’s (Hammerschmidt 1934) that under common operating conditions, gas hydrates can form and plug up the natural gas transmission lines. It was then realized that hydrates can represent a serious safety and operational challenge to the oil and gas industry as they could form inside the transmission pipelines and processing facilities, therefore blocking the flow of fluids 4 and potentially leading to explosions. The cost of flow assurance due to hydrate formation is estimated to be over 500 million USD/yr (Anderson et al. 2005). The oil and gas companies are still searching for effective hydrate formation inhibiting chemicals which are both economically and environmentally sound. A common inhibitor used in industry is methanol. As shown in Figure  1-2, methanol inhibits the formation of hydrates by shifting the equilibrium curve to more severe conditions (i.e. higher pressure and lower temperatures). 1.2.3 Other Applications of Gas Hydrates Gas hydrates have also the potential to be utilized in the following technologies: • energy storage and transportation (Hao et al. 2008; Kumar, Englezos, et al. 2009) • gas separation for pre- and post-combustion capture of CO2 (Linga et al. 2007; Kumar, Linga, et al. 2009) • CO2 sequestration in geological formation and depleted gas reservoirs (Zatsepina & Pooladi-Darvish 2010) • water desalination (Max & Pellenbarg 1999) 1.3 Hydrate Research Generally, practical challenges in hydrate research can be classified into two categories: control and prediction. Control relates to hydrate formation, decomposition and inhibition in connection with flow assurance or gas extraction from natural hydrate deposits. Prediction relates to where and how much natural gas is stored in the form of hydrate in natural formations around the globe. To overcome these challenges, hydrate research, as a multi-disciplinary activity, is carried out over a spectrum of length scales (see Figure  1-3).(Ripmeester 2000) Figure  1-3The molelength scamilliseconthe experof the simGas hydrdissertatioperspectivprocesseswhich is ofrom thisunderstanIt should (number restricted. The lengthcular mechales. Computd time-scaleiments. Interulations andates have been is to invese and condu. Molecular therwise ve work will adding of the be noted thaof molecules by the availa and time scnisms involvational sciens; whereas, eesting disco experimentn extensiveltigate the dict computesimulation isry challenginvance the mhydrate-relatt molecular  that can be ble computaales over whed in hydrattists are pusxperimentalveries can bes. y studied expssociation anr simulations an excellentg or not posolecular-scaled processesdynamics (Msimulated ustional poweich hydrate e-related prohing molecuists are attem made by stuerimentally d inhibition to determin technique thsible to capte knowledge (formation,D) simulatioing MD) andr. How well research is bcesses take plar simulatiopting to enhdying naturas well as th of gas hydrae the mechaat renders fure experime about hydr dissociationns have som the duratiothe forces co eing carriedlace at the nns on larger ance the timal phenomeneoretically. Ttes from thenisms goverundamental ntally. The iates and imp and inhibitie limitationn of the simnsidered in  out. ano-time anmicro- and e resolutiona at the ovehe aim of th molecular ning these information nsights arisinrove our on). s: the size ulation are these simula5 d  of rlap is g tions 6 represent the real interactions between the atoms and/or molecules is also a very important factor in capturing the system behavior correctly. Fortunately, with the help of the supercomputers, molecular simulations of methane hydrate systems with several tens to hundreds of thousands of particles (104-106) up to the micro-second time scales are now feasible. The simulations performed in this thesis are within the range of 20000-100000 particles (length of the simulation box size range: 6-22 nm) and a time scale of 2-200 nano-seconds. 1.4 Basics of the Molecular Dynamics Simulation From the time of Newton the deterministic mechanical interpretation of the Nature, where it is believed that the system behavior can be computed if we have the interaction forces plus a set of initial conditions for the system’s parts, has driven the science.(Haile 1997) Molecular dynamics is a computer simulations technique (in silico) where the time evolution of a set of interacting particles, i.e. their trajectory, is followed by integration of their Newtonian equations of motion.(Frenkel & Smit 2001) Through the application of statistical mechanics macroscopic properties such as temperature, pressure and free energies can be extracted from the trajectory.  The procedures followed in studying a system either theoretically or experimentally are similar. We adjust certain inputs, monitor the system behavior and then outputs are measured or computed. A molecular-scale simulation can be broken to 5 steps: (1) choosing an appropriate force model, (2) constructing an initial configuration of the system, (3) performing an equilibration phase during which the values of a set of monitored bulk properties become stable, and (4) running the production phase during which the trajectories are calculated, and (5) analysis of the trajectory.(Leach 2001) Molecular dynamics can be done in two ways: “Equilibrium” where the macroscopic property of interest, according to the ergodic hypothesis, can be calculated from the time average of that property during the simulation. On the other hand, in “Non-equilibrium” mode, the system is stimulated away from the equilibrium and the system response is followed.(Haile 1997) One of the common applications of MD is to predict the properties of materials. However, it could be used as a purely exploratory tool as well.(Frenkel & Smit 2001) The main goal of the simulations performed in 7 this thesis is not to quantitatively predict the properties of the system but is more in the nature of providing explanations by determining the mechanisms involved.  For an N-particle system the trajectory is obtained by simultaneously solving the differential equation of Newton’s second law for each particle: ݀ଶ࢘௜݀ݐଶ = ࢇ௜ =ࡲ௜݉௜ 		; ݅ = 1,2, … , ܰ Equation  1-1Here,	࢘௜ is the position of particle i in the space and ࡲ௜ is the vector sum of the forces acting on particle i arising from all other particles. An analytical solution of the equations of motion for an N-body system is not available and thus numerical techniques such as the finite difference method have to be used. The idea is to break the time evolution of the movement of the particles into many small time steps that are ݐߜ apart, during which the force is assumed to be constant. Taylor series expansion of the position yields: ݐ)࢘ + ݐߜ) = ݐ)࢘) + ݀ݐ݀(ݐ)࢘ᇣᇤᇥ	࢜(௧)ݐߜ + ଵଶ	݀ଶݐ)࢘)݀ݐଶᇣᇤᇥ	ࢇ(௧)ݐߜଶ + ଵ଺݀ଷݐ)࢘)݀ݐଷ ݐߜଷ + ݐߜ)ࣩସ) Equation  1-2Where ݐߜ)ࣩସ) represents the other terms of the expansion with the smallest power of time step being 4. One of the common algorithms for integrating the equations of motion is the Verlet algorithm.(Verlet 1967) Writing the series expansion for backward time step for positions and adding it to Equation  1-2 gives: ݐ)࢘ + ݐߜ) = 2ݐ)࢘) − ݐ)࢘ − ݐߜ) + ݐߜ(ݐ)ࢇଶ + ݐߜ)ࣩସ) ݐ)࢜) = ݐ)࢘ + ݐߜ) − ݐ)࢘ − ݐߜ2(ݐߜ + ݐߜ)ࣩଶ) Equation  1-3 Equation  1-4One of the disadvantages of the Verlet algorithm is that the velocity does not appear explicitly in Equation  1-3. The velocity at time t is calculated by dividing the position difference at t+δt and t-δt by 2δt. The other drawback is that the position is obtained by adding a small term (ݐߜ(ݐ)ࢇଶ) to the difference of two larger terms (2ݐ)࢘) − ݐ)࢘ − ݐߜ)) which can lead to a loss of precision. 8 A common variation of the Verlet algorithm to resolve the above drawback is the so-called leap-frog algorithm where the positions and velocities are calculated with half a time step lag between them:  ݐ)࢘ + ݐߜ) = ݐ)࢘) + ࢜ ቀݐ + ଵଶݐߜቁ ݐߜ ࢜ ቀݐ + ଵଶݐߜቁ = ࢜ ቀݐ −ଵଶݐߜቁ +ݐ)ࡲ)݉ ݐߜ  Equation  1-5 Equation  1-6These computationally expensive calculations are usually done by software packages. In this thesis, dl_poly (Smith & Forester 1996) and gromacs (Hess et al. 2008) software packages have been used to perform the computations.  1.4.1 Calculation of Forces In molecular dynamics simulations, it is assumed that the forces are pair-wise additive and thus must be evaluated between N(N-1)/2 pairs. This makes calculation of the forces the most computationally expensive part of an MD simulation which scales with N2. It should however be noted that in practice by applying the cut off radius the number of interactions that need to be calculated dramatically decreases (Leach 2001). Forces, ࡲ, are calculated by taking a negative derivative of the potential of the system, ܷ, with respect to the position.  ࡲ = −߲ܷ߲࢘  Equation  1-7In classical mechanics, the interactions between atoms and molecules are categorized into two groups, namely, the bonded and non-bonded interactions.  Therefore, depending on the type of the atoms and molecules of interest, different types of potential functions can contribute to the total potential of the system including bonded interactions such as bond stretching, angle bending, and dihedrals, and, non-bonded interactions such as electrostatic and the van der Waals. For simple molecules such as water and methane which are mostly used in this thesis, the internal degrees of freedom of the molecules are considered fixed, and the total potential of the system is defined as the sum of the electrostatic and Lennard-Jones (L-J) potentials. 9 It is usually common that for a class of molecules which are encountered in a special field of study, a particular force field is constructed. In principle, this force field is a sum of different types of potentials. For example, in biology, the AMBER family (Cornell et al. 1995) is a popular force field for the study of proteins: ஺ܷெ஻ாோ = ෍ ݇௕൫ݎ − ݎ௘௤൯ଶ௕௢௡ௗ௦+ ෍ ݇ఏ൫ߠ − ߠ௘௤൯ଶ௔௡௚௟௘௦+ ෍ ௡ܸ2 [1 + cos(݊߶ − ߛ)]ௗ௜௛௘ௗ௥௔௟௦ᇣᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇤᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇥ௕௢௡ௗ௘ௗ	௜௡௧௘௥௔௖௧௜௢௡௦+෍4ߝ௜௝ ቈ൬ఙ೔ೕ௥೔ೕ൰ଵଶ− ൬ఙ೔ೕ௥೔ೕ൰଺቉௜ழ௝+෍݇௘ݍ௜ݍ௝ݎ௜௝௜ழ௝ᇣᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇤᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇧᇥ௡௢௡ି௕௢௡ௗ௘ௗ ௜௡௧௘௥௔௖௧௜௢௡௦ Equation  1-8where ݇௕ and ݇ఏ	are the force constant for the harmonic oscillation of bond stretching and angle bending, ݎ௘௤ and ߠ௘௤ are the equilibrium bond length and angle, ߶ is the dihedral angle and ݊, ௡ܸ and ߛ are the dihedral potential parameters, i and j refer to the interacting atoms on different molecules, ߝ is the L-J potential depth, ߪ is the distance at which the inter-particle L-J potential is zero, ݎ is the distance between the particles, q is the partial charge of each particle and ݇௘ is the Coulomb’s constant. The parameters of these force fields, such as ݇௕, ݎ௘௤, ߝ and ߪ, are optimized to correctly reproduce molecular structure/properties and bulk material properties such as melting point and density. For instance, TIP4P-ice model of water fits the experimental phase diagram of water better than the other available 4-site models.(Abascal et al. 2005) 1.4.2 Calculation of Physical Properties Two main physical properties of interest in simulations as well as experiments are temperature and pressure. A brief description of how these properties are calculated in an MD simulation is given below. 1.4.2.1 Temperature The temperature is related to the kinetic energy, ܧ௞௜௡, of the system: 10 ܧ௞௜௡ = ଵଶ෍௠೔௩೔మಿ೔సభ Equation  1-9where mi and vi are the mass and velocity of the particle i. If needed, the temperature of a molecular simulation can be controlled via the operation of a computational “thermostat” by scaling the atomic velocities. From the classical theory of gases, it can be shown that the absolute temperature is related to the kinetic energy as below (Leach 2001): ଵଶ ௗܰ௙݇஻ܶ = ܧ௞௜௡ Equation  1-10where ௗܰ௙ is the degrees of freedom. For ideal gases the degrees of freedom is three times the total number of atoms due to having no degrees of freedom for vibration and rotation. In general: ௗܰ௙ = 3ܰ − ௖ܰ Equation  1-11where ௖ܰ is the number of constraints such as those used for fixing a bond length or a bending angle. An alternative way to keep the temperature constant in a simulation is to couple the system to an external hat bath that supplies or removes heat from the system as necessary. However, this strategy does not generate rigorous canonical averages (Leach 2001). To address this issue, the method of “extended system” was introduced by Nose (Nosé 1984) which was subsequently developed by Hoover (Hoover 1985). The Nose-Hoover method considers the heat bath as an integral part of the system and is represented by an additional degree of freedom. The potential and kinetic energy of the heat bath is added to the potential and kinetic energy of the system and the dynamics of this extended system will be monitored for adjusting the temperature. For more details on this topic please refer to (Frenkel & Smit 2001). 1.4.2.2 Pressure Pressure can be calculated using the following equation (Leach 2001): 11 ࡼ = ܸܰ ݇஻ܶ +13ܸ෍෍࢘௜௝. ࡲ௜௝௝ழ௜௜ Equation  1-12where V is the volume and ݇஻	is the Boltzmann constant. The first term is related to the kinetic energy and the second term is the virial contribution. Again, if needed, the pressure of a simulated system can be controlled via the operation of a computational barostat by regulating the simulation cell volume.  1.4.3 Treatment of the Long-Range Electrostatic Interactions The Coulomb potential decays very slowly with distance (it is proportional to 1/r). This leads to a significant loss of accuracy if the Coulomb interactions are to be truncated at a cut-off distance. A common reliable technique to treat the long-range electrostatic interactions properly is the Ewald summation (Ewald 1921). In this technique the electrostatic potential is not truncated and particle i is assumed to interact with all of the other charged particles as well as all its periodic images, but not itself.(Frenkel & Smit 2001) ܷ࢒࢛࢕࡯. =12෍ݍ௜෍ݍ௝ห࢘௜௝ + ܮ࢔หᇱ௝,࢔ே௜ୀଵ Equation  1-13Here the prime on the summation indicates that the summation runs over all periodic images n and over all particles j, except j = i when n = 0. L is the simulation box vector. In Ewald summation, each charge is assumed to be surrounded by a neutralizing distribution (commonly Gaussian) of a charge with equal magnitude but of opposite sign. A second charge distribution will be added to the system to keep the charge neutrality of the simulated system. This divides the summation into two parts that converge more rapidly compared to the original summation. The contribution from the particle charge and the neutralizing Gaussian charge distribution decays rapidly with respect to r and its value becomes negligible after some cut-off distance.  Using the Fourier transformation, the contribution from the second charge distribution is treated in the reciprocal space (the space in which the Fourier transform of a spatial function is 12 represented). The Ewald summation is computationally expensive and scales with ࣩ(N3/2). A number of modifications have been developed to reduce the computational cost of this technique. The particle-mesh Ewald (PME) summation (Darden et al. 1999) is an attractive alternative for an intermediate size system that scales with ࣩ(NlogN).  1.4.4 Ensembles Used Molecular dynamics simulations can sample the system phase space, i.e. space of atomic positions and momentums (velocities), in different ensembles where a specific set of parameters are controlled depending on the external conditions of interest. The ensembles used in this thesis are given in Table  1-1. Table  1-1. Different ensembles used in this thesis ensemble Constant entitiesNVE (micro-canonical) Number of particles, volume, total energyNVT (canonical) Number of particles, volume, temperatureNPT (isothermal isobaric) Number of particles, pressure, temperature 1.4.5 The General Algorithm for MD The general algorithm of an MD simulation is summarized below in Figure  1-4. Figure  1-4The time the convesystem wdictated bbond, vibfs. It is suof the sysAdditionaSHAKE . The generstep is usualrgence of thith a tolerancy the fastestrate at a freqggested thattem.(Leach 2lly, these fas(Ryckaert et al algorithm ly in the range total energe of 1 part i motion presuency of ab a reasonable001) Thereft vibrating bal. 1977) to bof moleculare of 0.5-2 fsy of the systn 102 -103. Tent in the syout 3000 cm time step foore, time steonds can bee able to ap dynamics s. The acceptem or other he appropriastem. Bonds-1. This is equr MD simulps of about  constrainedply larger tim imulations ance of the macroscopicte time step  containing ivalent to a ations is one1 fs are com using differe steps. result is asse average profor each simlight atoms, period of ap tenth of themonly used ent algorithmssed by checperties of thulation is such as C-Hproximately fastest motin the literatu such as 13 king e   10 ion re. 14 Chapter 2: Literature Review Over the last two decades, molecular dynamics has been shown to be an excellent tool that can provide useful insights into hydrate processes of nucleation (Moon et al. 2003; Walsh et al. 2009; Debenedetti & Sarupria 2009; Liang & Kusalik 2011), dissociation (Bagherzadeh et al. 2012; Bagherzadeh et al. 2013) and inhibition (Anderson et al. 2005), with results which correlate well with the experimental data.(Anderson & Jordon 2009) What follows is a summary of the key findings from the MD simulations in these areas. 2.1 Molecular Simulations of Hydrate Formation and Growth Due to the extremely low solubility of methane in water, many methane hydrate formation simulations that focus on the homogeneous nucleation mechanism, applied some type of methane delivery, usually in the form of artificially supersaturated methane solutions, to study the methane hydrate formation mechanism.(Razul et al. 2005; Vatamanu & Kusalik 2007; Walsh et al. 2009; Liang & Kusalik 2010; Vatamanu & Kusalik 2010) These simulations are useful in providing insight into the structures of hydrate nuclei that form in solution and the growth of hydrate after nucleation. However, supersaturation of the aqueous solution with respect to methane may alter the hydrate formation mechanism compared to the mechanism that naturally occurs.  Simulations of hydrate nucleation for heavily supersaturated aqueous solutions of methane for comparatively long time (micro-second time scale) reveal that the hydrate crystal nucleation is a cooperative process between the guest and water molecules and is induced by the close relative position of methane guests.(Walsh et al. 2009) However, it should be noted that the resulting structure was a mixture of sI and sII linked through unusual 51263 cages. Liang and Kusalik (Liang & Kusalik 2013) reported instantaneous nucleation of the H2S sI hydrate. Their NVE simulation did support the two-step mechanism proposed through constant temperature simulations, namely NVT and NPT, where firstly an amorphous structure forms which subsequently evolves to the stable solid crystalline hydrate structure under the temperature and pressure conditions of the simulation. However, in their simulation, the nucleation rate slowed down due to the heat release of the 15 nucleation event. A combined coarse-grained molecular dynamics and classical nucleation theory modeling showed that under realistic conditions the homogeneous nucleation is highly unlikely as the nucleation rate was calculated to be extremely low and the critical nuclei size was very large.(Knott et al. 2012) Other studies simulate the hydrate formation process in the presence of silica walls to mimic the natural conditions of hydrate formation in sediments. Such simulations have been performed for CO2 and CH4 hydrate. Methane hydrate formation, again under supersaturated methane solutions, showed that crystal growth in bulk solution (away from the silica surfaces) proceeded ahead of that adjacent the silica surface.(Liang et al. 2011) Bai and coworkers (Bai et al. 2011) proposed a three-step mechanism for CO2 hydrate formation in the presence of a silica surface: first, an ice-like layer of water molecules formed close to the solid surface. During the second stage, an intermediate layer formed which acted as nucleation seed for the hydrate formation on top of this layer. This buffer layer also alleviated the mismatch between the arrangement of water molecules in the ice-like layer and the hydrate phase. Upon completion of the intermediate layer, CO2 hydrate formed in a layer-by-layer fashion in the (100) and similar families of planes, provided that there is enough space for its growth. It was found that there should be a minimum distance of ~36 Å between the silica surfaces so that the formation of hydrate could occur. Again, it should be noted that the nucleated structure is different than the perfect sI hydrate as the ratio of the number of large to small cages was different than three. The authors focused on the regular small and large cages of sI hydrate and did not discuss formation of any irregular cages. 2.2 Molecular Simulations of Hydrate Dissociation Hydrate dissociation is an endothermic phase transition. So, in order to sustain the hydrate decomposition process, a source of heat must be continually provided: ܪܥସ. ݊ܪଶܱ	(ℎ݁ݐܽݎ݀ݕ) → ܪܥସ(݃) + ݊ܪଶܱ(ܽݍ. ) ; Δܪ > 0  Equation  2-1Here, n represents the hydration number. The heat of this reaction at 292 K and over a range of pressures up to 20 MPa is 54.4 ± 1.4 kJ/mol. (Gupta et al. 2008) 16 In principle, hydrate decomposition can be regarded as being similar to ice melting. The temperature of a melting cube of ice placed inside a liquid phase can reach a maximum of 273 K, the solid/aqueous equilibrium value at ambient pressure, while the surrounding can be at higher temperatures providing the heat necessary for melting via heat transfer. Similarly, while hydrate decomposes, one would then expect that the highest temperature that the decomposing hydrate crystals can reach will be the three-phase equilibrium temperature, at the pressure of experiment, and then heat would be transferred to the melting hydrate via the temperature gradient established through the produced liquid water to sustain the phase change. However, in the literature, the majority of the molecular dynamics simulations of hydrate decomposition are performed with the operation of a constant temperature thermostat, keeping the hydrate phase at a temperature higher than the three phase hydrate-water-gas equilibrium temperature, which is therefore not consistent with the physics of this phase transition. Kusalik and coworkers have used a special technique where a temperature pulse is introduced in their molecular simulations therefore providing a “localized” thermostat. (Razul et al. 2005; Vatamanu & Kusalik 2007) A more accurate and realistic approach is to dissociate the hydrate phase in an NVE ensemble where temperature gradients are able to establish.(Alavi & Ripmeester 2010) It should be noted that in isothermal MD simulations, the thermostats remove or provide “heat” to the system (by adjusting the velocities of the molecules to impose the constant temperature condition) and cannot sustain heat transfer, particularly for the size and time scales of the MD simulations of methane hydrate. The temperature and mass gradients play important roles in methane hydrate dissociation under experimental conditions. Even if methane hydrate decomposition occurs in a stirred reactor, nano-scale thermal gradients at the hydrate surface will still persist due to non-slip boundary conditions.  Molecular dynamics (MD) simulations of the dissociation of methane hydrate in contact with water phase have been widely performed.(Walsh et al. 2009; Alavi & Ripmeester 2010; Baez & Clancy 1994; Rodger 1994; Rodger 2000; Yasuoka & Murakoshi 2000; Tse & Klug 2002; C. Moon et al. 2003; English & MacElroy 2004; English et al. 2005; Ding et al. 2007; Hawtin et al. 2008; English & Phelan 2009; Myshakin et al. 2009; Iwai et al. 2010; Ripmeester & Alavi 2010; Sapna Sarupria & Debeneddestructio1994) Thsmaller anThe dissodecomponear the ithe cage-cage.(SaruFigure  2-in the midThe set-uFigure  2-“over-temdescribedvibrationdown theaggregateEnglish aneglectedmolecularhave beenCH4 hydrmechanisWateetti 2011) Mn and the gue melting temd the decomciation fronsition rate wnterface.(Myspecific occupria & Deb1. The commdle, is in cop of these si1). These simperature” h as a two-stes of water m lattice. Duri together (Ynd Phelan (E by running  dynamics si performedate in constam during mer elting of a naest moleculperature ofposition ratt was not shaas attributedshakin et al. pancy as welenedetti 201on configurntact with bumulations usulations areigher than thp process inolecules in thng the seconasuoka & Munglish & Phthe simulatiomulations o. Alavi and Rnt volume cthane hydraHydratno-crystal oes diffused a methane hye was slowerrp and exten to the reocc2009) The rl as the occu1) ation for hydlk water. ually include usually perfe hydrate-w which the ee hydrate lad stage, metrakoshi 200elan 2009), in in thermon methane hipmeester (Aonstant enerte decompose f methane hyway and formdrate was fou for higher cded over 5-upation of tate of CO2 dpancy of therate decomps the hydratormed underater-gas equinhanced diffttice) leads tohane molecu0; Ding et aln the above stated ensemydrate dissoclavi & Ripmgy (NVE) enition where Water drate proceed a gas-ricnd to decreage occupan7 Å and the he methane issociation w cages surro osition simue phase in co the operatiolibrium poinusive behavi cell distortiles escape fr. 2007; Iwai studies the eble (i.e. NViation undereester 2010)semble. Therows of hydeded throughh region.(Baase as the clucies.(Englishoscillation inmolecules inas found tounding a parlations. Thentact with bn of a MD t. Decompoor of host (eon which ulom these bret al. 2010). ffect of heatT and NPT) more realis studied decy observed rate cages pa individual cez & Clancyster size bec & Phelan 2 the  water caviti be dependenticular  hydrate phaulk water (sethermostat asition is xcessive timately breaoken cages aAs recognize transfer is . Recently tic conditionomposition a “concertedrallel to the 17 age  ame 009) es t on se, e t an ks nd d by s of ” interface phase andIn naturapermafrohydrates ipossible cof hydratengineersHong & PFigure  2-media. 2.3 MeAlthoughareas andnot fully substantiafrom aquwhereas, 0.148. Thsaturated and nucleWatedecomposed temperaturl formationsst and belown nature, it ionfiguratione decomposi to design a ooladi-Darv2. The configthane Wat there have b technologieknown. An il concentrateous solutionin methane his is equivaleliquid solutiation very inr  together. Me drops as hi, hydrates us sea floor ses essential to to study sution and gasmore efficienish 2005) uration for her Solubilieen significs, the fundamnteresting anion enhance. As an examydrate phasnt to about on (see Tablteresting subMineraHydratMineraoreover, temgh as 20 K wually occur indiments. He investigate ch systems is release fromt and controydrate decoty Issue ant advancementals of thd somewhatment necessple, CH4 soe with comp5000-fold coe  2-1 for simjects of sciel e l perature graere observe porous envnce, in orderthese proces shown in F the natural lled methanmposition sients in the e hydrate fo mysterious ary to form tlubility (molletely occupincentration ilar data for ntific and enWater dients were d due to hydironments w to understases in the prigure  2-2. Unhydrate resee recovery p  mulations inapplication ormation/deccharacteristiche solid hyde fraction) ined cages, theincrease of mCO2). This mgineering stuestablished irate phase dithin and bend the dissoesence of poderstandingrvoirs can asrocesses.(Mo the presencf gas hydratomposition  of gas hydrrates of wea water is 2.5 methane methane comakes both hdy, and as wn the liquid ecompositiolow the ciation of rous media. the mechansist field ridis 2003; e of porous es in differenprocesses arates is the kly soluble g5×10-5 at STole fraction pared to theydrate formell raises an18 n.  A ism t e still uests P; is  ation  importandecompo Table  2-1Guest SoCH4 2CO2 6 It has beereleased iKlug 200recently rincreasingHowevershould reCH4. In mThis impoand cannincluded out of thephase sim2006; ConFigure  2-reservoirsGas t question ofsition?  . Guest conclubility in liq.55×10-5 (Scha.09×10-4 (Carrn observed nto the liquid2; Vatamanueported that the dissocia, it is intuitivsult in bubblost of theseses an unreot escape to in the simula liquid phasulations (gade & Vega 3. The config to allow theWate how methaentration (muid water at Srlin & Battinoll et al. 1991that during m water phas & Kusalik 2 formation otion rate.(Yaely expectede formation  simulationsalistic conditthe gas phasted system ine and reach ts/water/hyd2010; Tung euration for h released gur ne evolves anole fraction)TP Co 1995) 0.) 0.ethane hydre and form m007; Bagherf bubbles regasaki et al.  that methandue to the li separate gasion where the reservoir. T order to prhe gas phaserate) for hydt al. 2010; Tydrate decoest reach theHydratd reaches th enhancemenoncentration 145 (Anderso139 (Udachinate decompethane-richzadeh et al. 2duces the dis2014) e hydrate diquid phase b spaces weree released mo address thovide a poss (see Figure rate growth ung et al. 20mposition si gas phase.e e gas phase t in hydrate n 2004)  et al. 2001)osition, the g regions.(Fo012; Uddin solved methssociation inecoming sup not initiallyethane mustis limitationibility for th 2-3). Only aand equilibri11; Tung et mulations inWater upon hydratRatio 5686 228 uest methanrrisdahl et al& Coombe ane concent a finite amoer-saturated included in  remain in th, gas spaces e released m few papers um curve pral. 2012) cluding sepGas e e molecules. 1996; Tse &2014) It wasration thus unt of water with respecthe system se liquid phashould be ethane to difconducted thediction.(Na arate gas 19  are   also  t to etup. se fuse ree-da 2.4 MoInspired bfield of gaopposed inhibitors(thermodto prevenlower andnucleationA green cincludingsample AinhibitionFigure  2-(B) Type AFP (D) blue; O, rthreonineThe wintantifreezelecular Simy naturally os hydrate floto the total p (KI’s).(Chriynamic) inhit hydrate for fall in the r and/or groategory of K fish, insectsFPs are show activities. 4. Sample AFIII, Ocean pYellow mealwed; S, yellow. er flounder ( protein (wf-ulations occurring loww assurancerevention ostiansen & Sbitors such amation, the ange of 0.1-1wth processI’s are antifr, plants, fungn in Figure Ps. (A) Typout (Zoarceorm beetle . The ice-binPseudopleuronAFP) whichf Hydrate-concentrat are being mf hydrate forloan 1994; Ss methanol wconcentratio.0 wt.%. Thes instead ofeeze proteini and bacteri 2-4. Despitee I, Winter fls americanus(Tenebrio mding surfaceectes americanu enables it to Inhibitionion antifreezade to managmation) by inloan 2005; Khere concens required ey work by k shifting thes (AFPs). AFa, to help th strikingly diounder (Pse): op-AFP, (olitor): mb-As of these prs) can surviv prevent the e substancese the risk ofjecting chemelland 2006ntrations of for KI’s to binetically slo equilibrium Ps have evoem survive sfferent strucudopleuroneC) Snow fleaFP. Color soteins are ree arctic tem growth of ic, in recent y hydrate plugicals known) Unlike conup to 50 wt.e effective awing the iceconditions. lved in livinubzero temptures they al ctes america (Hypogastrcheme: H, wlatively flat aperatures bye crystals beears efforts i formation  as kinetic ventional % may be nere significant/hydrate  g organismseratures. A l have shownnus): wf-AFura nivicolahite; C, cyannd rich in  producing ayond a size 20 n the (as eded ly , few  ice P, ): sf-; N, n that would leastructuresFigure  2-Threonin(i)THR rmodels. CThis protsurface, IEach TH(–CH3). TAFP adsosurface is1992) Figure  2-are emphchain whhydrogenIn early mthe ice suand cowoTHR withd to the cell  of the alani5. The primae (THR, T) esidue are poolor schemeein has four BS) of the hR has a hydrhe chemicalrbs to the { believed to 6. Chemical asized by traereas ALA on, white; nitroodels, the hrface was berkers (Haym valine (VAdeath. The pne (ALA or ry and threeresidues are inted by arr: carbon, cythreonine (Telix and repeophilic hydr structures o20-21} ice plbe a key fact structures ofnsparent volly has one mgen: blue. ydrogen bonlieved to drivet et al. 199L) residues irimary and tA)-rich 37-am-dimensionacircled. The ows. Other aan; oxygen, rHR or T) reat every elevoxyl (–OH) f the THR aane where por in the bin (A) the THRumes aroundethyl side cding betweee the adsorp8) carried oun the same lohree-dimensino acid wf l structures o(i+4) and (i+mino acid reed; hydrogesidues whichen residues wside chain annd ALA resirotein surfacding and ant residue an them. THRhain. Color sn the OH grtion of wf-At mutation stcations in thional (Siche-AFP are shof wf-AFP. A7)ALA residsidues are sn, white; nitr lie on one sith an apprd a pendantdues are shoe/structuralifreeze actiod (B) the AL has one hycheme: carboup of THRFP to the icudies on wf-e protein seri & Yang 19wn in Figurlanine (ALAues with reshown with bogen, blue. ide (known oximate dist hydrophobiwn in Figure complemenn.(Knight etA residue. Tdroxyl and oon, cyan; ox residues ane surface. HoAFP where quence to re95) helix e  2-5.  , A): white bpect to the all-and-stickas the ice binance of 1.65 c methyl gro  2-6. The wftarity with th al. 1991; Chhe side chaine methyl siygen, red; d water oxygwever, Haythey replacedmove the 21 ars.  ding nm. up -e ice ou ns de en at met  22 presumed hydrogen bonding of the THR groups to the ice surface. The mutation substitutes the pendant –OH group of THR with a pendant –CH3 group in VAL. The retention of antifreeze activity despite this mutation disproved the idea that hydrogen bonding of –OH groups of THR residues are solely responsible for the anti-freeze action of wf-AFP. Surprisingly, mutation of the THR residues to serine (SER) which only retains the pendant –OH group leads to the loss of antifreeze activity. In later studies, the role of hydrophobic ALA residues was investigated and it was concluded that the contribution of hydrogen bonding to the ice surface is secondary to the van der Waals interactions of the AFP with the surface.(Baardsnes et al. 1999; Davies et al. 2002) AFPs evolved in living organisms to deal with ice formation but some of them inhibit formation of gas hydrates which have very different microstructures from ice. This is an adaptation of a natural product for a different purpose than it was designed for, which could lead to green inhibition practices. For some AFP inhibitors, it was shown that nucleation and growth inhibition on ice and gas hydrates are independent processes and that good performance for ice inhibition does not necessarily lead to good hydrate inhibition activity.(Ohno et al. 2010) The mechanism of the action of antifreeze molecules involves interactions that occur on the molecular scale. This level of detail is difficult to obtain with readily accessible spectroscopic or other experimental techniques. Molecular simulation techniques are well established and render useful molecular insights on nanometer scales which are otherwise challenging or not possible to capture experimentally. Molecular simulation studies of AFPs in water and water/ice systems have been performed. It is reported that the ice binding surface (IBS) pre-organizes the surrounding solution water molecules, therefore maximizing compatibility with the ice surface and minimizing the entropy penalty of binding.(Nutt & Smith 2008) These simulations showed that water molecules around the IBS of the spruce budworm AFP arrange in structures that are commensurate with a specific surface of the ice lattice and thus water molecules at the IBS bind to the grooves of the rough ice surface. On the other hand, water molecules in the vicinity of the non-IBS of the AFP take on arrangements that are incompatible with the ice lattice, and this further contributes to the antifreeze action. Protein 23 surfaces containing the THR-ALA-ALA sequence are also found to interact more favorably with the ice/water interface but no explanation is given about why this particular combination is favored.(Wierzbicki et al. 2007)  By investigating the interactions of sample AFPs with the hydrate surface and determining the structural/functional relationship for the inhibitory action it should become possible to rationally design a more effective inhibitor molecule which can potentially outperform the currently available kinetic inhibitors or modify them for improving their inhibition activity. Molecular dynamics simulations of the inhibition mechanism of the non-biological antifreeze substances on gas hydrates have also been studied by Rodger and coworkers (Storr et al. 2004) and Anderson et al.(Anderson et al. 2005) They studied inhibition mechanisms of polyvinyl pyrollidone (PVP) and polyvinyl caprolactam (PVCap) in contact with methane hydrate surfaces. A challenge in simulating these polymeric inhibitors is that their structure in solution is not known. Anderson et al. performed simulations of single vinylcaprolactam monomer molecules in the presence of a structure II (sII) clathrate hydrate with methane molecules in the cages (the usual methane hydrate structure is sI). They concluded that charge distribution at the edge of the inhibitor molecule which coincides with charge separation of water molecules at the surface of the hydrate led to strong binding between the inhibitor and the hydrate surface. However, it is not clear if this same mechanism will hold for the binding of these molecules on the empty cages of the structure I (sI) methane clathrate hydrate surface.  2.5 Thesis Objectives The following knowledge gaps are identified. • The phenomenon of hydrate formation and decomposition becomes more complex when the water and gas phases are in contact with a solid surface or porous environments as in their natural occurrences in sediment phases.(Linga et al. 2009; Seo et al. 2005; Bagherzadeh et al. 2011) The nature of the interactions between the water, gas-phase guest, and solid substrate phases must be understood in order to gain more insight about the fundamentals 24 of this process. Since these interactions mostly occur at the interface between the phases and usually extend over nanoscale lengths, the investigation of such systems using molecular simulation approaches is feasible and can give rise to useful insights on the process of natural gas hydrate formation and decomposition. Since sandy sediments have the highest saturation for hydrates, silica surfaces are proposed to be included in the simulated system to investigate their effect on the dissociation process of methane hydrate. An MD setup for the simulations of methane hydrate dissociation in the presence of the silica surface that is closer to the experimental conditions would be to have two large thermostated water reservoirs (both at the same target temperature and pressure) in contact with non-thermostated water reservoirs (adiabatic NVE) at a temperature Twater, which are then in contact with the non-thermostated (adiabatic NVE) hydrate phase at an initial equilibrium temperature of Thydrate< Twater. The heat conductance between the water and hydrate phase in this set-up would more accurately represent the experimental set-up for the hydrate melting. Performing the simulation in this manner would be feasible; however, the simpler NVE simulations are useful enough for illustrating the general features of the hydrate dissociation process.  • Understanding the fundamentals of methane transfer from the hydrate phase to the gas phase will be useful in optimizing the extraction process and maximizing the gas yield from hydrate dissociation. Gas reservoirs should be added to the simulated system to account for the transfer of released methane to the gas phase instead of artificially limiting the released methane to remain in the liquid phase. • Although, KI’s and specifically the antifreeze proteins are experimentally proven to effectively inhibit hydrate growth (Zeng et al. 2006; Daraboina et al. 2011; Sharifi et al. 2014; Perfeldt et al. 2014) their mechanism of action is still not well understood. The elucidation of the physical and chemical interactions between the inhibitors and the gas hydrate crystal is critical to design green and more effective inhibitors with potentially superior performance.  25 The objectives of the current research work are to advance the fundamental knowledge and give rise to the understanding of the hydrate decomposition and inhibition processes by employing molecular dynamics to study the mechanisms involved in these processes in more realistic simulated systems. More specifically: 1. To delineate the influence of silica surface on the decomposition behaviour of methane hydrate and the interaction between different phases (silica, water, hydrate, gas) 2. To assess the possibility and conditions required for the formation of methane nano-bubbles due to the decomposition of hydrates and their effect on the growth of hydrate 3. To determine the functional mechanism of binding of the winter flounder AFP to the methane hydrate surface  2.6 Thesis Organization The rest of this dissertation is organized as the following. Chapter 3 studies the influence of the silica surface on the water characteristics and the interaction among the gas hydrate former gases (CH4 and CO2), water phase and the silica surfaces. In Chapter 4, the decomposition of methane hydrate in the presence of silica surface is investigated. In Chapter 5, by adding the gas spaces to the simulated system, the possibility of the formation of methane nano-bubbles as a result of hydrate decomposition and their effect on the growth rate is examined. Chapter 6 discusses the mechanism with which the wf-AFP is bound to the hydrate surface. Finally, Chapter 7 summarizes the key findings and provides recommendations for future work. 26 Chapter 3: State of Water in the Presence of Silica 3.1 Introduction The formation mechanism of methane clathrate hydrate in natural deposits is important from both industrial and fundamental scientific points of view. No conclusive mechanism is as yet available for the formation of vast naturally occurring methane hydrate reserves. The industrial implications of understanding and controlling hydrate processes include the synthesis of methane hydrate from methane gas and water for storage and transportation purposes,(Englezos 1993; Gudmundsson 1996) extraction of methane from natural methane hydrate formations in permafrost and ocean sediments,(Collett 2002) and sequestration of carbon dioxide in depleted gas reservoirs.(Zatsepina & Pooladi-Darvish 2010)  As sandy reservoirs have the highest saturation of hydrates (Riedel et al. 2006) hydroxylated silica is added to the simulated system. A study of the interaction between water and the silica surface, especially when a separate reservoir of the gas phase guest molecules is present, can help in understanding the state of the system before the onset of hydrate crystal nucleation. Pure silica (taken to be the quartz phase) has the molecular formula SiO2, but adjacent to water surfaces, the outer layer of Si atoms is considered to be hydroxylated. Argyris and coworkers(Argyris et al. 2008) analyzed the structure and dynamics of interfacial water on silica with different densities of hydroxyl group on the surface. They found that up to about ~14 Å away from the surface, the local water has different structural and dynamical characteristics from bulk water. Experimental techniques including the IR spectroscopy,(Aines & Rossman 1984) differential scanning calorimetry (DSC),(Konrad 1990) thermoporometry,(Handa & Stupin 1992) DSC and FTIR (Ping et al. 2001) have been utilized to demonstrate the existence of a special non-freezable state of water of thickness from 3 to 15 Å, on different types of substrate such as quartz, clayey silt, silica gel and hydrophilic polymers.   To gain insight into the heterogeneous hydrate nucleation process in the presence of silica surfaces, the structure of the water meniscus between two hydroxylated silica surfaces (a slit-like pore), under 27 vacuum, methane, or carbon dioxide atmospheres is investigated. The structure of the silica/water/gas interfaces is studied and determined whether these sites have potential as nucleation sites for hydrate formation. 3.2 Computational Methods The molecular dynamics (MD) simulations were run with the DL_POLY (Smith & Forester 1996) program version 2.20. The intermolecular potential for atoms on different molecules is considered as a sum of Lennard-Jones and Coulombic interactions with force field parameters given in Table  3-1. A united-atom potential was used for methane,(Tse et al. 1983) the EPM potential (Harris & Yung 1995) for CO2, the TIP4P-ice potential (Abascal et al. 2005) for water, and the potential of Lopes et al. (Lopes et al. 2006) was employed to represent the silica atomic interactions. Standard combination rules were used to calculate the Lennard-Jones parameters for the unlike pairs of atoms, ߝ௜௝ = ඥߝ௜ߝ௝  and ߪ௜௝ = ൫ߪ௜ + ߪ௝൯ 2⁄  . These combinations of potentials have been used in simulations of methane and carbon dioxide clathrate hydrate nucleation and decomposition (Walsh et al. 2009; Bai et al. 2011) and are successful in capturing physical features of systems containing hydrates. Table  3-1. Parameters of the force field applied in this chapter Atom ε (kJ·mol-1) σ (Å) q (e) C (CH4) 1.3650 3.6400 0 C (CO2) 0.2411 2.7850 0.6645 O (CO2) 0.6901 3.0640 -0.33225O (H2O) 0.8822 3.1668 0 E (H2O) 0 0 -1.1794 H (H2O) 0 0 0.5897 Si (SiO2-bulk) 2.5104 3.9200 1.08 O (SiO2-bulk) 0.6364 3.1538 -0.53 O (SiO2-surface) 0.6364 3.1538 -0.64 H (SiO2-surface) 0.1925 0.4000 0.32  The cutoff radius for all of interatomic interactions in the simulation box was set to 15 Å and the long-range Coulombic interactions were treated using the Ewald summation method.(Ewald 1921) 28 All simulations were performed in the NVT ensemble with the thermostat constant of 0.2 ps and initial pressure of 1 bar. The equations of motion were integrated every 1 fs with the total time of 2 ns. Periodic boundary conditions were applied in all three dimensions. Thus replicated thin films of water between silica layers were simulated. The crystallographic unit cell of silica (quartz) (Levien et al. 1980) is used to construct the silica layers. Surface silica atoms potentially exposed to water were hydroxylated using an O-H bond length of 1.1 Å, according to set-up of Lopes et al.(Lopes et al. 2006) The details of the construction of silica surfaces are given in  Appendix A. The two silica slabs are placed in the simulation box in the xz-plane (38.6 Å apart along the y-axis). Variable numbers of water molecules (from 432 to 2400) were placed between the silica slabs. The maximum number of 1584 water molecules initially occupied about half of the simulation cell length along the z-direction. In simulations with methane and carbon dioxide, gas molecules were placed in the remaining empty space of the simulation box, surrounding the water phase from both sides along the z-direction. Table  3-2. The temperature, number of molecules in the water and gas phases and the calculated contact angles (degrees). The simulations all have a total run time of 2 ns. Simulation No. T(K) Number of molecules Contact angle H2O Gas 1 275 432 0 38.3 ± 2.8 2 1584 0 27.3 ± 3.6 3 1584 306 (CH4) 36.7 ± 4.5 3a 1920 375 (CH4) 36.8 ± 7.3 3b 2400 461 (CH4) 39.8 ± 1.2 4 1584 612 (CH4) 51.1 ± 7.2 5 1584 306 (CO2) 39.0 ± 7.0 6 1584 459 (CO2) 52.4 ± 5.7 7 300 864 0 32.6 ± 4.2 8 1584 0 22.5 ± 4.2 9 1584 306 (CH4) 28.5 ± 6.0 10 1584 918 (CH4) 34.1 ± 4.1 11 1584 306 (CO2) 29.9 ± 1.0 12 1584 459 (CO2) 31.9 ± 5.2  Two sets in contacy-directionumbers simulatioof the simphase areThe atomconfigurashape of Figure  3-moleculemeniscusregions dto calcula3.3 ReSnapshot300 K areuniform gThe concsurfaces isurface inof simulatiot with the wan) to study thof moleculesns are given ulation box considered s of the siliction of the sthe meniscus1. (a) Initial s between th illustrates thesignated aste average gsults and Ds of the simu shown in Fas distributientration ens striking and contact withns were perfter meniscue effect of s in each simin Table  3-2 at 275 K is afor studying a slabs are kimulation ce shows that set-up and (be silica slabse hydrophil G (in gas phas concentraiscussionlation with tigure  3-2. Thon, are showhancement o is quantifie methane anormed at 275s are performlab separatioulation and t. The ideal gapproximatelthe effect ofept frozen dull at 275 K wwater has we) final confi at the top anic character ase), I (at thtion in the si he water mee initial set-n in  Appendf methane ad below. Fod carbon di and 300 Ked with silicn on the mehe numberins pressure oy 20 MPa. L pressure. ring the simith 1584 wattted the silicguration for d bottom ofof the silica se water-gas mulation boniscus in conup configuraix B. nd carbon dir reference, aoxide gas are. Two additioa slab separaniscus and wg scheme uf 306 methaarger numbeulations. Ther moleculea surface indsimulation ( the simulaturfaces usedinterface) anx. tact with CHtions for theoxide on the set of simu performed nal simulatitions of 50 ater propersed to represne moleculers of gas moe initial setus is shown inicating that  2) at 275 K wion cell. The in this workd C (in wate4 and CO2 se simulatio water menilations of a funder similaons for methand 62 Å (inties. The ent the diffes in the gas plecules in thp and final  Figure  3-1. it is hydrophith 1584 wat formation o. Three distr phase) are gas at 275 anns, which hascus and siliclat bulk water simulation29 ane  the rent hase e gas The ilic. er f the inct used d ve a r  conditionB.  Figure  3-K  (simulillustratedgas near tindicated3.3.1 WThe strucsilica surfrelative nFigure  3-thick layeadjacent tcenter of close to es. A sample 2. The final (ations 9 and  as cyan sphhe interface  in Table  3-2ater and Gture of the waces) is deterumber of wa3 along with rs of water mo the silica sthe simulatioach silica sursnapshot of t = 2 ns) sna12) in the preres and COis evident. T. The asterisas Distribuater phase amined by cater moleculea typical snaolecules. Thurface (L1), an cell far froface (roughlthe water/gapshot of theesence of  CH2 as linear cyhe numbersk marks a potion long the y-dilculating thes as we movpshot from ese layers, s water layer m the silicay within L1 as distributio simulations 4 and CO2 gan-red mole in parenthestential nuclrection of th water molece perpendicuthe simulatiohown in Fignext to the f surface (C). nd L2) shown in these sy at T = 275 (aseous phacules. The ees corresponeation site foe simulationule distributlar to the siln. Analyses ure  3-3, reprirst layer (L2The peaks o the formatistems is givesimulations 3ses. CH4 monhanced cond to simular hydrate cry cell (perpenion. The proica surfaces are performeesent the bo) and the waf enhanced won of two din in  Append and 5) and lecules are centration otion numberstals. dicular to thfile of the is shown in d on three 3und water ter layer neaater populastinct hydrat30 ix 300 f the s e  Å r the tion ion layers neacenter wiFigure  3-of enhancis shown.Three-dimsimulatiowater mothe watersurface cosurface mlocationsr the surfaceth little indic3. Water disted water den ensional den cell are callecules at the and silica inrrespond toake strong h. . Beyond thation of furtribution alonsity are presnsity plots foculated and s silica surfacterfaces are c the locationydrogen bonese layers, thher water layg the distanent near ther the concenamples are se and the enlearly visibles of surface ds with these number ofer structurin ce between t surface. Fortration of thhown in Fighanced conc in these figO-H groupse O-H grou water molecg. he two silica clarity, onlye species inure  3-4. Theentration ofures. The hig of silica. Waps and are pules decreas surfaces. Tw half of the s different pa enhanced c methane gah density loter moleculreferentially es towards to major peymmetric prrts of the oncentrations molecules ci near the sies near the silocated at th31 he aks ofile  of at lica lica ese Figure  3-275 K (simgas interfIn Figureenhancemexperimewith liquithe adsorwas obserLehmkühformationformationIn Figuresimulatio4. Concentraulation #3)aces can be   3-4(a), it is sent is even ntally observd water undeption of gueved for gas ler et al.(Leh of a thin fil.    3-1, three bn box dimention surfaces. The enhanobserved. Ween that megreater on seed that for ur high pressst moleculesphase Xe. Hmku ̈hler et am of liquid-loxes with 7×sion in the x for (a) methced concentrater density ithane concengments of thnstirred liquures, spontaon the wateowever, X-Rl. 2008) for ike CO2 on 7 Å2 cross se-direction ar ane and (b)ations of mes also enhantration is ene silica surfaid CO2, propneous macror-guest interfay diffractiogaseous COthe water intctions in the selected. T water densitthane at theced at the sihanced at thce not coveane, and isoscopic hydraace.(Boewern and reflect2-liquid wateerface but the yz-plane anhese distincty projected  silica surfaclica surface. e water surfred with watbutane phaste formation et al. 2012) ivity measurr system shoe absence od a depth eq regions are on the yz-plae-gas and waace, but er. It is es in contact is triggeredSimilar behaements of wed the f hydrate ual to the shown as G32 ne at ter-  by vior  (in 33 gas phase), I (at the water-gas interface) and C (in water phase). The 7 Å distance is twice the distance of the first minimum of OW-OW radial distribution function. The average number of water and gas molecules in the boxes during the final nanosecond of the simulation trajectories were counted and used to obtain the density of gases in each box. The ratio (I/G) is given in Table  3-3 for the different simulations. A notable density enhancement (I/G > 1.0) of methane and carbon dioxide at the site of the meniscus is notable. Carbon dioxide accumulation is enhanced on the surface to a greater degree than methane. Table  3-3. The density of CH4 and CO2 molecules at the interface and gas phase of the simulation cell for different simulation set-ups. All values are reported with three significant digits. The standard deviation of the mean for gas density is approximately in the order of 10-7 mol.cm-3 for all simulations. A sample value is given for simulation 3 in brackets. Simulation (gas) Gas density (mol.cm-3) G (gas phase) I (water-gas interface) Ratio (I/G) 3 (CH4) 2.49E-03 [2.16E-07] 5.57E-03 [2.82E-07] 2.24 4 (CH4) 9.46E-03 1.35E-02  1.43 5 (CO2) 4.75E-04 6.55E-03  13.8 6 (CO2) 8.69E-04 1.08E-2 12.4 9 (CH4) 3.22E-03 4.73E-03  1.47 10 (CH4) 2.00E-02 1.08E-02 0.542 11 (CO2) 9.36E-04 6.00E-03 6.41 12 (CO2) 1.96E-03 8.88E-03 4.53  For reference, the density profile of bulk water at 275 K in contact with methane/carbon dioxide phases through a flat interface are given in Figure  3-5. A similar density enhancement of the gas molecules at the water interface is seen in these simulations as well. The number of water and gas molecules is similar to that in simulations 3 and 5. 34  Figure  3-5. Distribution of hydrate forming gases adjacent to a flat water interface along the length of simulation box (no silica surface) for (a) the methane – water system and (b) the carbon dioxide – water system. The behavior of the methane-water system is similar to that of Figure  3-2. 3.3.2 Contact Angle To characterize the distribution of water molecules, the space between the silica surfaces is binned into volume elements of 1 Å length in yz-plane and a depth equal to the length of the simulation box along the x-axis (i.e., 39.328 Å). The silica-water and water-gas (water-vacuum) boundaries are found by determining the average water densities in the bins over the trajectory of the simulation during the final nanosecond where the shape of the meniscus has stabilized. The distribution of water molecules in the form of a meniscus is shown in Figure  3-6. The outer boundary of the water phase is not sharp and the density of water decreases gradually from ~1 g/cm3 in the water phase to 0 as we move away from the water phase. To calculate the contact angle, α, between the water and hydroxylated silica phase, the outer boundary of water phase is chosen as sites where the density has the mean value of 0.5 g/cm3. In all cases the section of the arc of the meniscus near the silica surface (~2.5 Å perpendicular distance) was fit to a straight tangent line, as illustrated in Figure  3-6, which yielded an r-squared value of greater than 0.98 in all cases. The average of the tangents at the four points of contact of the meniscus with the silica surface is used to report the contact angle. This method is similar to the method of Shi.(Shi 2006)  The average values of the contact angle between the water and silica phases in different simulations are given in Table  3-2. Snapshots of the final configurations of the water/CH4 and water/CO2 systems at two temperatures are shown in Figure  3-2. Comparing the snapshots shows that at the lower temperature, the curvature of the meniscus decreases and water spreads less on the (a) (b) 35 hydrophilic surface and therefore the contact angle increases. This is consistent with the fact that the surface tension of SPC/E water increases from 61.3 to 64.5 mN/m as temperature decreases from 300 to 275K.(Chen & Smith 2007) The presence of methane and CO2 molecules on the gas phase affects the contact angle. For each gas, greater gas pressure leads to larger water-silica contact angles. The CO2 gas with its higher solubility leads to smaller contact angles than similar simulation with CH4 gas.    Figure  3-6. A typical distribution of water molecules between the silica surfaces obtained by binning of the areas in the yz-plane of the simulation. The higher local concentration of water molecules near the silica surfaces is apparent. The outer boundary of the water meniscus is used to estimate the contact angle. Four tangent lines were drawn on four points of contact of water to the silica surfaces and the contact angles are reported in Table  3-2. A schematic of the method for drawing tangent lines on the meniscus to determine the contact angle is shown in the bottom panel. An arc of approximate length of 2.5 Å is used to draw the tangent to the meniscus near the silica surface. It should be noted that the contact angle is a property of a bulk sample of water. In simulations where the number of water molecules is small (simulations #1 and #7), the amount of water in the 36 simulation is not sufficient to observe the bulk spread of water on the surface. As the number of water molecules in the simulation decreases further, the layer connecting the silica surfaces becomes thinner and ultimately tears. The contact angles of the water phase with hydroxylated silica in the vacuum and in the presence of gas are related to the surface tensions between various phases through Young’s equations for the force balance at the meniscus. For the direction perpendicular to the silica surface, ஺݂ = ௪݂௚ sin ߙ Equation  3-1where fA is the adhesive force of the liquid water to the silica surface and fwg is the tension force between the water and gas. For the direction parallel to the silica surface,  ௦݂௪ + ௪݂௚ cos ߙ = ௦݂௚ Equation  3-2where fsw and fsg are the tension forces between the solid silica-water and solid silica-gas phases, respectively. The latter equation is usually written in terms of the surface tensions, ߛ௦௪ − ߛ௦௚ = −ߛ௪௚ cosߙ Equation  3-3where the subscripts on the surface tension terms is similar to those of the tension forces. The silica-gas and water-gas interfacial tensions depend on the pressure of the gas in the simulation. In all simulations we observe that the contact angle α < 90° which implies that γwg > 0 > γsw - γsg. The contact angle increases as gas molecules are added to the water-silica systems and this angle increases further at higher gas pressures as seen in Table  3-2. As the gas pressure in the simulation increases, the interfacial tension force between the water and gas phases, fwg, decreases (Nielsen et al. 2012) requiring an increase in the contact angle α to retain the force balance condition in Equation  3-1.  The ranges of contact angles at 300 K for the pure water–silica simulations (between 29 and 34°) are in good agreement with the experimental value of 32.79 ± 1.12° reported for distilled water on quartz at 23°C.(Sklodowska & Matlakowska 1997) 3.3.3 WTo undersurface, tand the nThe dipoperpendicFrom the0.2) are o0.5). The extends awith the ffrom therpore radiifrom theiwhere theFigure  3-surface. Tnormal tomoleculelines are tater Orienstand the orihe average oormal vectorle moment vular (i.e., <c panel (b) ofriented diffeconclusion it least 6 Å ininding of Hamoporomet in the ranger 1H NMR ry reported t7. (a) Schemhe blue dott the silica sus and normahe L2 layer vtation entational chf the cosine  of the silicaector of watosφ> = 0) to this figure, irently than ths that the orito the waternda et al.(Hric analysis o 23-70 Å anelaxation timhe thicknessatic of the ored arrow shorface, n. (b) l vector to thalues and tharacteristicsof the angle, surface, n iner molecules the silica sut is also evide water moentational ef phase. This anda & Stupf ice formatid with the ste experimen of the affectientation of ws the directhe average e silica surfae blue lines  of water mo φ, between t Figure  3-7( in the L1 layrface. This ient that the lecules in thefect of the sthickness foin 1992) whon in porouudy of Totlatal results ined water laywater molecution of dipolangle φ betwce. The blacare the C laylecules in thhe dipole ma), is calculater appears tos schematicawater molec central, bulilica surface r the bound o reported a s Vycor glasnd and cowo silica sampler to be at leles in the lae moment veen the polak lines are ther values. e layers relatoment of waed and plott be approxilly shown inules in layer k-like C layeon the waterlayer is in govalue betwes and variourkers (Totlaes (average dast 8 Å.  yer adjacentector M and r moment vee L1 layer vaive to the silter moleculeed in panel (mately  Figure  3-7. L2 (<cosφ>r (<cosφ>  ~ molecules od agreemenen 4.5 to 5.5s silica gels wnd et al. 201iameter of 4 to the silica the vector ctor of watelues, the red37 ica s b).  ~  t  Å ith 1) 0 Å) r  38 3.3.4 Residence Correlation Function A residence correlation function, CR,(Argyris et al. 2008) is defined to study the exchange and dynamics of water molecules in each layer, ܥோ(ݐ) =〈 ௪ܱ(ݐ). ௪ܱ(0)〉〈 ௪ܱ(0). ௪ܱ(0)〉 Equation  3-4Here ௪ܱ(ݐ) is a variable assigned values of 0 if the water molecule oxygen atom has moved outside the layer it occupied at time t = 0, and 1 if the water molecule oxygen atom has remained inside its corresponding t = 0 layer. The value of CR decays from 1 to 0 as water molecules move between layers during the course of simulation. The residence correlation function is defined to investigate the dynamics and displacement of water molecules in the layer adjacent to the silica surface (L1), next water layer away from the silica surface (L2) and the water layer in the center of the meniscus (C) and is plotted for different simulations in Figure  3-8. The water molecules in layer L1 are noticeably less mobile than those in layers L2 and C. The water molecules in the L2 are also less mobile than those of the C layer. As the temperature increases from 275 K (top row) to 300 K (bottom row), water molecules diffuse more quickly between horizontal layers and CR decays faster. Water is also spread out over a larger cross section of the surface at the higher temperature which could influence its dynamics. In the presence of gases in contact with the silica surface and water meniscus, the mobility of the water molecules in L1 increases (compare the first column of Figure  3-8 with the second and third columns.). This suggests that the presence of gas molecules in the gas space affects the dynamic behavior of water molecules in layer L1. This could be due to the dissolution of gas molecules in water phase as well as gas adsorption on the silica surface.  Figure  3-numbers L1 water lshown in blue. Topwith 306 mpanel. 3.3.5 HAs the wathe silica the adjaceChandler If two wof less thbonding ois countedstep. 8. The residein parenthesayer adjacenred, and the row is at 27ethane molydrogen Bter moleculesurface and ent water mo 1996) is defiܥு(ݐ) =ater moleculan 30º) at timf each parti and the avence correlates correspont to silica su correlations5 K and bottecules in theond Analyss in the hydspecially thelecules. Thened to inspe〈ݐ)ܪ). 0)ܪ)〉〈0)ܪ). 0)ܪ)〉 es are hydroge t, H(t) = 1cular pair is crage value ision parameted to the simrface are sho for the wateom row at 30 middle panis rate structur surface O-H hydrogen boct the rate oen-bonded  for that paiontinued ov calculated fr for the watulation numwn in black,r layer in the0 K. The gasel and 306 cae are hydrog groups affnd correlatif hydrogen-b(O-O distanr, otherwise er the trajecor different  er layers in dbers in Table the correlati center of th phase spacrbon dioxiden bonded, iect the hydroon functionond formatce of less thaH(t) = 0. Thtory and thewater layers ifferent sim  3-2. The coons for L2 we meniscus Ce is empty foe molecules t is importangen bondin(HBCF), CHion and decan 3.5 Å ande checking f number of h(L1, L2 and Culations. Thrrelations foater layer are are shown r left panel, in the right t to study hog dynamics o(t),(Luzar & y, Equatio O-H…O anor the hydroydrogen bo) at each tim39 e r the  in filled w f n  3-5gle gen nds e 40 Water molecules in L1 form hydrogen bonds with other water molecules as well as the O-H groups of the hydroxylated silica surface. These two contributions are separately computed and are labelled as CWW and CSW, respectively. Figure  3-9(a) shows the decay of HBCF for simulation 4. The surface-water (SW) hydrogen bonding in L1 has the slowest decay. This indicates that hydrogen bond formed between the silica surface and water molecules in L1 is the most stable. The water-water (WW) HBCF in L1 shows slower decay when compared to waters in L2 which shows that the silica surface indirectly affects the WW hydrogen bonding dynamics in L1. It also appears that this influence is extended to L2, although to a lesser extent. Figure  3-9(b) compares two similar simulations (5 and 11) at different temperatures. As expected, the simulation at higher temperature shows faster decay for HBCF as water molecules have greater kinetic energy and hydrogen bonds break more quickly.  Figure  3-9. The hydrogen bond correlation functions for water molecules in the L1, L2, and C layers (each 3 Å thick). (a) hydrogen bonding dynamics in layers near the silica surface (L1 and L2) and the central layer (C) for simulation 4. The hydrogen bonding in L1 is separated into two contributions: hydrogen bonding between water molecules (WW) and hydrogen bonding between water molecules and OH of the silica surface (SW). (b) hydrogen bonding dynamics in layers L1 and L2 for simulations at 275 K (5) and 300 K (11). 3.3.6 Effect of the Separation between the Silica Surfaces To assess whether the separation distance between the silica surfaces affects system properties, two simulations (3a and 3b) with separations of 50 and 62 Å are performed (see Figure  3-10). These simulations were prepared such that their water/methane ratio is roughly the same as that in simulation 3 (separation distance of 38 Å). As reported in Table  3-2, the contact angles in these (a) (b)41 simulations are identical within the error margin of the measurements. Additionally, the CR(t) orientation parameter, residence time and hydrogen bonding show similar trends for L1, L2 and C which is assigned to be in the middle of the simulation box to represent the bulk water phase. This indicates that the effect of silica surface extends to approximately 6 Å away from the surface and water molecules beyond this distance behave similar to bulk water. New phenomena related to confinement can happen in silica pores with diameter less than 12 Å where bulk-like water would not be present. Figure  3-and meth38 Å (bot 10. Snapshotane gas phatom). The co from Simulase for systemntact angle  tions 3, 3a, as with silica remains relatnd 3b of Talayer separaively the samble  3-2 after tions of y = 6e in these t2 ns showing2 Å (top), 50hree cases.  the menisc Å (middle)42 us , and 43 3.3.7 Discussion The simulations show that the mobility of water molecules adjacent to the silica surface is significantly decreased and as we move away from the surface, the behavior of the water molecules becomes more bulk-like. Methane and carbon dioxide molecules bind extensively to the liquid interface and silica surfaces. These results suggest a heterogeneous nucleation mechanism for hydrate formation in porous media. As relatively high guest/water molecule ratios are needed to form hydrates, the interfacial regions can deliver greater concentrations than delivered by gas solubility in bulk water alone. The interface is an appropriate location for hydrate nucleation due to the enhanced relative concentration of gaseous guest molecules. The mobilities of water molecules in layers up to 6 Å away from the solid silica interface are also smaller than in the bulk liquid phase and this can further enhance the nucleation of hydrate phase. The entropy penalty for forming the solid hydrate from the less mobilized water near the surface is smaller than forming hydrate from bulk water molecules which are fully mobile in the liquid phase. A potential nucleation site is marked by an asterisk in Figure  3-2. From the results of Table  3-3, it is clear that the local gas concentration at the interface (I) is considerably larger than the local gas concentration inside the gas phase (G). Therefore, as expected, the interface is more prone to nucleation of hydrates as compared to the bulk phase. As the number of gas molecules increases, further surface concentration enhancement becomes less pronounced and eventually disappears. The interface has a finite capacity to adsorb gas molecules and after reaching this limit additional gas molecules remain in the gas phase. Moreover, the concentration enhancement decreases at higher temperature which shows that adsorption capacity of the interface becomes less with increasing temperature. 3.4 Conclusions The distribution and dynamics of water molecules in a slit-like pore made up of fully hydroxylated silica surfaces in the absence and presence of hydrate forming gas molecules are studied. Water molecules initially placed between the silica surfaces are allowed to evolve in an NVT ensemble simulation for 2 ns and formed a meniscus (curved interface), confirming that hydroxylated silica is 44 hydrophilic. The water – silica contact angle was determined by drawing tangents to the meniscus and the results were in good agreement with the experimental values. The water molecules at the silica surface are partially immobilized and the strength of the hydrogen bonding network is different from that of bulk water. Water molecules are preferentially populated near the O-H groups of silica surface which indicates formation of hydrogen bond between interfacial water molecules and the fully hydroxylated silica surfaces. Methane and carbon dioxide are preferentially adsorbed on the water interface at the meniscus and at the surface of the exposed silica not wetted by water. In accordance with experimental X-ray reflectivity measurements of Lehmkühler et al., it was observed that the hydrate forming gases accumulate adjacent to the interface and the calculated gas concentration in interfacial regions was two orders of magnitude greater than that attained solely through solubility. It is proposed that a possible nucleation site for heterogeneous nucleation of hydrates in water wetting porous media is at the curved interfacial region between water phase and the hydrate forming gas, where abundant amount of water and hydrate former molecules are available, and a few Ǻngstroms away from the solid surface immediately after the bound water layer (e.g. L2 in our simulation) where the mobility of water molecules is slightly less than the bulk phase.  45 Chapter 4: Methane Hydrate Dissociation in the Presence of Silica 4.1 Introduction Understanding the details of methane hydrate dissociation, especially under conditions similar to those of the hydrate sediments found on continental margins, under the ocean floor, and under the permafrost, i.e. in porous media and under restricted geometry, can be beneficial in designing optimal and controlled methane recovery processes from natural methane hydrate reservoirs.(MacDonald 1990; Milkov et al. 2003; Hong & Pooladi-Darvish 2005) Since it is known that sandy environments contain highest hydrate saturation,(Riedel et al. 2006) silica is incorporated in the simulations.  In this chapter, the dissociation of methane hydrate in contact with a hydroxylated silica surface with NVE MD simulations starting at different initial temperatures, above the equilibrium temperature, is studied. Three sets of simulations are performed. In the first set, methane hydrate decomposition in contact with water is simulated as the reference system. In the next two series, simulations of methane hydrate in contact with water on two sides (in the z-direction) and having interfacial contact with hydroxylated silica on two other surfaces (on the xz-plane) are performed. One set of simulations has the solid methane hydrate phase in direct contact with the hydroxylated silica surfaces and the other set of simulations have a relatively thin layer of water (3 Å) acting as a buffer between the hydroxylated silica and methane hydrate surface to relieve the structural mismatch between these two solid lattices. The effect of the silica surface on the decomposition of methane hydrate is examined as a first step in simulating hydrate dissociation under conditions closer to the real conditions. By performing NVE simulations, the endothermic nature of the hydrate dissociation is accounted for and the formation of temperature gradients between the decomposing hydrate surface and the bulk water phase is allowed. As discussed previously, NPT or NVT simulations which “equilibrate” the temperature between the hydrate phase and the aqueous phase at regular intervals in the trajectory through the operation of the thermostat do not correctly simulate the hydrate decomposition process. The NVE simulations are still one step away from reality since the initial temperature of the hydrate in the simulation is still set to be that of the water and still not 46 equal to the equilibrium temperature of the hydrate on the phase diagram. The present simulations do however demonstrate the establishment of temperature gradients and are a step towards greater realism in modeling the hydrate decomposition and its rate of occurrence. The goal is to gain a better phenomenological and semi-quantitative understanding of the hydrate decomposition process starting with molecular scale concepts. 4.2 Computational Methods A similar computational method as Chapter 3 is applied. As a reference system, a 3×3×6 unit cell replica of the cubic sI hydrate (a = 11.94 Å) was prepared with methane molecules placed at the center of the small and large cages. A row of water and methane molecules were eliminated from one end of the hydrate phase in the z-direction to render the hydrate phase symmetric along this direction of water contact. A total of 2322 water and 357 methane molecules were present in the initial methane hydrate phase. Two reservoirs of water (a total of 5406 water molecules) were added to both sides of the hydrate phase along the z-direction such that the lengths of these water reservoirs on each side were approximately half the size of the hydrate phase. This arrangement is called the hydrate-water (HW) scenario. The initial configuration of the system is shown in Figure  4-1(a). To equilibrate the system and relax the water molecules at the hydrate/water interface, 400 ps constant volume-constant temperature (NVT) runs were performed followed by 400 ps NPT simulations for pressure relaxation. Thermostat and barostat constants of 0.2 and 0.5 ps were used in these simulations. The hydrate phase was allowed to dissociate in a final 2 ns NVE simulation. Using these protocols three sets of simulations were conducted with initial temperatures of 283, 293 and 303 K.  Figure  4-scenario; and 6 andmethane z-directiobox alongFor simulscenariosin Figure simulatio2670 watIn the secthe silica between tinitial liqu“bound wchange.(H1. The initial(b) SHW sce 8 slabs perphydrate dissn in HW sce the z- and xations of me are consider 4-1(b), the mns) is placed er moleculesond case, calayers are mohe methane id phase in tater”, whichanda 1986;  configurationario and (cendicular toociation are nario and al-direction arthane hydraed. In the firethane hydrin between t are placed olled the silicaved further hydrate and his set-up. I sticks to poClennell et ans for simul) SWHW sce the y-directalso shown. Fong y-directie Lz = 145 Åte decomposst scenario, ate phase (owo hydroxylutside the m-water-hydrout along ththe silica surt is known thre walls andl. 1999) Usin ation of metnario. The 1ion used to cor clarity, son in SWHW and Lx = 36ition in the pcalled the silf the same siated silica suethane hydraate-water (SWe y-directionfaces. Thereat water in p does not ung thermopohane hydrate0 slabs perpharacterize ampling laye scenario. T Å, respectivresence of tica-hydrate-wze and strucrfaces separte phase, inHW) scena and a ~3 Å are total of orous medidergo a meltrometric ana dissociationendicular to the system brs are only she length ofely. he silica surater (SHW)ture as used ated by 38.6 between therio, shown i layer of wat3413 water ma forms a thiing/freezinglysis, Handa; (a) HW the z-directiehavior durihown along  the simulatface, two  scenario, shin the HW  Å. A total o silica surfacn Figure  4-1er is placed iolecules in n layer, calle phase  et al. (Hand47 on ng the ion own f es. (c), n the d the a 48 1986) obtained a thickness of 4.5 to 5.5 Å for the bound water layer in porous silica gels with pore radii range of 23 to 70 Å. The SWHW scenario is meant to mimic these conditions. Similar equilibration stages and production run as in the HW scenario are performed. The silica phase is frozen throughout these simulations.  The F3 order parameter of Baez and Clancy (Baez & Clancy 1994) is used to characterize the local geometrical arrangement of the water molecules during the hydrate decomposition simulation. F3 measures the deviation of the water molecule network from the near tetrahedral arrangement observed in hydrate phase: ܨଷ,௜ = 	 〈ൣcos ߠ௝௜௞ หcos ߠ௝௜௞ห + cosଶ 109.47൧ଶ〉௝,௞ = ቄ~0.0 ℎ0.1~݁ݐܽݎ݀ݕ ݎ݁ݐܽݓ Equation  4-1Here jikθ represents the angle formed by a triplet of water oxygen atoms where i is in the center and j and k indices correspond to neighboring water oxygen atoms which lie within a spherical shell of radius 0.35 nm around the central oxygen i.  To analyze the structure of the phases, the simulation cell is subdivided into 10 layers in the z-direction, each ~12 Å in thickness, as shown in Figure  4-1(a). The order parameter can be separately averaged over all water molecules in the simulation cell, designated as F3, over all water molecules initially in the hydrate phase, F3,hydrate, or only over water molecules within the ith layer perpendicular to the z-axis, F3,Li. The decomposition of methane hydrate is associated with the increase of the F3 value from ~0.0 to ~0.1. The variation of the F3 order parameters is also monitored in eight layers in the y-direction to study the methane hydrate decomposition parallel to the silica surface, see Figure  4-1(c). The potentials used in the simulations in Chapter 3 and 4 have not been calibrated against the methane hydrate dissociation curve and so the temperature/pressure values of the simulations may be shifted with respect to the experimental phase diagram. Although the absolute values of the pressure and temperature of the phase transformations may differ from the experimental values, the trends are quantitatively reasonable. 4.3 ReThe confof the NVdecompoFigure  4-scenario. bonds, anThe evolushown inperpendicand L8, grphases, hliquid watF3,Li incresimulatioFor the Tappreciaboutmost Lsults and Digurations ofE simulatiosed to an ap2. The configMethane mod, red dottedtion of the l Figure  4-3 bular to the zeen and violave order paer (F3 ≈ 0.09ases slightly n, the outmoinit = 293 K le extent. Th3 and L8 layiscussion the water-hns for the Hpreciable exturations at tlecules, liqu hydrogen bocal structury calculating-direction. Iet lines) whirameters, F3 5) values. Twith time. Inst L3 and L8 simulation, Le stepwise ders at the me ydrate systemW scenario aent in the simhe end of thid water andonds, respeces of the ph the local wan Figure  4-3ch contain w≈ 0.04, betwhere is a sma the HW scelayers of the3 and L8 are issociation othane hydras with threere given in Fulations wi e 2 ns runs w hydrate are tively. ases in the Hter order pa, the layers atater moleculeen that of tll amount onario for Tin methane hythe only twof the hydratte/water inte different inigure  4-2. Tth Tinit = 283ith Tinit = 28representedW simulatiorameters F3,L the methanes from bothe solid metf hydrate disit = 293 anddrate phase  layers that he phase at Tirface dissocitial temperahe hydrate h K. 3, 293, and 3 by cyan sphn with Tinit =i for layers ine hydrate/wh hydrate phhane hydratsociation at t 303 K, by thhave dissociave dissocianit = 303 K iiate first. Aftures at the eas not 03 K in the Heres, red-wh 283 K are  the simulatater interfacease and liquie (F3 ≈ 0.02)hese layers se end of theated (F3 ≈ 0.ted to an s apparent. Tter these laye49 nd W ite ion  (L3 d  and ince  2 ns 1). he rs have subslayers (L4Figure  4-different lrow), 293dissociatiand 8), shFor the sindicate tThe chanand 303 Kthe SWHadjacent tthe decomrelatively tantially diss and L7, blue3. The F3,Li oayers perpen(middle rowon of the meow F3,Li valuet of potentihat the methges in the lo are also givW and SHWo the silica sposition offlat, a phenoociated, as c and black, rrder paramedicular to th), and 303 Kthane hydraes of 0.04 whal energy funane hydrate/cal structure en in Figure scenarios aturface. It is c rows of hydmenon repoharacterizedespectively) ter for the we z-direction (bottom rowte phase occich are betwctions chosewater systemof phases fo  4-3. The hig t = 0 are dulear that decrate cages parted in other by F3,Li ≈ 0.begins. ater and hyd (see Figure). In the HWurs and the ieen the hydn for this wo equilibriumr the SHW aher starting e to the quicomposition rallel to the  studies as w08, the dissorate phases a  4-1) for sim scenario atncomplete crate-like andrk, the F3 v temperaturnd SWHWvalues of thk decomposunder HW cinterface andell.(English ciation of th s a functionulations with 283 K, no laages at the i liquid-like salues for thee is near 283scenarios foe F3 parametition of the onditions p the dissociet al. 2005; Ae next inner  of time for  Tinit = 283 (rge scale nterface (slabtates of wate different lay K.  r Tinit = 283, er (F3 ≈ 0.05hydrate phasroceeds throation front islavi & 50 top s 3 r. ers 293 ) for e ugh  51 Ripmeester 2010) The F3 curves in Figure  4-3 (particularly for the SWHW scenario) show that the dissociation of methane hydrate in the presence of silica interfaces still proceeds in a stepwise manner with outer layers decomposing before inner layers. However, due to the effects of the additional boundaries, the methane hydrate also dissociates from the sides in contact with the silica surface and therefore decomposition of each layer perpendicular to the z-direction does not proceed as clearly defined as in the HW case. The corresponding plots of the variation of the F3 order parameter in the slabs along the y-direction showed a similar trend and are given in  Appendix C. During methane hydrate decomposition in an NVE ensemble, the kinetic energy is expended to overcome the potential energy holding the hydrate structure together and therefore the temperature of the system decreases. The time variation of the temperature profiles are shown in Figure  4-4. Dissociation of hydrates is also associated with an increase in pressure which is a result of the release of methane gas (see Figure  4-5). Due to the greater extent of methane hydrate decomposition, the temperature drop and pressure increase for the SWHW and SHW simulations are higher than those of the HW scenario. Temperature drops as large as ~20 K are observed in the SWHW scenario at 303 K. For the HW scenario at 303 K, as a result of the endothermic hydrate decomposition, the temperature drops to ~290 K at 1.3 ns and after this time hydrate decomposition effectively ceases. Since the simulations are adiabatic, as a consequence of the decrease in temperature, the rate of decomposition of the remaining methane hydrate phase and the rate of further temperature drop decreases. This is known to be a problem in naturally occurring methane hydrate reservoirs where the depressurization technique is employed to extract methane from the hydrate deposits without supplying the heat required for the endothermic reaction of hydrate decomposition. Some of the temperature plots in Figure  4-4 do not have clear temperature plateaus. This shows that the corresponding simulations have not reached equilibrium. The goal in this work is not to continue these non-equilibrium simulations until they reached a final equilibrium state, but rather to use the evolving temperature as a measure of the rate of hydrate dissociation. The simulations would 52 all either reach equilibrium at the hydrate melting point between hydrate and water phases or proceed until all the hydrate is dissociated and the system consists of water and methane gas.  Figure  4-4. Temperature profiles for the simulations at 283 K (top), 293 K (middle) and 303 K (bottom) in the NVE ensemble for three different scenarios introduced in the text, HW (left), SHW (middle) and SWHW (right). The drop in temperature is an indication of hydrate dissociation. There is an associated increase in pressure with hydrate dissociation as shown in Figure  4-5. Due to keeping the silica surfaces frozen during the simulations, the absolute values of the calculated pressure in the SHW and SWHW simulations may not be accurate. Nevertheless, changes in pressure as the hydrate dissociates and methane is released in the simulation cell are accurately characterized. Since the methane hydrate phase decomposes to a larger extent, the pressure increase for the simulations in the presence of silica is larger than those of the HW scenario.  53  Figure  4-5. Pressure profiles for the decomposition of methane hydrate. The first, the second and the third rows are at 283, 293 and 303 K. The first, second and third columns correspond to the HW, SHW and SWHW scenarios, respectively. Interesting microscopic details revealed by studying the hydrate/water interface for the Tinit = 283 K HW simulation at different times as shown in Figure  4-6. At the start of the simulation, the blue color-coded water molecules represent water in the liquid phase and red color-coded water molecules are in the methane hydrate phase. The same color-coding is retained for each water molecule throughout the simulation. In the three snapshots it is observed that although hydrate cages remain mainly intact, there are exchanges of water molecules between the cages and the liquid phase. The exchange of water molecules between the liquid and hydrate phases, also observed in the simulations of Walsh et al. (Walsh et al. 2009) and Vatamanu et al. (Vatamanu & Kusalik 2006) illustrates the dynamic nature of the equilibrium between the water and hydrate phases. It is speculated that temporary opening and reforming of cages that accompanies the exchange of water molecules may provide a mechanism for the guest exchange between the surrounding liquid and methane hydrate phase. Figure  4-moleculeat 283 K. moleculeSnapshotsurface wdifferent simulatioFigure  4-Tinit = 3036. Snapshotss colored witSome distorts between ths of the confith the bounTinit values arn conditions3. After 2 ns K and met of the NVEh red and whion in the ine liquid phasigurations od water layere shown in F is much fast, the methanhane gas mo simulation oite) in contaterface of thee and intactf the water-ms (SWHW sigure  4-7. Mer in these se hydrate phlecules have f the HW scct with liqui methane hy methane hyethane hydrcenario) at thethane hydrimulations thase has comaggregated ienario with d water (watdrate phasedrate cages iate systems e end of theate dissociatan the HWpletely dissonto nano-bu methane hyder molecule and exchans observed. in the presen 2 ns simulaion under idsimulations ciated in thebbles. rate (water s colored as bge of water ce of silica tions with entical as shown in  simulation w54 lue) ith Figure  4-hydrate inThe decoby showinmethane the restricdirection.contours distance o7. Simulation the SWHWmpositions og F3 contouhydrate decoted mobility Moreover, iwhere F3 chf 5 to 8 Å in configurati scenario at f the methars in the simmposes fast of the watet is interestinanges betwee the z-directons at the eninitial tempene hydrate pulation yz-pler along the r layer betweg to note thn ~0.04 to ~ion. This pre d of the 2 nsratures of 28hases in the ane in Figurz-direction ten the solid at the dissoc0.1) is neithdicted range NVE runs f3, 293, and 3SWHW ande  4-8. In the han the y-dirsilica and hyiation front er straight n of 5 – 8 Å or decompo03 K.  HW scenariSWHW simection. This drate phases(shown in Fior sharp andover which dsition of metos are compulation, the could be du in the y-gure  4-8 by  extends ovecompositio55 hane ared e to er a n occurs is 1994),(MyFigure  4-(left) andIn the SWdissociaticonditionspecified Stupin 19(Uchida etemperatunatural poalso a diswater andinterestinin agreemenshakin et al8. Contours o HW scenariHW scenaron is faster ts in small poP and T, the92; Henry ett al. 2002) pre in poroures may be tribution of p the silica sug to study ott with the rep. 2009) for df the local Fo (right) at Tio, hydrate dhan the HWres shift to h driving forc al. 1999; Seredicted a des glass with adifferent thaore sizes. Hrface and thher pore shaorted simulecomposition3 order paraminit = 303 K.issociation o scenario. Thigher pressue for decomshadri et al. 2crease of 12verage pore n the slit-likeowever, the e gas phase)pes such as ation values  of methaneeter showinccurs more sis is consisteres and lowposition of m001; Wilder.3 and 0.5 K radii of 4 an pore studiegeneral beharemain similcylindrical poranging from hydrate in g the dissoclowly than tnt with the er temperatuethane hyd et al. 2001)for methaned 100 nm, red in this chavior of the sar for hydrores in mole 5 to 15 Å (contact with iation trendhe SHW casobservation res and therrate is higherFor example hydrate equspectively. Tpter and in rystem (interphilic mineracular simulatRodger  water.  under SWHe, but the that equilibrefore at a .(Handa & , Uchida et ailibrium he shape ofeality there iactions betwls. It would ions as well.56 W ium l.  the s een be  The formfeature ofunder difFigure  4-with wateFigure  4-303 K aftescenario. hydrate aTo quantz-directiot = 0 is inof cages omoleculedecompoation of nan the methanferent simula10. The methr and hydrop9. Bubble forr 2 ns simulMethane mond liquid waify the forman at differendicative of tf the sI hyds diffuse intoses faster, tho-bubbles fre hydrate detion conditiane bubbleshilic hydroxmation duriation times. Tlecules are ster phases artion of nanot times durinhe hydrate strate. As time the surroune methane nom the releacompositionons are show form to minylated silica ng the decomop: HW scehown by cyae represente-bubbles in g the simularucture whe passes and ding liquid wano-bubblesse of methan in these simn qualitativeimize the susurface.   position of nario, middn spheres. Td by hydrogthe solutiontion are coure methane mthe methaneater. In the  merge to a ge gas in theulations. Thly in Figure rface contacmethane hydle: SHW scehe water moen bonds in , the numbernted. The peolecules are hydrate phaSHW scenarreater exten liquid watere methane b 4-9 and quant of the metrate at initianario, bottomlecules initired and blue of methaneriodic metha entrapped ise decompoio where thet. As seen in phase is a ubbles formtitatively inhane molecul temperatu: SWHW ally in the , respectively molecules inne distributinside the roses and meth hydrate  Figure  4-1057 ed les re of .  the on at ws ane , at t 58 = 1.0 ns there are three small nano-bubbles at z = −40, −8 and 15 Å. After 2 ns simulation time, these bubbles merge into two larger nano-bubbles at z = −30 and 11 Å with approximate diameters (peak widths) of 17 and 30 Å. The presence of these bubbles in the liquid phase may explain the “memory effect” as it has been shown experimentally that nano-bubbles with radius of 50 nm constitute 600 cm3 of gas per 1 dm3 of water with internal pressures of up to 6 MPa. These nano-bubbles persist for more than two weeks at atmospheric pressure and room temperature.(Ohgaki et al. 2010) Previous simulation works (Rodger 2000; Hawtin et al. 2008; Vatamanu & Kusalik 2010) also concluded that the memory effect is related to the locally high concentration of aqueous methane in the melt rather than the presence of remaining hydrate structural precursors. Figure  4-simulatio4.4 CoMethane surfaces ~temperatuconditiontemperatu10. Profile ofns are all at Tnclusions hydrate deco40 Å apart res of 283, 2s. As expectres. The dis the number init = 303 K.mposition uwas simulate93 and 303 ed, the rate osociation ratof methane  nder a confid with moleK to investigf decomposes in the pre molecules alned geometrcular dynamiate the dissoition was lowsence of the ong the z-axy between twcs in an NVciation procer at lower silica surfaceis of the simo fully hydrE ensemble ess closer toinitial simulas were founulation box. oxylated siliat initial  reservoir tion d to be faste59 The ca r 60 than the non-confined water-methane hydrate simulations. During hydrate dissociation in the presence of silica, the hydrate core shrank in a stepwise manner with a curved decomposition front which extended over 5 to 8 Å. The presence of a water layer between the hydroxylated silica and methane hydrate alleviates the mismatch between the methane hydrate/hydroxylated silica lattices at the interface and therefore hydrate decomposed slower in the silica-water-hydrate-water (SWHW) scenario than the silica-water-hydrate (SWH) scenario. In both SWHW and SWH scenarios, the additional exposed surfaces made the methane hydrate decomposition faster than the hydrate-water (HW) scenario. Labeling of water molecules according to their origin in the hydrate or liquid water phase revealed that there is exchange between water molecules from the surrounding liquid with the cages of the hydrate phase. This demonstrates the dynamic nature of the hydrate structure when it is in equilibrium with an aqueous phase and provides a mechanism for guest exchange between the liquid and hydrate phases.   61 Chapter 5: Nano-bubble Formation during Methane Hydrate Dissociation 5.1 Introduction The novelty of this chapter is that gas reservoirs are added to the simulation box with methane hydrate and water, providing a possibility for methane released from the decomposing hydrate to migrate to the gas phase instead of artificially remaining in the liquid phase due to having no distinct gas phase to escape. Hydrate dissociation was found to lead to formation of methane nano-bubbles in the liquid phase even in the presence of gas reservoirs. The process of nano-bubble formation and the stability of these tiny bubbles are studied in detail, and the minimum required concentration of dissolved methane for the formation of bubbles as well as their effect on the hydrate growth rate are also determined. Nano-bubbles are known to have unique properties such as an unexpectedly very long life time (Zhang et al. 2007; Ohgaki et al. 2010) and a high adsorption of impurities.(Uchida et al. 2011) With the advancement of nanotechnology they are gaining industrial applications in areas such as water treatment,(Agarwal et al. 2011) nanoscopic cleaning (Seddon et al. 2012) and froth floatation.(Fan et al. 2010a; Fan et al. 2010b) The existence of surface nano-bubbles, also in the shape of nano-pancakes, has been proven by the atomic force microscopy technique (AFM).(Ishida et al. 2000; Lou et al. 2000) The freeze-fractured replica TEM technique (Uchida et al. 2011) and dynamic light scattering (Hampton & Nguyen 2010) have been used to capture the existence of nano-bubbles in aqueous solution. 5.2 Simulation System and Methods A 3×3×6 super cell of methane hydrate is constructed based on the crystallographic data (McMullan & Jeffrey 1965) of its unit cell and was placed between two bodies of water along the z-axis (approximately 3.5 nm thick). Additionally, the water phases are sandwiched by gas reservoirs along the z-axis. Figure  5-1 shows the initial configuration used for all dissociation simulations in this chapter.  To make the hydrate phase symmetrical along the z-direction, a row of water and methane were removed from one end. The resulting hydrate block is composed of 2430 water and 405 methane and 531, they corre0.148). Figure  5-the yz-plais sandwilayers, eainterface.shown in section ofSimulatiomotion arboundarydistance wcoupling with a timShake algWater anGuissani combinatthe availaphase equmolecules. Trespectively. spond to th1. Projectionne. A 3×3×6ched by gas ch 6 Å thick  Methane mred and whi the simulatns are perfore integrated conditions aas set to 1.5constant of 0e constant oorithm. d methane ar1993) force ion of force ble force fielilibrium temhe number oThe numbere stoichiome of the initial super cell ofreservoirs (Galong the z-dolecules are te. Hydrate wion box. med using G according tore applied in nm. Tempe.2 ps. Parrinf 0.5 ps. The modelled ufield, respectfields gives td models. Aperature at tf water and s of water antry of the fu configuratio hydrate pha). For analyirection sucshown as cyaater is showROMACS  the leap-fro all three dimrature was cello-Rahmane bond lengtsing TIP4Pively. Condehe best predll of the prodhis pressuremethane mod gas moleclly occupied n of G-W-Hse (H) is plasis, liquid wah that layersn spheres. On as hydrogsoftware (Heg algorithm ensions anontrolled usi coupling wh and HOH-ice (Abascal and Vega (Ciction for thuction runs is 302 K forlecules in thules in the shydrate (i.e. -W-G three-ced betweenter and hydr 5 and 16 lie xygen and hen bonds in ss et al. 200with a time d the short-rng Nosé-Hoas used to k angle of wat et al. 2005)onde & Vege methane h were perfor this force fie liquid and ystem are seoverall CH4 phase simul two bodiesate phase arat the water/ydrogen of red evident t8) version 4.step of 1 fs. ange potentiover algoriteep the preser are kept rand a uniteda 2010) shoydrate phasemed at 400 beld. gas phase is lected such t mole fractioation system of water (We sliced intohydrate liquid water he middle 5.5. EquationPeriodic al cut-off hm with a sure constanigid using th-atom (Guilwed that thi diagram amar. The thre62 3055 hat n =  on ) and  20 are s of t e lot & s ong e-63 The hydrate-water-gas system was first equilibrated at 302 K with a 500 ps NVT run. Next, a 500 ps NPT run at 302 K and 400 bar was performed to equilibrate the system under pressure-temperature conditions on the hydrate equilibrium curve. To study hydrate decomposition, first the temperature of the liquid water and methane reservoir phases were raised during a 250 ps simulation to the desired temperature (320-360 K). During all of these runs, the hydrate phase was frozen by position-restraining the oxygens of water and carbons of methane molecules. Finally, the system was allowed to evolve in an NPT ensemble at 400 bar and temperature of interest for a few tens of nano-seconds with no position-restraint on the hydrate phase. 5.3 Results and Discussion Hydrate dissociation is an endothermic process and energy in the form of heat must be provided to break the crystalline order of the hydrate (hydrogen bonds). Therefore, hydrate dissociation is associated with an increase in the potential energy of the system. This trend can be clearly seen in Figure  5-2 which shows the variation of the potential energy of the system with time for simulations at different temperatures. As expected, the hydrate dissociation rate is faster at higher temperatures. The block of hydrate is completely dissociated after about 1.9, 4.3, 14 and 41 ns at 360, 350, 340 and 330 K, respectively. Although the equilibrium temperature at 400 bar is 302 K, the dissociation process is relatively slow at 320 K, thus higher temperatures are applied to speed up the process and observe decomposition within a reasonable amount of time for the computations. 64   Figure  5-2. Potential energy of the system during hydrate decomposition. The increase in potential energy indicates the dissociation of hydrate. As expected, rate of hydrate decomposition is faster at higher temperatures. It takes approximately 1.9, 4.3, 14 and 41 ns at 360, 350, 340 and 330 K for the hydrate block to melt completely. 5.3.1 Dynamics of Hydrate Dissociation To probe the dynamics of hydrate dissociation the local order parameter, F3 (Rodger et al. 1996) defined in Equation  4-1 of Chapter 4 is used. It distinguishes water environments in liquid and solid (ice and hydrate) phases by measuring the deviation of triplets of oxygen atoms of neighboring water molecules from the ideal tetrahedral arrangement observed in the solid ice phases. The F3 parameter takes a value of ~0.09 for liquid water and ~ 0.015 for the hydrate phase. Since, the temperatures studied in this paper are all above the melting point of ice for the TIP4P-ice water model, ~272 K, certainly no ice is present in the system. To extract spatial information, the yz-projection of liquid water and hydrate phases is partitioned into 20 layers along the z-direction, each 6 Å thick. This gives each water pool four layers and the hydrate block 12 layers such that layers 5 and 16 fall on the water/hydrate interface at the smaller and larger z-values, respectively (see Figure  5-1). As reported previously,(Alavi & Ripmeester 2010; Bagherzadeh et al. 2012) the rows of the hydrate phase parallel to the exposed surface decompose in a concerted and stepwise fashion where each exposed row of cages parallel to the hydrate/water interface dissociate together. The decomposition front then moves inwards into the hydrate phase. Figure  5-3, which shows the time variation of the 65 F3 parameter for the layers in the simulation, demonstrates the step-wise behavior of hydrate decomposition at 330 K. The behavior of the system is similar at other temperatures except that the decomposition rate is faster at higher temperatures. Layers 5 and 16 are initially at the water/hydrate interface and therefore their corresponding F3 value is ~0.05 at time 0. Layers 6 through 15 are initially in the hydrate thus their F3 values is ~0.01 at time zero. Layer 4 and 17 are in the liquid phase and their F3 values remain unchanged at ~0.09. As dissociation progresses, the outermost layers of the hydrate phase are destabilized and after approximately 15 ns, they collapse completely. It is worth noting that after the F3 of the outermost layers (decomposition front) reaches ~0.05, corresponding to about 50% destruction of hydrate structure, the F3 of the next inner layer begins to rise which indicates a destabilization of this inner layer of hydrate cages. This sequence of cage collapse continues until the whole hydrate phase dissociates. An interesting trend observed in Figure  5-3 is that, after only three layers of hydrate remain, specifically, L9, L10 and L11, the hydrate phase collapses collectively and does not follow the stepwise decomposition. This is also evident in the potential energy profile of Figure  5-2, where near the end of decomposition, a sharp increase in the potential energy is visible.   Figure  5-3. The local order parameter, F3, as a function of time for different layers in z-direction during the simulation at 330 K. For layer assignment refer to Figure  5-1. Hydrate decomposition progresses in a concerted manner where rows of cages parallel to the hydrate/water interface collapse sfrom ~ 0.5.3.2 NBy followliquid phahomogenagglomerdecomponano-bubhowever,along thebubbles dFigure  5-hydrate pin water aFigure  5-projectionshown in Methane liquid phshow the has a roumethane imultaneous01 to ~ 0.09. ano-bubbling the trajese and formeously in theate and nuclesition continbles causing after the dia x- or y-direcistort into cy4 illustrates xhase after 10nd present a4. A sample s on the xz-pink and cymolecules frase as sphericross sectionghly sphericaremains dissly within a 1e Size ctories of the nano-bubbl liquid phaseate into tinyues and mor them to grometer of thetion, due to lindrical-shaz- and yz –p ns at 350 Ks bubbles (ssnapshot of t and yz-planan spheres, rom the decocal and cylins of the nanl shape wheolved in wat5 ns time pe system it ises. Initially, . After, reac gas spherese CH4 molew. The nano bubbles grothe applicatiped bubbleslane snapsho. In this statpherical and he decompoes. Methaneespectively. mposing hyddrical nano-o-bubbles onreas the left-er. riod. This is  found that ras the hydrathing the satu in the vicinicules enter th-bubbles at ws to becomon of period in order to ts of CH4 me, CH4 moleccylindrical).sing hydrate molecules inWater molecrate enter inbubbles. Th these planehand bubbleassociated weleased CH4e is dissociatration limit ty of hydratee liquid watfirst grow rae equal to thic boundaryminimize tholecules whules are see -water-gas sitially in theules are showto the gas ree projectionss. As seen, t is roughly cith the incre molecules aing, CH4 moin water, CH/water interer phase, thedially in all de simulation conditions, e interfacial tich were initn in the gas pystem at 350 hydrate andn by hydroservoirs and on the xz anhe right-hanylindrical. Aase of F3 valgglomerate ilecules disp4 molecules face. As y diffuse intirections;  box length the sphericalension. ially in the hase, dissol K shown in  gas phases gen bonds in remain in thd yz planesd nano-bubb small amou66 ue n the erse o the  ved are  red. e  le nt of 67 To quantify the size of these nano-bubbles, the water density map of the system was calculated by binning the yz projection of the simulation box into 1 Å × 1 Å squares, similar to that performed in Chapter 3. The boundary of the bubbles was defined as the locus of points with water density of 0.5 g/cm3; about half of its bulk value. A typical water density map of the system is shown in Figure  5-5 for a simulation at 350 K. Once, the boundary of the bubble is determined, the radius was estimated by fitting the profile to a circle.  Figure  5-5. Water density distribution for simulation at 350 K at 14 ns, after the dissociation of the hydrate phase. Blue area represents the gas phase (water density of zero) and green area corresponds to liquid water with average density of ~ 1 g·cm-3. Only simulations at 340 K and 350 K resulted in formation of nano-bubbles. Hydrate decomposition at lower temperatures was not fast enough to generate methane supersaturation in the liquid phase which would lead to the formation of methane-rich regions. At the lower temperatures, the methane molecules released from hydrate decomposition have enough time to diffuse into the gas reservoirs. The bubble radius profile for simulations at 340 K and 350 K as a function of the simulation time are plotted in Figure  5-6. 68  Figure  5-6. (A) yz-snapshot of the mass density distribution of methane molecules for the simulation at 350 K at 5 ns where the hydrate phase has decomposed entirely. Nano-bubble 1 (nb1), 2 (nb2) and 3 (nb3) are shown by circles. (B) The bubble radius profiles for simulations at 340 K and 350 K. At 350 K, initially three bubbles formed. The smaller bubbles, nb2 and nb3, gradually shrank by dissolving in water. The larger nano-bubble remained stable for about 20 ns and it also eventually began to dissolve in water and collapsed. At 340 K a bubble nucleates, at about 14.5 ns, grows to a cylindrical shape with a radius of ~11 Å, and remains stable for ~15 ns duration of the simulation. After formation of the nanobubble, the methane molecules inside the bubble start to diffuse out of the bubble into the liquid phase causing the bubble to shrink.  A similar phenomenon occurred for the simulation at 350 K. Initially, three bubbles formed. The smallest bubble (nb2) did not grow and collapsed after 1.5 ns (green curve in Figure  5-6). The middle size bubble (nb3) grew to a radius of ~8 Å (blue curve in Figure  5-6); however, it also eventually dissolved in liquid water after about 12 ns. The largest bubble (nb1), red plot in Figure  5-6, grew to a stable size of ~11 Å. It remained stable for about 25 ns; however, it also began to decrease in size and slowly dissolved in water and the methane molecules diffused into the methane gas reservoirs at the two ends of the system. 69 Figure  5-7 illustrates the distribution of methane molecules in different phases in the simulation during the hydrate decomposition. The number of methane molecules in the gas phase, which corresponds to the gas uptake curve usually measured in experiments, increases as a result of hydrate decomposition. The number of methane molecules in the aqueous phase decreases except when bubble 1 collapses at about 40 ns.  Figure  5-7. Number of methane molecules in different phases during simulation at 350 K. Overall, the number of molecules in the gas phase increases, indicating release of methane to the gas reservoir as a result of hydrate decomposition.  5.3.3 Critical Concentration of Dissolved Methane Required for Bubble Formation The accumulation of methane in liquid water depends on two factors: (1) the rate of hydrate dissociation which governs the rate with which methane molecules are released into the liquid phase and (2) the rate of diffusion of methane out of the liquid phase and into the gas phase reservoirs. In a saturated solution, this second rate should counterbalance the accumulation of methane in the aqueous phase. Temperature increases the rate of both of these processes. Ultimately, it is the balance of the two rates that determines the accumulation rate of methane in the liquid water phase. Hydrate dissociation is slowest at 320 K. Therefore, the rate of release of methane in liquid water is low. This provides methane molecules with enough time to diffuse through the liquid phase and 70 encounter the gas phase instead of accumulating in the liquid phase and forming methane-rich regions. On the other hand, at 350 K, the rate of methane release from hydrate decomposition is greater than the rate of methane diffusion into the gas phase. This leads to the accumulation of methane in the aqueous phase. The diffusion coefficient of methane in water-methane mixture (xCH4=0.06) at 320 K and 350 K are estimated to be ~3×10-5 and ~4×10-5 cm2/s for the same united atom model for methane as used in this work but with the SPC/E water model.(Shvab & Sadus 2014) The decomposition rate is also slow at low temperatures. However, the dissociation rate increases faster than the diffusion rate at higher temperatures and this causes accumulation of methane in water and thus formation of the nano-bubbles at higher temperatures. Although the accumulation of methane molecules in the aqueous phase (supersaturation) is necessary, it is not sufficient for the formation of nano-bubbles. A minimum degree of supersaturation is required to provide the non-equilibrium driving force for the dissolved methane to nucleate into a nano-bubble spontaneously (Rebata-Landa & Santamarina 2012). In other words, if the concentration of dissolved methane molecules (xg: mole fraction), surpasses a critical limit (xgc), methane bubbles can form. Additionally, the tiny bubbles that form after nucleation will remain stable only after reaching a critical size. The critical size is related to the supersaturation via the following equation (Lubetkin 2003): ݎ∗ = 2ܲݏߛ Equation  5-1Here, ݎ∗ is the critical size, ߛ is the surface tension, ݏ is the supersaturation and ܲ is the pressure at which the bubbles nucleate. For instance, methane nucleation at atmospheric pressure and 25 ºC requires a supersaturation of 80 (Bowers et al. 1996), i.e. methane mole fraction of ~ 0.002, which then yields an approximate critical radius of 18 nm. For a cylindrical bubble, the radius would be 9 nm. It is worth noting that despite the limitation of the system size in the simulations performed, the radius of the bubbles estimated in the previous section (~ 1.1 nm) are within the same order of magnitude and in acceptable agreement with the radius calculated from Equation 5-1. 71 To find the critical dissolved methane mole fraction (CDMM), a box of water molecules surrounded by methane gas reservoirs was simulated (see Figure  5-8). To achieve the desired dissolved methane mole fraction, a number of water molecules were randomly replaced with methane. To keep these simulations analogous to the hydrate dissociation simulations, the size of the simulation box and the overall composition of methane and water were kept the same and similar equilibration stages were applied as in the hydrate decomposition simulations. It should be noted that simulation box size may affect the bubble size and bubble shape. In a larger simulation box, the released methane molecules may not coalesce to create a large bubble with a diameter equal to the simulation box length along the x- or z-directions. Consequently, the bubble(s) may remain spherical and the maximum diameter (size) may depend on the simulation box size. Limitations of scope and computational resources prevented resolving of this issue in a comprehensive manner in this thesis. Despite these limitations of the simulation methodology, the phenomenon of nanobubble formation is a general feature of the hydrate dissociation process under certain conditions. Initially, a simulation at 350 K in which the methane concentration was half of that in the fully occupied hydrate phase (xCH4 = 0.074) was tested (Figure  5-8(A)). As seen in Figure  5-8(B), bubbles of methane formed in this simulation. Afterwards, the methane concentration was further halved (0.074)/2 to test bubble formation. Subsequently, the bisection method was followed to find the CDMM. Table  5-1 summarizes all the simulations performed to determine the CDMM. Figure  5-8(C) and (D) show the snapshots of the simulation with dissolved methane mole fraction of 0.042 at time 0 and 5 ns, respectively.  Table  5-1. Number of water molecules, nW, dissolved methane, nCH4 (aq), and gaseous methane, nCH4(g), for different simulations performed to determine the CDMM. simulation 1 2 3 4 5 nW 5744 5974 5859 5917 5945nCH4 (aq) 460 230 345 287 259 nCH4 (g) 460 808 674 742 775 xCH4 (aq) 0.074 0.037 0.056 0.046 0.042Bubble formed? yes no yes yes no  The simuand did nthe criticaand 400 b(xCH4 = 0Figure  5-in a waterand 400 bsubsequewith dissosimulatiothis simuThe profiand 400 bare given pressure a lations with ot lead to thl dissolved mar is 0.044 w.148).  8. (A) yz-sna-methane mar with the int depletionlved methann at 350 K anlation and thles of the zzar and the sin  Appendixround 400 binitial dissolve formation ethane molhich is ~30pshot of the ixture. (B) ynitial configu of methane e mole fractd 400 bar wie methane m-componentimulation 5 o D. With thear. ed methaneof bubbles. He fraction wh% of the moinitial state oz-snapshot oration showin the liquid ion of 0.042.th the initialolecules dif of the pressf this sectio operation o mole fractioence, it is eich is necesle fraction of the simulaf the simulan in (A). Thphase are ev (B) yz-snap configuratiofused into thure tensor fon (initial dissf the barostn of 0.046 astimated thasary for the f methane in tion with metion after 1 ne formation oident. (C) yzshot of the sin shown in e gas reservor the decomolved methaat both profind 0.042, rest for the prebubble form a fully occuthane mole s of NPT simf bubbles an-snapshot omulation aft(C). No bubir. position simne mole fracles show a flpectively, disent simulatiation at 350 pied sI hydrafraction of 0ulation at 3d the f the simulater 5 ns of NPbles formed ulation at 35tion of 0.04uctuating 72 d ons K te .074 50 K ion T in 0 K 2) 73 5.3.4 Effect of Nano-Bubbles on the Growth Rate of the Hydrate Phase The methane nano-bubbles in our simulations are not stable beyond a few ten nanoseconds but their lifetimes under the real conditions where the bubble nucleation is heterogeneous (due to the presence of impurities in the liquid water and the surface of the reactor or the porous media) could be longer. Surface bubbles (nano-pancakes) are known to be stable for long times, of order of seconds. If the nano-bubbles remain stable in the system, they may act as nucleation sites for the heterogeneous nucleation of hydrate crystals. Conde and Vega (Conde & Vega 2010) observed that in some of their methane hydrate growth simulations, methane gas formed bubbles in liquid water before leading to the growth of hydrate phase. The persistence of these nano-bubbles either in the solution or on the reactor surface may be the cause of, or may contribute to the so-called memory effect, introduced at the end of section  1.1. The bubbles can provide a readily available source of methane at higher pressures relative to the gas phase due to the young-Laplace equation and thus facilitate nucleation of hydrate crystals. Here the effect of the methane nanobubbles in the aqueous phase on the methane hydrate growth rate is investigated. After 4 ns of hydrate decomposition at 350 K the simulation is stopped and the configuration of the system was used to start a new simulation by cooling the liquid water to 280 K, a temperature below the hydrate-water-gas equilibrium temperature of 302 K at 400 bar for the applied force field.  Figure  5-9(A) illustrates the initial configuration for the cooled system. Panel (B) shows the configuration after 50 ns where hydrate growth, represented by the ordered hydrogen bonding network of water molecules, is evident. Figure  5-initially inFor claritthe hydraof the hydsides of thwhere ini AdditionFigure  5-parameteboth sidemethane the otherto the prefacilitatesmethane methane Laplace emethane there is a the growtmoleculegrowth o9. (A) yz-sna the hydratey, only methte phase is arate growthe initial hydtially a bubbal hydrate gr9(A)) which r for layers ps of the initiagas. Three laside only twsence of rea mass transfeat relatively hgas inside thquation.(Gomolecules cobubble, bothh. The numbs where founbserved on tpshot at 4 ns and liquid pane moleculelso shown as simulation arate templatle was presenowth on botis visually aparallel to thel hydrate coyers (18 Å ino layers havedily availabler of methanigh pressuree bubble is inldman 2009)ntribute to t the dissolver of methad to be incohe bubble sid of the partiahase are shos in the nan a reference t 280 K aftere. However,t. h sides of thparent in Fig hydrate surre. However total) of hy formed (12 methane soe to the hyd, which favoversely prop On the sidehe growth oed methane ane moleculerporated in te where hydlly decompown by hydroo-bubble areof position i 50 ns. As se more hydrate initial hydrure  5-9(B), wface. As seen, more growdrate cages a Å in total). urce (the burate interfacers the hydraortional to  where initiaf the hydratend methanes initially in the hydrate prate grew ansed hydrategen-bonds i shown. Onenside the hyden, hydrate ce cages are fate phase (sas monitor in Figure  5-th occurred re formed oFurther growbble) near th by providinte formationthe bubble silly no bubbl. However,  molecules fhe bubble ishase. This n extra 6 Å a at 350 K. Wan red and bl of the methrate phase. ages have gormed on thhown by theed quantitati10, hydrate on the side on the bubbleth on the be interface. Tg abundant /growth. Thze accordinge is present, on the side wrom the bub 99. 36 of thumber agreelong the z-diter moleculeue, respectivane molecu(B) yz-snapsrown on bothe right hand vertical linevely by the Fcages formedf the bubble side while oubble side is he bubble source of e pressure o to the Youonly dissolvehere initiallyble contribuese methanes with the furection (abo74 s ely. les in hot   side s in 3  on  of n due f ng-d  te to  rther ut 75 half of the length of the sI unit cell). The cross section on the xy-plane is a 3×3 replica and thus about 36 methane molecules should be captured inside the hydrate cages based on 8 methane molecules per unit cell (3×3×0.5×8 = 36). The bubble remains in the vicinity of the hydrate interface and supplies methane for the hydrate growth. Once the bubble becomes so small that it is no longer stable, it collapses and creates a large local concentration of methane in the aqueous phase and further facilitates the growth of hydrate.   Figure  5-10. F3 profile of the hydrate growth simulation at 280 K and 400 bar. L4 to L10 and L11 to L17 refer to the slices on the non-bubble side and bubble side, respectively. As can be seen, three layers of hydrate cages have grown on the bubble side whereas there are only two layers of hydrate growth on the other side of the initial hydrate template.  It is noted that the thickness of water layers (mass transfer time of methane to the gas phase) on the sides of hydrate phase affects the stability and the conditions required for the formation of these bubbles.(Weijs et al. 2012) Studying this effect is beyond the scope of this work; however, I believe that the behavior of the system leading to bubble formation is properly/qualitatively captured in these simulations. 5.4 Conclusions Molecular dynamics simulations show that methane hydrate decomposition may lead to the formation of nano-bubbles. Nano-bubbles form as a result of the release of methane in the liquid water phase which super-saturates the aqueous solution. Super-saturation is countered by the 76 diffusion of methane molecules out of the liquid phase (release to the gas phase). Ultimately, the balance between these phenomena determines the accumulated level of methane in the liquid phase. It is found that at low temperatures, the hydrate dissociation rate is slow and this allows the methane molecules to escape to the gas phase without bubble formation. However, as the temperature increases, the rate of release of hydrate methane to the liquid phase becomes greater than that of the diffusion rate and therefore methane accumulation level in the liquid phase increases. If dissolved methane concentration passes a threshold, then a bubble forms. For the specific system size in these simulations, these bubbles grow to a diameter of about 11 Å, remain stable for about 45 nano seconds and eventually re-dissolve in the liquid phase and diffuse to the gas phase reservoir. It is also estimated that a critical dissolved methane mole fraction which is required to form bubbles at 350 K and 400 bar is 0.044, corresponding to about 30% of the methane mole fraction in a fully occupied sI hydrate. Bubbles may form under different supersaturation conditions at different temperatures and pressures. In the micro-second simulation of methane hydrate nucleation by Walsh and coworkers (Walsh et al. 2009) at 250 K and 500 bar  where they observed bubble formation before the nucleation of hydrates, the average methane mole fraction in bulk water after bubble formation was 0.039.   77 Chapter 6: Mechanism of Binding of wf-AFP to the Hydrate Surface 6.1 Introduction The detailed description of the wf-AFP molecule is given in section  2.4 of Chapter 2. The current chapter discusses the first molecular dynamics simulation studying the interactions of the winter flounder antifreeze protein (wf-AFP) with the surface of a methane hydrate crystal. The aim is to obtain direct molecular insights on the mechanism of action of wf-AFP in inhibiting hydrate growth and to determine the responsible structural features of the AFP. Such studies lead to a more comprehensive understanding of the function of AFPs and will be a first step in the rational design of new green and more efficient kinetic hydrate inhibitor molecules with potential superior activity. 6.2 Computational Methods Three sets of simulations at 275 K and 800 bar were performed: APF/Water to equilibrate the protein in water under the simulation conditions and two production runs where the interactions between the wf-AFP and the hydrate surface are closely monitored. A 7×7×3 unit cell replica block of methane hydrate (unit cell taken form the crystallographic data (McMullan & Jeffrey 1965)) with all cages occupied by methane is prepared and one (case I) and two (case II) AFP molecules (to test the effect of AFP concentration) are placed ~0.1 nm away from the hydrate surface with threonine residues (which constitute the ice binding surface of the wf-AFP) pointing toward the hydrate surface (see Figure  6-1). The 7×7×3 replica is selected such that the protein will not interact with its own image due to periodic boundary conditions. The length of the wf-AFP molecules is ~ 5 nm and the cut off radius is 1.5 nm. Therefore, a simulation box length of at least 8 nm is required. The unit cell length of the cubic sI hydrate is ~ 1.2 nm and a replication of at least 7 times will satisfy the above requirement. The hydrate is in contact with water on both sides with a thickness of 3.6 nm along the z-axis and the system is sandwiched between two reservoirs of methane gas such that the overall composition of water and methane in the simulation cell is kept approximately equal to that of the stoichiometry of a fully occupied sI hydrate (i.e. methane mole fraction = ~0.148). The number oresulting Figure  6-The AFP phase andrespectivemiddle ofnm, respesimulatioused for tthe threonLeft pane Table  6-1Simulati[total timProtein Water (lWater (hMethaneMethane f molecules simulation b1. yz-projectiis initially al water in thly. Hydrate  the box. Simctively. (C) n. The protehe case II siine residues: xy view; top. Number ofon set e] (wf-AFP) iquid) ydrate)  (hydrate)  (gas) A B present in eaox is 8.42×8on of the iniigned with te liquid phaswater and mulation boxRelative posiin is roughlymulation. Th are shown.  right: xz vi molecules pAFP/Wat[5 ns] 1 19669 - - - ch set of sim.42×15.64 ntial configurhe long axis e are represeethane are s sizes along tion of wf-A 0.1 nm awaye cyan helixLiquid wateew; bottom rresent in eacer AFP/HCase I [200 ns1 165526468 1029 2879 ulations is sm.        ations of thealong the x-dnted by a cyhown as red the x, y and FP with resp from the hy represents tr molecules aight: yz viewh simulationydrate/Wat] Case[1352 167664681029301Cummarized  case I (A) anirection. wfan ribbon, chydrogen boz-directions ect to the hydrate surfache wf-AFP. Fnd gas phas.  set er/Gas  II  ns] 2   1 in Table  6-1d case II (B-AFP, methayan spheres nds and cyaare 8.421, 8.4drate block e. Same relaor clarity one methane a. The size of ) simulationne in the gaand blue linn spheres in21 and 15.63in case I tive distancely side chainre not shown78 the s. s es,  the 9  is s of . 79 The GROMACS 4.6.1 MD code (Hess et al. 2008) was used to run all simulations with the leap-frog algorithm and a time step of 2 fs. Effect of time step (0.5, 1.0 and 2.0 fs) is tested for a small system with no inhibitor in  Appendix F. Water and protein bonds were constrained using the LINCS algorithm and short range interactions were truncated at 1.5 nm. Coulombic interactions were treated using the particle mesh Ewald (PME) method. The temperature and pressure were controlled by Nosé-Hoover thermostat with a time constant of 0.1 ps and by Parrinello-Rahman barostat with a time constant of 2 ps. Periodic boundary conditions were applied in three dimensions. Water was modelled by TIP4P-ice,(Abascal et al. 2005) methane by a united atom potential,(Guillot & Guissani 1993) and the amber99sb-ildn force field (Lindorff-Larsen et al. 2010) was used for the protein. In each set of simulations, to equilibrate the system, first, an energy minimization of water around the wf-AFP was run. This was followed by a 0.5 ns constant volume – constant temperature (NVT) simulation for temperature equilibration and a subsequent 0.5 ns constant pressure – constant temperature (NPT) run to equilibrate the pressure and temperature to the final desired values. For the case I and case II simulations, the hydrate phase was kept frozen during the equilibration stages. Finally, the hydrate is unfrozen and hydrate-AFP-water-methane gas system was allowed to naturally evolve in an NPT ensemble at 275 K and 800 bar for 135-200 nanoseconds. This thermodynamic condition is well within the hydrate phase stable region for the TIP4P-ice methane hydrate force field where the equilibrium temperature for a pressure of 400 bar is estimated at 302 K (Conde & Vega 2010) which is reasonably close to the experimental value of 297.23 K.(Sloan & Koh 2008) These conditions guarantee that there will be no large scale dissociation of the hydrate phase during the simulations. 6.3 Results and Discussion 6.3.1 AFP/Water Simulation Root Mean Square Deviation (RMSD) of the backbone of the protein is plotted in Figure  6-2. It initially increases and after ~0.3 ns, gradually levels off at about 0.2 nm, suggesting that there is no substantial change in protein structure at the conditions of simulation (275 K and 800 bar).  80  Figure  6-2. Root mean square deviation of the backbone of the protein during the wf-AFP/Water simulation. RMSD remains below 0.3 nm suggesting no significant change in the α-helical structure of wf-AFP under simulation conditions of 275 K and 800 bar. The initial and final configurations for the wf-AFP/Water simulation are shown in Figure  6-3(A) and (B). It is evident that the wf-AFP remains stable under the function of the amber99sb-ildn force field at simulation conditions of 275 K and 800 bar.  The radial distribution function (RDF) of water oxygen (OW) around hydroxyl oxygen (OG1) and methyl carbon (CG2) of THR residues during the wf-AFP/Water simulation are plotted in Figure  6-3(C). The first peak for the OW-OG1 RDF occurs at 0.282 nm and the first minimum is at 0.356 nm which is indicative of hydrogen bonding between the hydroxyl group of the THR residue and water molecules. The first peak for the OW-CG2 RDF occurs at a larger distance, 0.37 nm, which is expected, due to the hydrophobicity of the methyl group. For reference, the radial dimensions from the center of the elliptical structure I hydrate large cages are 0.254 and 0.306 nm. Water structuring around the hydrophobic methyl group is more pronounced and the OW-CG2 RDF shows a clear second peak at 0.632 nm. However, this is not evident in the OW-OG1 profile which does not show evidence of water structuring at intermediate distances. The influence of hydrophobic solutes on local water order is a well-known phenomenon.(Galamba 2014) Figure  6-helical stris 8.421 nmblue, resp(OW-OGOW, bluepeak, whimethyl grThe averarespectiveal. 1963) within th 6.3.2 AThe final segment nfor the reare approroughly presidue 26interface.3. (A) Initial ucture of wf and blue lectively. (C)1, black), the). The first ple the first poup. ge length anly. Furthermfor rotation e range for aFP/Waterconfiguratioear the headst of the simximately parerpendicular have unrav  (0ns) and (B-AFP is evidines represen Radial distr THR termieak of the Oeak of OW-Cd radius of tore, averagearound the Nn α-helical st/Hydrate Sin of the case  of the AFPulation. At 2allel to the h to the hydraeled (denatu) final (5 ns)ent in the mit water moleibution funcnal methyl cW-OG1 RDFG2 occurs ahe helix duri Φ and Ψ, th-Cα and Cαructure (Φ=mulation (cI simulation became anc00 ns, parts ydrate xy-surte surface bred), and the configuratioddle of the scules. x, y ations for watarbon (OW-C occurs at 0t 0.37 nm dung 5 ns of sie Ramachan-C bonds are-60°±10° anase I)  is shown in hored to theof the wf-AFface. The pretween resid protein tail ns of the wfimulation bond z axes areer oxygen wG2, red) an.282 nm, sime to hydrophmulation aredran dihedr -67.2° and d Ψ=-45°±Figure  6-4(A hydrate surP corresponotein bends ues 15 and 2has deforme-AFP/Waterx. The size  colored as ith the THRd other wateilar to the Oobic charac 4.77 nm anal angles (Ra-37.5°, respe10°). ). In case I, aface and remding to residat residue 145. Loops of d near the g simulation.of the cubic red, green an hydroxyl oxr oxygens (OW-OW first teristics of Td 0.242 nm, machandranctively, whict about 55 nained robusues 1 throug such that itthe helix aftas/water 81   α–box d ygen W-HR  et h are s, a tly so h 13  is er Figure  6-ns. The bbars in paresidues ilines. In Figurein the texsurface. T9ALA res6.3.3 PFigure  6-protein none can s4. (A) yz- andinding residnel (C). Liqn the empty  6-4(B) one ct), shown in he zoomed idues in the rotein Defo5(A) illustratear the wateree that altho (B) xy-projues, 2THR, uid water is nhalf-cages isan see that tyellow, are ain view of Fihalf-cages atrmation dues snapshots/gas interfacugh betweenection and (6ALA and 9Aot shown in evident in phe methyl pligned approgure  6-4(C) the hydratering the ca of the protee has denatu 120 and 14C) zoomed inLA, are sho panel (B) foanel (C). Hyendant groupximately abodepicts the c/water interfse I Simulain structure red and lost0 ns a substa yz-view of wn in yellowr clarity. Thdrogen bons of the binve the emptapture of thace.  tion at different t its helical stntial portion the case I sim in panel (B)e entrapmends are shownding residuey half-cages e binding 2Times. The taructure. Mo of the proteulation at 2 and as thickt of binding  as red dashs (discussed of the hydraHR, 6ALA, il segment ore importantin has unfol82 00  ed later te and f the ly ded, 83 it refolds to its helical structure afterwards. Nevertheless, the section of the wf-AFP which is immersed in the liquid water phase, especially the part near the hydrate block, retains its helical structure for the duration of the simulation. The thickness of liquid water layer plays a role in the observation of this behavior. If larger amount of water is placed on both sides of the hydrate phase such that the thickness of the water layers is longer than the length of the protein (~ 5 nm), the protein may not migrate to the gas phase, becomes unsolvated and thus deforms.   Figure  6-5. (A) The structure of wf-AFP at different times during the case I simulation. Cyan, red, black, orange and transparent yellow correspond to 5, 80, 120, 140 and 170 ns, respectively. The blue line on the left shows the approximate position of the hydrate/water interface while the one on the right shows that of the water/gas interface. Other molecules are not shown for clarity. (B) RMSD of the wf-AFP during the case I simulation. There is a significant change in the backbone structure of the protein after 60 ns. This substantial structural change can be seen as unfolding of the section of the protein which is not adjacent to the water/hydrate interface. Figure  6-5(B) shows the RMSD of the backbone of the wf-APF for case I during the 200 ns of the simulation. The RMSD slowly increases until about 75 ns, after which substantial structural deformation occurs leading to RMSD values as high as 1 nm. This significant structural change is related to the end of protein diffusing into the gas/water interface. 6.3.4 AFP/Water/Hydrate Simulation (case II) The final configuration of the case II simulation is shown in Figure  6-6(A). In the case II simulation, one protein has bound to the methane hydrate surface. Neither of the proteins in this simulation have deformed to a significant extent and both maintained their helical structure throughout the 135 84 ns simulation (the RMSD of both of the proteins remained well below 0.4 nm). Therefore, the protein distortion behavior observed in the case I simulation is not a general feature of the AFP-hydrate-water system. Figure  6-6 shows that the binding pattern of the wf-AFP to the methane hydrate surface in case II is similar to that of the case I. Similar to Figure  6-4(B), in Figure  6-6(B), it can be seen that the binding residues of the wf-AFP are positioned approximately above adjacent empty half-cages on the surface of the hydrate phase. Panel (C) of Figure  6-6 also depicts the entrapment of these residues in the empty half-cages at the hydrate/water interface. Note that in both case I and case II simulations, the protein molecules are not aligned parallel to the x- and y-directions (rows of cages) of the hydrate surface. Any set of neighboring empty half-cages can act as binding sites for the wf-AFP.  6.3.5 The F3 Order Parameter  To probe the local water environment around the wf-AFP residues, the F3 order parameter (Baez & Clancy 1994) around the side chains (pendant groups) of all protein amino acid residues is calculated. F3 measures the deviation of water network from the ideal tetrahedral arrangement observed in ice and hydrate phases. F3 values of 0.015 correspond to water molecules arranged in hydrate-like environments whereas values of 0.09 correspond to a liquid water environment in these simulations. The F3 for water molecules near the hydrate/water interface is expected to have the average value of ~0.05 since water molecules at the interface have characteristics of both the ordered hydrate phase and the disordered liquid water phase. Here the F3 values were averaged for water molecules within a spherical shell of radius 0.6 nm around each terminal heavy atom of the pendant groups of all 37 amino acid residues of wf-AFP (see  Appendix F). To reduce fluctuations, 125 evenly spaced sample configurations within each 5 ns time frame of the simulation (i.e., one sample every 40 ps) are taken and then the F3 is averaged.  Figure  6-simulatio(B) and aof bindinred dashe 6. (A) yz- andn at 100 ns. Ts thick bars g residues ind lines.  (B) xy-projhe binding in panel (C). the empty hection and (residues, 24T Liquid watealf-cages is C) zoomed inHR, 28ALAr is not showevident in pa yz-view of  and 31ALAn in panel (nel (C). Hyd the protein 1, are shown iB) for clarityrogen bond in the case n yellow in p. The entraps are shown 85 II anel ment as 86  Figure  6-7. (A) F3 values for the side chains of all residues of wf-AFP for the case I simulation at 172.5 ns. 2THR, 6ALA and 9 ALA have F3 values below 0.05 suggesting they are trapped within the half-cages at the hydrate/water interface. Refer to  Appendix F for the naming scheme of the atoms. (B) F3 vs. time profiles of the methyl side chain of the four THR residues and (C) ALA 3, 6,7 and 9 of the wf-AFP in case I simulation. The adsorption of wf-AFP to the hydrate surface is a collaborative action between THR and ALA residues. In this case: 2THR, 6ALA and 9 ALA. (D) A snapshot of the case I simulation at 195 ns which shows that the methyl group of 2THR is trapped inside a half-cage at the hydrate surface and its hydroxyl side chain forms two hydrogen bonds, shown as green lines (emphasized by arrows), with local water molecules. 87 A typical F3 profile of all residues of wf-AFP in case I simulation at 172.5 ns is plotted in Figure  6-7(A). Interestingly, only 2THR, 6ALA and 9ALA residues display F3 values below the expected interfacial values, suggesting that these residues have an ordered neighborhood of water molecules and are effectively incorporated into the empty half-cages at the hydrate surface.  The time variation of the order parameter, F3(t), around the methyl side chain of all THR residues and ALA residues number 3, 6, 7 and 9 during the case I simulation are plotted in Figure  6-7(B) and (C), respectively. The carbon atoms of the THR and ALA methyl group side chains are designated as CG2 and CB. Figure  6-7(B) illustrates that the F3 values of THR residues number 13, 24 and 35 remain above the interfacial value of ~0.05. Therefore, these residues are not anchored to the hydrate surface. However, F3 of 2THR begins to decrease at ~55 ns and stays below the interfacial value afterwards suggesting that 2THR is trapped inside a half-cage at the hydrate/water interface as evident in Figure  6-7(D). The F3 profile of ALA residues in Figure  6-7(C) displays an interesting trend. F3 of 9ALA starts to decrease at about the same time as 2THR methyl group and ~100 ns later 6ALA shows a sudden drop, which implies its entry into an empty hydrate half-cage as well. The F3 profiles for the two protein chains for the case II simulation are given in Figure  6-8. Only protein 1 was effectively attached to the hydrate surface and 24THR, 28ALA, 31ALA and 35THR displayed F3 values noticeably lower than the interfacial values. Pendent groups of other wf-AFP residues are seen in Figure  6-7(A) and Figure  6-8 to lead to partial ordering of water molecules around the protein. Even residues which do not bind to the hydrate surface have F3 values less than the solution value of 0.09.  88  Figure  6-8. F3 values of all side chains of wf-AFP residues for case II simulation at 67.5 ns. As seen 24 THR, 28 ALA and 31 ALA (designated by red circles) as well as 35 THR for protein 1 have F3 values well below 0.0525 (red broken line). This suggests these residues have an ordered water neighborhood and are trapped within the half cages at the hydrate surface. Protein 2 is not effectively anchored to the hydrate surface at this time.  The F3(t) profiles of the two protein chains during the case II simulation are given in Figure  6-9. For protein 1 as seen in panels (A) and (B), after about 35 ns, THR residues number 24 and 35 as well as ALA residues number 28 and 31 show F3 values well below the interfacial value. Within the time frame of the simulation, protein 2 did not adsorb to the hydrate surface and its THR methyl F3 values remain above the interfacial value except for the 35THR residue where it is fluctuating below and above the interfacial value. The ALA residues of protein 2 are also oscillating near the interfacial 89 values. 3ALA and 32ALA residues seem to be adsorbed to the hydrate surface. However, there is no consistent anchoring of protein 2 to the hydrate surface due to the absence of stable (durable) docking of the THR residues to the hydrate surface.   Figure  6-9. F3 vs. time profiles of the methyl side chain of (A) all of the THR residues and (B) ALA 28, 31, 32 and 36 of protein 1 in case II simulation. The adsorption of wf-AFP to the hydrate surface is a collaborative action between THR residues and ALA residues, in this case 24THR, 28ALA and 31 ALA. F3(t) profiles of the methyl side chain of (C) all of the THR residues and (D) ALA 3, 7, 14 and 32 of protein 2. During the 135 ns of the simulation, this protein is not effectively bound to the hydrate surface. The protein 2 THR residues do not illustrate F3 values which stably remain below the interfacial value of ~0.05. Based on 13C solid-state NMR spin lattice relaxation time analysis (Mao & Ba 2006) and mutations studies (Baardsnes et al. 1999) of the AFP – ice binding, it is suggested that the IBS of wf-AFP consists of the THR side chains and the (i+4)ALA and (i+8)ALA residues. However, it is not microscopically determined why this particular combination of THR and ALA residues contributes to the ice antifreeze action. Our simulation results indicate that the wf-AFP adsorbs to the methane 90 hydrate surface by the cooperative entrapment of the pendant methyl group of THR residue at location i and two pendant methyl groups of (i+4)ALA and (i+7)ALA residues further down the protein sequence. These groups constitute the hydrate-binding surface (HBS) of the wf-AFP (i = 2, 13 or 24). These hydrophobic methyl pendant groups tend to enter into the empty half cages at the methane hydrate interface and effectively act as a guest. This stabilizes the hydrate surface by optimizing the hydrophobic interaction of methyl groups with surrounding polar water molecules, similar to what occurs in the formation of gas hydrates themselves from liquid water and methane gas. The size of the methyl group very closely resembles that of CH4 thus making it an ideal guest to stabilize hydrate cages. An interesting feature observed in Figure  6-7(B) and Figure  6-9(A) is that for wf-AFP to adsorb on the hydrate surface, all four THR residues may not be necessary and wf-AFP can bind to the hydrate surface via only one or two of its four THR residues. Longer simulation times may show binding along the entire length of the AFP. 6.4 Hydrogen Bonding Analysis The (i)THR hydroxyl group can form hydrogen bonds with nearby water molecules at the hydrate surface, leading to further stabilization of the wf-AFP binding at the hydrate surface. Figure  6-7(D) is a snapshot at 195 ns of the case I simulation illustrating how the pendant methyl side chain of 2THR residue is trapped inside one of the half cages at hydrate surface while the THR hydroxyl group forms two hydrogen bonds (green lines) with the local water molecules.  Figure  6-195.608 anwater/hythe hydroFigure  6-195.608 acage at thnearby wThe averaresidues awf-AFP inhydrogenhydrate suobservati10. The yz-prd (C) 195.80drate interfaxyl group an10 shows thrnd 195.800 ne water/hydater moleculge number ond neighbor water and t bonding in rface and laon agrees wiojection of t0 ns. The 2Tce while therd nearby waee more snas where it crate interfacees, respectivef hydrogen ing water mhe case I whethe first 5 nsst 5 ns wherth the findinhree snapshHR methyl e are 1, 0 andter moleculepshots of than be seen th and the hydly. bonds betweolecules, givere wf-AFP in of the case Ie the proteings that no sigots of the casgroup is trap 2 hydrogens. e 2THR residat its methyroxyl groupen the pendn in Table  6teracts with simulation w was adsorbnificant gaine I simulatioped in the c bonds (showue in the cal group is tra forms 0, 1 oant OH grou-2, are not sthe hydrate here the wf-ed to the hyd of hydrogen at (A) 195enter of the hn in green se I simulatiopped in the r 2 hydrogeps of indiviignificantly dsurface. MorAFP was norate surfacen bonding is .418 ns, (B) alf-cages atlines) between at 195.418center of a hn bonds withdual THR ifferent betweover, the t adsorbed t are similar. T obtained fo91  the n , alf- the een o the his r the 92 wf-AFP/ice complex in the interfacial region compared to the solution.(Dalal et al. 2001) This explains why hydrogen bonding between THR residues and water molecules at the hydrate surface is not the main driving force for the adsorption of wf-AFP to the hydrate surface. Hydrogen bonding between the antifreeze molecule and the ice or hydrate surface leads to more efficient inhibition, but it is not sufficient for antifreeze effects. A mutant of wf-AFP where hydroxyl groups of THR were substituted by pendant methyl groups still showed inhibitory effects, retaining 30% to 50% of the activity of the original AFP.(Haymet et al. 1998; Zhang & Laursen 1998) The partial loss of inhibition activity is attributed to the absence of hydrogen bonding. But the fact that this mutant is still active for inhibition demonstrates that hydrogen bonding is not the structural feature solely responsible for inhibition activity. Table  6-2. Average numbers of hydrogen bonds between individual THR residues and neighboring water molecules. System THR  2 13 24 35 wf-AFP/Water (0-5 ns) 1.78 1.48 1.53 1.32 Case I (5-10 ns, not bound to hydrate surface)  1.84 1.57 1.61 1.75 Case I (195-200 ns, bound to hydrate surface) 1.77 1.92 1.54 1.22  6.5 Conclusions It is found that for hydrate inhibition, wf-AFP is bound to the hydrate surface by cooperative anchoring of the pendant methyl side chain of (i)THR, (i+4)ALA and (i+7)ALA residues to the half-cages at the hydrate/water interface. However, compared to the ice inhibition where AFP adsorption should occur along a very particular direction on the ice surface, in hydrate inhibition the wf-AFP is free to align on the hydrate surface and position its methyl side chains from THR and ALA residues to optimally bind to the half-cages at the hydrate surface. The average distance of the centers of the neighbor cages of methane hydrate is between 0.602 and 0.672 nm. In wf-AFP, (i)THR - (i+4)ALA residues, and, (i+4)ALA - (i+7)ALA residues are spaced similarly for i = 2, 13 and 24 at (0.541, 0.558), (0.538, 0.569), and (0.551, 0.600) nm, respectively. The protein structure is flexible and can complement the hydrate surface by matching the half-cage spacing.  93 The major consequence of the adsorption of wf-AFP on the hydrate surface is that it partially covers the hydrate surface and adds resistance to mass transfer of water and methane to the hydrate surface, which leads to growth inhibition. In addition, the side chain functional groups on the protein residues facing the solution structure surrounding waters differently than the hydrate, thus adding entropic resistance for the hydrate lattice growth. Hydrate growth simulations under methane pressure in the presence of the wf-AFP are planned to perform to verify these predictions. The formation of cages of water molecules which surround hydrophobic groups may be a general aspect of the function of antifreeze proteins on ice as well as hydrates. The presence of hydrate half-cages facilitates the AFP binding on gas hydrates, but the nature of the interaction is general.  The simulations identify the involved antifreeze groups and other structural aspects of the protein and hydrate surface which are needed to incorporate the binding.  In principle, with additional information on the pendant groups of other active AFP molecules, one should be able to design new molecules that integrate different residues with required spacing between them to attain superior performance. Simulations indicate that the non-binding sites of the wf-AFP protein are not essential to the operation of the protein as a hydrate inhibitor, except to give it the overall helical structure and enhancing its water solubility.  94 Chapter 7: Conclusions and Recommendations 7.1 Contributions to the Hydrate Knowledge The following list summarizes the key findings of the work undertaken in this thesis: 1. In order to understand the fundamentals of formation and decomposition of natural gas hydrates it is essential to include solid surfaces that can represent the porous media. Silica as the main constituent of sandy sediments is incorporated into the simulated systems. Prior to investigating hydrate dissociation in the presence of silica, it is important to characterize the interaction between solid surface and water. It is found that the water layer adjacent to the silica surface has different dynamic and hydrogen-bonding characteristics than the bulk water. This structural perturbation of water can be detected up to about 6 Å away from the solid surface. The water phase formed a concave meniscus reflecting the hydrophilicity of the surface. The calculated water-silica contact angle was in good agreement with the experimental values. Gas molecules (CH4 and CO2) were preferentially accumulated on the meniscus and the silica surface. It is proposed that a nucleation site for heterogeneous nucleation of hydrate is located on the meniscus (interface between water phase and gas phase) a few Å away from the solid surface where guest molecules are abundantly available and the mobility of host molecules (water) are restricted as compared to the bulk.  2. By studying the decomposition of hydrates between two fully hydroxylated silica surfaces ~40 Å apart it is found that the partially immobilized buffer layers that form on the surface of silica (bound water) are essential to alleviates the mismatch between the crystalline structure of silica and hydrate. In fact a hydrate which was in direct contact with silica surface was less stable and decomposed faster compared to a case where the partially immobilized buffering water layer exists. The decomposition of the hydrate phase progressed faster near the silica surface and led to a concave dissociation front. As expected, the dissociation rate was faster in comparison with simulations of decomposition of hydrate in the absence of a silica surface. Similar to the decomposition of the hydrate only in contact 95 with bulk water, the hydrate block melted in a step-wise fashion where rows of hydrate at the decomposition front collapsed sequentially. 3. Molecular dynamics simulations show that methane hydrate decomposition can lead to the formation of nano-bubbles. Nano-bubbles of methane will form in the liquid phase if (1) the dissociation rate is fast enough to generate a super-saturation that would lead to bubble nucleation and (2) the diffusion rate of methane is slow enough to prevent fast migration of methane to the gas phase and (3) the concentration of dissolved methane is above a threshold. It is suggested that hydrate decomposition can be used as a method for forming nano-bubbles in a solution. This could be useful in water treatment technologies where nano- and micro-bubbles of oxygens are used. 4. This study indicates that the mechanism of binding of wf-AFP to hydrates is different from ice and involves specific engagement of the methyl pendant groups of the THR and ALA (four and seven residues down the protein sequence) residues of the AFP molecule with the empty half-cages at the hydrate surface. The end results of understanding the mechanism of wf-AFP binding to the hydrate surface could be to design a synthetic, organic analog which functions based on a similar mechanism or modify the existing commercial inhibitors to improve their activity.   7.2 Recommendations for Future Work 1. Only sI methane hydrate is studied in this work. It would be beneficial to experiment with the sII hydrates of natural gas (C1/C2/C3) to better represent the natural gas hydrate samples.   2. The formation of nano-bubbles was studied where there was no silica surface present in the system. It would be interesting to investigate the formation of nano-bubbles under conditions where heterogeneous nucleation for the bubbles is possible. The stability of the surface nano-bubbles may be different than the bulk nano-bubbles. It would also be informative to test other minerals such as clays in the simulations. 96 3. When methane hydrate crystals are formed inside the hydrocarbon transportation/processing facility, industry needs to remove the blockage by dissociating the hydrate. Since AFPs work by adsorbing to the hydrate surface, it would be of practical interest to determine the effect of AFP on the dissociation behavior of the hydrate plug. Sharifi and coworkers (Sharifi et al. 2014) observed that in the presence of AFP I and III hydrate decomposition occurred at higher temperatures compared to the control solution and an additional peak was present in their differential scanning calorimetry profiles. 4. The inhibition simulations can become comprehensive by testing other types of hydrate-active AFPs such those given in Figure  2-4 to determine their functional mechanisms. This can lead to the design of new inhibitor molecules with potential superior activity or improvement of the performance of the existing hydrate inhibitors by modifying their molecular structure. 5. Differential scanning calorimetry (DSC) experiments can be performed to estimate the average induction time (as a measure of the effectiveness of the inhibitor) for methane hydrate formation from a solution containing the mutated wf-AFP to confirm the mechanism proposed in Chapter 6. The methyl group of hydrate binding residue sets (iTHR, (i+4)ALA and (i+7)ALA) can be replaced by a hydroxyl group (i.e. substituting the THR and ALA residues by SER residues) to remove the contribution of the hydrophobic interactions.    97 Bibliography Abascal, J.L.F. et al., 2005. A potential model for the study of ices and amorphous water: TIP4P/Ice. The Journal of chemical physics, 122, p.234511. Agarwal, A., Ng, W.J. & Liu, Y., 2011. Principle and applications of microbubble and nanobubble technology for water treatment. Chemosphere, 84(9), pp.1175–1180. Aines, R.D. & Rossman, G.R., 1984. Water in minerals? A peak in the infrared. Journal of Geophysical Research, 89(B6), pp.4059–4071. Alavi, S. & Ripmeester, J.A., 2010. Nonequilibrium adiabatic molecular dynamics simulations of methane clathrate hydrate decomposition. The Journal of Chemical Physics, 132, p.144703. Anderson, B.J. et al., 2005. Properties of inhibitors of methane hydrate formation via molecular dynamics simulations. Journal of the American Chemical Society, 127(50), pp.17852–17862. Anderson, B. & Jordon, W., 2009. Molecular Level Modeling in Gas Hydrate Studies, The National Energy and Technology Laboratory. Available at: http://www.netl.doe.gov/technologies/oil-gas/FutureSupply/MethaneHydrates/newsletter/newsletter.htm. Anderson, G.K., 2004. Enthalpy of dissociation and hydration number of methane hydrate from the Clapeyron equation. The Journal of Chemical Thermodynamics, 36(12), pp.1119–1127. Anon, International Energy Statistics - EIA. Available at: http://www.eia.gov/cfapps/ipdbproject/IEDIndex3.cfm?tid=3&pid=26&aid=2 [Accessed April 29, 2014]. Argyris, D. et al., 2008. Molecular structure and dynamics in thin water films at the silica and graphite surfaces. The Journal of Physical Chemistry C, 112(35), pp.13587–13599. Baardsnes, J. et al., 1999. New ice-binding face for type I antifreeze protein. FEBS letters, 463(1), pp.87–91. Baez, L.A. & Clancy, P., 1994. Computer Simulation of the Crystal Growth and Dissolution of Natural Gas Hydratesa. Annals of the New York Academy of Sciences, 715(Natural Gas Hydrates), pp.177–186. Bagherzadeh, S.A. et al., 2013. Evolution of methane during gas hydrate dissociation. Fluid Phase Equilibria, 358, pp.114–120. Bagherzadeh, S.A. et al., 2011. Magnetic Resonance Imaging of Gas Hydrate Formation in a Bed of Silica Sand Particles. Energy & Fuels, 25(7), pp.3083–3092. Bagherzadeh, S.A. et al., 2012. Molecular Modeling of the Dissociation of Methane Hydrate in Contact with a Silica Surface. The Journal of Physical Chemistry B, 116(10), pp.3188–3197. Bagherzadeh, S.A. et al., 2012. Molecular simulation of non-equilibrium methane hydrate decomposition process. The Journal of Chemical Thermodynamics, 44(1), pp.13–19. 98 Bai, D. et al., 2011. Microsecond Molecular Dynamics Simulations of the Kinetic Pathways of Gas Hydrate Formation from Solid Surfaces. Langmuir. Baldwin, B.A. et al., 2009. Using magnetic resonance imaging to monitor CH4 hydrate formation and spontaneous conversion of CH4 hydrate to CO2 hydrate in porous media. Magnetic resonance imaging, 27(5), pp.720–726. Boewer, L. et al., 2012. On the Spontaneous Formation of Clathrate Hydrates at Water–Guest Interfaces. J. Phys. Chem. C, 116(15), pp.8548–8553. Bowers, P.G., Bar-Eli, K. & Noyes, R.M., 1996. Unstable supersaturated solutions of gases in liquids and nucleation theory. J. Chem. Soc., Faraday Trans., 92(16), pp.2843–2849. Carroll, J.J., Slupsky, J.D. & Mather, A.E., 1991. The solubility of carbon dioxide in water at low pressure. J. Phys. Chem. Ref. Data, 20(6), pp.1201–1209. Chen, F. & Smith, P.E., 2007. Simulated surface tensions of common water models. The Journal of chemical physics, 126, p.221101. Chou, K.-C., 1992. Energy-optimized structure of antifreeze protein and its binding mechanism. Journal of molecular biology, 223(2), pp.509–517. Christiansen, R.L. & Sloan, E.D., 1994. Mechanisms and kinetics of hydrate formation. Annals of the New York Academy of Sciences, 715(1), pp.283–305. Clennell, M.B. et al., 1999. Formation of natural gas hydrates in marine sediments 1. Conceptual model of gas hydrate growth conditioned by host sediment properties. Journal of Geophysical Research, 104(B10), pp.22985–23. Collett, T.S., 2002. Energy resource potential of natural gas hydrates. AAPG bulletin, 86(11), p.1971. Committee on Assessment of the Department of Energy, Methane Hydrate Research and Development Program: Evaluating Methane Hydrate as a Future Energy Resource, 2010. Realizing the Energy Potential of Methane Hydrate for the United States, Washington, D.C.: The National Academies Press. Conde, M. & Vega, C., 2010. Determining the three-phase coexistence line in methane hydrates using computer simulations. The Journal of chemical physics, 133, p.064507. Cornell, W.D. et al., 1995. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. Journal of the American Chemical Society, 117(19), pp.5179–5197. Dalal, P. et al., 2001. Hydrogen bond analysis of Type 1 antifreeze protein in water and the ice/water interface. PhysChemComm, 4(7), pp.32–36. Daraboina, N. et al., 2011. Natural gas hydrate formation and decomposition in the presence of kinetic inhibitors. 2. Stirred reactor experiments. Energy & Fuels, 25(10), pp.4384–4391. 99 Darden, T. et al., 1999. New tricks for modelers from the crystallography toolkit: the particle mesh Ewald algorithm and its use in nucleic acid simulations. Structure, 7(3), pp.R55–R60. Davidson, D.W., 1973. Clathrate Hydrates. In F. Franks, ed. Water: A Comprehensive Treatise. New York: Plenum, pp. 115–234. Davies, P.L. et al., 2002. Structure and function of antifreeze proteins. Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences, 357(1423), pp.927–935. Davy, H., 1811. On a Combination of Oxymuriatic Acid and Oxygen Gas. Phil Trans, 101, pp.155–162. Debenedetti, P.G. & Sarupria, S., 2009. Hydrate molecular ballet. Science, 326(5956), pp.1070–1071. Ding, L.Y. et al., 2007. Molecular dynamics simulation on the dissociation process of methane hydrates. Molecular Simulation, 33(12), pp.1005–1016. Englezos, P., 1993. Clathrate hydrates. Industrial & Engineering Chemistry Research, 32(7), pp.1251–1274. English, N.J., Johnson, J.K. & Taylor, C.E., 2005. Molecular-dynamics simulations of methane hydrate dissociation. The Journal of chemical physics, 123, p.244503. English, N.J. & MacElroy, J.M.D., 2004. Theoretical studies of the kinetics of methane hydrate crystallization in external electromagnetic fields. The Journal of chemical physics, 120, p.10247. English, N.J. & Phelan, G.M., 2009. Molecular dynamics study of thermal-driven methane hydrate dissociation. The Journal of chemical physics, 131, p.074704. Ewald, P.P., 1921. Die Berechnung optischer und elektrostatischer Gitterpotentiale. Annalen der Physik, 369(3), pp.253–287. Fan, M. et al., 2010a. Nanobubble generation and its application in froth flotation (part I): nanobubble generation and its effects on properties of microbubble and millimeter scale bubble solutions. Mining Science and Technology (China), 20(1), pp.1–19. Fan, M. et al., 2010b. Nanobubble generation and its applications in froth flotation (part II): fundamental study and theoretical analysis. Mining Science and Technology (China), 20(2), pp.159–177. Forrisdahl, O.K., Kvamme, B. & Haymet, A.D.J., 1996. Methane clathrate hydrates: melting, supercooling and phase separation from molecular dynamics computer simulations. Molecular Physics, 89(3), pp.819–834. Frenkel, D. & Smit, B., 2001. Understanding molecular simulation: from algorithms to applications 2nd edition., San Diego, USA: Academic press. 100 Galamba, N., 2014. Water Tetrahedrons, Hydrogen-Bond Dynamics, and the Orientational Mobility of Water around Hydrophobic Solutes. The Journal of Physical Chemistry B, 118(15), pp.4169–4176. Gale, P.A. & Steed, J.W. eds., 2012. Supramolecular Chemistry: From Molecules to Nanomaterials, Goldman, S., 2009. Generalizations of the Young–Laplace equation for the pressure of a mechanically stable gas bubble in a soft elastic material. The Journal of chemical physics, 131(18), p.184502. Gudmundsson, J.S., 1996. Method for production of gas hydrates for transportation and storage. Guillot, B. & Guissani, Y., 1993. A computer simulation study of the temperature dependence of the hydrophobic hydration. The Journal of chemical physics, 99(10), pp.8075–8094. Gupta, A. et al., 2008. Measurements of methane hydrate heat of dissociation using high pressure differential scanning calorimetry. Chemical Engineering Science, 63(24), pp.5848–5853. Haile, J.M., 1997. Molecular Dynamics Simulation: Elementary Methods 1st edition., New York: Wiley-Interscience. Hammerschmidt, E.G., 1934. Formation of gas hydrates in natural gas transmission lines. Industrial & Engineering Chemistry, 26(8), pp.851–855. Hampton, M.A. & Nguyen, A.V., 2010. Nanobubbles and the nanobubble bridging capillary force. Advances in colloid and interface science, 154(1), pp.30–55. Handa, Y.P., 1986. Compositions, enthalpies of dissociation, and heat capacities in the range 85 to 270 K for clathrate hydrates of methane, ethane, and propane, and enthalpy of dissociation of isobutane hydrate, as determined by a heat-flow calorimeter. The Journal of Chemical Thermodynamics, 18(10), pp.915–921. Handa, Y.P. & Stupin, D.Y., 1992. Thermodynamic properties and dissociation characteristics of methane and propane hydrates in 70-. ANG.-radius silica gel pores. The Journal of Physical Chemistry, 96(21), pp.8599–8603. Hao, W. et al., 2008. Evaluation and analysis method for natural gas hydrate storage and transportation processes. Energy conversion and management, 49(10), pp.2546–2553. Harris, J.G. & Yung, K.H., 1995. Carbon dioxide’s liquid-vapor coexistence curve and critical properties as predicted by a simple molecular model. The Journal of Physical Chemistry, 99(31), pp.12021–12024. Hawtin, R.W., Quigley, D. & Rodger, P.M., 2008. Gas hydrate nucleation and cage formation at a water/methane interface. Phys. Chem. Chem. Phys., 10(32), pp.4853–4864. Haymet, A.D.J. et al., 1998. Valine substituted winter flounderantifreeze’: preservation of ice growth hysteresis. FEBS letters, 430(3), pp.301–306. 101 Henry, P., Thomas, M. & Clennell, M.B., 1999. Formation of natural gas hydrates in marine sediments, 2, Thermodynamic calculations of stability conditions in porous sediments. JOURNAL OF GEOPHYSICAL RESEARCH-ALL SERIES-, 104, pp.23–23. Hess, B. et al., 2008. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. Journal of chemical theory and computation, 4(3), pp.435–447. Hester, K.C. & Brewer, P.G., 2009. Clathrate hydrates in nature. Annual Review of Marine Science, 1, pp.303–327. Holder, G.D., Kamath, V.A. & Godbole, S.P., 1984. The Potential of Natural Gas Hydrates as an Energy Resource. Annual Review of Energy, 9(1), pp.427–445. Hong, H. & Pooladi-Darvish, M., 2005. Simulation of depressurization for gas production from gas hydrate reservoirs. Journal of Canadian Petroleum Technology, 44(11). Hoover, W.G., 1985. Canonical dynamics: equilibrium phase-space distributions. Physical Review A, 31(3), p.1695. Ishida, N. et al., 2000. Nano Bubbles on a Hydrophobic Surface in Water Observed by Tapping-Mode Atomic Force Microscopy. Langmuir, 16(16), pp.6377–6380. Iwai, Y. et al., 2010. Analysis of dissociation process for gas hydrates by molecular dynamics simulation. Molecular Simulation, 36(3), pp.246–253. Kelland, M.A., 2006. History of the development of low dosage hydrate inhibitors. Energy & Fuels, 20(3), pp.825–847. Knight, C.A., Cheng, C.C. & DeVries, A.L., 1991. Adsorption of alpha-helical antifreeze peptides on specific ice crystal surface planes. Biophysical Journal, 59(2), pp.409–418. Knott, B.C. et al., 2012. Homogeneous Nucleation of Methane Hydrates: Unrealistic under Realistic Conditions. Journal of the American Chemical Society, 134(48), pp.19544–19547. Konrad, J.-M., 1990. Unfrozen water as a function of void ratio in a clayey silt. Cold Regions Science and Technology, 18(1), pp.49–55. Kumar, R., Englezos, P., et al., 2009. Structure and composition of CO2/H2 and CO2/H2/C3H8 hydrate in relation to simultaneous CO2 capture and H2 production. AIChE Journal, 55(6), pp.1584–1594. Kumar, R., Linga, P., et al., 2009. Two-Stage Clathrate Hydrate/Membrane Process for Precombustion Capture of Carbon Dioxide and Hydrogen. Journal of Environmental Engineering, 135(6), pp.411–417. Leach, A., 2001. Molecular Modelling: Principles and Applications 2nd edition., Harlow, England: Prentice Hall. 102 Lehmkühler, F. et al., 2008. The Carbon Dioxide- Water Interface at Conditions of Gas Hydrate Formation. Journal of the American Chemical Society, 131(2), pp.585–589. Levien, L., Prewitt, C.T. & Weidner, D.J., 1980. Structure and elastic properties of quartz at pressure. American Mineralogist, 65(9-10), p.920. Liang, S. & Kusalik, P.G., 2010. Explorations of gas hydrate crystal growth by molecular simulations. Chemical Physics Letters, 494(4), pp.123–133. Liang, S. & Kusalik, P.G., 2011. Exploring nucleation of H2S hydrates. Chemical Science, 2(7), pp.1286–1292. Liang, S. & Kusalik, P.G., 2013. Nucleation of Gas Hydrates within Constant Energy Systems. The Journal of Physical Chemistry B, 117(5), pp.1403–1410. Liang, S., Rozmanov, D. & Kusalik, P.G., 2011. Crystal growth simulations of methane hydrates in the presence of silica surfaces. Phys. Chem. Chem. Phys., 13(44), pp.19856–19864. Lindorff-Larsen, K. et al., 2010. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins: Structure, Function, and Bioinformatics, 78(8), pp.1950–1958. Linga, P. et al., 2009. Gas Hydrate Formation in a Variable Volume Bed of Silica Sand Particles. Energy & Fuels, 23(11), pp.5496–5507. Linga, P., Kumar, R. & Englezos, P., 2007. The clathrate hydrate process for post and pre-combustion capture of carbon dioxide. Journal of Hazardous Materials, 149(3), pp.625–629. Lopes, P.E.. et al., 2006. Development of an empirical force field for silica. Application to the quartz-water interface. The Journal of Physical Chemistry B, 110(6), pp.2782–2792. Lou, S.-T. et al., 2000. Nanobubbles on solid surface imaged by atomic force microscopy. Journal of Vacuum Science & Technology B, 18(5), pp.2573–2575. Lubetkin, S.D., 2003. Why is it much easier to nucleate gas bubbles than theory predicts? Langmuir, 19(7), pp.2575–2587. Luzar, A. & Chandler, D., 1996. Hydrogen-bond kinetics in liquid water. Nature, 379(6560), pp.55–57. MacDonald, G.J., 1990. The Future of Methane as an Energy Resource. Annual Review of Energy, 15(1), pp.53–83. Makogon, Y.F. et al., 1971. Finding of a Pool of Gas in the Hydrate State. Moscow, DAN SSSR, 196(1), pp.197–206. Mao, Y. & Ba, Y., 2006. Insight into the Binding of Antifreeze Proteins to Ice Surfaces via 13C Spin Lattice Relaxation Solid-State NMR. Biophysical journal, 91(3), pp.1059–1068. 103 Max, M.D. & Pellenbarg, R.E., 1999. Desalination through methane hydrate, Google Patents. Available at: http://www.google.com/patents/US5873262 [Accessed November 28, 2014]. McMullan, R.K. & Jeffrey, G.A., 1965. Polyhedral Clathrate Hydrates. IX. Structure of Ethylene Oxide Hydrate. The Journal of Chemical Physics, 42(8), pp.2725–2732. Milkov, A.V. et al., 2003. In situ methane concentrations at Hydrate Ridge, offshore Oregon: New constraints on the global gas hydrate inventory from an active margin. Geology, 31(10), pp.833–836. Moon, C., Taylor, P.C. & Rodger, P.M., 2003. Clathrate nucleation and inhibition from a molecular perspective. Canadian journal of physics, 81(1-2), pp.451–457. Moon, C., Taylor, P.C. & Rodger, P.M., 2003. Molecular dynamics study of gas hydrate formation. Journal of the American Chemical Society, 125(16), pp.4706–4707. Moridis, G.J., 2003. Numerical studies of gas production from methane hydrates. Spe Journal, 8(4), pp.359–370. Myshakin, E.M. et al., 2009. Molecular Dynamics Simulations of Methane Hydrate Decomposition†. The Journal of Physical Chemistry A, 113(10), pp.1913–1921. Nada, H., 2006. Growth mechanism of a gas clathrate hydrate from a dilute aqueous gas solution: A molecular dynamics simulation of a three-phase system. The Journal of Physical Chemistry B, 110(33), pp.16526–16534. Nielsen, L.C., Bourg, I.C. & Sposito, G., 2012. Predicting CO2–water interfacial tension under pressure and temperature conditions of geologic CO2 storage. Geochimica et Cosmochimica Acta, 81, pp.28–38. Nosé, S., 1984. A unified formulation of the constant temperature molecular dynamics methods. The Journal of chemical physics, 81(1), pp.511–519. Nutt, D.R. & Smith, J.C., 2008. Dual function of the hydration layer around an antifreeze protein revealed by atomistic molecular dynamics simulations. Journal of the American Chemical Society, 130(39), pp.13066–13073. Ohgaki, K. et al., 1996. Methane Exploitation by Carbon Dioxide from Gas Hydrates—Phase Equilibria for CO2-CH4 Mixed Hydrate System. Journal of chemical engineering of Japan, 29(3), pp.478–483. Ohgaki, K. et al., 2010. Physicochemical approach to nanobubble solutions. Chemical Engineering Science, 65(3), pp.1296–1300. Ohno, H. et al., 2010. Interaction of antifreeze proteins with hydrocarbon hydrates. Chemistry-a European Journal, 16(34), pp.10409–10417. Park, Y. et al., 2006. Sequestering carbon dioxide into complex structures of naturally occurring gas hydrates. Proceedings of the National Academy of Sciences, 103(34), p.12690. 104 Perfeldt, C.M. et al., 2014. Inhibition of Gas Hydrate Nucleation and Growth: Efficacy of an Antifreeze Protein from the Longhorn Beetle Rhagium mordax. Energy & Fuels, 28(6), pp.3666–3672. Ping, Z.H. et al., 2001. States of water in different hydrophilic polymers—DSC and FTIR studies. Polymer, 42(20), pp.8461–8467. Ramachandran, G.N., Ramakrishnan, C. & Sasisekharan, V., 1963. Stereochemistry of polypeptide chain configurations. Journal of Molecular Biology, 7, pp.95–99. Razul, M.S.. et al., 2005. Computer simulations of heterogeneous crystal growth of atomic systems. Molecular Physics, 103(14), pp.1929–1943. Rebata-Landa, V. & Santamarina, J.C., 2012. Mechanical effects of biogenic nitrogen gas bubbles in soils. Journal of Geotechnical and Geoenvironmental Engineering, 138(2), pp.128–137. Riedel, M. et al., 2006. Gas hydrate transect across Northern Cascadia Margin. Eos, Transactions American Geophysical Union, 87(33), pp.325–332. Ripmeester, J.A., 2000. Hydrate Research-From Correlations to a Knowledge-based Discipline: The Importance of Structure. Annals of the New York Academy of Sciences, 912(1), pp.1–16. Ripmeester, J.A. & Alavi, S., 2010. Molecular Simulations of Methane Hydrate Nucleation. ChemPhysChem, 11(5), pp.978–980. Rodger, P., 2000. Methane hydrate: melting and memory. Annals of the New York Academy of Sciences, 912(1), pp.474–482. Rodger, P.M., 1994. Towards a Microscopic Understanding of Clathrate Hydrates. Annals of the New York Academy of Sciences, 715(1), pp.207–211. Rodger, P.M., Forester, T.R. & Smith, W., 1996. Simulations of the methane hydrate/methane gas interface near hydrate forming conditions conditions. Fluid phase equilibria, 116(1-2), pp.326–332. Ryckaert, J.-P., Ciccotti, G. & Berendsen, H.J., 1977. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. Journal of Computational Physics, 23(3), pp.327–341. Sarupria, S. & Debenedetti, P.G., 2011. Molecular dynamics study of carbon dioxide hydrate dissociation. The Journal of Physical Chemistry A, 115(23), pp.6102–6111. Sarupria, S. & Debenedetti, P.G., 2011. Molecular Dynamics Study of Carbon Dioxide Hydrate Dissociation. The Journal of Physical Chemistry A. Scharlin, P. & Battino, R., 1995. Solubility of CCl2F2, CClF3, CF4, and CH4 in Water and Seawater at 288.15-303.15 K and 101.325 kPa. Journal of Chemical & Engineering Data, 40(1), pp.167–169. 105 Seddon, J.R. et al., 2012. A Deliberation on Nanobubbles at Surfaces and in Bulk. ChemPhysChem, 13(8), pp.2179–2187. Seo, Y.T. et al., 2005. Efficient recovery of CO2 from flue gas by clathrate hydrate formation in porous silica gels. Environ. Sci. Technol, 39(7), pp.2315–2319. Seshadri, K., Wilder, J.W. & Smith, D.H., 2001. Measurements of equilibrium pressures and temperatures for propane hydrate in silica gels with different pore-size distributions. J. Phys. Chem. B, 105(13), pp.2627–2631. Sharifi, H. et al., 2014. Insights into the Behavior of Biological Clathrate Hydrate Inhibitors in Aqueous Saline Solutions. Crystal Growth & Design, 14(6), pp.2923–2930. Shi, B., 2006. Molecular Dynamics Simulation of the Surface Tension and Contact Angle of Argon and Water, University of California, Los Angeles. Shvab, I. & Sadus, R.J., 2014. Thermodynamic properties and diffusion of water + methane binary mixtures. The Journal of Chemical Physics, 140(10), p.104505. Sicheri, F. & Yang, D.S., 1995. Ice-binding structure and mechanism of an antifreeze protein from winter flounder. Nature, 375, pp.427–431. Sklodowska, A. & Matlakowska, R., 1997. Influence of exopolymers produced by bacterial cells on hydrophobicity of substrate surface. Biotechnology techniques, 11(11), pp.837–840. Sloan, E.D., 2005. A changing hydrate paradigm—from apprehension to avoidance to risk management. Fluid Phase Equilibria, 228, pp.67–74. Sloan, E.D., 2003. Fundamental principles and applications of natural gas hydrates. Nature, 426(6964), pp.353–363. Sloan, E.D. & Koh, C.A., 2008. Clathrate Hydrates of Natural Gases 3rd ed., Boca Raton, Fl.: CRC Press. Smith, W. & Forester, T.R., 1996. DL_POLY_2. 0: A general-purpose parallel molecular dynamics simulation package. Journal of Molecular Graphics, 14(3), pp.136–141. Storr, M.T. et al., 2004. Kinetic inhibitor of hydrate crystallization. Journal of the American Chemical Society, 126(5), pp.1569–1576. Tabuchi, H., 2013. An Energy Coup for Japan: “Flammable Ice.” The New York Times. Available at: http://www.nytimes.com/2013/03/13/business/global/japan-says-it-is-first-to-tap-methane-hydrate-deposit.html [Accessed April 15, 2013]. The Council of Canadian Academies, 2008. Energy from Gas Hydrates: Assessing the Opportunities and Challenges for Canada, Canada. Available at: http://www.scienceadvice.ca/en/assessments/completed/gas-hydrates.aspx. 106 Totland, C. et al., 2011. Water Structure and Dynamics at a Silica Surface: Pake Doublets in 1H NMR Spectra. Langmuir, 27(8), pp.4690–4699. Tse, J.S., Klein, M.L. & McDonald, I.R., 1983. Molecular dynamics studies of ice Ic and the structure I clathrate hydrate of methane. The Journal of Physical Chemistry, 87(21), pp.4198–4203. Tse, J.S. & Klug, D.D., 2002. Formation and decomposition mechanisms for clathrate hydrates. Journal of Supramolecular Chemistry, 2(4-5), pp.467–472. Tung, Y.T. et al., 2011. Growth of Structure I Carbon Dioxide Hydrate from Molecular Dynamics Simulations. The Journal of Physical Chemistry C. Tung, Y.-T. et al., 2012. Molecular Dynamics Study on the Growth of Structure I Methane Hydrate in Aqueous Solution of Sodium Chloride. The Journal of Physical Chemistry B, 116(48), pp.14115–14125. Tung, Y.-T. et al., 2010. The growth of structure I methane hydrate from molecular dynamics simulations. The Journal of Physical Chemistry B, 114(33), pp.10804–10813. Uchida, T. et al., 2002. Effects of pore sizes on dissociation temperatures and pressures of methane, carbon dioxide, and propane hydrates in porous media. J. Phys. Chem. B, 106(4), pp.820–826. Uchida, T. et al., 2011. Transmission electron microscopic observations of nanobubbles and their capture of impurities in wastewater. Nanoscale research letters, 6(1), pp.1–9. Udachin, K.A., Ratcliffe, C.I. & Ripmeester, J.A., 2001. Structure, composition, and thermal expansion of CO2 hydrate from single crystal X-ray diffraction measurements. The Journal of Physical Chemistry B, 105(19), pp.4200–4204. Uddin, M. & Coombe, D., 2014. Kinetics of CH4 and CO2 Hydrate Dissociation and Gas Bubble Evolution via MD Simulation. The Journal of Physical Chemistry A. Available at: http://pubs.acs.org/doi/abs/10.1021/jp410789j [Accessed May 2, 2014]. Vatamanu, J. & Kusalik, P.G., 2007. Molecular dynamics methodology to investigate steady-state heterogeneous crystal growth. Journal of Chemical Physics, 126(12), pp.124703–124703. Vatamanu, J. & Kusalik, P.G., 2006. Molecular insights into the heterogeneous crystal growth of sI methane hydrate. The Journal of Physical Chemistry B, 110(32), pp.15896–15904. Vatamanu, J. & Kusalik, P.G., 2010. Observation of two-step nucleation in methane hydrates. Phys. Chem. Chem. Phys., 12(45), pp.15065–15072. Verlet, L., 1967. Computer“ experiments” on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Physical review, 159(1), p.98. Walsh, M.R. et al., 2009. Microsecond simulations of spontaneous methane hydrate nucleation and growth. Science, 326(5956), p.1095. 107 Weijs, J.H., Seddon, J.R. & Lohse, D., 2012. Diffusive shielding stabilizes bulk nanobubble clusters. ChemPhysChem, 13(8), pp.2197–2204. Wierzbicki, A. et al., 2007. Antifreeze proteins at the ice/water interface: three calculated discriminating properties for orientation of type I proteins. Biophysical journal, 93(5), pp.1442–1451. Wilder, J.W. et al., 2001. Modeling hydrate formation in media with broad pore size distributions. Langmuir, 17(21), pp.6729–6735. Yagasaki, T. et al., 2014. Effect of Bubble Formation on the Dissociation of Methane Hydrate in Water: A Molecular Dynamics Study. The Journal of Physical Chemistry B, 118(7), pp.1900–1906. Yasuoka, K. & Murakoshi, S., 2000. Molecular Dynamics Simulation of Dissociation Process for Methane Hydrate. Annals of the New York Academy of Sciences, 912(1), pp.678–684. Zatsepina, O. & Pooladi-Darvish, M., 2010. CO2 Storage as Hydrate in Depleted Gas Reservoirs. In Canadian Unconventional Resources and International Petroleum Conference. Zeng, H. et al., 2006. Effect of Antifreeze Proteins on the Nucleation, Growth, and the Memory Effect during Tetrahydrofuran Clathrate Hydrate Formation. Journal of the American Chemical Society, 128(9), pp.2844–2850. Zhang, W. & Laursen, R.A., 1998. Structure-function relationships in a type I antifreeze polypeptide the role of threonine methyl and hydroxyl groups in antifreeze activity. Journal of Biological Chemistry, 273(52), pp.34806–34812. Zhang, X.H., Khan, A. & Ducker, W.A., 2007. A nanoscale gas state. Physical review letters, 98(13), p.136101.  AppendAppendiThe unit Databasefrom the propertieTable  A-1HexagonOrthorh  Figure  A-simulatioFigure  A-plane is hices x A cell for silica.(Levien et ahexagonal us of the orth. Unit cell pal unit cell ombic simul1. The originn cell (yellow 2. The hydroydroxylated was downlol. 1980) An onit cell (paralorhombic unarameters foration cell al hexagona rectangle).xylated cub. aded from thrthorhombilelogram) ofit cell are giv the originalܽ = ܾܽ = 4.916 Հl unit cell (p Si (blue), O ic unit cell oe Americanc simulation  silica. The uen in Table  and the new= 4.916 Հ; ܿ =; ܾ = 8.5148arallelogram(red). f silica. Si (bl Mineralogiscell was chonit cell para A-1 and sho silica unit c5.4054 Հ Հ; ܿ = 5.4054) and the neue), O (red)t Crystal Strsen (yellow rmeters for siwn in Figurell ߙ = ߚՀ ߙ =w selected o, H (white). ucture ectangle) stalica and e  A-1. = 90°; ߛ = 12ߚ = ߛ = 90°rthorhombicThe (010) ba108 rting 0°   sal 109 All the silica atoms on the surface were hydroxylated using a bond length of 1.1 Å for the O-H bond, according to Lopes et al.(Lopes et al. 2006) A 8×1×14 supercell of this new unit cell was used for the construction of silica slabs. The simulation box contains two silica slabs on the xz-plane (38.6 Å apart in the y-direction). The atoms of the silica slabs were kept frozen during the simulations.   Appendi Figure  B-gas) and  Figure  B-concentra      x B 1. The initia5 (water and2. Snapshotstion enhancl setup confi carbon diox of CH4 andement near tgurations foride gas).  CO2 gas dishe interface  the simulattribution adjis evident. ion 3 of Tablacent to a flae  3-2 (water t water surfa and methance. The gas110 e   111 Appendix C The F3 order parameters for the SHW and SWHW scenarios along the y-axis are plotted in Figure  C-1. As expected, the outer layers decompose faster. Along the y-axis, all the layers start to dissociate almost at the same time, although at different rates. This is due to the fact that the outer edges of the hydrate along the z-direction begin decomposing immediately. Decomposition along the y-axis also shows a layer-by-layer nature.  Figure  C-1. The F3 order parameter along the layers in the y-direction, see Figure  4-1.     AppendiFigure  D(G-W-H- x D -1. zz-compoW-G) at 350 nent of the pK and 400 baressure tensr. Red line ror during meepresents ththane hydrae 20-points m te decompooving averasition simulage trend line112 tion . Figure  Dfraction o   -2. zz-compof 0.042 at 350nent of the p K and 400 b ressure tensar. Red lineor during th represents the simulatione 20-points   with dissolvmoving avered methane age trend lin113 mole e. AppendiThe methexamine t(TIP4P-ic800 bar athe hydraA 2×2×2 methane 0.24 and Figure  E-simulatiowhite (hyalso showsimulatioA 0.2 ns Nkept frozperformeOther runIn all of tacceptabl2 fs, it wi x E od of Condhe effect of e for water and 305 K to te stable regireplica of thmolecules. T0.665 nm, re1. Initial (0 nn is 800 bar adrogen) bonn by red dasn box at t=0VT simulaten, are perfod in an aniso parametershe three runy within a coll not significe and Vega (time step. Tnd united atmake sure thon. e unit cell ofhe size of thspectively, ans) and final nd hydrate ids. The hydrhed lines an ns along theion followedrmed for temtropic NPT  are set simils, hydrate grmparable raantly affect Conde & Vehe equilibriuom model fat 800 bar a sI methane e resulting sd is shown (50 ns) confis grown at 3ogen bondind methane m x, y and z d by a 0.2 ns perature anensemble foar to those gow and the pnge. This cothe results oga 2010) is fm temperatuor methane) nd 280 K (chydrate is plimulation boin Figure  E-gurations of05 K. Water g between wolecules areirections areNPT simulad pressure er 50 ns and tiven in Chapotential enenfirms that af the simulatollowed to ore at 400 bais 302 K. A onditions usaced in contx along the x1.  the system. molecules arater molecu shown by c 0.24, 0.24 antion, during quilibration. ime steps ofter 6.  rgy profiles, s long as theions performbserve hydrr for the forcsimulation ised in Chapteact with 368, y and z dir The pressure shown by les in the hyyan spheres.d 0.665 nm. which the hyThe produc 0.5, 1.0 andshown in Fig time step ised in this thate growth ae field appli performed r 6) fall with water and 6ections are 0e of the red (oxygen)drate phase  The size of drate phasetion runs are 2.0 fs are teure F-2, fall between 0.5esis. 114 nd ed at in 4 .24,  and is  is  sted.   and 115  Figure  E-2. Potential energy profile for hydrate growth simulation at 305 K and 800 bar with different time steps. The energy trends show similar behavior and the end state at 50 ns is essentially the same.    116 Appendix F Table  F-1. Heavy atoms of the side chain of all residues of wf-AFP around which the F3 is calculated. The wf-AFP molecular structure is given in Figure  2-5. Residue Side chain Atom Residue Side  chain Atom 1ASP carboxyl OD1 20ALA methyl CB OD2 21ALA methyl CB 2THR methyl CG2 22GLU carboxyl OE1 hydroxyl OG1 OE2 3ALA methyl CB 23LEU methyl CD1 4SER hydroxyl OG methyl CD2 5ASP carboxyl OD1 24THR methyl CG2 OD2 hydroxyl OG1 6ALA methyl CB 25ALA methyl CB 7ALA methyl CB 26ALA methyl CB 8ALA methyl CB 27ASN carbonyl OD1 9ALA methyl CB amine ND2 10ALA methyl CB 28ALA methyl CB 11ALA methyl CB 29ALA methyl CB 12LEU methyl CD1 30ALA methyl CB CD2 31ALA methyl CB 13THR methyl CG2 32ALA methyl CB hydroxyl OG1 33ALA methyl CB 14ALA methyl CB 34ALA methyl CB 15ALA methyl CB 35THR methyl CG2 16ASN carbonyl OD1 hydroxyl OG1 amine ND2 36ALA methyl CB 17ALA methyl CB 37ARG amine NH1 18LYS ammonium NZ NH2 19ALA methyl CB   

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0166109/manifest

Comment

Related Items