Monte Carlo and molecular dynamics simulations of methane in potassium montmorillonite clay hydrates at elevated pressures and temperatures

Journal of Colloid and Interface Science 282 (2005) 422–427 www.elsevier.com/locate/jcis

Monte Carlo and molecular dynamics simulations of methane in potassium montmorillonite clay hydrates at elevated pressures and temperatures James O. Titiloye a,∗ , Neal T. Skipper b a Chemical Engineering and Applied Chemistry, SEAS, Aston University, Birmingham B4 7ET, UK b Department of Physics and Astronomy, University College, Gower Street, London WC1E 6BT, UK

Received 3 June 2004; accepted 17 August 2004 Available online 22 September 2004

Abstract The structure and dynamics of methane in hydrated potassium montmorillonite clay have been studied under conditions encountered in sedimentary basin and compared to those of hydrated sodium montmorillonite clay using computer simulation techniques. The simulated systems contain two molecular layers of water and followed gradients of 150 bar km−1 and 30 K km−1 up to a maximum burial depth of 6 km. Methane particle is coordinated to about 19 oxygen atoms, with 6 of these coming from the clay surface oxygen. Potassium ions tend to move away from the center towards the clay surface, in contrast to the behavior observed with the hydrated sodium form. The clay surface affinity for methane was found to be higher in the hydrated K-form. Methane diffusion in the two-layer hydrated K-montmorillonite increases from 0.39 × 10−9 m2 s−1 at 280 K to 3.27 × 10−9 m2 s−1 at 460 K compared to 0.36 × 10−9 m2 s−1 at 280 K to 4.26 × 10−9 m2 s−1 at 460 K in Na-montmorillonite hydrate. The distributions of the potassium ions were found to vary in the hydrates when compared to those of sodium form. Water molecules were also found to be very mobile in the potassium clay hydrates compared to sodium clay hydrates.  2004 Elsevier Inc. All rights reserved. Keywords: Methane; K; Na; Montmorillonite; Surface complex; Clay; Simulations

1. Introduction The structure and mobility of aqueous and organic fluids in swelling clay minerals are crucial to the location and migration of petroleum hydrocarbons in the soil. Clays are layered aluminosilicates with structures and interlayer pores that are ideal for studying fundamental properties of confined fluids. They are found in terrestrial environments and marine sediments, and form a basic component of soil. 2:1 clays are made up of sheets of octahedrally coordinated cations sandwiched between two sheets of tetrahedrally coordinated cations. The internal pores are saturated by hydroxyl groups and several molecular layers of absorbed water. Negative layer charges arise from isomorphic cation * Corresponding author. Fax: +44(0)121-359-4094.

E-mail address: [email protected] (J.O. Titiloye). 0021-9797/$ – see front matter  2004 Elsevier Inc. All rights reserved. doi:10.1016/j.jcis.2004.08.131

substitution of the tetrahedral and octahedral cations. These layer charges are balanced by charge compensating counterions such as Na+ and K+ and tend to form ionic solutions with water molecules found in hydrated clays thereby forcing the clay layers apart and resulting in swelling. A number of studies have been performed to highlight and understand how swelling clay minerals play a crucial role in oil and gas production and exploration [1,2]. In addition, there have been several simulations and experimental studies directed at understanding the structure and dynamics of interlayer species in clays [3–6]. Sposito et al. [4,5] reported the simulation of radial distribution function for interlayer water in lithium, sodium and potassium montmorillonite. Sutton and Sposito [7] also used molecular simulation to investigate the interlayer water structure of Cs-smectite hydrates. Understanding the structure and behavior of these interlayer species will help alleviate the problem of borehole instability

J.O. Titiloye, N.T. Skipper / Journal of Colloid and Interface Science 282 (2005) 422–427

commonly associated with the uptake of water by smectite from drilling fluids. Previous work has shown that sorption of water by different smectite clays is dependent on the size and charge of the hydrating counterion which is largely responsible for the swelling behavior of clays [8–10]. However, most of these previous works have focused largely on clay–water-cation systems, with few comparable studies existing on the interaction of hydrocarbons and organics fluids with interlayer species. In a recent study, Park and Sposito [11] used simulation techniques to investigate the effect of montmorillonite surfaces on methane hydrate formation and proposed the occurrence of methane hydrate in natural sediment in presence of clay surfaces at ambient condition. One of our aims in this work is to use simulation techniques to investigate the underlying mechanism of smectite clay swelling in the presence of organic fluids. We determine the structure and mobility of methane in the interlayer region of potassium montmorillonite hydrated clay and monitor the impact of hydrophobic molecules on the microscopic properties of other interlayer species. We compared our results with previous studies on sodium montmorillonite hydrated clay [12,13] and investigate the difference between the role of K+ and Na+ ions as clay swelling inhibitors in the presence of organic molecules.

2. Methods A combination of constant stress Monte Carlo (MC) and constant volume molecular dynamics (MD) methods using the codes MONTE and DLPOLY was used in the simulation [14,15]. Skipper and Sposito et al. [10,14] have developed methodology for investigating clay minerals while the general principles of these simulations techniques have been described by Frenkel and Smit [16]. The model clay mineral is a “Wyoming”-type smectite clay, of composition [Al28 Mg4 ][Si62 Al2 ]O160(OH)6− 32 , with 6 charge compensating potassium cations, 4 methane particles, and 64 water molecules corresponding to a two-layer hydrate system. The simulation cell in the Monte Carlo runs is composed of opposing 21.12 × 18.28 Å2 clay mineral particles. This simulation cell was doubled in the direction normal to the clay sheets for the molecular dynamics runs, thus allowing us a cutoff distance of 9 Å in the minimum image convention used in the MD program DLPOLY. Temperature and pressure conditions corresponding to specific burial depth up to 6 km were introduced in the simulation to reproduce the pressures and temperature encountered in a typical petroleum-rich sedimentary basin. A compacting stress gradient of 150 bar km−1 and a geothermal gradient of 30 K km−1 were used in the simulations [2]. In the simulation cell, water and methane molecules were initially arranged randomly, while the potassium counterions were assigned positions opposite the clay negative charge

423

sites. The interlayer species were then selected randomly to move for over 5 million MC steps. This gave information on total potential energy of the system, interlayer species positions including orientations, and clay layer spacing. The equilibrium configurations obtained from MC calculations were then used as starting configurations in the MD simulations. The MD simulations were performed with a timestep of 0.5 fs for up to 800 ps using an NVT ensemble. Interlayer molecular configurations were stored every 1000 steps and used for results analysis. The simulation used Ewald summation technique with 3-dimensional periodic boundary conditions. The potential model used for the water molecules is based on the TIP4P model of Jorgensen at al. [17] where water was treated as a rigid molecule with four intermolecular interaction sites. Methane was represented by a single Lennard– Jones in the OPLS model [18], and the interactions involving the clay mineral with interlayer species follow the methods described by Boek et al. [19]. These interactions are pairwise additive involving Coulomb and Lennard–Jones 6–12 potential of the form 
 qi qj Aij Bij − 6 + 12 , U (r) = (1) rij rij rij where U is the potential energy of the system, the indices i and j indicate charge sites on each of the two interacting species, q is the effective charge on a site, rij is the intermolecular site separation, and A and B represent the van der Waals attractive and short-range repulsive terms, respectively. The potential parameters employed in the simulation follow those used in previous studies for Namontmorillonite [13] while the additional interaction parameters used are K–O/Ow : A = 1.17 × 103 , B = 6.33 × 105 ; K–K: A = 2.25 × 103 , B = 6.68 × 105; K–CH4 : A = 2.67 × 103 , B = 2.39 × 106 , with the units of kcal Å6 mol−1 and kcal Å12 mol−1 for A and B, respectively. The self-diffusion coefficients for the interlayer species were calculated using 3D Einstein relation:  2  D = 1/6d/dt r(t − t0 ) , (2) where |r(t − t0 )|2  is the mean-square displacement, t is the simulation elapsed time and t0 is arbitrary starting point. The slope of plots of the mean-square displacement against simulation time 20 < (t − t0 ) < 200 after equilibration was used in obtaining the diffusion coefficients.

3. Results and discussion Two-layer hydrate K-montmorillonite simulations were performed with 64 water molecules, 4 methane, and 6 potassium ions as interlayer species. Five different sets of simulation were carried out with each set simulated at different temperature and pressure that correspond to burial depths of 0, 1, 2, 3 and 6 km under gradients of 150 bar km−1

424

J.O. Titiloye, N.T. Skipper / Journal of Colloid and Interface Science 282 (2005) 422–427

Table 1 Calculated clay layer spacing for 2-layer hydrated K-montmorillonite with methane particles as a function of depthsa compared to hydrated Na-form Depth (km)

2-layer K-clay layer spacing (Å)

2-layer Na-clay layer spacing (Å) [12,13]

0.0 1.0 2.0 3.0 6.0

15.37 (0.06) 15.60 15.45 15.67 15.99

15.36 (0.06) 15.45 15.40 15.40 15.82

a The uncertainties in the calculated layer spacing are given in brackets.

Fig. 1. Density profiles for CH4 molecules in (a) K- and (b) Na-montmorillonite clays as a function of distance from clay surface at different depths. The midplane is taken as the origin in the z direction. The Na curves are displaced by 0.015 units.

and 30 K km−1 . The Monte Carlo simulation on each set resulted in an equilibrium layer spacing as shown in Table 1. The calculated layer spacing was compared to that of hydrated Na-form and found to be consistently higher in the hydrated K-montmorillonite, and also found to increase with burial depths. These are in line with both experimental (15.6 ± 0.5 Å) and simulated (15.5 ± 0.2 Å) observations of 2-layer hydrates of K-montmorillonite [5] and the current calculated values reflecting the presence of organic particles in the clays. The distribution of methane particles in the interlayer region is shown in Fig. 1 compared with methane distribution in the interlayer of sodium montmorillonite hydrate under similar conditions. This shows the number density profile of methane particles as a function of distance z measured from the center of the interlayer region. It was observed that in these two layer hydrate systems, the hydrocarbon molecule does not discriminate between the two hydrated forms. The radial distribution functions (RDFs) for CH4 –O and CH4 – Ow for hydrated K-form shown in Fig. 2 indicate strong correlation peaks at around 3.7 Å (for CH4 –Ow ) and 3.6 Å (for CH4 –O). The relative peak height in the figure for the first hydration shell is a function of coordination number calculated for methane particles in the system. The coordination number was calculated by integration up to the first minimum following the main correlation peaks. Methane parti-

Fig. 2. Radial distribution functions for (a) CH4 –O and (b) CH4 –Ow at different depth for hydrated K-montmorillonite.

Table 2 Coordination number for methane particle in hydrated K-montmorillonite as a function of depthsa compared with hydrated Na-formb Depth (km)

2-layer K-montmorillonite

2-layer Na-montmorillonite

nCH4 –Ow

nCH4 –Os

nt

nCH4 –Ow

nCH4 –Os

nt

0.0 1.0 2.0 3.0 6.0

12.28 13.64 12.98 13.21 12.89

6.34 5.95 6.28 6.20 5.56

18.62 19.59 19.26 19.41 18.45

12.71 12.17 12.09 12.47 12.33

6.23 5.84 5.78 5.88 5.67

18.94 18.01 17.87 18.35 18.00

a O and O denote water and clay surface oxygen, respectively. n is the w s t total number of oxygen around methane. b Refs. [12,13].

cles are seen to be surrounded by about 13 water molecules and 6 clay surface oxygen (Table 2), with a clear preference to move closer to the surface and away from the center. The clay surface affinity for methane particles appears to be higher in the K-form than in the Na-form, and this trend was noticeable as burial depth increases up to 3 km. However, the dynamics and mobility of methane were found to be similar in both forms especially at lower depths. The self diffusion coefficients were calculated from MSD versus time curves. At 3 km depth methane was noted to move faster in K-montmorillonite with Dmethane = 2.09 × 10−9 m2 s−1 compared to 1.83 × 10−9 m2 s−1 found with the Na-form (Table 3). However, the rate of movement at 6 km in the

J.O. Titiloye, N.T. Skipper / Journal of Colloid and Interface Science 282 (2005) 422–427

K-form was slower (Dmethane = 3.27 × 10−9 m2 s−1 ) than observed with Na-form (Dmethane = 4.26 × 10−9 m2 s−1 ). At higher temperatures of 370 and 460 K corresponding to burial depths of 3 and 6 km, respectively, the evidence of increased hydrophobic interaction observed with hydrated Na-montmorillonite [13] persists in the K-form (Fig. 3). As there is a reduction in clay porosity with depth and as we concluded in earlier studies that maximum diffusion of methane through hydrated clays occurs at around 3 km depth, our current results on K-hydrated forms supports this conclusion. Furthermore, the increase in the diffusion coefficients noticed for methane particles as depth increases reinforces the suggestion that two-layer clay hydrates are unlikely to trap methane effectively during migration.

425

3.1. Water and potassium counterions distribution and mobility Two populations of potassium ions were identified in the simulation, similar to those found for hydrated sodium form [12,13]. However, it is interesting to see how these two populated sites have switched roles in the two hydrated forms (Fig. 4). The first population sites are those adjacent to the tetrahedral negative charge sites in close proximity to the clay surface. This site has attracted more potassium ions than sodium ions as previously observed. This is due to the reluctance of potassium ions to fully associate with water molecules in the interlayer pores. A secondary peak located

Table 3 Self-diffusion coefficients (m2 s−1 ) for interlayer water, K+ and CH4 in hydrated K-montmorillonite as a function of burial depthsa compared to hydrated sodium formb T, P, “depth”

280 K, 1 bar_0 km 310 K, 150 bar_1 km 340 K, 300 bar_2 km 370 K, 450 bar_ 3.0 km 460 K, 900 bar_ 6.0 km

2-layer K-clay

2-layer Na-clay

DH2 O DK (×109 )

DCH4

DH2 O DNa (×109 )

DCH4

0.77 1.56 2.57 3.40 7.91

0.39 1.20 1.25 2.09 3.27

0.57 1.50 1.77 3.27 5.82

0.36 0.96 1.20 1.83 4.26

0.12 0.28 0.72 0.89 2.06

0.12 0.30 0.18 0.89 1.34

a The uncertainties for the calculated D values are D H2 O (±0.03), DK/Na (±0.02), DCH4 (±0.05). b Refs. [12,13].

Fig. 3. Radial distribution functions for methane–methane interactions at different depths for 2-layer hydrated K-montmorillonite.

Fig. 4. Interlayer density profiles for (a) K+ and (b) Na+ counterions as a function of distance from clay surface at different burial depths. The midplane is taken as the origin in the z direction.

426

J.O. Titiloye, N.T. Skipper / Journal of Colloid and Interface Science 282 (2005) 422–427

Fig. 5. Comparing the density profiles of interlayer water in hydrated forms of (a) K- and (b) Na-montmorillonite as a function of distance from clay surface at different depths. The midplane of the interlayer region is taken as the origin in the z direction.

in the K-form (Fig. 4a) at around 1.7 to 1.8 Å appears to be part of migrating K ions from the midplane towards the surface sites. The second population site located near the midplane of the interlayer region contains fewer potassium ions and this contrasts with observed behavior in the Naform where sodium ions are more in populations and are fully solvated with outer sphere complex of six water molecules. This contrasting behavior between the two ions in the different hydrated forms confirms the effect of weak solvation of potassium ions earlier observed by Radnai et al. [20] in clay minerals. The calculated self-diffusion coefficients are of the same order in both hydrated forms except at 2 km depth, where it was noticed that K ions were moving considerably faster (DK = 0.72 × 10−9 m2 s−1 ) than Na ions (DNa = 0.18 × 10−9 m2 s−1 ). Water molecules are also noticed to move faster in the K-form than in the hydrated sodium form which is a consequence of reduced tendency to solvate K ions. The water molecules in both K- and Na-forms were also observed to bond strongly with the clay surface oxygen. This is evident from Fig. 5 showing the density profile for water molecules in the two hydrated forms as a function of burial depths. The interlayer density profiles also confirm that water is layered in both systems. The relative position of layered water molecules to the K-counterions and methane particles is depicted in the simulation snapshot shown in Fig. 6.

4. Summary Computer simulation has been used to study the structure and dynamics of methane in hydrated K-montmorillonite

(a)

(b)

Fig. 6. Molecular dynamics snapshots at 800 ps showing layered water molecules position, K+ ions (colored purple and coordinated to oxygen atoms) and methane particles (colored yellow) at (a) 0 and (b) 6 km.

clay and compared to structure and dynamics of methane in hydrated Na-montmorillonite. The simulation revealed that K+ ions preferred sites adjacent to the surface of the clay sheets while Na+ ions have the tendency to stay at the center of clay sheets where they can be fully solvated. The total coordination number calculated for the K-ions from the nearest neigbor oxygens ranges from 9.7 to 10.7 (depending on the burial depth) compared to a range of 8.6 to 10.8 found for the hydrated Na-ions. Although the clay surface affinity for methane particles is higher in the hydrated

J.O. Titiloye, N.T. Skipper / Journal of Colloid and Interface Science 282 (2005) 422–427

K-montmorillonite than in the sodium form, methane molecules do not discriminate between a Na- and a K-hydrated form in their distribution, and their dynamics and mobility are similar in both forms especially at lower depth. The calculated layer spacing was found to be consistently higher in the hydrated K-montmorillonite than in the sodium form and in general consistent with experimental values for Khydrated clays. The observed increase in the diffusion coefficients of methane with depth strongly suggests rapid migration of hydrocarbons within the clay layer and this requires an optimum depth for effective molecular trapping of hydrocarbons in the soil. Further work is in progress investigating the effect of increasing methane concentration on migrations and on the behavior of counterions in hydrated clay systems.

[3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]

References

[17] [18]

[1] G. Sposito, N.T. Skipper, R. Sutton, S.H. Park, F.R. Chang, A.K. Soper, J.A. Greathouse, Proc. Natl. Acad. Sci. USA 96 (1999) 3358. [2] F.K. North, Petroleum Geology, second ed., Unwin Hyman, Boston, 1990.

[19] [20]

427

E.S. Boek, M. Sprik, J. Phys. Chem. B 107 (2003) 3251. S.H. Park, G. Sposito, J. Phys. Chem. B 104 (2000) 4642. G. Sposito, S.H. Park, R. Sutton, Clays Clay Miner. 47 (1999) 192. J. Greathouse, G. Sposito, J. Phys. Chem. B 102 (1998) 2406. R. Sutton, G. Sposito, J. Colloid Interface Sci. 237 (2001) 174. E.S. Boek, P.V. Coveney, N.T. Skipper, J. Am. Chem. Soc. 117 (1995) 12608. F.C. Chang, N.T. Skipper, G. Sposito, Langmuir 11 (1995) 2734. A.S. Bains, E.S. Boek, P.V. Coveney, S.J. Williams, M.V. Akbar, Mol. Sim. 26 (2001) 101. S.H. Park, G. Sposito, J. Phys. Chem. B 107 (2003) 2281. J.O. Titiloye, N.T. Skipper, Chem. Phys. Lett. 329 (2000) 23. J.O. Titiloye, N.T. Skipper, Mol. Phys. 99 (2001) 899. N.T. Skipper, G. Sposito, F.-H. Chang, Clays Clay Miner. 43 (1995) 294. T.R. Forrester, W. Smith, DLPOLY Package, CCLRC, Daresbury Laboratory, Daresbury, 1993. D. Frenkel, B. Smit, Understanding Molecular Simulation, Academic Press, London, 1996. W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. W.L. Jorgensen, J.D. Madura, J.D. Swenson, J. Am. Chem. Soc. 106 (1984) 6638. E.S. Boek, P.V. Coveney, N.T. Skipper, J. Am. Chem. Soc. 117 (1995) 12608. H. Ohtaki, T. Radnai, Chem. Rev. 93 (1993) 1157.