How Properties of Solid Surfaces Modulate the Nucleation of Gas Hydrate

Molecular dynamics simulations were performed for CO2 dissolved in water near silica surfaces to investigate how the hydrophilicity and crystallinity of solid surfaces modulate the local structure of adjacent molecules and the nucleation of CO2 hydrates. Our simulations reveal that the hydrophilicity of solid surfaces can change the local structure of water molecules and gas distribution near liquid-solid interfaces, and thus alter the mechanism and dynamics of gas hydrate nucleation. Interestingly, we find that hydrate nucleation tends to occur more easily on relatively less hydrophilic surfaces. Different from surface hydrophilicity, surface crystallinity shows a weak effect on the local structure of adjacent water molecules and on gas hydrate nucleation. At the initial stage of gas hydrate growth, however, the structuring of molecules induced by crystalline surfaces are more ordered than that induced by amorphous solid surfaces.

In nature, gas hydrates can form in different geological environments, for example, in sedimentary rock or unconsolidated clay. Thus, different solid surfaces from crystalline to amorphous and those with hydrophilic to hydrophobic characteristics should be considered in the study of hydrate nucleation. As such, the important fundamental question on how surface properties affect the hydrate formation must be understood.
In this work, a number of solid surfaces with different hydrophilicity and crystallinity are considered to study the nucleation mechanism and dynamics of hydrate formed from those surfaces. We choose the formation of CO 2 hydrate from silica surfaces for this study. The choice of CO 2 hydrate comes from the interest in sequestering CO 2 , a greenhouse gas, by ocean storage technologies, such as the direct injection of CO 2 into the ocean [17][18][19][20][21] . In such case, the injection of CO 2 into sediments could sequester CO 2 in the hydrate form 1,[22][23][24][25] , and also create a CO 2 hydrate layer that would prevent liquid CO 2 from dissolving into the seawater. The successful application of this approach requires a fundamental understanding of the mechanism of CO 2 hydrate nucleation in the presence of silica surfaces. The consideration of different types of surfaces (crystalline/amorphous and hydrophilicity) are important since most actual surfaces, including sedimentary rock or unconsolidated clay, are rough and non-uniform in terms of hydrophilicity.

Results
Series of microsecond scale molecular dynamics simulations were performed to investigate the nucleation mechanisms of hydrate with silica surfaces have different properties. A typical initial configuration is shown in Fig. 1, and the hydrate formation processes for several systems are shown in Fig. 2. To simplify the description, systems with different silica surfaces are marked as crys-n or amor-n, for which the 'crys' or 'amor' corresponds to a crystalline or amorphous solid surface, and the number of n indicates the percentage of -OH groups saturating the surface (see Fig. 3 in Methods section for details). Therefore, from the crys-100 to the crys-0 systems, the -OH groups saturating the crystalline surfaces are increasingly changed into -H groups, decreasing the solid surface hydrophilicity.   28 . Since the simulation systems in this study are not a bulk water system, the structure of water close to the surface might be distorted due to the existence of silica. Hence, we note that a solid-like structure of water with a minus value of ϕ F 4 will be classified into an ice-like structure. Figure 4a shows the obtained values of ϕ F 4 for various systems with different surface hydrophilicity. In general, CO 2 hydrate was formed within 4 μ s for all the silica surfaces studied, indicating the presence of hydrophilic surfaces prompts hydrate nucleation.
In our previous study 15 , we considered a super-saturated CO 2 solution system with 100% -OH groups and crystalline silica surface (denoted as crys-100, see Fig. 3), and it shows a three-step nucleation mechanism: an ice-like layer closest to the silica surface formed first, followed by an intermediate thin layer subsequently formed to compensate for the structure mismatch between the ice-like layer and the final stable CO 2 hydrate, and finally a CO 2 hydrate layer generated from the intermediate structure acting as  Fig. 4) are shown in xy-plane. For clarity, CO 2 molecules are omitted. Color code for the snapshots is the same as in Fig. 1, and the hydrate cages are shown with red wire-frame bonds.
Scientific RepoRts | 5:12747 | DOi: 10.1038/srep12747 Figure 3. Top view of various silica surfaces studied in this work. The red and white spheres represent -OH and -H groups on silica surfaces, respectively. In those systems, the type of -OH or -H group is assigned randomly while keeping the percentage of -OH groups to a given value.
Scientific RepoRts | 5:12747 | DOi: 10.1038/srep12747 a nucleation zone. In this work, although the concentration of CO 2 is decreased, the three-step nucleation mechanism is also existed for crys-100 system. However, as we will discussed below, the nucleation mechanism will be changed when hydrophilicity of the surface changes.
For the solid surfaces with weak hydrophilicity (crys-25 and crys-0), the amorphous ice-like layer closest to the silica surface vanishes, only a single liquid-like intermediate layer is formed directly on the surface, and as a result the nucleation mechanism changes. The diffusion coefficients of water molecules in the region closest to the solid surface are calculated based on the mean square displacement, and the results are listed in Table 1. One can see clearly that for the system of crys/amor-100, 75 and 50, the diffusion coefficients are small enough, indicating a solid-like layer formed. For the system of crys/amor-25 and 0, however, the diffusion coefficients increased by one order of magnitude, indicating a liquid-like layer close to the solid surface. Figure 4a also shows that for weakly hydrophilic surfaces, the ϕ F 4 in region A is close to that for liquid water, indicating that the CO 2 hydrate directly nucleates through an intermediate (liquid-like) water layer coating on the surfaces, rather than from an ice-like water layer as seen for strongly hydrophilic surfaces (crys-100). As a result, the three-step nucleation mechanism 15 changes into a two-step nucleation mechanism.  To determine the reason for the variation of nucleation mechanism as a function of surface hydrophilicity, we investigated the local water structure near the silica surfaces. The results show that the water layer close to the solid surfaces becomes less structured when the surface hydrophilicity decreases. For example, for the crys-0 system (the solid surfaces with the lowest degree of hydrophilicity), the value of ϕ F 4 close to the silica surface increases to 0.05, significantly deviating from − 0.2 for the crys-100 interface (Fig. 4a). This observation again proves that for the solid surfaces with low hydrophilicity the amorphous ice-like interfacial layer vanishes, and a liquid-like interfacial water layer appears.
An interesting observation found is that the amorphous ice-like water layer close to the silica surfaces can persist even for the case with 50% -OH groups. As is shown in Fig. 4a for the crys-75 and crys-50 systems, the ϕ F 4 value for the interfacial layer (region A) fluctuates around − 0.2, similar to that in the crys-100 system. However, for the crys-25 system, the small number of -OH groups is insufficient to stabilize the ice-like structure of the water layer, and as a result, the ice-like layer disappears. This is confirmed by the observation that near the silica surface, both the value of ϕ F 4 increases to 0.05 (Fig. 4a), and the diffusion coefficient increases to ~10 −9 m 2 /s, similarly to that for the crys-0 system.
Weakening surface hydrophilicity not only changes the local structure of water molecules at the interface, but also weakens the tendency for CO 2 molecules to leave from the solid surface. Figure 5 gives the evolution of CO 2 density in the region close to the silica surfaces (within a distance of 0.75 nm). The figure shows that the final CO 2 density near a silica surface depends strongly on the surface hydrophilicity: the more hydrophilic the surface is, the lesser the number of CO 2 molecules are remained. For weakly hydrophilic surfaces, a relatively large amount of CO 2 molecules are remained on the surfaces, and inevitably, the existence of CO 2 molecules near the silica surface would affect the local structure of the water layer closest to the surface. The existence of CO 2 in the water layer can hinder the full development of the hydrogen-bond network, and should be considered as another reason for the disappearance of amorphous ice-like layer (region A), especially for the system with surfaces weakly hydrophilic.
The effect of surface hydrophilicity on the distribution of CO 2 molecules is consistent with that of local structure of the water layer. From the crys-75 to crys-50 systems, the corresponding CO 2 densities of those systems are just slightly increased compared with the crys-100 system, whereas the hydrophilicity of the silica surface is significantly decreased (Fig. 5). For the crys-25 and crys-0 systems, however, there are less -OH groups to hold the water molecules, and instead the surface can adsorbed more CO 2 molecules. As such, one could suggest that the -H-modified silica surface can be considered CO 2 -philic.
Moreover, the calculation results of the density of water molecules in region A (Table 1) show that the amorphous ice-like layer is a low-density "ice". With the weakening of the surface hydrophilicity, one can see that the density of water is decreased. This is partly because the adsorption capacity of solid surfaces decreased when -OH groups are changed to -H groups (see Discussion section for details), and partly because the amount of CO 2 molecules remained within the layer becomes larger.
The presence of hydrophilic surfaces also affects the dynamics of hydrate nucleation. By using the ring perception 29 and the cage identification 30 algorithms, we calculated the time evolution of the number of hydrate cages, as shown in Fig. 6. From the figure it is surprising to find that the induction time for the crystal nucleation is reduced when less -OH groups exist on the silica surface. In other words, CO 2 hydrates can be more easily formed from less hydrophilic solid surfaces.
Effect of surface crystallinity. The study of the effect of the silica surfaces crystallinity on hydrate formation process is started by comparing the crystalline and amorphous surfaces. Figure 4b shows the cules near amorphous and crystalline surfaces have about the same value as long as the surface hydrophilicity is the same. As such, the water and hydrate structures near silica surfaces mainly depend on the surface hydrophilicity (i.e., number of -OH groups), rather than surface crystallinity. We note that the roughness introduced in our model amorphous surface is in the atomic scale, different from the larger surface roughness on amorphous surface, such as pits, wedges, and pore. Sear and coworkers 31 demonstrated that the roughness in larger scale would cause a stronger influence on the liquid-solid nucleation. However, simulation on the amorphous surfaces gave the same trend on the effect of surface crystallinity, and showed again that the effect of the atomic scale roughness on hydrate nucleation is rather weak.
On the other hand, in comparison with the systems with crystalline silica surfaces (Fig. 5a), the CO 2 density on amorphous silica surface is kept slightly higher (Fig. 5b). More CO 2 molecules remained on amorphous surfaces in turn perturb the formation of the hydrogen bond network, and consequently, more free water molecules are present in the adjacent water layer. Figure 6 shows that the induction time for the hydrate nucleation is nearly independent on surface crystallinity because the local water structure is weakly affected by the surface crystalline. The growth rate of the hydrate after nucleation is also roughly the same for different systems (see Fig. 6). This is because the subsequent hydrate growth after nucleation is far from silica surfaces, and any surface effect is shielded by the initial hydrate layer formed. As a result, the growth rate of hydrate is mainly controlled by temperature, pressure and local concentration of CO 2 , rather than the properties of solid surfaces.
Although different systems have similar nucleation time and growth rate as long as the hydrophilicity of these surfaces is similar (Fig. 6), the degree of order for those hydrate crystals does depend on the surface crystallinity. Figure 7 shows the time evolution of the order parameter p n n n n sI sII sI sII 16 which was designed to describe the order of the hydrate crystal. Figure 7 shows that all hydrate crystals are dominated by the sI-like component, as indicated by the positive value of p. However, it is found that the hydrate induced by a crystalline silica surface is more ordered than by an amorphous one, even though the growth rates of those hydrates are similar. This observation follows from the fact that amorphous surfaces lead to a more disordered hydrogen bond network (Fig. 8), and hence the intermediate layer contains more free water molecules, which may reduce the order of the hydrate crystals.

Discussion
The change of nucleation mechanism (i.e. three-step nucleation or two-step nucleation) is mainly caused by the local structure of water close to the solid surface. To address the dependence of local water structure on surface hydrophilicity, we further analyzed in detail the structure of water molecules closest to the silica surface. Figure 8 shows that for the crys-100 system, an -OH group can adsorb on average about 2.8 water molecules, and most of these adsorbed water molecules are part of hexagonal or pentagonal water rings (identified by the ring perception algorithm 29 ). However, we note that although this layer contains several hexagonal rings, it not contains single-crystal motifs of any ice lattice; it is an amorphous layer with solid state. This is also the reason that the adsorbed water layer is called amorphous "ice-like".
When the hydrophilicity of the silica surfaces is continuously weakened from crys-100 to crys-0 systems, the average number of water molecules hydrogen bonded with silica (per -OH or -H group) decreases from ~3 to ~1, in which the fraction of water molecules forming 5-and 6-membered water rings decreases accordingly and that of free water molecules increases (Fig. 8). For example, for the Scientific RepoRts | 5:12747 | DOi: 10.1038/srep12747 crys-0 system (the least hydrophilic system), the number of water molecules hydrogen bonding to the silica surface significantly decreases. Only 1.1 water molecules are adsorbed to each -H group, and the percentage of free water molecules can be as large as 42% (Fig. 8). In this case, therefore, the mobility of water molecules increases, and consequently the ice-like water layer disappears (region A in Fig. 4a).
Consistent with the change of ϕ F 4 as a function. of surface hydrophilicity, the structure of the water layer near the surface is still relatively unchanged even though the number of -OH groups decreases to less than 50%. This is because even though the number of water molecules hydrogen bonding with the silica surface directly decrease with the decrease of -OH groups, the structure can still be stabilized alternatively by forming hydrogen bonds with its neighboring water molecules which are stabilized by hydrogen bonding with the -OH or -H groups. The MD simulation study on confined water 32 indicated  The adsorbed water molecules may come from hexagonal water rings, pentagonal rings, and free ones. Note that free water molecules here are the water molecules that do not belong to hexagons and pentagons, and they might appear as dimers, trimers, tretamers, or even rarely observed heptamers.
Scientific RepoRts | 5:12747 | DOi: 10.1038/srep12747 that the confining of water molecules on translation is stronger than that on rotation. It might be an explanation for the issue that why the local structure of water in the region close to the surface in the crys-50 system is similar to that in crys-100 system, rather than that in crys-0 system. This is because the rotation of water molecules can hold the network of hydrogen bonds undestroyed. But when the -OH groups decreases to 25% (crys-25) or less (crys-0), the ability of the silica surface to form hydrogen bond with water molecules significantly decreases so that the surface cannot stabilize the ice-like structure in the interfacial layer. Thus, the percentage of free water molecules significantly increases (Fig. 8). Similar results are obtained in confined silica-water binary system 33 . They found that the structure of water layer closest to the silica is more ordered and highly confined when silica surface has larger -OH density, i.e. more hydrophilic. In general, the surface hydrophilicity can change the local water structure, determining the presence or absence of the amorphous ice-like layer and inducing a change in the nucleation mechanism. Figure 9 gives the lifetime of 5 12 and 5 12 6 2 cages during the nucleation stages. The figure shows that hydrophilic surfaces prolong the lifetime of cages, indicating that the hydrophilic surfaces act as stabilizer for the adjacent cages, as suggested in our previous simulation results 16 . Although a cage hydrogen bonding with solid surfaces indeed increases its lifetime, a larger inaccessible region for hydrate cages near the solid surface, however, is existed for the cages formed on strong hydrophilic surfaces. We ascribe it to the structural mismatch between the hydrate lattice and the local structure of adjacent water layer, which leads to an increase of nucleation time (Fig. 6) and an ice-like water layer is formed firstly on strong hydrophilic surfaces.
The effect of structural mismatch on hydrate nucleation depends on the surface hydrophilicity. For strong hydrophilic surfaces, the structural mismatch between the hydrate lattice and amorphous ice-like layer inhibits the adsorbed hydrate cages growth into fully developed critical nucleus, and hence an intermediate layer has to form to bridge the structure mismatch. On the other hand, in the ice-like layer, the mobility of water molecules is strongly restricted by the silica surfaces, which also inhibits the hydrate nucleation and slows down its dynamics. But for weak hydrophilic surfaces, a liquid-like water layer closest to surfaces appears instead, and the restrictions on the adsorbed water molecules are substantially weakened. This can be confirmed from the diffusion coefficients of water listed in Table 1, which indicating that less hydrophilic surfaces pose less restriction in the mobility of the water molecules, and in consequence, accelerate the nucleation process. We should note that the diffusion coefficients of the ice-like or liquid-like layer close to the surface are heterogeneous: the diffusivity of water molecules in the xy-direction is stronger than that in the z-direction with respect to the solid surfaces. The anisotropy of water in confined space is a common phenomenon, and similar results are obtained for water between graphite surfaces 34 . In general, the hydrate nucleation tends to occur more easily on less hydrophilic silica surfaces, partly because of the stabilization effect of silica on hydrate cages, and partly because of the availability of free or less-constrained water molecules which are needed for the formation of hydrate cages. Moreover, the diffusion coefficients of the hydrate crystals formed in region C are in the range of 0.18 ~ 0.33 × 10 −9 m 2 /s for all of the systems, showing a solid feature.
The structure of water molecules at solid-liquid interfaces was also analyzed in detail, as shown in Fig. 8. In comparison with crystalline surfaces, amorphous surfaces induce an increase in the amount of adjacent free water molecules as well as a decrease in pentagonal water rings, while the amount of hexagonal water rings remains unchanged. This is partly due to the amorphous silica surfaces, as surface roughness induces a stronger structural perturbation on adjacent water rings and cages, leading to a more disordered hydrogen bond network. However, the total number of water molecules hydrogen bonding to -OH or -H groups does not significantly change. Consistently, Fig. 4b shows that there is no significant difference between ϕ F 4 for water molecules near crystalline and amorphous surfaces. In general, our simulations demonstrate that enhanced structural perturbation of amorphous surfaces on hydrate rings and cages leads to a slightly disordered hydrogen bond network when the surface hydrophilicity is similar.
The molar fraction of CO 2 in water we used is 0.042 in our simulation systems, which is a slightly supersaturated solution (see Method section for details). Since that some of the CO 2 molecules are remained near the solid surface, the effective degree of supersaturation in different systems will be different, and it is another important factor to affect the hydrate formation rate and the mechanism. At ~0.8 μs for all of the systems, the expulsion of CO 2 molecules from the region close to the silica surfaces is finished (Fig. 5) and it is ready to start the formation of hydrate (Fig. 6). Hence, we calculated the CO 2 density in the border region (region A and B) and the central region (region C and D) at 0.8 μs, respectively. Take crys-0 and crys-100 systems for example, the density in border region is 0.775 molecules/nm 3 for crys-0 system and 0.042 molecules/nm 3 for another, with a difference of 18.45 times (0.775 divided by 0.042). The density in central region, however, is 1.15 and 1.47 molecules/nm 3 for the two systems respectively, only have a difference of 1.28 times. Compared with the former, the density difference of CO 2 in central region between different systems is negligible. Therefore, we point out that the obvious difference in the induction time for hydrate nucleation occurred near the surface (Fig. 6) is owing to both the local structure of water and the density of CO 2 are significantly different for systems with different silica surfaces. However, the hydrate growth is mainly occurred in the central region, and the rate for all of the systems is roughly the same. This is caused by the similar CO 2 density in central region and the weakened surface effect. We note that although it is a superstaturated system, the main factor to affect the hydrate nucleating is the properties of water and CO 2 close to the solid surface, rather than that in the central region of the system. Since the CO 2 molecules are all expulsed to some extent in different systems, the local degree of supersaturation of CO 2 in border region well be further reduced. So, the main results we obtained can reflect the real situation of the system in a certain extent. In real system, perhaps the local structuring of water and the expulsing of CO 2 near the surface become less obvious owing to the low CO 2 concentration, but the nucleation mechanism will not change.
In summary, we performed microsecond MD simulations to investigate how the presence of hydrophilic silica surfaces modulates the hydrate nucleation. The hydrophilicity of solid surfaces can change the local structure of adjacent water layers: the nucleation mechanism varies from three steps into two steps when surface hydrophilicity becomes weak. It also affects the nucleation dynamics of the hydrate formation: the induction time for nucleation is reduced when the surface is less hydrophilic. While for strong hydrophilic surfaces, the structure mismatch between the hydrate lattice and amorphous ice-like layer inhibits the nucleation and then the nucleation dynamics is slowed down. The crystallinity of solid surfaces affects the hydrate formation process in a rather weak manner. The hydrates induced by crystalline silica are more ordered than by amorphous silica. Owing to the fact that hydrate growth after nucleation is far from the silica surface, the growth rates of the hydrate seem to be independent on both surface hydrophilicity and surface crystallinity.

Methods
Molecular dynamics simulations were performed by using LAMMPS 35 . The simulations were implemented in a two-phase system that contains a pair of silica layers (solid phase) and a mixed H 2 O/CO 2 fluid (liquid phase) with 200 CO 2 and 4600 H 2 O molecules, as is shown in Fig. 1. In the initial configurations, the CO 2 and H 2 O molecules were randomly placed in a simulation box of 6.22 nm × 6.22 nm × 4.98 nm (roughly equals to 5 × 5 × 4 unit cells of sI hydrate). At the top and bottom of the simulation box, two silica layers were added to represent a slit pore of sedimentary rock. Due to the lack of the solubility of CO 2 solution within a confined space in our simulation conditions, we extrapolated the solubility of bulk CO 2 solution to our conditions based on the experimental data [36][37][38] , which is ~0.037 of CO 2 molar fraction in water. In our simulation systems, the molar fraction of CO 2 is 0.042, showing a supersaturated solution. Considering the confining effect, we estimate that the degree of supersaturation of CO 2 used in the simulations is not higher than 4 times of the real system at the same conditions. In all of the simulations, the positions of SiO 2 molecules were fixed and periodic boundary conditions were imposed in all three Cartesian directions.
The CO 2 molecules were represented by the EPM2 model 39 , which has three Lennard-Jones sites with charges centered at each atom and with rigid bond lengths but a harmonic bond angle. The H 2 O molecules were described by the extensively used TIP4P model 40 , in which the rigidity of H 2 O molecules was restricted by the SHAKE algorithm 41 . The silanols/silanes model 42 was adopted for the SiO 2 layers. The Lennard-Jones interaction parameters for molecules are reported elsewhere 15 . Note that the force fields used were successfully applied to study the nucleation of CO 2 hydrates in both two-phase 15 and three-phase systems 16 . A cutoff radius of 12.0 Å was employed for the short-ranged interactions, and the PPPM algorithm 35 was used for long-ranged electrostatic interactions. Constant temperature and pressure simulations were maintained at 265 K and 150 bar, respectively, with the Nosé-Hoover algorithm [43][44][45] . According to experimental phase diagram 1 , under the conditions the hydrate phase is thermodynamic stable. The Newtonian equations of motion were integrated based on the velocity Verlet algorithm 46 with a time step of 2 fs.
In order to investigate the effect of surface properties on the hydrate nucleation process, we introduced several model solid surfaces with different hydrophilicity and crystallinity, as summarized in Fig. 3. The crystalline silica layers were taken from the (1 1 1) plane of the β-cristobalite 47 , while the amorphous silica layers with a roughness down to atomic scale were built by using in total 1060 silicon and 2120 oxygen atoms with a density closing to the experimental value of 2.2 g/cm 3 . For different model SiO 2 surfaces, dangling bonds on the surfaces were saturated by -OH or -H groups. The -OH group has no rotational freedom around the Si-O bond. Obviously, a silica surface saturated by an -OH group alone (hydroxylated model) is more hydrophilic than that saturated by an -H group (hydrogenized model) because an -OH group covering the silica surface can be associated in hydrogen bonds with three other water molecules. While an -H group can hydrogen bond to a single water molecule at most.
To generate different surface hydrophilicity, the usual method controlling the surface density of the -OH group is to select two locations on the (1 1 1) plane of the β-cristobalite 33 or to select different crystallographic faces of the crystal. When using this method to change the hydrophilicity of solid surfaces, however, the surface structure is also changed. To separate the effects of surface hydrophilicity and surface structure, we introduced another method in which a number of dangling bonds on silica surface are chosen randomly to be saturated with -H groups, to control the surface hydrophilicity. Hence, in this method the variation of the hydrophobicity of a model silica surface does not necessarily change the surface structure.
A typical MD simulation was divided into two steps. First, an NpT relaxation process of 2 ns was performed at 265 K and 150 bar to eliminate the effect of the initial configuration. Within the first 1 ns, the oxygen atoms of water molecules were fixed to minimize the total dipole moment of water molecules; then the restriction on the oxygen atoms was gradually removed and the system continued to relax for another 1 ns. A typical configuration after the relaxation process is shown in Fig. 1. The configuration is used as the initial configuration of next step. In the second step, hydrate nucleation and growth simulations were performed at the same NpT conditions to 4 μs in total time.