Methane Hydrate Nucleation within Elastic Confined Spaces: Suitable

Subscriber access provided by University of South Dakota

Interfaces: Adsorption, Reactions, Films, Forces, Measurement Techniques, Charge Transfer, Electrochemistry, Electrocatalysis, Energy Production and Storage

Methane Hydrate Nucleation within Elastic Confined Spaces: Suitable Spacing and Elasticity Can Accelerate the Nucleation Jingpeng Hou, Dongsheng Bai, and Wei Zhou Langmuir, Just Accepted Manuscript • DOI: 10.1021/acs.langmuir.8b02387 • Publication Date (Web): 30 Aug 2018 Downloaded from http://pubs.acs.org on August 31, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Methane Hydrate Nucleation within Elastic Confined Spaces: Suitable Spacing and Elasticity Can Accelerate the Nucleation

Jingpeng Hou1, Dongsheng Bai1*, Wei Zhou1

1

Department of Chemistry, School of Science / Key Laboratory of Cosmetic, China National Light Industry, Beijing Technology and Business University, Beijing 100048, PR China * [email protected]

ABSTRACT Elastic materials are candidates for process intensification of gas storage by forming gas hydrate. In this work, molecular dynamics simulations of hydrate nucleation in elastic silica double layers were performed to study the effect of elastic confined spaces on hydrate formation. It is found that in narrow confined spaces, hexagonal rings dominated hydrogen-bond network of water molecules established rapidly by a multi-site nucleation mechanism. With molecules added, a bilayer water structure was formed finally because elastic space can adapt the volume expansion. In medium and wide confined spaces, hydrates were formed from a series of “pseudo cages” which is considered as precursors of complete hydrate cages. Moreover, the induction time for nucleation has a minimum when elasticity of the silica layer changes: nucleation is fastest in the weak-elastic system. When elasticity increases, it becomes hard to adapt the volume expansion during nucleation, and it also difficult to nucleate in very weak-elastic systems due to the fluctuation of the layers.

Keywords: methane hydrate, nucleation, elastic confined space, molecular dynamics simulation

ACS Paragon Plus Environment

Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

INTRODUCTION The consumption of traditional fossil energies is greatly accelerated by the industrialization process in developed countries. So far, all the countries dominated by the consumption of fossil energies in the world are facing with the challenge of the gradual depletion of energy resource. Natural gas hydrates, with its characteristics of wide distribution, high reserves, high energy density, and cleaning, are unanimously considered as a new potential alternative energy resource in 21st century. Natural gas hydrates are ice-like crystalline clathrate compounds consisting of water and methane molecules, in which methane molecules are enclosed in cavities formed by hydrogen bonded water molecules.1 Several types of hydrate lattice structures have been found, including sI (cubic), sII (cubic), sH (hexagonal),2 and HS-1(hexagonal) structure.3,4 Natural gas hydrates are usually sI structure, including six 51262 cages and two 512 cages in the unit cell. But when a small amount of larger guest molecules existed, it may become sII structure.5,6 As a new energy resource, natural gas hydrates are a geologically important class of materials and are abundant in the arctic region and marine sediment. Compared to the traditional energy storage technology, the energy density is expected to be higher in hydrate form. However, the gas hydrates prepared in laboratory usually have a certain amount of empty cages, which reduced its stability, leading to a discount on the potential advantages of energy storage. This is mainly caused by the low solubility of CH4 during the nucleation and growth of hydrates, which is much lower than the concentration of it in the prefect sI hydrate. Hence, improving the enrichment degree of CH4 in water is likely to decrease the nucleation barrier and increase the energy density of gas hydrate, so as to further enhancing the energy storage by using hydrate form. Porous materials, such as activated carbon, silica, zeolite, metal organic frameworks (MOF), can be introduced to enhance the adsorption of methane.7,8 In hydrate nucleation stage, using porous materials to enrich CH4 molecules might be beneficial to form the hydrates with high energy density. We should note that the hydrate formation is an expansion process (liquid water with density of ~1.0 g·cm-3 transformed into solid hydrate with ~0.9 g·cm-3). Most of porous materials, however, have rigid pores or channels, which shows a restraining effect on hydrate growth. Therefore, a class of elastic or flexible porous materials are considered, such as flexible MOF, covalent organic frameworks (COF), dry water,9 and so on. When elastic porous materials are absorbed with gas molecules, they usually show a breathing effect like a sponge.10 It can adapt to the expansion of volume during the hydrate formation and might be an ideal adsorption material.


Page 2 of 23

Page 3 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

So, the influence of the material elasticity on the dynamics of hydrate formation within elastic medium is an important topic, but few studies focus on it. Since the temporal and spatial limitations of the monitoring techniques, molecular simulation becomes a powerful method of providing molecular details of hydrate formation. Plenty of studies about hydrate nucleation and growth mechanism on solid surfaces were investigated by performing molecular dynamics (MD) simulations.11-14 Chen et al.12 investigated the CH4 hydrate formation in amended montmorillonite interlayer with leonardite humic acid sediment. They found a kinetic inhibition effect of the sediment on hydrate formation when CH4 adsorbed on the aggregated sediment. Cygan et al.15 found that the structure of CH4 hydrate in clays is different than in aqueous solution and bulk hydrate. Yan et al.13 performed MD simulations of CH4 hydrate nucleation in a system with clay layer. In their study, the diffusion of CH4 to the clay surface promotes the formation of semi-cages, and with the hydrate formed, it blocks the entrance of the pore, resulting a slow diffusion of CH4 molecules. He et al.14 studied the effect of substrates on CH4 formation in microsecond MD simulations by introducing a hydrophilic silica surface and a hydrophobic graphite surface. Water molecules are H-bonded with silanol groups in silica, while CH4 is preferred to adsorbed on the graphite surface. Kusalik et al.16 presented a MD simulation of CH4 hydrate growth in a partial rigid silica layers. They found the hydrate is preferred to grow towards the bulk solution, and the silica surfaces can serve as a CH4 source to promote the hydrate growth. In our previous studies, we investigated the nucleation process of CO2 hydrate induced by hydroxylated silica surfaces17,18 and hydrogenized silica surfaces.19 We found the existence of silica surfaces can accelerate hydrate nucleation, and the hydrophilicity of solid surfaces can alter the mechanism and dynamics of gas hydrate nucleation.19 In this work, we perform a series of MD simulations to investigate the hydrate nucleation process within elastic silica layers, a simplified model system for elastic confined space. And we will compare the results with rigid systems, to study the effect of elasticity of the silica surfaces on dynamics of hydrate nucleation. Model and Simulation Method The simulations were performed by using LAMMPS,20 and a typical initial configuration is shown in Figure 1. To obtain the initial configuration, we firstly added a pair of silica layers at the top and bottom of the simulation box, where the surfaces were saturated by –OH groups. The hydroxylated silica surfaces are suitable to act as confined space for the study of hydrate formation


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 23

and dissociation.14,16-19,21-23 The size of box is 6.25 × 6.25 × (L0+1.2) nm3, in which L0 equals to 0.5, 1.25, and 2.5 nm for different simulation systems, indicating the difference of the restriction in the confined space. For different system, a certain amount of water and methane molecules were added between the two silica layers randomly. The amount of the molecules is summarized in Table 1.

Figure 1. Typical initial configuration of the simulation system. In the figure, silica layers and water molecules are represented in stick style, methane molecules are shown in purple spheres, and orange dashed lines represent hydrogen bonds.

Table 1. The amount of water and methane molecules added in the initial configurations.

system

amount of water

amount of methane

L0 = 0.5 nm

460

60

L0 = 1.25 nm

1150

150

L0 = 2.5 nm

2300

300

In our simulations, the extensively studied TIP4P/ice water model24 was used and the rigidity was restricted with the SHAKE algorithm,25 while a single point model26 was employed to describe methane molecules. The intermolecular interactions of all these components are composed of a Lennard-Jones and a Coulomb part, which can be written as:

 σ Eij = 4ε ij  ij  rij 

12

  σ ij  −    rij

  

6

 qq + i j .  rij 

The parameters are listed in Table 2 for convenience. The force field can be successfully used for the simulation of hydrate formation.27,28 The silanols model29 was adopted for the silica layers,


Page 5 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

which is CHARMM style and the parameters can be found in ref. 29. The unlike parameters of Lennard-Jones potentials were obtained by the Lorentz-Berthelot mixing rule. A cutoff radius of 1.2 nm was utilized for the short-ranged interactions, and the long-ranged coulomb interactions were evaluated by using the pppm algorithm.30 As far as we know, the elastic restriction of silica layers has not been previously reported. Completely rigid silica models are commonly used in most of the simulations. In the work by Kusalik et al.,16 partial rigid model is proposed to allow the movement between the structural units of the silica, but the spacing between the layers cannot be changed obviously. To obtain the elastic of the confined space, we added a harmonic bond interactions between the pair of silica layers. Therefore, it seems to be a connection between the two layers by a dummy spring. The introduction of the spring brings an extra potential energy of U = k ( L − L0 ) 2 , showing the elasticity of the confining environment. For a system with its

equilibrium distance of L0, we described it as strong-elastic (SE), mid-elastic (ME), or weak-elastic (WE) if the relative deformation

∆L = 0.5, 1.0, or 1.5 when the potential increment L0

reaches 10 kcal, respectively. Take the system with L0 = 0.5 nm as an example, the SE, ME and WE system is equivalent to that the elastic coefficient k equals to 160, 40, and 17.8 kcal·nm-2, respectively. Furthermore, we also investigated the system of k = ∞ and k = 0, which are denoted as rigid (R) and flexible (F) system. The elastic coefficient k for different systems is listed in Table 3. To simplify the description, systems will be marked as “X-L0” in the following sections, for which X = R, SE, ME, WE, or F.

Table 2. Potential parameters for water and methane.

q (e)

ε (kcal/mol)

σ (Å)

water-O

-1.1794

0.21084

3.1668

water-H

0.5897

0

0

methane

0

0.294

3.730


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 23

Table 3. The elastic coefficient k used for different systems. Note: the unit of k in the table is kcal·nm-2.

L0 (nm)

0.5

1.25

2.5

rigid (R)

∞

∞

∞

strong-elastic (SE)

160

25.6

6.4

mid-elastic (ME)

40

6.4

1.6

weak-elastic (WE)

17.8

2.84

0.71

flexile (F)

0

0

0

elasticity

The whole simulation process was divided into three steps. First, an NVT relaxation process of 1 ns was performed at 260 K to eliminate the effect of the initial configuration. During the process, the position of silica layers was fixed. Then, a NpxyT simulation at 260 K and 25 MPa was performed on timescale of 50 ns (for systems with L0 = 0.5 nm) or 200 ns (for L0 =1.25 and 2.5 nm). Constant temperature and lateral xy-pressure were maintained by using the Nosé-Hoover algorithm.31 Finally, we added 1 methane and 8 water molecules randomly at the central region of the simulation box every 1 ns (for L0 = 0.5 nm) or 2 ns (for L0 =1.25 and 2.5 nm), to study the elastic response properties of the systems. The temperature and lateral pressure were maintained as the same as the previous step. This finally step runs for 50 ns (for L0 = 0.5 nm), 500 ns (for L0 = 1.25 nm) and 1500 ns (for L0 = 2.5 nm). The timestep we used is 0.5 fs. During the simulations, periodic boundary conditions were imposed in x- and y-direction, while in z-direction, shrink-wrapped non-periodic boundary was used. RESULTS AND DISCUSSION For elastic silica layers, there are two quantities describing its elastic properties: equilibrium distance L0, and elastic coefficient k. In this section, we will discuss the influence of confined spaces on hydrate nucleation from these two aspects. Self-organization of the system in narrow confined space. In a narrow-confined system with L0 = 0.5 nm, the layer spacing is roughly the same as the diameter of hydrate cages. Hence, the structure formed of the system is no longer the traditional gas hydrate structures. Since L0 is comparable to the thickness of water mono layer,32 the guest molecules will be enclosed by the water rings rather than water cages.


Page 7 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

The initial numerical density of water molecules is 23.55 nm-3, which is a low-density state compare to the density of bulk liquid water of 31.32 nm-3 (0.937 g·cm-3 in mass density33). At the NpxyT simulation stage (step 2), we found that the water molecules initially adjust their positions by translational motions between the layers, and form hydrogen bonds with –OH group of the silica surface and/or adjacent water molecules. As the simulation proceeded, a linear-branch-like hydrogen bond network forms rapidly, and the diffusion of liquid water is reduced. Subsequently, more hydrogen bonds are formed by frequent rotation and vibration motions of water molecules. And finally, water molecules form a structure dominated by six-membered rings, and methane molecules are packed inside.

Figure 2. Time evolution of order parameter ξ for different systems. The words “Geometric” and “HBond” marked in the legend are different criterions used to determine hydrogen bonds. We note that the curves of SE and WE are omitted for clarity since they are basically the same as those of ME.

To illustrate the formation of hydrogen bonds and six-membered rings, we define an order parameter ξ to represent the order of water molecules. In a quasi-planar water cluster, ξ equals to the number of edges (hydrogen bonds) shared by two rings divided by the number of rings. For example, ξ = 0 for a cluster consisted only by several isolated rings, and ξ = 0.5 when it contains a pair of adjacent rings sharing one edge. With the development of clusters (i.e. the increase in the number of hydrogen bonds), the value of ξ will gradually increase. For a full-developed system with planar six-membered rings, like graphite layer structure, ξ equals to 3. We note that a cluster


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

contains all interconnected water molecules, and different connectivity criterion will lead to different cluster sizes. Figure 2 shows the time evolution of ξ. The dashed lines in the figure indicates that the connectivity of two water molecules is determined by geometric criterion: the distance of oxygen atom of two water molecules is within 2.76 Å.1 While, the solid line means the connectivity is determined by hydrogen-bond criterion: both the distance of two oxygen atoms is within 2.76 Å and the H－O…O angle is less than 30°.34 In Figure 2, one can see clearly that at the initial stage, the value of ξ based on geometric criterion is higher than that based on hydrogen-bond criterion for all of the simulation systems. This indicates that the water molecules adjust its position quickly, forming a hexagonal structure within ~12 ns. However, the hydrogen bond network is not fully established at that moment, which can be seen in Figure 3a. Then, water molecules rotate themselves to form more hydrogen bonds. After ~18 ns, the ξ value becomes consistent by the two criteria, indicating the finally formation of a complete hydrogen bond network. Figure 3 shows the structure change of system R-0.5 from 12 ns to 18 ns. Obviously, the hydrogen bond network is established in two steps. We note that this asynchronous phenomenon is also found in hydrate nucleation process.18 Interestingly, one can be found from Figure 3 that the formation of hydrogen bond network is simultaneous in many sites, rather than spreading from one “nucleus”. This probably because the formation of hydrogen bonds requires only a slight “in situ” rotation of water molecules after ~12 ns. In confined space, the energy barrier of rotation is much smaller than that of translation. Therefore, we think it is possible to simultaneously nucleate in multiple positions.


Page 8 of 23

Page 9 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Figure 3. Structure change of water layer in system R-0.5. (a) 12 ns, (b) 14 ns, (c) 16 ns, and (d) 18 ns. The gray filled parts in the figure are closed hexagonal rings formed by water molecules. The arrows indicate the direction of the movement of CH4 and the location after 2 ns.

After the simulations performed for 50 ns, additional water and methane molecules were added to the system, as described in the previous section. For the rigid system R-0.5, the increase of molecules leads to structural disorder very quickly, showing a significant decrease in ξ order parameter in Figure 2. The snapshot in Figure 4 also demonstrated the structural disorder clearly. For elastic and flexible systems, however, the order parameter has only slight decrease because of the adaptive variation of the silica layers. As seen from the snapshot showing in Figure 4, water molecules form a double-layered six-membered rings structure. The root-mean-square deviation of those rings (marked in Figure 4) indicates that most of the rings are non-planar six-membered rings with chair conformation, and methane molecules are adsorbed on one side of the rings. This structure is very similar to the precursor that forms hydrate cages. Jacobson et al.35 suggested that the adsorbed guest molecules should be considered as a part of the crystal nucleus. While, before the molecules added, the rings formed in each system are almost all planar six-membered rings, and methane is filled in the center of the ring, like the structure reported by Zeng et al.32


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4. Structure evolution of water layers in different systems. The numbers marked at the bottom-right corner of each panel are root-mean-square deviation (RMS) of hexagonal water rings, and we note that the RMS value of the chair conformation of water in ice Ih lattice is ~1.277.


Page 10 of 23

Page 11 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

The time evolution of the water density and the distance of silica layers is shown in Figure 5 for each system. With the increase of the molecules added to the system, the expansion of the space between silica layers becomes slow for SE-0.5 system. The bilayer water structure began to break, and methane molecules aggregated (Figure 4). A high-density ice phase without methane filled in is formed within a local region at ~87 ns, and the local density of water reaches ~33.2 nm-3. At this moment, the size of the confined space fluctuated in a certain degree. The ME-0.5, WE-0.5 and F-0.5 systems, however, are still keeping the double-layered structure until the end of the simulations.

Figure 5. Time evolution of water density and the distance of silica layers.

Hydrate nucleation in medium confined space. In this section, the initial distance between the silica layers is 1.25 nm, which is a weak restriction on the motion of individual water molecules but is still a strongly restricted space for hydrate formation since it can only accommodate one unit cell of the hydrate in thickness. In addition, it is note that almost all water molecules in the system are interfacial water, because the effect of the solid surface on water molecules can reach ~1.4 nm.36 Considering the induction time of hydrate nucleation within confined spaces is ~500 ns in our previous studies,17 we extended the simulation time and slowed down the rate of molecule adding, as is described in section 2.


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 6. Time evolution of the number of rings and cages.

Firstly, by using the ring perception algorithm37 and the cage identification algorithm,38 we calculated the number of rings (pentagon and hexagon) and cages (512 and 51262), and the results are shown in Figure 6. In R-1.25 and ME-1.25 systems, only pentagon and hexagonal rings formed, and there are almost no hydrate cages. While in F-1.25 systems, the hydrate cages are formed obviously after ~420 ns, indicating the induction time of nucleation in these systems is ~420 ns. Comparing to the hydrate formation in rigid confined spaces,17 the nucleation has been intensified. This may be because the continuous addition of water and methane molecules has a promoting effect on nucleation. We guess that the dynamic evolution of local hydrogen bond network might be enhanced when the molecules artificially added, which is beneficial for the structure transform from liquid water to hydrate. However, more detailed research needs to be further developed. Changing the rate and amount of the added molecules may be a good starting point to validate our guess, and we will discuss it in our subsequent studies.


Page 12 of 23

Page 13 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Figure 7. The F4 order parameter profile along z-direction at 650 ns. To facilitate comparison, we normalized the thickness of each system in z-direction.

Besides, we calculated the F4 order parameter profile of water molecules39,40 along z-direction, and the results for the systems at 650 ns are shown in Figure 7. We found that an ice-like layer is appeared near the silica surfaces, reflecting in a negative F4 value marked in Figure 7. The results are consistent with our earlier studies.17,19 In central area, the F4 order parameter increases for F-1.25 system, indicating the tendency to form hydrates. For R-1.25 and ME-1.25 systems, to our surprise, although the F4 order parameter is near 0 (which means the water in system should have been liquid state41), a certain amount of “pseudo cage” structures were formed in all systems. A 512 pseudo cage is defined as a structure that each edge of a central pentagon ring is connected with at least one five-membered rings. While a 51262 pseudo cage is a structure similar to 512 pseudo cage, and the difference is that the central ring is a hexagon. Pseudo cages can be considered as precursors of hydrate cages. In the relaxation stage of nucleation, plenty of pseudo cages appear frequently in the system, as shown in Figure 8. For 51262 pseudo cages, since the central ring is mostly chair type, the five-membered rings connected are usually located on both sides of the central ring, forming several structures listed in Figure 8. In our simulations, individual half 51262 cages were never observed, while double-half cages were observed only when the full 51262 cage immediately formed. The transition of the central hexagonal ring from nonplanar to planar might be the resistance to form full 51262 cages. We found that only when the connection number of the central ring (i.e. the number of


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

five-membered rings connected with the central ring) is high enough, like 5961, 51261 pseudo cages shown in Figure 8, can the central ring becomes planar.

Figure 8. Typical pseudo cages appeared during the nucleation stage of the simulations.

Figure 9. Typical migration of central ring of pseudo cages during the nucleation stage of the simulations.

For 512 pseudo cages, almost all structures listed in Figure 8 can be observed, and the migration of central ring occurs frequently. A typical migration process is shown in Figure 9. With the breaking and generation of hydrogen bonds, the original central ring marked as A (with low connection number of five) at the left side becomes a side ring, while its neighbor pentagon ring marked as B on the right side becomes a new central ring which has a high connection number of seven. For a central ring, the larger the connection number it has, generally speaking, the more stable it is. Actually, we need to note that the stability is also related to the presence or not of guest molecules.


Page 14 of 23

Page 15 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Figure 10. Time evolution of the number of central rings in the pseudo cages during the nucleation process.

In the initial stage of hydrate nucleation, by the manner of hydrogen bond formation, water rings coalescence, and central rings migration, more stable structures with high connection number are gradually generated in the system. Then, clusters with single central ring developed into that with dual or even multi central rings, and eventually, form complete hydrate cages. Figure 10 shows the time evolution of the average number of central rings in each cluster for the system of ME-1.25. Based on our definition of pseudo cages, we note that in 512 cages, each ring is central ring, while in 51262 cages, only the two hexagonal rings are central rings. This is caused by the difference of symmetry of the two kinds of cages. The higher symmetry of 512 cages is facilitate to the migration of central rings, which may explain why 512 cages are always formed firstly in the nucleation stage, even it is sometimes empty.35,42

Hydrate nucleation in wide confined space. The initial distance of silica layers is 2.5 nm we used in this section, which provides sufficient space for hydrate nucleation, but it still has certain restrictions on the growth of hydrate. By calculating the time evolution of cage amounts, we can obtain the induction time of hydrate nucleation roughly. The induction time of the system with different elastic is shown in Figure 11. To make it clear, we performed several additional simulations for the systems with other elastic coefficients, and also extended the simulation time in order to observe the nucleation events.


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 11. The hydrate nucleation induction time for systems with different elastic.

One can see that for weak-elastic system, the induction time is the shortest. With the increase of elasticity, the induction time increased slightly. At the end of the induction time (i.e. the time hydrate began to growth after the nucleation stage), these systems (with k from 0.5 to 3.2 kcal· nm-2) have a similar spacing of silica layers, which is roughly equal to 2.43 (±0.08) unit cells of sI methane hydrate. This might be the smallest size of confined spaces that required for hydrate growth continuously from its nucleus. When k is larger enough, the confined space becomes almost rigid, and the induction time increases sharply, indicating the difficult for hydrate growth. For the system of R-2.5, our previous study at similar simulation conditions shows that the induction time is ~1.9 µs.17 The rigid system has a longer induction time, because the hydrate formation is usually accompanied by an increase in volume, and rigid confined space cannot adapt to the change of volume, inhibiting the phase transition of hydrate. To our surprise, we found that as the elasticity becomes weaker (k < 0.5 kcal·nm-2), instead of the decreasing of the induction time, it increases dramatically. This is because the smaller elasticity leads to an obvious fluctuation of the silica layers during the whole simulation process, and the distance of the two layers can be very far when water and methane molecules added, as is shown in Figure 12. We note that although the silica layer has a stabilizing effect on hydrate cages,18 the dramatic fluctuation of the position of silica surfaces may cause the damage of hydrogen bond network close to it. As a result, the nucleation of hydrates close to the silica surface is inhibited. For a very weak-elastic system with k = 0.1 kcal·nm-2, nucleation was not observed in


Page 16 of 23

Page 17 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

our simulation time scale. Considering our previous work,17 excessive spacing between the pair of silica layers is not the reason of nucleation difficulties. The main reason should be the influence of the fluctuation of the layers on the hydrogen bond network. In this situation, quasi-homogeneous nucleation at the central region of the system may be easier than heterogeneous nucleation induced by solid surfaces.

Figure 12. Time evolution of the distance of silica layers. Note that the curves are interrupted at ~910 ns for R-2.5 and ~1450 ns for ME-2.5, because the energy of the systems is high enough at that time which caused the failure of the simulation.

We noted that the CH4 concentration used (0.115 in mole fraction) is particularly high, which may lead to phase separation due to the spinodal decomposition in bulk system. To prove this point, we performed an additional simulation of hydrate nucleation in bulk phase without silica surfaces. The simulation box with the size of 6.25 × 6.25 × 2.5 nm3 is used, which is the same as the wide confined space. Due to the spinodal decomposition of the system at ~7.2 ns, the concentration of CH4 in water solution becomes very low, obviously increasing the induction time for hydrate nucleation (~3.71 µs marked in Figure 11). Within confined space, we think it should be phase-separated too in the final configuration, but it not occurred spontaneously. It may be ascribed to two aspects: 1) when the solid hydrate crystals formed, the free space becomes narrow between the silica layers, which is difficult for CH4 to form a bubble; and 2) because the pentagonal and hexagonal rings of the hydrate cages which exposed to the liquid mixture can adsorb guest molecules,43,44 it can be expected that some CH4 molecules in liquid phase will


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

adsorbed on the faces of the cages, which facilitated the dispersing of methane. Both the two aspects are unfavorable to form a CH4 bubble. Besides, the confining can significantly affect the thermodynamic properties (e.g. decrease on melting temperature) of the system. So, it is hard to give the exact degree of supersaturation of our composition within confined space. As a result, the phase separation by spinodal decomposition in confined space seems to be more difficult. Moreover, when CH4 concentration is close to the concentration of spinodal decomposition, hydrate nucleation can also be accelerated. In this situation, unlike the classical nucleation theory, the energy barrier of phase transition is very small. Nucleation is so fast that it is almost no induction time, even it is occurred in multiple sites at the same time. At 240 K and 20 MPa, Sarupria et al.45 found a fast nucleation of hydrates when methane mole fraction xm is ~0.07, and Jimenez et al.46 found that the spinodal decomposition occurred when xm is higher than 0.09 at 285 K and 50 MPa. However, we should note that those studies were performed in bulk system. It will become more complex in the confined space, which is still an open question.

CONCLUSIONS In this work, by using molecular dynamics simulation method, we investigated the methane hydrate nucleation process within elastic confined silica layers systematically. Our simulation results reveal that the elastic strength and the size of the confined space has obvious influence on hydrate formation. In a narrow confined space, it can form a six-membered ring network structure, and methane molecules located on one side of the ring. The formation of the network follows a multi-site nucleation mechanism. When molecules added, elastic confined space can adapt the change of volume, and finally a bilayer water structure was formed. In medium confined systems, hydrate cages were only formed in weak-elastic and fixable systems. While in the initial stage, a certain number of pseudo cages, precursor of complete hydrate cages, were formed in all systems. For wide confined space, we found that there is a minimum for induction time of hydrate nucleation: it is fastest in the weak-elastic system. With the elasticity increases, the confined space becomes more and more unable to adapt the increase in volume during phase transition. While for very weak-elastic system, it also difficult to nucleate owing to the fluctuation of silica layers. Therefore, only the system with proper elasticity and interlayer spacing can accelerate the hydrate formation process.

ACKNOWLEDGEMENT This work is supported by National Natural Science Foundation of China (No. 21403009).


Page 18 of 23

Page 19 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

REFERENCES 1.

Sloan, E. D.; Koh, C. A. Clathrate Hydrates of Natural Gases, 3rd ed.; CRC Press: Boca Raton, FL, 2008.

2.

Buffett, B. A. Clathrate Hydrates. Annu. Rev. Earth Planet. Sci. 2000, 28, 477-507.

3.

Vatamanu, J.; Kusalik, P. G. Unusual Crystalline and Polycrystalline Structures in Methane Hydrates. J. Am. Chem. Soc. 2006, 128, 15588-15589.

4.

Yang, L.; Tulk, C. A.; Klug, D. D.; Moudrakovski, I. L.; Ratcliffe, C. I.; Ripmeester, J. A.; Chakoumakos, B. C.; Ehm, L.; Martin, C. D.; Parise, J. B. Synthesis and Characterization of a New Structure of Gas Hydrate. Proc. Natl Acad. Sci. USA 2009, 106, 6060-6064.

5.

Hendriks, E. M.; Edmonds, B.; Moorwood, R. A. S.; Szczepanski, R. Hydrate Structure Stability in Simple and Mixed Hydrates. Fluid Phase Equilib. 1996, 117, 193-200.

6.

Sloan, E. D.; Subramanian, S.; Matthews, P. N.; Lederhos, J. P.; Khokhar, A. A. Quantifying Hydrate Formation and Kineticinhibition. Ind. Eng. Chem. Res. 1998, 37, 3124-3132.

7.

Billemont, P.; Coasne, B.; Weireld, G. D. Adsorption of Carbon Dioxide, Methane, and Their Mixtures in Porous Carbons: Effect of Surface Chemistry, Water Connect, and Pore Disorder.

Langmuir 2013, 29, 3328-3338. 8.

Pantatosaki, E.; Pazzona, F. G.; Megariotis, G.; Papadopoulos, G. K. Atomistic Simulation Studies on the Dynamics and Thermodynamics of Nonpolar Molecules within the Zeolite Imidazolate Framework-8. J. Phys. Chem. B 2010, 114, 2493-2503.

9.

Li, Y.; Zhang, D.; Bai, D.; Li, S.; Wang, X.; Zhou, W. Size Effect of Silica Shell on Gas Uptake Kinetics in Dry Water. Langmuir 2016, 32, 7365-7371.

10. Serre, C.; Mellotdraznieks, C.; Surblé, S.; Audebrand, N.; Filinchuk, Y.; Férey, G. Role of Solvent-Host Interactions That Lead to Very Large Swelling of Hybrid Frameworks. Science

2007, 315, 1828-1831. 11. Park, S. H.; Sposito, G. Do Montmorillonite Surfaces Promote Methane Hydrate Formation? Monte Carlo and Molecular Dynamics Simulations. J. Phys. Chem. B 2003, 107, 2281-2290. 12. Ji, H.; Wu, G.; Zi, M.; Chen, D., Microsecond Molecular Dynamics Simulation of Methane Hydrate Formation in Humic-Acid-Amended Sodium Montmorillonite. Energy Fuels 2016,

30, 7206-7213.


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

13. Yan, K. F.; Li, X. S.; Chen, Z. Y.; Xia, Z. M.; Xu, C. G.; Zhang, Z. Molecular Dynamics Simulation of the Crystal Nucleation and Growth Behavior of Methane Hydrate in the Presence of the Surface and Nanopores of Porous Sediment. Langmuir 2016, 32, 7975-7984. 14. He, Z.; Linga, P.; Jiang, J. CH4 Hydrate Formation between Silica and Graphite Surfaces: Insights from Microsecond Molecular Dynamics Simulations. Langmuir 2017, 33, 11956-11967. 15. Cygan, R. T.; Guggenheim, S.; Koster van Groos, A. F. Molecular Models for the Intercalation of Methane Hydrate Complexes in Montmorillonite Clay. J. Phys. Chem. B

2004, 108, 15141-15149. 16. Liang, S.; Rozmanov, D.; Kusalik, P. G. Crystal Growth Simulations of Methane Hydrates in the Presence of Silica Surfaces. Phys. Chem. Chem. Phys. 2011, 13, 19856-19864. 17. Bai, D.; Chen, G.; Zhang, X.; Wang, W. Microsecond Molecular Dynamics Simulations of the Kinetic Pathways of Gas Hydrate Formation from Solid Surfaces. Langmuir 2011, 27, 5961-5967. 18. Bai, D.; Chen, G.; Zhang, X.; Wang, W. Nucleation of the CO2 Hydrate from Three-Phase Contact Lines. Langmuir 2012, 28, 7730-7736. 19. Bai, D.; Chen, G.; Zhang, X.; Sum, A. K.; Wang, W. How Properties of Solid Surfaces Modulate the Nucleation of Gas Hydrate. Sci. Rep. 2015, 5, 12747. 20. Plimpton, S. J. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput.

Phys. 1995, 117, 1-19. 21

Bagherzadeh, S. A.; Englezos, P.; Alavi, S.; Ripmeester, J. A. Molecular Modeling of the Dissociation of Methane Hydrate in Contact with a Silica Surface. J. Phys. Chem. B 2012,

116, 3188-3197. 22

Bagherzadeh, S. A.; Englezos, P.; Alavi, S.; Ripmeester, J. A. Influence of Hydrated Silica Surfaces on Interfacial Water in the Presence of Clathrate Hydrate Forming Gases. J. Phys.

Chem. C 2012, 116, 24907-24915. 23

Liang, S.; Kusalik, P. G. The Nucleation of Gas Hydrates Near Silica Surfaces. Can. J. Chem.

2015, 93, 791-798. 24. Abascal, J. L.; Sanz, E.; García, F. R.; Vega, C. A Potential Model for the Study of Ices and Amorphous Water: Tip4p/Ice. J. Chem. Phys. 2005, 122, 234511-234519. 25. Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian


Page 20 of 23

Page 21 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J.

Comput. Phys. 1977, 23, 327-341. 26. Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Optimized intermolecular potential functions for liquid hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638-6646. 27. Walsh, M. R.; Beckham, G. T.; Koh, C. A.; Sloan, E. D.; Wu, D. T.; Sum, A. K. Methane Hydrate Nucleation Rates from Molecular Dynamics Simulations: Effects of Aqueous Methane Concentration, Interfacial Curvature, and System Size. J. Phys. Chem. C 2011, 115, 21241-21248. 28. Walsh, M. R.; Rainey, J. D.; Lafond, P. G.; Park, D.; Beckham, G. T.; Jones, M. D.; Lee, K.; Koh, C. A.; Sloan, E. D.; Wu, D. T.; Sum, A. K. The Cages, Dynamics, and Structuring of Incipient Methane Clathrate Hydrates. Phys. Chem. Chem. Phys. 2011, 13, 19951-19959. 29. Lopes, P. E.; Murashov, V.; Tazi, M.; Demchuk, E.; Jr, M. A. Development of an Empirical Force Field for Silica. Application to the Quartz-Water Interface. J. Phys. Chem. B 2006, 110, 2782-2792. 30. Hockney, R. W.; Eastwood, J. W. Computer Simulation Using Particles; New York, Taylor & Francis, 1989. 31. Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A

1985, 31, 1695-1697. 32. Bai, J.; Angell, C. A.; Zeng, X. C. Guest-Free Monolayer Clathrate and Its Coexistence with Two-Dimensional High-Density Ice. Proc. Natl Acad. Sci. USA 2010, 107, 5718-5722. 33. Vega, C.; McBride, C.; Sanz, E.; Abascal, J. L. F. Radial Distribution Functions and Densities for the SPC/E, TIP4P and TIP5P models for liquid water and ices Ih, Ic, II, III, IV, V, VI, VII, VIII, IX, XI and XII. Phys. Chem. Chem. Phys. 2005, 7, 1450-1456. 34. Zhang, J.; Hawtin, R. W.; Yang, Y.; Nakagava, E.; Rivero, M.; Choi, S. K.; Rodger, P. M. Molecular Dynamics Study of Methane Hydrate Formation at a Water/Methane Interface. J.

Phys. Chem. B 2008, 112, 10608-10618. 35. Jacobson, L. C.; Hujo, W.; Molinero, V. Amorphous Precursors in the Nucleation of Clathrate Hydrates. J. Am. Chem. Soc. 2010, 132, 11806-11811. 36. Argyris, D.; Tummala, N. R.; Striolo, A.; Cole, D. R. Molecular Structure and Dynamics in Thin Water Films at the Silica and Graphite Surfaces. J. Phys. Chem. C 2008, 112, 13587-13599.


Langmuir 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

37. Matsumoto, M.; Baba, A.; Ohmine, I. Topological Building Blocks of Hydrogen Bond Network in Water. J. Chem. Phys. 2007, 127, 134504-134512. 38. Jacobson, L. C.; Hujo, W.; Molinero, V. Thermodynamic Stability and Growth of Guest-Free Clathrate Hydrates: A Low-Density Crystal Phase of Water. J. Phys. Chem. B 2009, 113, 10298-10307. 39. Rodger, P. M.; Forester, T. R.; Smith, W. Simulations of the Methane Hydrate/Methane Gas Interface near Hydrate Forming Conditions Conditions. Fluid Phase Equilib. 1996, 116, 326-332. 40. Báez, L. A.; Clancy, P. Computer Simulation of the Crystal Growth and Dissolution of Natural Gas Hydrates. Ann. N.Y. Acad. Sci. 1994, 715, 177-186. 41. Moon, C.; Hawtin, R. W.; Rodger, P. M. Nucleation and Control of Clathrate Hydrates: Insights from Simulation. Faraday Discuss. 2007, 136, 367-382. 42. Jacobson, L. C.; Hujo, W.; Molinero, V. Nucleation Pathways of Clathrate Hydrates: Effect of Guest Size and Solubility. J. Phys. Chem. B 2010, 114, 13796-13807. 43. Walsh, M. R.; Koh, C. A.; Sloan, E. D.; Sum, A. K.; Wu, D. T. Microsecond Simulations of Spontaneous Methane Hydrate Nucleation and Growth. Science 2009, 326, 1095-1098. 44. Guo, G.; Zhang, Y.; Liu, H. Effect of Methane Adsorption on the Lifetime of a Dodecahedral Water Cluster Immersed in Liquid Water: A Molecular Dynamics Study on the Hydrate Nucleation Mechanisms. J. Phys. Chem. C 2007, 111, 2595-2606. 45. Sarupria, S.; Debenedetti, P. G. Homogeneous Nucleation of Methane Hydrate in Microsecond Molecular Dynamics Simulations. J. Phys. Chem. Lett. 2012, 3, 2942-2947. 46. Jiménez-Ángeles, F.; Firoozabadi, A. Nucleation of Methane Hydrates at Moderate Subcooling by Molecular Dynamics Simulations. J. Phys. Chem. C 2014, 118, 11310-11318.


Page 22 of 23

Page 23 of 23 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

TOC Image


Methane Hydrate Nucleation within Elastic Confined Spaces: Suitable

Recommend Documents