Molecular Simulations on the Structure and Dynamics of Water

Jun 12, 2013 - The structure and dynamics of water–methane fluids between clay surfaces are investigated through the grand-canonical Monte Carlo (GC...
0 downloads 7 Views 4MB Size
Download PDF
Article pubs.acs.org/JPCC

Molecular Simulations on the Structure and Dynamics of Water− Methane Fluids between Na-Montmorillonite Clay Surfaces at Elevated Temperature and Pressure Qi Rao, Yuan Xiang, and Yongsheng Leng* Department of Mechanical and Aerospace Engineering, The George Washington University, Washington, D.C. 20052, United States ABSTRACT: The structure and dynamics of water−methane fluids between clay surfaces are investigated through the grand-canonical Monte Carlo (GCMC) and molecular dynamics (MD) simulations. The chemical potentials of water and methane at the temperatures and pressures corresponding to different burial depths are calculated. These chemical potentials are used in the GCMC simulations to determine the water and methane contents in the clay interlayer at a burial depth of 6 km. The results are used as initial inputs for further MD simulations to investigate the static and dynamic properties of the confined fluid. Simulation results show that initial clay swelling is dominated by water adsorption into clay interlayer, followed by the formation of methane hydrates as the basal spacing increases. Methane content in clay is found to increase in a step fashion, from initial inner-sphere complex to both inner- and outersphere structures. It is found that methane is not fully coordinated by water molecules due to the low density of water content in Na-montmorillonite clay.

1. INTRODUCTION The properties of aqueous and aqueous−hydrocarbon fluids between charged surfaces, particularly between aluminum− silicate clay minerals, are of considerable interest and importance in surface and colloid science1,2 and are directly related to many practical industrial and geological applications.3 The underlying mechanisms of the surface chemistry process involve solvation of interlayer cations by water molecules, ion exchange between different species, and, under the circumstances of organic molecules being intercalated into pore spaces,4,5 the competition of surface sorption between solvated ions, water molecules, and hydrocarbon species. These properties are particularly crucial to the oil and gas production in petroleum industry. Examples include fundamental understanding of the location and primary and tertiary migration of petroleum hydrocarbon,6,7 clay swelling,8,9 and well bore stability in oil-well drilling industries.10 Aqueous and aqueous−organic interlayer fluids between clay minerals are usually subjected to a variety of thermodynamic conditions, ranging from ambient to high temperature and pressure (T/P) sedimentary basin conditions. In oil recovery industry, many oil reservoirs contain a high fraction of clay minerals with swelling properties. The uptake of water molecules into clay interlayer space from surrounding micropores11 results in clay swelling and the formation of damage of oil well bores.9,12 This situation becomes even more severe if improper stimulation operations apply, such as hydraulic stimulation for naturally fractured reservoirs.13 A fundamental question concerns how the variations of thermodynamic conditions and stimulations impact on the structure and dynamics of aqueous and aqueous−hydrocarbon pore fluids between charged clay surfaces.14−17 © XXXX American Chemical Society

The knowledge of structure and dynamics of aqueous hydrates between charged clay surfaces has been advanced substantially in the past decade through molecular simulations.3,9,11,12,15,18−26 Experimental measurements using different approaches provide first-hand information on the adsorption structure and ion exchange mechanism in aqueous systems.27−31 Recent experiments4,5 also revealed that methane and other small organic hydrocarbon molecules can penetrate into clay interlayer to form complex adsorption structures, which attracted considerable interests in theoretical modeling.3,16,17,32−34 The latter is particularly related to the oil and gas recovery since the uptake and retention of hydrocarbon molecules in clay hydrates may significantly alter the properties of the clay matrix. Molecular modeling for the structure and dynamics of the water−methane mixture between charged clay surfaces has just begun recently, owing to the importance of the role of organic molecules intercalated into clay hydrates. The structure and dynamics of methane in hydrated sodium montmorillonite clays under different T/P conditions were investigated by Titiloye and Skipper.16 Park and Sposito32 found clay mineral surfaces have “thermodynamic promotion effect” on methane hydrate formation. Cygan et al.33 investigated the structure and dynamics of methane and water in the interlayer of montmorillonite clays through molecular dynamics simulations. Different cations such as Na+, K+, Ca2+, and Mg2+ were considered in the simulations. The structure and dynamics of methane in hydrated potassium montmorillonite clays under Received: April 4, 2013 Revised: June 3, 2013

A

dx.doi.org/10.1021/jp403349p | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

different T/P conditions were investigated by Titiloye and Skipper.17 Zhang and Choi35 also investigated the structure of methane, water, and cation in hydrated potassium montmorillonite clays via molecular dynamics simulation. Zhou et al.36 used molecular dynamics simulation to investigate the behavior of methane hydration in Na-smectite interlayer with different layer-charge distribution. However, in previous studies it is not clear how the numbers of water and methane molecules in the clay slit pore were determined under different T/P conditions, and in many cases a priori determined sets of methane and water molecules were used in the simulation. Different results may be obtained due to the different initial conditions. Therefore, it is important to find a reasonable way to determine the numbers of methane and water molecules in clay interlayer. The present paper will focus on the clay swelling mechanism in water−methane mixture and the adsorption structure of water, ion, and methane confined between charged clay surfaces. Specifically, the grand-canonical ensemble Monte Carlo (GCMC) and molecular dynamics (MD) simulations are performed to understand the clay swelling mechanism in an aqueous−hydrocarbon environment from a dehydrated state.

f pure (Tmix , Ppure) = fmix (Tmix , Pmix , x)

(1)

where f pure is the fugacity for pure methane or water, f mix is the fugacity for methane or water in the mixture, Tmix is the temperature of the mixture, Pmix is the pressure of the mixture, and x is the mole fraction for methane or water in the mixture. The fugacity is a function of temperature, pressure, and mole fraction and can be determined from EOS. Since μpure = μ° + RT ln f pure

(2)

μmix = μ° + RT ln fmix

(3)

where μpure is the chemical potential for pure methane or water, μmix is the chemical potential for methane or water in the mixture, and μ° is the standard chemical potential for methane or water. From eqs 1−3, we have μpure (Tmix , Ppure) = μmix (Tmix , Pmix , x)

(4)

The chemical potentials of the species in the mixture can be obtained by calculating the chemical potentials at the corresponding pure component pressures Ppure. These pressures Ppure can be used as input data in the Widom’s insertion Monte Carlo simulations to determine the chemical potentials. In order to calculate Ppure, the Peng−Robinson EOS equation48 is used in an iterative manner to calculate the mole fractions and fugacities of the species. The EOS equation takes the form

2. MOLECULAR MODELS AND SIMULATION METHODS The GCMC and constant-NVT MD simulations are performed using the computational packages Towhee 37,38 and LAMMPS,39 respectively. The model of clay mineral is the sodium-saturated Wyoming-type montmorillonite that has the chemical formula Na0.75[Si7.75Al0.25](Al3.5Mg0.5)O20(OH)4. The clay surfaces contain eight clay unit cells forming a clay patch of 20.77 × 18.03 Å2. The realistic CLAYFF20 force field for clay minerals is used in this work. The potential model used for water molecules is the SPC40 model and methane is represented by the Jorgensen’s OPLS all atom (OPLS-AA)41 model. The GCMC run usually takes at least 4 million (4M) moves for rigid water and methane molecules. In MD simulations the intramolecular bond stretch and bond angle terms are introduced into the SPC water model,42 as well as into the methane molecule,41 to ensure full flexibility of these molecules and hydroxide components on a clay surface.20 A cutoff distance of 9.0 Å is used for the dispersive interactions, while the long-range electrostatic interactions are treated by the Ewald summation method.43 Prior to the GCMC simulations, the chemical potentials of water and methane are computed using the Widom’s insertion method in the constant-NPT ensemble Monte Carlo simulations.44,45 In this ensemble, the temperature was set to the temperature of mixture. Here we only consider a simple situation where the mixture contains only water and methane, and the mixture reaches a vapor−liquid equilibrium. A simple method (method I) to determine the pressure for each component is to set the pressure of water equal to its saturated vapor pressure at the temperature of mixture, while the pressure of methane is equal to the rest of the total pressure.46 However, a more suitable method (method II) to determine the pressures for each component in the mixture is to use an equation of state (EOS) to compute the fugacities of species in the mixture and to find the pressures with which pure methane or water has the same fugacity47 with those in the mixture, i.e., find the pressure Ppure for pure methane or water to let

P=

RT aα − 2 Vm − b Vm + 2bVm − b2

(5)

where a, b, and α are constants in the equation, R is the universal gas constant, and Vm is the molar volume. In the present study, we use method I to calculate the chemical potentials due to its simplicity. We also use method II to calculate the chemical potentials under specific T/P conditions to compare with the results obtained from method I. The Widom’s insertion method 44 is based on the thermodynamics definition of chemical potential, which includes insertions of a test particle into the N-particle simulation system at the desired temperature and pressure. The interaction energy of the virtual test particle with the Nparticle system during the insertion moves is used to calculate the chemical potential of the fluid. The chemical potential can be computed by45 48

βμ = ln(Λ3βP) − ln

βPV N+1

∫ dsN +1 exp(−βΔU ) (6)

where μ is the chemical potential, β is the inverse of kT, Λ is the de Broglie wavelength, V is the volume of simulation box, sN is the scaled coordinates, and ΔU is the system energy change induced by the insertion of test particle. The swelling behavior of Na-montmorillonite was investigated by GCMC simulations where the temperature, volume, and chemical potential were fixed. This ensemble allows us to determine actually how many water and methane molecules that can enter the clay interlayer at specific T/P conditions. During the GCMC simulations, the fluid confined between clay mineral surfaces can exchange water and methane molecules with a reservoir of water−methane mixture, which has fixed chemical potentials for water and methane. To investigate the dynamics behavior of confined fluid, the interlayer molecular configurations obtained from GCMC B

dx.doi.org/10.1021/jp403349p | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

different T/P conditions. The results are compared with those from a semiempirical equation,51 as shown in Figure 2.

simulations were chosen as the initial configurations for the subsequent MD simulation runs. The simulations were performed for a total of 5 ns production runs for the data collection to obtain the interlayer structural and dynamic properties.

3. RESULTS AND DISCUSSION 3.1. Chemical Potential. The chemical potential of water at 298 K and 1 bar pressure is found to be −45.56 kJ/mol. This value is approximately equal to the value for the TIP4P water (−44.5 kJ/mol).49 Using this chemical potential value, we first performed GCMC simulation to study the swelling behavior of Na-montmorillonite under saturated vapor pressure of water. As shown in Figure 1, the swelling curve obtained from the GCMC simulations is in agreement with experimental results50 as well as with other MD simulation results.20 Figure 2. Variations of mole fractions of water and methane in the mixture versus T/P in gas phase.

The mole fraction of methane in our calculation is slightly larger than that from Duan’s semiempirical equation, and correspondingly, the mole fraction of water in our calculation is smaller than that from Duan’s equation. The difference increases as the temperature increases. However, the trends of the variations of mole fractions for both species versus temperature are the same. Therefore, the chemical potentials obtained from the constant-NPT Monte Carlo simulations are reasonable for the subsequent GCMC simulations. 3.2. Grand-Canonical Monte Carlo Simulation. With the confidence of the obtained appropriate chemical potentials for water and methane, we have performed the GCMC simulations to study the clay swelling in water−methane mixture at T = 460 K and P = 900 bar, which correspond to the temperature and pressure at the burial depth of 6 km. This case has been studied previously through MD simulations by Titiloye and Skipper.16 The chemical potentials of water and methane are given in Table 1. A series of GCMC simulations involving at least 4 million moves of molecules are performed with various basal spacing from 10.5 to 24 Å. These GCMC runs are followed by 12−16 million continuing runs to get statistic values for water and methane contents. The average numbers of water molecules and methane molecules entering the clay interlayer are shown in Figure 3. The swelling curves in Figure 1 are also displayed in Figure 3 for comparison. Moreover, we have also calculated the average numbers of water molecules entering the clay interlayer without methane at the same T/P condition, as shown in the figure. Evidently, the water content at this high T/P condition is significantly reduced. When methane is present, we find that water

Figure 1. Swelling curves for the hydration of Na-montmorillonite obtained from different methods.

The constant-NPT ensemble Monte Carlo simulation results for the chemical potentials of water and methane under different T/P conditions are listed in Table 1. The data in the parentheses are chemical potentials for water and methane obtained from method II. It is seen that water pressures Ppure obtained from the Peng−Robinson EOS equation48 are higher than the saturated vapor pressure of water except at burial depth of 1 km. The difference reaches the largest at 6 km. However, the difference of chemical potentials at this T/P condition is only about 1%. Therefore, we employ the chemical potentials calculated from method I, i.e., assuming water pressure reaches its saturated pressure. To further evaluate the suitability of chemical potentials that will be used in GCMC simulations, we have calculated the mole fractions of water−methane mixture in gas phase under

Table 1. Chemical Potentials of Water and Methane in the Mixture under Different T/P Conditions (Data in Parentheses Are Calculated from Method II) T, P, and depth 310 340 370 400 430 460 a

K, K, K, K, K, K,

150 300 450 600 750 900

bar, bar, bar, bar, bar, bar,

1 2 3 4 5 6

km km km km km km

saturated pressure of water (kPa)a 6.23 (5.93) 27.19 (30.51) 90.54 (116.55) 245.78 (356.25) 570.27 (921.81) 1170.90 (2110.50)

chemical potential of water (kJ/mol) −45.9 −46.8 −47.7 −48.8 −50.2 −51.8

± ± ± ± ± ±

0.04 0.1 0.3 0.7 0.8 2.7 (−52.4 ± 3.2)

chemical potential of methane (kJ/mol) −24.0 −24.9 −26.2 −27.6 −29.0 −30.7

± ± ± ± ± ±

0.07 0.1 0.2 0.2 0.2 0.2 (−31.0 ± 0.2)

See ref 52. C

dx.doi.org/10.1021/jp403349p | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

Figure 4. Density profiles of methane molecules (represented by the carbon atom Cm in methane) in the interlayer region. The origin corresponds to the clay surface oxygen. Figure 3. Average numbers of water and methane molecules entering the clay interlayer.

molecules enter clay interlayer first, while methane molecules cannot enter clay interlayer until the basal spacing reaches 12 Å. The number of water molecules entering clay interlayer at this burial depth is smaller than that at ground condition. This difference is becoming larger as the basal spacing is beyond 12 Å. Further, methane content increases in a step fashion. As the basal spacing increases, more methane hydrates are formed during clay swelling. In Titiloye and Skipper’s work,16 four methane molecules were put into the interlayer hydration film at different clay spacings, resulting in high methane contents at smaller spacings (such as one-, two-, and three-layer hydrate systems) and low methane contents at larger spacings (four-layer hydrate and above). The present GCMC simulation results show that the increase in basal spacing generally allows more methane hydrates to form in the clay interlayer. 3.3. Molecular Dynamics Simulation. Molecular dynamics simulations based on the GCMC equilibrium simulation results at burial depth of 6 km are further performed to investigate the structure and dynamics of methane hydrates. Five basal spacings are chosen in the simulations. The basal spacings and the corresponding GCMC results for water and methane molecules are listed in Table 2.

Figure 5. Density profiles of water molecules (represented by water oxygen Ow) in the interlayer region. The dashed lines show the density of water hydrogen atoms in different spacings. The origin corresponds to the clay surface oxygen.

surfaces. As the clay spacing further increases to 20 and 24 Å, some methane molecules begin to migrate into the midplane area to form outer-sphere hydrates. The overall density profiles of methane are similar to those in Titiloye and Skipper’s work,16 but all of the peaks of the profiles except 24 Å spacing are lower than their results, simply because less methane molecules are present in the clay interlayer. The density profiles of water in the interlayer region are displayed in Figure 5. When the basal spacing is 13 Å, there are essentially two small oxygen peaks close to each other, namely the one-to-two-layer-transition hydrate confined between clay surfaces. This is also indicated by the hydrogen density distribution where two small peaks are around the central large peak. When the basal spacing increases to 16 Å, two layers of hydration water are clearly seen. There are four peaks indicating four layers of water forming when the basal spacing reaches 20 Å. When the basal spacing increases to 24 Å, water density in the central area becomes flattened, with a value lower than the bulk value of water (0.0311 Å−3)52 at 460 K and 900 bar. This is partly due to the volume exclusive effect of methane in the central region. Figure 6 shows the density profiles of sodium cations in the interlayer region. The distribution curves are quite asymmetric due to the unevenly distributed surface charges on clay surfaces. The sodium ions are in the inner-sphere3 coordination at the basal spacing of 13 Å, since all the Na+ peaks are close to the

Table 2. Numbers of Water and Methane Molecules in Different Basal-Spacing Methane Hydrates at T = 460 K and P = 900 bar basal spacing (Å)

CH4

H2O

13 16 18 20 24

2 2 3 3 7

36 62 74 96 128

Methane Hydrate Structure. The density profiles of methane in the interlayer region are presented in Figure 4. The one-layer central peak shows that the methane molecules are located at the midplane confined tightly by clay surfaces. Two peaks appear in the 16 and 18 Å spacings, corresponding to two and three water hydration layers shown in Figure 5. This indicates that methane molecules have affinity with the clay D

dx.doi.org/10.1021/jp403349p | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

coordination numbers are listed in Table 3, which are calculated by integrating the RDF to the first minimum. The snapshots of Table 3. Coordination Number of Water around Methane at T = 460 K and P = 900 bar

a

Figure 6. Density profiles of sodium ions in the interlayer region. The origin corresponds to the clay surface oxygen.

basal spacing (Å)

nCH4−Ow

nCH4−Os

nt a

13 16 18 20 24 bulk water

6.7 11.2 11.3 12.4 12.3 19.9

11.7 5.2 3.9 3.3 2.2

18.4 16.4 15.2 15.7 14.5 19.9

nt is the total coordination number of oxygen around methane.

the structure of the methane hydrate intercalate with Namontmorillonite are displayed in Figure 8. We also calculate the RDF for CH4−Ow and the coordination number in bulk water for comparison, as shown in Figure 7a and Table 3. The calculated coordination number of bulk water molecules around methane at T = 460 K and P = 900 bar is found around 19.9. This value is consistent with previous results under a variety of simulation conditions (i.e., between 19 and 23).53 When the basal spacing is 13 Å, the CH4−Ow RDF peak is at a smaller distance of 3.5 Å, and the first minimum is located at 5.5 Å. These two values are typical distances associated with the RDF for methane in bulk water,53,54 as shown in Figure 7a. At the same basal spacing, the CH4−Os RDF peak is located at 3.8

clay surface. For all other methane hydrates, the nearest sodium peak close to the clay surface shows the inner-sphere complex, while the peaks in the central regions show the fully hydrated outer-sphere3 complexes. These distributions at different methane hydrates are similar to the work by Titiloye and Skipper.16 In order to further understand the detailed methane hydrate structure, we have calculated the radial distribution functions (RDFs) for methane−water (CH4−Ow) and methane−surface oxygen (CH4−Os) interaction pairs. The results are presented in Figure 7a,b. The partial CH4−Ow and CH4−Os and the total

Figure 7. Radial distribution functions for (a) CH4−Ow, (b) CH4−Os, (c) Na+−Ow, and (d) Na+−Os interaction pairs at T = 460 K and P = 900 bar. The functions are normalized by the bulk water density at the same T/P condition. E

dx.doi.org/10.1021/jp403349p | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

Figure 8. Snapshots of the methane hydrates in Na-montmorillonite interlayer at different basal spacing of (a) 13 Å (36H2O, 2CH4), (b) 18 Å (74H2O, 3CH4), and (c) 24 Å (128H2O, 7CH4). Color scheme: oxygen = red, hydrogen = white, carbon = brown, silicon = blue, aluminum = cyan, magnesium = orange, and sodium = yellow.

Å and the first minimum is at 5.0 Å. The total coordination number of oxygen around methane molecules is 18.4, with 6.7 of these arising from water oxygen and 11.7 from surface oxygen (see Table 3). These results are similar to those given by Titiloye and Skipper.16 It is clear that the methane hydrate is an inner-sphere complex.3 As the thickness of methane hydrates increases, the CH4−Ow RDF peaks gradually shift from 3.5 to 3.7 Å, while the first minima of these RDFs are essentially at the same distance of 5.5 Å. These RDFs gradually approach to the bulk RDF, indicating that more water molecules come into the methane first hydration shell, but nevertheless, they are still much lower than the bulk curve. This can be also seen from the gradual decreasing of RDFs for the CH4−Os pair correlations as the interlayer spacing increases. This trend also shows that the structure of methane hydrate has a tendency to change from an inner-sphere complex to an outer-sphere complex. Figure 7b also shows that the first RDF peaks remain at 3.8 Å, but the minima gradually shrink from 5.0 Å to a distance of 4.5 Å. In the two-layer hydrate, the total coordination number of oxygen around methane molecule is reduced to 16.4. Compared to the coordination number of methane in bulk water (N = 19.9, see Table 3), it is clear that methane molecules are not fully coordinated. This trend becomes even more severe when the clay spacing is further increased. Table 3 shows that the total coordination numbers of oxygen around methane molecules are further decreased from 16.4 to 14.5. These results, which are derived from combined GCMC and MD simulations, are 2−3 water molecules less than the coordination numbers given by Titiloye and Skipper.16 It suggests that methane under this T/P condition is not fully hydrated. The main reason is that the actual number of water molecules given by GCMC simulation or the density of water in methane hydrates is significantly reduced due to the high T/ P condition and the presence of methane in the clay interlayer. This is in contrast to the results shown in Figure 1, where water molecules can easily come into the interlayer at the ground level and when methane is not present.

As shown in Figure 3, water content in clay interlayer is strongly influenced by the T/P conditions. Our simulation results show that water density in the interlayer decreases with the increase of temperature and pressure. This trend is consistent with previous simulation work.55 Also, the T/P condition has stronger influence on the water density than the effect of methane. For example, at basal spacing of 16 Å, the water number drops from 84.1 at T = 298 K and P = 1 bar to 66.9 at T = 460 K and P = 900 bar. If methane is present, the water number further drops to 63.9. The differences caused by both T/P condition and the presence of methane increase as the basal spacing is becoming larger. The RDFs for sodium−water (Na+−Ow) and sodium− surface oxygen (Na+−Os) interaction pairs are presented in Figure 7c,d. For comparison, the RDF for Na+ in bulk water under the same T/P condition is also shown in Figure 7c. The positions of the first peak for Na+−Ow and Na+−Os are the same in all the cases. The total oxygen coordination numbers and the partition of Ow and Os around Na+ are shown in Table 4. It is seen that the total O coordination number for sodium is very close to the bulk coordination number, indicating that sodium ions are fully coordinated. Diffusion Behaviors of Different Species in Methane Hydrates. To understand the dynamics behavior of methane hydrate in clay interlayer, we calculate the self-diffusion Table 4. Coordination Number of Water around Sodium at T = 460 K and P = 900 bar

a

F

basal spacing (Å)

nNa+−Ow

nNa+−Os

nta

13 16 18 20 24 bulk water

3.4 4.3 4.2 4.3 4.2 5.6

2.4 1.7 1.8 1.7 1.7

5.8 6.0 6.0 6.0 5.9 5.6

nt is the total coordination number of oxygen around sodium. dx.doi.org/10.1021/jp403349p | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

4. CONCLUSIONS This work investigates the structure and dynamics of water− methane fluid between clay surfaces through GCMC and MD simulations. The chemical potentials of methane and water in mixture at different T/P conditions are calculated by the Widom’s insertion method and applied to the GCMC simulations. We specifically focus on the GCMC simulations to determine the water and methane contents in Namontmorillonite clay at a burial depth of 6 km. Simulation results show that initial clay swelling is dominated by water adsorption into clay interlayer. As the basal spacing increases, methane molecules begin to enter the interlayer to form methane hydrates. This process continues in a step fashion as the clay spacing increases. Simulation results show that water density in clay interlayer is lower than the bulk value at the T/P condition studied, while methane molecules are not fully hydrated. The overall water coordination numbers around methane at different spacings are less than the results given by Titiloye and Skipper’s studies,16 in which a fixed number of methane particles were assumed at different methane hydrates. As the basal spacing increases, the diffusion of interlayer species increases and the hydrated methane complexes begin to include inner-sphere and outer-sphere models. The present study may provide a reasonable way to understand the clay swelling mechanism in aqueous−hydrocarbon fluids.

coefficients D for different species, which are determined from the Einstein relation ⟨|r(t ) − r(0)|2 ⟩ = 2dDt

(7)

where r(t) is the position of the particle at the time t and d is the dimensionality of space in which the diffusion is considered. Here d = 2 for interlayer since we only consider the diffusion in the plane parallel to the clay surface, and ⟨|r(t) − r(0)|2⟩ is the mean-square displacement (MSD). The diffusion coefficients of methane, water, and sodium ions at 460 K and 900 bar are listed in Table 5. As the basal spacing Table 5. The 2D-Diffusion Coefficients of Methane, Water and Sodium Ions at T = 460 K and P = 900 bara basal spacing (Å) 13 16 18 20 24

(13.08) (15.82) (19.1) (21.4)

DCH4 (×10−9 m2/s)

DH2O (×10−9 m2/s)

2.32 (2.12) 9.74 (4.26) 11.38 (8.41) 11.46 (8.64) 12.82

3.72 (2.92) 9.92 (5.82) 13.31 (12.65) 14.87 (12.01) 15.38

DNa (×10−9 m2/s) 0.64 1.75 2.36 2.79 3.23

(0.10) (1.34) (3.82) (3.10)

a

The data in the parentheses are 3D-diffusion constants calculated by Titiloye and Skipper16 through MD simulations under the same T/P condition.

■

increases, the overall diffusion coefficients increase. These results are larger than those given by Park and Sposito at ground level32 due to the high T/P in this study. For example, at basal spacing of 18 Å the diffusion coefficients in our work for methane, water, and sodium species are 4.7, 2.7, and 16.8 times larger than those at ground conditions.32 Our calculation results have roughly the same orders of magnitude as those in Titiloye and Skipper’s work under the same T/P condition.16 The detailed number differences (see Table 5) are attributed to the methods the diffusion coefficients are calculated. In our work we use eq 7 to calculate 2D diffusion coefficients, but in Titiloye and Skipper’s work, they calculated 3D diffusion coefficients.16 Since the diffusion of molecular species in the normal direction of the clay surface is much smaller than those in the lateral directions, the 3D diffusion coefficients are generally less than the 2D coefficients. Table 5 also shows that sodium ions have lower mobility compared with water and methane, and its overall diffusion coefficients at different interlayer spacings increase at a slower rate than those of other two species. This is because the strong electrostatic interaction between Na+ and the charged clay surface significantly slows down the diffusion of metal ions. This slowing down behavior was also observed for K+ ions confined between mica surfaces in aqueous solution.56 A final check on the diffusion coefficients shown in Table 5 finds that at high T/P condition studied in this work, as the basal spacing is increased from 13 to 24 Å (corresponding to one-to-two-layer-transition state to a very thick hydrate), the mobility of methane and water molecules and Na+ ions in Namontmorillonite is increased by only about 5 times. However, in highly charged mica clay at ambient conditions, the difference of diffusion for different species in similar spacings could be as high as 2 orders of magnitude.56 Consequently, the T/P condition and charge density of clay minerals may have profound effects on the diffusion of different species in interlayer hydrate.

AUTHOR INFORMATION

Corresponding Author

*Tel 202-994-5964; e-mail [email protected] (Y.L.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is supported by the American Chemical Society Petroleum Research Fund (PRF# 49596-DNI 5) and the National Energy Research Scientific Computing Center (NERSC).

■

REFERENCES

(1) Israelachvili, J. N. Intermolecular and Surface Forces, 2nd ed.; Academic Press: New York, 1991. (2) Hunter, R. J. Foundations of Colloid Science, 2nd ed.; Oxford University Press: New York, 2001. (3) Sposito, G.; Skipper, N. T.; Sutton, R.; Park, S. H.; Soper, A. K.; Greathouse, J. A. Surface geochemistry of the clay minerals. Proc. Natl. Acad. Sci. U. S. A. 1999, 96 (7), 3358−3364. (4) Guggenheim, S.; van Groos, A. F. K. New gas-hydrate phase: Synthesis and stability of clay-methane hydrate intercalate. Geology 2003, 31 (7), 653−656. (5) Cosultchi, A.; Cordova, I.; Valenzuela, M. A.; Acosta, D. R.; Bosch, P.; Lara, V. H. Adsorption of crude oil on Na+-montmorillonite. Energy Fuels 2005, 19 (4), 1417−1424. (6) Newman, A. C. D. Chemistry of Clays and Clay Minerals; Mineralogical Society: London, 1987. (7) North, F. K. Petroleum Geology; Unwin Hyman: Boston, 1990. (8) Hensen, E. J. M.; Smit, B. Why clays swell. J. Phys. Chem. B 2002, 106 (49), 12664−12667. (9) Skipper, N. T.; Lock, P. A.; Titiloye, J. O.; Swenson, J.; Mirza, Z. A.; Howells, W. S.; Fernandez-Alonso, F. The structure and dynamics of 2-dimensional fluids in swelling clays. Chem. Geol. 2006, 230 (3−4), 182−196. (10) Cuss, R. J.; Rutter, E. H.; Holloway, R. F. Experimental observations of the mechanics of borehole failure in porous sandstone. Int. J. Rock Mech. Min. Sci. 2003, 40 (5), 747−761. G

dx.doi.org/10.1021/jp403349p | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

(32) Park, S. H.; Sposito, G. Do montmorillonite surfaces promote methane hydrate formation? Monte Carlo and molecular dynamics simulations. J. Phys. Chem. B 2003, 107 (10), 2281−2290. (33) Cygan, R. T.; Guggenheim, S.; van Groos, A. F. K. Molecular models for the intercalation of methane hydrate complexes in montmorillonite clay. J. Phys. Chem. B 2004, 108 (39), 15141−15149. (34) Odriozola, G.; Aguilar, J. F.; Lopez-Lemus, J. Na-montmorillonite hydrates under ethane rich reservoirs: NPzzT and mu PzzT simulations. J. Chem. Phys. 2004, 121 (9), 4266−4275. (35) Zhang, J. F.; Choi, S. K. Molecular dynamics simulation of methane in potassium montmorillonite clay hydrates. J. Phys. B: At., Mol. Opt. Phys. 2006, 39 (18), 3839−3848. (36) Zhou, Q.; Lu, X. C.; Liu, X. D.; Zhang, L. H.; He, H. P.; Zhu, J. X.; Yuan, P. Hydration of methane intercalated in Na-smectites with distinct layer charge: Insights from molecular simulations. J. Colloid Interface Sci. 2011, 355 (1), 237−242. (37) Martin, M. G.; Siepmann, J. I. Transferable potentials for phase equilibria. 1. United-atom description of n-alkanes. J. Phys. Chem. B 1998, 102 (14), 2569−2577. (38) Martin, M. G.; Siepmann, J. I. Novel configurational-bias Monte Carlo method for branched molecules. Transferable potentials for phase equilibria. 2. United-atom description of branched alkanes. J. Phys. Chem. B 1999, 103 (21), 4508−4517. (39) Plimpton, S. Fast parallel algorithms for short-range moleculardynamics. J. Comput. Phys. 1995, 117 (1), 1−19. (40) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction models for water in relation to protein hydration. In Intermolecular Forces; Pullman, B., Ed.; D. Reidel: Amsterdam, 1981; p 331. (41) Jorgensen, W. L.; Maxwell, D. S.; TiradoRives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 1996, 118 (45), 11225−11236. (42) Teleman, O.; Jonsson, B.; Engstrom, S. A molecular-dynamics simulation of a water model with intramolecular degrees of freedom. Mol. Phys. 1987, 60 (1), 193−203. (43) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (44) Widom, B. Some topics in the theory of fluids. J. Chem. Phys. 1963, 39 (11), 2802−2812. (45) Frenkel, D.; Smit, B. Understanding Molecular Simulation; Academic Press: New York, 1996. (46) Duan, Z. H.; Moller, N.; Greenberg, J.; Weare, J. H. The prediction of methane solubility in natural-waters to high ionicstrength from 0 °C to 250 °C and from 0 to 1600 bar. Geochim. Cosmochim. Acta 1992, 56 (4), 1451−1460. (47) Ott, J.; Boerio-Goates, J. Chemical Thermodynamics: Principles and Applications; Academic Press: New York, 2000. (48) Peng, D. Y.; Robinson, D. B. A new two-constant equation of state. Ind. Eng. Chem. Fundam. 1976, 15 (1), 59−64. (49) Tambach, T. J.; Bolhuis, P. G.; Hensen, E. J. M.; Smit, B. Hysteresis in clay swelling induced by hydrogen bonding: Accurate prediction of swelling states. Langmuir 2006, 22 (3), 1223−1234. (50) Fu, M. H.; Zhang, Z. Z.; Low, P. F. Changes in the Properties of a montmorillonite-water system during the adsorption and desorption of water-hysteresis. Clays Clay Miner. 1990, 38 (5), 485−492. (51) Duan, Z. H.; Mao, S. D. A thermodynamic model for calculating methane solubility, density and gas phase composition of methanebearing aqueous fluids from 273 to 523 K and from 1 to 2000 bar. Geochim. Cosmochim. Acta 2006, 70 (13), 3369−3386. (52) Harvey, A. H. Thermodynamic Properties of Water: Tabulation from the IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use, 1998. (53) Meng, E. C.; Kollman, P. A. Molecular dynamics studies of the properties of water around simple organic solutes. J. Phys. Chem. 1996, 100 (27), 11460−11470. (54) Koh, C. A.; Wisbey, R. P.; Wu, X. P.; Westacott, R. E.; Soper, A. K. Water ordering around methane during hydrate formation. J. Chem. Phys. 2000, 113 (15), 6390−6397.

(11) Rotenberg, B.; Marry, V.; Vuilleumier, R.; Malikova, N.; Simon, C.; Turq, P. Water and ions in clays: Unraveling the interlayer/ micropore exchange using molecular dynamics. Geochim. Cosmochim. Acta 2007, 71 (21), 5089−5101. (12) Boek, E. S.; Coveney, P. V.; Skipper, N. T. Monte Carlo molecular modeling studies of hydrated Li-, Na-, and K-smectites: Understanding the role of potassium as a clay swelling inhibitor. J. Am. Chem. Soc. 1995, 117 (50), 12608−12617. (13) Rahman, M. K.; Hossain, M. M.; Rahman, S. S. A shear-dilationbased model for evaluation of hydraulically stimulated naturally fractured reservoirs. Int. J. Numer. Anal. Methods Geomech. 2002, 26 (5), 469−497. (14) Chavez, M. D.; de Pablo, L.; de Pablo, J. J. Monte Carlo molecular simulation of the hydration of K-montmorillonite at 353 K and 625 bar. Langmuir 2004, 20 (24), 10764−10770. (15) de Pablo, L.; Chavez, M. L.; Sum, A. K.; de Pablo, J. J. Monte Carlo molecular simulation of the hydration of Na-montmorillonite at reservoir conditions. J. Chem. Phys. 2004, 120 (2), 939−946. (16) Titiloye, J. O.; Skipper, N. T. Molecular dynamics simulation of methane in sodium montmorillonite clay hydrates at elevated pressures and temperatures. Mol. Phys. 2001, 99 (10), 899−906. (17) Titiloye, J. O.; Skipper, N. T. Monte Carlo and molecular dynamics simulations of methane in potassium montmorillonite clay hydrates at elevated pressures and temperatures. J. Colloid Interface Sci. 2005, 282 (2), 422−427. (18) Chang, F. R. C.; Skipper, N. T.; Sposito, G. Monte Carlo and molecular dynamics simulations of interfacial structure in lithiummontmorillonite hydrates. Langmuir 1997, 13 (7), 2074−2082. (19) Chang, F. R. C.; Skipper, N. T.; Sposito, G. Monte Carlo and molecular dynamics simulations of electrical double-layer structure in potassium-montmorillonite hydrates. Langmuir 1998, 14 (5), 1201− 1207. (20) Cygan, R. T.; Liang, J. J.; Kalinichev, A. G. Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field. J. Phys. Chem. B 2004, 108 (4), 1255−1266. (21) Karaborni, S.; Smit, B.; Heidug, W.; Urai, J.; vanOort, E. The swelling of clays: Molecular simulations of the hydration of montmorillonite. Science 1996, 271 (5252), 1102−1104. (22) Chang, F. R. C.; Skipper, N. T.; Sposito, G. Computersimulation of interlayer molecular-structure in sodium montmorillonite hydrates. Langmuir 1995, 11 (7), 2734−2741. (23) Skipper, N. T.; Chang, F. R. C.; Sposito, G. Monte-Carlo Simulation of interlayer molecular-structure in swelling clay-minerals. 1. Methodology. Clays Clay Miner. 1995, 43 (3), 285−293. (24) Delville, A. Modeling the clay-water interface. Langmuir 1991, 7 (3), 547−555. (25) Delville, A. Structure of liquids at a solid interface - an application to the swelling of clay by water. Langmuir 1992, 8 (7), 1796−1805. (26) Delville, A. Structure and properties of confined liquids - a molecular-model of the clay water interface. J. Phys. Chem. 1993, 97 (38), 9703−9712. (27) Cheng, L.; Fenter, P.; Nagy, K. L.; Schlegel, M. L.; Sturchio, N. C. Molecular-scale density oscillations in water adjacent to a mica surface. Phys. Rev. Lett. 2001, 87 (15), 156103. (28) Park, C.; Fenter, P. A.; Nagy, K. L.; Sturchio, N. C. Hydration and distribution of ions at the mica-water interface. Phys. Rev. Lett. 2006, 97, (1), 016101. (29) Schlegel, M. L.; Nagy, K. L.; Fenter, P.; Cheng, L.; Sturchio, N. C.; Jacobsen, S. D. Cation sorption on the muscovite (001) surface in chloride solutions using high-resolution X-ray reflectivity. Geochim. Cosmochim. Acta 2006, 70 (14), 3549−3565. (30) Xu, L.; Salmeron, M. An XPS and scanning polarization force microscopy study of the exchange and mobility of surface ions on mica. Langmuir 1998, 14 (20), 5841−5844. (31) Bowers, G. M.; Bish, D. L.; Kirkpatrick, R. J. Cation exchange at the mineral-water interface: H3O+/K+ competition at the surface of nano-muscovite. Langmuir 2008, 24 (18), 10240−10244. H

dx.doi.org/10.1021/jp403349p | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

(55) de Pablo, L.; Chavez, M. L.; de Pablo, J. J. Stability of Na-, K-, and Ca-montmorillonite at high temperatures and pressures: A Monte Carlo simulation. Langmuir 2005, 21 (23), 10874−10884. (56) Leng, Y. S. Hydration force between mica surfaces in aqueous KCl electrolyte solution. Langmuir 2012, 28 (12), 5339−5349.

I

dx.doi.org/10.1021/jp403349p | J. Phys. Chem. C XXXX, XXX, XXX−XXX