Clathrate hydrates simulations

Oil and water do not mix. This we have learned from an early age and have been recently reminded of in the Gulf of Mexico oil spill in 2010. Oil and other non-polar compounds are thus said to be “hydrophobic.” Similarly, small hydrophobic, gaseous molecules, such as methane, have low solubility in water (typically on the order of one part in a thousand). It can strike us as counter-intuitive, then, that at low temperatures and/or high pressures, water and gas crystallize together into an “icy” solid with high gas composition (water to gas ratio of about 6:1). The explanation lies in the fact that hydrophobicity is not merely the tendency of dissimilar substances to demix, but also involves the tendency of water to form hydrogen-bond networks on a molecular scale. When the hydrophobic molecule is of the same size as the cells of the hydrogen-bond network, they can become trapped or “enclathrated” into a regular lattice of hydrogen-bond cages, forming a crystalline structure known as clathrate hydrate. (See the Rath article in the MRS Bulletin April 2008 issue.)

These clathrate hydrates are often incorrectly termed “methane ice” or “icy crystals” as in, for example, popular media coverage of the Gulf of Mexico oil spill. In fact, the formation of gas hydrates has long been the foremost problem to the oil and gas industry because of their tendency to plug pipelines. Large capital and financial resources are usually in place as remediation strategies to prevent, thermodynamically, or at least control kinetically the formation of hydrates in pipelines, such as the injection of inhibitors, insulation, and/or heating (to shift or avoid the phase boundary).

Beyond their role in pipeline plugs, clathrate hydrates have been brought into focus in several areas of energy science. Broadly speaking, one can identify one category of areas as being concerned when the energy is contained in the guest molecules, which are thus a fuel. This is the case for the vast reserves of naturally occurring methane hydrate in the seafloor or permafrost (amount of methane in these hydrate deposits are on the order of 700,000 TCF (trillion cubic feet);Reference Boswell¹ as a base, the U.S. annual natural gas consumption is ~23 TCFReference Boswell¹), where the methane is extracted and possibly replaced by its greenhouse gas byproduct, carbon dioxide.Reference Graue, Kvamme, Baldwin, Stevens, Howard, Aspenes, Ersland, Husebø and Zornes² It is also the case for exploratory strategies envisioning molecular hydrogen as a fuel, where the clathrate water lattice acts as a storage medium to reversibly capture, concentrate, store, and release hydrogen fuel.Reference Sugahara, Haag, Prasad, Warntjes, Sloan, Sum and Koh^3, Reference Sugahara, Haag, Warntjes, Prasad, Sloan, Koh and Sum⁴ More concrete is the storage and transportation of natural gas in clathrate hydrates, which is moving from laboratory to pilot scale to full industrial reality.Reference Nogami and Watanabe⁵ Another major category of areas focuses on the energy contained in the structuring of the clathrate itself (e.g., as heat). Examples include the use of hydrates to extract and release the heat of formation in a refrigeration cycle.Reference Ogawa, Ito, Watanabe, Tahara, Hiraoka, Ochiai, Ohmura and Mori^6, Reference Fournaison, Delahaye, Chatti and Petitet⁷

The broad relevance of clathrate hydrates reflects the ubiquity of the constituents (water and gas), as well as the simplicity of the conditions for formation (low temperature/high pressure). The simplicity of the formation conditions, however, belies the complex rearrangements of molecules and hydrogen bonds that must occur to create the hydrate lattice-guest structure. A fundamental understanding of this complex process would be helpful in the design of strategies to control hydrate formation in any of the areas of energy science mentioned previously. Specifically, there are three types of information that are required: (1) thermodynamics (i.e., knowledge of phase boundaries and equilibrium guest compositions and thus energy content), (2) kinetics, which is important in process design, whether it be for plug formation/removal, storage, or refrigeration, and (3) mechanism, knowledge of which helps guide the design of additives to affect rates or to understand dependencies of formation on phase variables and composition. For each of these areas, high-performance computing (HPC) can help to answer questions that are difficult to address experimentally. We give a few examples in the next sections for what can be gained from HPC and how it might be exploited in the design of materials.

Clathrate hydrates basics

The crystalline structure formed by clathrate hydrates is based on a hydrogen-bonded network of water molecules (host), forming regular polyhedra cages that typically enclose one gas molecule (guest). The name derives from the Latin clathratus, from clathri “lattice bars.” Figure 1 illustrates the cages formed and their connectivity to form the commonly found structure I (sI) and structure II (sII), both of which are cubic with unit cell dimensions of 12.0 Å and 17.2 Å, respectively.Reference Sloan and Koh⁸ The structure formed largely depends on the size of the guest molecules enclosed in the cages; these guest molecules are unbounded to the cage framework and are able to freely rotate and translate to some degree within the cage. Structure I is formed with the smallest guests, such as methane, ethane, or carbon dioxide, while sII can enclose larger guests, such as propane or i-butane. Interestingly, very small molecules, such as N₂ and H₂, also form sII hydrates; for these molecules, size alone would favor sI, but because sII has a larger number of the small cages and N₂ and H₂ can better stabilize the small cages, sII is the equilibrium structure formed. Hydrates formed with N₂ and H₂ are also unusual because unlike all other known hydrates, multiple molecules are known to occupy the large cages (up to two N₂ and up to four H₂ in the 5¹²6⁴ cages of sII). A unique property of clathrate hydrates is that they form non-stoichiometric compounds, that is, not all cages must be occupied by guest molecules, leaving some cages empty. However, there must be a minimum number of cages filled in order to stabilize the structure. The occupancy of the structure and distribution of guests (for multi-guest systems) among the cages depend on the formation temperature, pressure, and composition of the gas phase.

Figure 1.

Water molecules self-assemble around small hydrophobic molecules to form cavities that are the building blocks of clathrate hydrate. The right combination of cages leads to the crystal structure, the two most common being structure I and structure II, both of which are cubic. The stability of each structure is determined by the guest molecules occupying the cavities. The notation for the cages is as follows: 5¹², 12 pentagonal faces; 6², 2 hexagonal faces; and 6⁴, 4 hexagonal faces.

Understanding of the thermodynamics of clathrate hydrates is an area well established, where the formation temperature and pressure can be reasonably well predicted by the solid solution theory developed by van der Waals and Platteeuw (vdWP).Reference Sloan and Koh^8, Reference van der Waals and Platteeuw⁹ This model is an excellent example of how statistical thermodynamics can be elegantly applied to describe a complex system, which can be analogously viewed as a Langmuir adsorption problem (where guest molecules can be only placed in specific sites [i.e., cages]).

Thermodynamics from simulations

When performing stability and kinetic studies of clathrate hydrates, one essential piece of information that often is lacking is the hydrate stability boundary for the given molecular models used. In laboratory measurements, hydrate experiments are planned such that the formation conditions (temperature, T; pressure, P; and composition, x) are known to be within the hydrate stability region (e.g., as shown in Figure 2). This is certainly true for molecular simulations as well. Clearly, one must also know the stability boundary for clathrate hydrates to correctly determine the driving force and mechanism for hydrate formation. However, determination of the equilibrium phase boundary for clathrate hydrates from molecular simulations is a major challenge, and is thus, with a few exceptions, neglected despite its critical importance. Because every molecular model, whether for water or the guest species, will give a different phase boundary, this determination is typically missing and difficult to obtain. In the most common systems, the three-phase equilibrium of liquid water, hydrate, and gas vapor must be determined by satisfying the equality of chemical potential for each species (water and guest(s)) at a given temperature and pressure. These data must be rigorously obtained either by performing computationally intensive simulations of the phases in equilibrium (e.g., liquid water, hydrate, and vapor)Reference Conde and Vega¹⁰ or by thermodynamic integrationReference Jensen, Thomsen, von Solms, Wierzchowski, Walsh, Koh, Sloan, Wu and Sum¹¹ from known standard states. Due to the difficulty in obtaining such data, most studies published on the dynamics and stability of clathrate hydrates have assumed a temperature, pressure, and hydrate composition (typically at full occupancy) for the simulations. However, without a more physical basis for determining this information, one is left to make inconclusive hypotheses about the mechanism, since the thermodynamic state for the system has not been defined.

Figure 2.

Phase diagram for methane hydrates showing the liquid water (L_w) + hydrate (H) + methane vapor (V) boundary. The open circles correspond to experimental data;Reference Sloan and Koh⁸ the solid squaresReference Jensen, Thomsen, von Solms, Wierzchowski, Walsh, Koh, Sloan, Wu and Sum¹¹ and solid diamondsReference Conde and Vega¹⁰ correspond to points along the phase boundary obtained from molecular simulations using molecular models for water (TIP4P/ice) and methane (united atom LJ).

One of the most important data from thermodynamics studies on clathrate hydrates is the occupancy of the hydrate. Because hydrates are non-stoichiometric compounds, the filling of the cages depends on the formation conditions (T, P, and x). The composition of hydrates is not readily obtained from experimental measurements, as it is difficult to directly measure the composition of the guest species in the hydrate lattice. The composition and occupancy of the hydrates is particularly important in the energy storage of hydrates and has a direct impact on the amount of gas stored in the cages of the crystal structure. For example, in order to maximize the amount of H₂ stored in hydrates, it is critically important to know which cages H₂ is trapped inside and how many H₂ molecules there are per cage.

Clathrate hydrates metastability

For almost two decades, since the first molecular simulation study of clathrate hydrate nucleation was attempted,Reference Baez and Clancy¹² the major difficulty in achieving a complete description of the molecular events leading to the crystallization of a hydrate structure from a disordered state has been to overcome the metastability. Specifically, either liquid water or ice in contact with a gaseous vapor can persist for a stochastic and possibly extremely long time before transitioning to the ordered hydrate state. This metastability, which is also commonly encountered in laboratory measurements, is only attenuated in molecular simulations, since there are far greater limitations in simulating the nucleation process, such as limited system size (tens of nanometers) and timescale (several microseconds). Nonetheless, simulations are invaluable and provide an insightful molecular level description of clathrate hydrate nucleation that is unattainable by laboratory measurements. The metastability of the system, reflected as the rare event in forming a critical nucleus, invariably introduces daunting challenges in the simulation of hydrate nucleation.

Figure 3 illustrates the process for hydrate nucleation from a melt state. A metastable state of gas and water must overcome a significant free energy barrier for spontaneous nucleation to occur. Furthermore, the low solubility of hydrocarbon gases in water presents an extra layer of difficulty to the investigation of hydrate nucleation. Therefore, for a mixture consisting of liquid water and methane vapor, the most probable onset of hydrate nucleation is at or near the vapor-liquid interface, further complicating the investigation of nucleation by potentially introducing interfacial and diffusive transport phenomena. Finally, once initiated, hydrate nucleation can take place in nanoseconds, a time scale not directly accessible experimentally.

Figure 3.

Schematic of a free energy landscape for clathrate hydrates nucleation from gas and water. A transition state, associated with a high free energy barrier, separates the melt and crystal states.

The period of time that the system exists in a metastable state corresponds to the induction time which, for hydrates, depends on the driving force (T, P, and x). Experimentally, this induction time can be observed to last for hours to days, depending on the setup and system conditions. Therefore, the induction time for hydrate nucleation poses a severe challenge if the simulations are to be performed at similar conditions to those in laboratory measurements. As evidenced from our own recent studies, a successful hydrate nucleation was only captured after over a microsecond of simulation time,Reference Walsh, Koh, Sloan, Sum and Wu¹³ a timescale that is about two orders of magnitude longer than any other hydrate simulation reported previously. Moreover, microsecond simulations for any type of system are still uncommon due to computational demands in total time required to complete the simulations. However, long time simulations may be an effective route to obtaining molecular level detail for hydrate nucleation (the ice nucleation simulation by Matsumoto et al.Reference Matsumoto, Saito and Ohmine¹⁴ required on the order of 400 ns).

Clathrate hydrates nucleation

Even though a molecular level mechanism for clathrate hydrate formation has been sought for decades, it was not until recently that a breakthrough was achieved that gave a first glimpse into the complex interactions among water and gas molecules from a completely disordered state to an “ordered” state representative of a clathrate hydrate structure,Reference Walsh, Koh, Sloan, Sum and Wu¹³ as illustrated by the progression of structuring and ordering in Figure 4. This was accomplished through a direct microseconds molecular dynamics (MD) simulation at high-driving force conditions; the timescale required for the simulations demonstrates the intensive computational demand needed to capture a rare event phenomena that is associated with the nucleation of the hydrate nuclei. As in any major findings, more questions than answers have resulted from the first computational realization of hydrate nucleation at the molecular level. Other prior simulations studies used a similar approach to observe the direct hydrate nucleation, but within the timescale probed, only pre-nucleation molecular events were captured.Reference Zhang, Hawtin, Yang, Nakagava, Rivero, Choi and Rodger^15–Reference Radhakrishnan and Trout¹⁸

Figure 4.

Snapshots from the complex multi-step mechanism of hydrate nucleation revealed from a microseconds MD simulation.Reference Walsh, Koh, Sloan, Sum and Wu¹³ The sequence in the panels shows the evolution of the structuring of the methane and water molecules in the early stages of nucleation (a) to (c) and the arrangement of the cages after hydrate growth (d). Spheres and gray specks are methane and water molecules, respectively. Red lines are hydrogen bonds between water molecules. In panels (a), (b), and (c), the gold methane spheres are the first two molecules to induce the initial clustering of water into a cage-like structure. The purple methane sphere is the molecule occupying the central cage showed in panel (b). In panel (d), the different colors of the cages correspond to the 5¹² (pink), 5¹²6² (red), 5¹²6³ (brown), and 5¹²6⁴ (light blue).

More recently, other studies have also reported the nucleation and growth of hydrates and exposed additional insight into the transition from a disordered to ordered (hydrate-like) structuring.Reference Jacobson, Hujo and Molinero^19–Reference Liang and Kusalik²¹ A common characteristic of all hydrate nucleation simulations is the seemingly high local concentration of gas (i.e., methane) in solution leading to the initial ordering of hydrate-like structures. This observation on the concentration illustrates how simulations reveal a key aspect preceding the hydrate nucleation that would be nearly impossible to obtain from experimental measurements.

From the collection of studies on hydrate nucleation, it is clear that an outstanding question remains as to how to define and identify a critical nucleus. Simulations provide a detailed account of the molecular events leading up to the identification of hydrate-like structuring and cages formation. However, it is not immediately obvious what is the rate-limiting step before spontaneous growth of the crystalline structure. See Walsh et al.Reference Walsh, Koh, Sloan, Sum and Wu¹³ and Jacobson et al.Reference Jacobson, Hujo and Molinero¹⁹ for a more detailed description of these early molecular events.

Another intriguing aspect of the simulations on nucleation and growth is the hydrate structure formed upon nucleation. Unlike the expected equilibrium structure for methane hydrates, the crystal formed results in an amorphous structureReference Walsh, Koh, Sloan, Sum and Wu^13, Reference Jacobson, Hujo and Molinero¹⁹ that contains several other types of cages, besides the 5¹² and 5¹²6² expected in sI hydrates. Figure 5 illustrates the five main cages encountered in a collection of simulation trajectories of hydrate nucleation. These cages, with regular polyhedra having square, pentagonal, and hexagonal faces, may precede the crystallization into the known sI for methane hydrates. Because of the abundance of these cages in the simulations, one can also speculate how the transition from one cage type to another may be particularly important in the understanding of structural transition (solid-solid rearrangement).

Figure 5.

Sample images of water cages formed during hydrate nucleation. The red circles correspond to the oxygen atoms in the water molecules; hydrogen atoms have been omitted for clarity. In parentheses are the number of water molecules in the cage.

These simulations for probing the nucleation and growth mechanism of clathrate hydrates are computationally demanding but provide significant new insights that further our understanding of the rate-limiting steps and controlling factors impacting and determining the structure and properties of clathrate hydrates. However, direct MD simulations are limited to the exploration of single trajectories, which may or may not lead to hydrate nucleation and growth. As such, conventional MD has been mostly unsuitable for studying activated processes.Reference Frenkel and Smit²² This difficulty has led to the development of rare-event simulation techniques, often combining Monte Carlo (MC)/MD, equipped to overcome the difficulty in sampling a wider domain of phase space,Reference Chandler²³ within which homogeneous hydrate nucleation is one such rare event. Transition path sampling (TPS) and other rare-event simulation techniques are among the most promising methods to illuminate on the dynamics of hydrate nucleation.Reference Chandler^23–Reference Peters, Zimmerman, Beckham, Tester and Trout²⁵ These techniques harvest an ensemble of unbiased pathways of activated processes in a manner akin to a MC walk in trajectory space without assuming a mechanism a priori, thus overcoming the limitations in sampling on rough energy landscapes.

The application of rare-event simulations to study hydrate nucleation requires significant computational resources, as a large number of trajectories must be generated in order to sample the phase space containing a representative set of paths connecting the states along the free energy landscape, as shown in Figure 3. It is estimated that upward of several hundred trial trajectories are necessary to determine the mechanism and rate of hydrate nucleation for a single set of conditions, and that each path will require about a few hundred nanoseconds of simulation time for both forward and reverse paths. Thus, it is estimated that altogether, many 10s of microseconds of non-sequential paths must be simulated to arrive at a complete nucleation mechanism for a given temperature and pressure condition.

A fundamental understanding of hydrate nucleation and growth at the molecular level can potentially have a transformative impact in every aspect of the energy science of hydrates. For instance, in flow assurance, hydrates must either be avoided or managed; knowing how hydrates are formed will give better insight into how to prevent or control hydrate formation. The development and use of a hydrate kinetic inhibitor is a prime example of where the molecular details can greatly contribute in the design of better chemicals to prevent or delay the formation of hydrates. Specifically, a classical strategy would target the structure of the critical nucleus for binding by inhibitor molecules, and so knowledge of the critical nucleus would help in the rational design of the inhibitors. In the storage area for gases (e.g., natural gas, H₂), that process must be efficient; knowing how hydrates are formed will help to design a system that overcomes typical limitations of kinetics, mass, and heat transfer. In CO₂ sequestration and CH₄ recovery, the dissociation of CH₄ and formation of CO₂ hydrates must be controlled; knowing how hydrates are formed will determine the limitations and the viable approaches to sequester and recover the resources. Scientifically, many questions still exist on the metastability and phase transition of hydrates; knowing how hydrates are formed will give the insight into how to probe these phenomena. In a broader context, further studies of hydrate nucleation will enhance our knowledge toward nucleation theory and rare events, fundamental understanding arising from the hydrophobic effect, and new approaches for structural quantification, which ultimately would lead to clathrate hydrate material design for energy science.

HPC simulations of clathrate hydrates

The molecular level description of the nucleation and growth of clathrate hydrates is typically studied using MD simulations, so that the time evolution of the molecular events prior, during, and after nucleation can be captured. Because of the metastability and stochastic nature of hydrate nucleation, these types of simulations must be run for long times in order to fully observe the transition from a completely disordered to ordered state. As such, MD simulations into the microseconds timescale are often required on system sizes that may range from several to tens of nanometers (size of cubic box dimension) containing tens of thousands of molecules (water + guest). Simulations of this scale are most effectively run on high-performance computer clusters, taking advantage of parallel code architecture to efficiently generate the molecular trajectories. A number of open-source MD codes are available (e.g., GROningen MAchine for Chemical Simulations, GROMACS; Not (just) Another Molecular Dynamics program, NAMD; Large-scale Atomic/Molecular Massively Parallel Simulator, LAMMPS; DL_POLY), which can be readily compiled to run in most parallel computing environments. For the system sizes considered for hydrate simulations, it is typical to use upward of 50 computing cores. In terms of actual timescale for the simulations into the microsecond timescale would require several weeks to possibly months of dedicated computing time.

Summary and conclusions

Clathrate hydrates are unique crystalline structures formed from the encaging of small hydrophobic gas molecules in a hydrogen-bonded water framework. These structures are increasingly important in energy science, as they store vast amounts of natural gas in natural arctic and oceanic deposits, provide a medium for storage and transportation of natural gas and H₂, and pose a central problem in the management of flow assurance in oil and gas lines. Moreover, they form unique and complex structures that are scientifically challenging to study and understand due to the hydrophobic/hydrophilic interactions between the water and gas molecules, their non-stoichiometry, the multi-component nucleation, the metastability, and the solid-solid phase transition.

High-performance computing (HPC) will continue to play an important role in advancing the energy science of clathrate hydrates. Many fundamental questions are still unanswered on the formation and properties of clathrate hydrates and the associated structural transitions. Given this knowledge, how can we control/inhibit the formation of hydrates, and what are the limits of metastability for hydrates? How can we control the hydrate composition, or what new crystal structures are to be found? To begin to answer these and many other fundamental but practical questions, HPC can clearly be of benefit. Advanced computing platforms, perhaps designed to be particularly well-suited for a class of calculations, need to be coupled with advanced computational methods, which are capable of overcoming spatial and temporal challenges related to the dynamics of formation, structural rearrangement, and compositional differences.