Understanding methane/carbon dioxide partitioning in clay nano- and meso-pores with constant reservoir composition molecular dynamics modeling

The interactions among fluid species such as H2O, CO2, and CH4 confined in nanoand meso-pores in shales and other rocks is of central concern to understanding the chemical behavior and transport properties of these species in the earth’s subsurface and is of special concern to geological C-sequestration and enhanced production of oil and natural gas. The behavior of CO2, and CH4 is less well understood than that of H2O. This paper presents the results of a computational modeling study of the partitioning of CO2 and CH4 between bulk fluid and nanoand meso-pores bounded by the common clay mineral montmorillonite. The calculations were done at 323 K and a total fluid pressure of 124 bars using a novel approach (constant reservoir composition molecular dynamics, CRC-MD) that uses bias forces to maintain a constant composition in the fluid external to the pore. This purely MD approach overcomes the difficulties in making stochastic particle insertion–deletion moves in dense fluids encountered in grand canonical Monte Carlo and related hybrid approaches. The results show that both the basal siloxane surfaces and protonated broken edge surfaces of montmorillonite both prefer CO2 relative to CH4 suggesting that methods of enhanced oil and gas production using CO2 will readily displace CH4 from such pores. This preference for CO2 is due to its preferred interaction with the surfaces and extends to approximately 20 Å from them.


Introduction
2][3][4][5][6] The behavior of water and cations at mineral surfaces and in the interlayer galleries of expandable clays (smectites) has been studied experimentally and computationally for many years, [7][8][9][10][11][12][13][14][15][16][17][18][19] and more recently there has been increased interest in carbon dioxide and hydrocarbons.  The tual interactions between CO 2 and CH 4 in nano-and meso-scale confinement, however, remain poorly understood.Most simulation studies have focused on single component adsorption in the interlayers of clay minerals.8][49][50][51][52] In shales and other sedimentary rocks, pores have a wide range of sizes, and the effects of pore size on CO 2 / CH 4 partitioning are poorly understood.This paper presents the results of a computational molecular dynamics (MD) modeling study of the partitioning of CO 2 and CH 4 between bulk fluid and slit-like nano-and meso pores from 4 to 73 Å thick and the structure of the fluid in those pores.In contrast to most computational studies of clays, the montmorillonite substrate model here is of finite size (Fig. 1), allowing study of the fluid interaction with the siloxane basal surfaces and the protonated sites on broken edges of the clay layers.The results show that the partitioning of CO 2 and CH 4 between the bulk fluid (representing fluid in large pores in shales and other rocks) and the nano-and meso-pores is substantially influenced by the pore thickness, with CO 2 being largely preferred in narrow confined pores.CO 2 is also greatly preferred at the protonated edges of broken clay layers.
The MD simulations were done using a novel approach that uses self-adjusting bias forces to maintain a constant composition in the fluid external to the pore (constant reservoir composition molecular dynamics, CRC-MD).Such self-adjusting bias forces have been used previously to maintain a constant chemical potential of the fluid species in modelling the growth of urea crystals from solution (constant chemical potential MD) 53,54 and to create a concentration gradient across a membrane for modelling the gas transport and separation of mixtures in membranes (concentration gradient driven MD). 55The advantages of using a purely MD approach to maintain the composition of dense fluids rather than using the stochastic particle insertiondeletion moves used in grand canonical Monte Carlo and related hybrid approaches, such as dual control volume grand canonical molecular dynamics, 56 have been discussed in detail in ref. 55.Inspired by these studies, [53][54][55] we have implemented the CRC-MD method, which has proven to be an efficient way to evaluate partitioning of fluid species between bulk fluid and nano-and meso-pores and also onto the surfaces of finite size solid substrates.

Methods
In the CRC-MD method, each simulation cell consists of a bulk fluid reservoir, two bias force regions, two composition control regions, two transition regions, and the substrate (clay) + pore assemblage in the center (Fig. 1a).The concentrations of the fluid species in the control regions are maintained at constant values by forces in the bias force regions.These bias forces act in such a way that if the concentration of a given species in a control region is different than the target concentration molecules of that species are moved into or out of the control region from or to the reservoir.A detailed explanation of the functional forms of the bias forces and how they work can be found elsewhere. 55he substrate used in the simulations was the expandable smectite clay mineral, montmorillonite, which develops a permanent negative structural charge by the isomorphic substitution of Al 3+ for Si 4+ in the tetrahedral sheet and Mg 2+ for Al 3+ in the octahedral sheet.The model here has a structural formula of M + 0.75 (Si 7.75 Al 0.25 )(Al 3.5 Mg 0.5 )O 20 (OH) 4 . 12The distribution of the isomorphic substitutions has a quasi-disordered pattern and was prepared using the ''Supercell'' program, 57 which creates substitutions in accordance with an extension of Lowenstein's rule, 58 which forbids tetrahedral Al-O-Al and octahedral Mg-O-Mg linkages.The simulated montmorillonite particles have lateral dimensions of B73.0 Å Â 41.4 Å along the a and b crystallographic axes, respectively, with an orthorhombic unit cell.
To create a finite size particle, the montmorillonite structure was cleaved along (010), creating broken edge sites on that surface.On the broken edges, the dangling tetrahedral Al 3+ and Si 4+ sites are saturated by single OH À groups, and most of the dangling octahedral Al 3+ sites are saturated with 1 OH À and 1 H 2 O molecule. 59Meanwhile, according to ab initio calculations, 60 when the octahedral Mg 2+ sites are saturated with 1 H 2 O molecule and 1 OH À group, proton transfer reactions occur with the neighboring octahedral Al 3+ sites and tetrahedral Si 4+ sites to generate octahedral Mg 2+ sites saturated with 2 H 2 O molecules (as predicted by the crystal growth theory), 61 and a neighboring Al 3+ site with 2 OH À groups.These coordination environments were used in the clay substrate.Furthermore, the oxygen atoms of the OH À groups of the broken edge sites were assigned a slightly higher negative charge (À0.9659 |e|) compared to the original CLAYFF 62 value (À0.95 |e|).This was done not only to ensure that the broken edges of the montmorillonite particle are electrostatically neutral, but also because the OH À groups of the edge sites have fewer bond connections compared to those in the interior of the clay structure.Importantly, the two (010) surfaces of this model have different compositions.One contains only unsubstituted tetrahedral Si 4+ and octahedral Al 3+ sites (left side of the particle in Fig. 1b), whereas the other contains two substituted tetrahedral Al 3+ and two octahedral Mg 2+ sites (right side of Fig. 1b).
The modeled montmorillonite particle consists of three T-O-T layers which encompass two anhydrous interlayers and expose two external basal surfaces to the pore (Fig. 1a).Because of the 3-dimensional periodic boundary conditions used, only the broken (010) surfaces and the basal surfaces bounding the pore are exposed to the fluid phase.The lateral dimensions of the simulation cell are 300 Å Â 41.4 Å (Fig. 1a).The thickness of the slit-like pore between the external basal surfaces (Fig. 1b) were varied from 4 to 73 Å.For pore thicknesses less than 10 Å, at the start of each simulation the charge compensating Na + ions were placed at the midplane of the pore far from the substitution sites on the external surfaces, thereby allowing the Na + ions to choose their preferred adsorption sites during the course of the simulation.For models with pore thicknesses 410 Å, the Na + ions were initially placed 5 Å from the external basal surface.The CO 2 and CH 4 molecules were initially placed external to the clay particle at distances 4B30 Å from the two broken (010) surfaces.The number of CO 2 and CH 4 molecules in the full simulation cell varies from 3456 to 8063 and increases with increasing pore thickness in order to maintain the desired target densities in the control region.We analyze the CO 2 /CH 4 partitioning and density profiles near the clay particle in three different dimensions: (i) normal to the basal surfaces bounding the pore (z direction); (ii) parallel to the broken edges in the transition region (z direction), and (iii) normal to the broken edges in the transition region (x direction).We define the pore thickness as the distance between the centers of the oxygen atoms of the basal surfaces bounding the pore.Similarly, we define the surfaces of the broken edges as the plane containing the centers of the protonated oxygen atoms on the surface.][33] We performed the molecular dynamics simulations in the canonical NVT ensemble at 323 K using the LAMMPS simulation package. 63A Nose ´-Hoover thermostat was used to control the temperature. 64The interatomic interactions for the montmorillonite were obtained from the CLAYFF 62 force field, and the parameters for the broken edge sites were obtained from newly developed metal-O-H bending potentials consistent with ClayFF. 65,66The CO 2 and CH 4 molecules were represented by the EPM2 67 and TraPPE-UA 68 models, respectively.Three-dimensional periodic boundary conditions were employed with a cutoff of 14.0 Å for short-range non-electrostatic interactions, and the long-range electrostatic interactions were computed using the particleparticle particle-mesh (PPPM) summation algorithm 69 with an accuracy of 10 À6 .A time step of 1 fs was employed to integrate the equations of motion.Each system was equilibrated for 10 ns followed by another 2 ns of data production with the atomic coordinates recorded every 10 fs.Importantly, we fixed the positions of 18 octahedral atoms (Al 3+ or Mg 2+ ) in the simulated montmorillonite structure (2 left: 2 middle: 2 right per TOT layer) in order to prevent the clay layers from moving with respect to each other.
A modified version of the PLUMED 2.3.0 plugin 70 was used to apply the bias forces to maintain the fluid composition in the control regions.The target CO 2 and CH 4 densities in the control regions were set to 2 CO 2 and 2 CH 4 molecules per nm 3 .At 323 K, these densities yield a total fluid pressure of 124 bars, as calculated by using the Peng-Robinson equation of state. 71These temperature and pressure conditions are relevant to the upper part of the Earth's crust and petroleum reservoirs.The width of the transition regions, control regions and the bias force regions were set to 25 Å, 25 Å and 2.5 Å, respectively.The target composition of the fluid species, i, in the control regions was maintained using the force constants k L i = k R i = 5000 kJ nm 3 mol À1 which are placed at the center of the bias force regions located on the left, L, and right, R, side of the clay particle.The compositions in the control regions were monitored at intervals of 0.5 ps during the entire simulation.

Results and discussion
The mole fractions and number densities of CO 2 and CH 4 molecules in the pores vary greatly depending on the pore thickness (Fig. 2a and b).At the smallest thickness studied (4.0 Å), neither CO 2 nor CH 4 enters the pore.At a thickness of 6.0 Å, CO 2 fills the pore and there is only a negligible fraction of CH 4 molecules near the entry points.These results are in good agreement with recent experimental and simulation studies that show that the basal spacing required for CO 2 intercalation in the interlayer galleries of smectite clays is less than for CH 4  This journal is © the Owner Societies 2019 due to steric effects. 6,31,32,38,39At a pore thicknesses of 7.5 Å some CH 4 enters the pore, and the CO 2 /CH 4 ratio decreases rapidly with increasing pore thickness up to 33.0 Å (Fig. 2a).At larger pore thickness, the CO 2 /CH 4 ratio decreases less rapidly and slowly approaches the composition in the control regions (2 molecules per nm 3 ), although for the thicknesses studied here it never reaches the control region composition.The number density of CO 2 molecules in the pores is always larger than in the control region, with the highest value of B8.2 molecules per nm 3 at 9.0 Å and progressively smaller values at larger thicknesses.In contrast, the CH 4 density increases with increasing thickness to 23.0 Å and then decreases with increasing thickness.However, it is never much greater than the value in the control region of 2 molecules per nm 3 .(Fig. S1, ESI, † shows that the fluid compositions in the control regions were kept at the target values with great accuracy throughout the production runs.)These changes highlight the important conclusion that pores with dimensions of 1 to a few nm bounded by basal clay surfaces preferentially incorporate CO 2 relative to CH 4 , with the smallest pores showing the greatest preference (Fig. 2b).These results are in agreement with recent GCMC simulation studies of CO 2 /CH 4 adsorption in the interlayers of Na-montmorillonite at 318 K and P fluid = 0 to 300 bars. 48Importantly, the CO 2 /CH 4 ratios at small pore sizes (o23.0Å) are in excellent agreement with the selectivity parameters obtained from recent GCMC simulation studies at similar thermodynamic conditions (323 K and P fluid = 100 bar). 49he probability density profiles (PDPs) of CO 2 and CH 4 in the pores perpendicular to their surfaces show that the changes in their concentrations with pore thickness are due to preferential sorption of CO 2 and structuring of the fluid near the pore surfaces (Fig. 3).Significant structuring of the fluid into three discernable layers extends to B15.0 Å from each pore surface, with CO 2 /CH 4 ratios larger than in the composition control regions extending to B20 Å from the surfaces.At a pore thickness of 6.0 Å, the CO 2 molecules are located at the mid-plane of the pore (3.0 Å from each surface), as expected since this thickness corresponds to a clay interlayer with 1 layer of intercalated molecules. 32At a pore thickness of 7.5 Å, the CO 2 begins to form two layers, and the small amount of CH 4 is at the middle of the pore.At 9.0 Å there are two well developed layers of CO 2 and CH 4 .At 12.3 Å a third layer of CO 2 and CH 4 begins to develop, and the layer of CO 2 nearest the surfaces develops a shoulder.At 15.6 Å the central layer begins to split into two, and at 20.5 Å a fifth layer begins to develop.By 43.0 Å the fluid structuring near the surfaces is fully developed, and the CO 2 and CH 4 concentrations in the middle of the pore are essentially equal.At 73.0 Å the concentrations in the central region of the pore are close to those in the control volume, 2 molecules per nm 3 .For all pore thicknesses greater than 33.0 Å the fluid at each surface is structured into three layers of CO 2 and CH 4 with maxima near 3.0, 6.8, and 10.4 Å and peak densities that decrease with increasing distance from the surface.The PDP's are in good agreement with those from GCMC simulation studies of single component adsorption in slit nanopores in Na-montmorillonite at 323 K and P fluid = 100 bars with a thickness of B21.0 Å 49 and also in interlayers of Na-montmorillonite at 298 K and P fluid = 40 bars. 45,46Such three-layer structures commonly occur near surfaces for many fluids, including H 2 O. 72 The differences in the structure of the CO 2 and CH 4 layers near the pore surfaces are due to the differences in their interaction with the basal surface of the clay.For the CO 2 layer nearest to the surface, the peak at 3.0 Å is due to CO 2 molecules adsorbed with one O CO 2 located above the center of a ditrigonal cavity on the basal surface and the other O CO 2 located above a Si tetrahedron, as observed in earlier simulation studies of smectite interlayers. 25,31,32The CO 2 molecules at 3.9 Å are located with their C CO 2 closer to a tetrahedral site.25][26][27][28][31][32][33][34]37,38,40 In contrast, the lack of structuring in the CH 4 peak nearest the surface and the equal distances between the three layers suggest that these molecules are interacting much more weakly with the surfaces and behave like hard spheres with the layers packed on one top of each other.This conclusion is consistent with the incorporation of CH 4 into the interlayer galleries of expandable clays occurring by a passive, space filling mechanism with basal spacings similar to 7.5 Å. 39 Increasing CO 2 and CH 4 incorporation with increasing pore thickness also greatly effects the coordination of the Na + ions to the pore surfaces (Fig. 3).At 4.0 Å, most of the Na + ions are located at the mid-plane of the pore (near 2.0 Å from each surface) with broad shoulders at 0.5 Å from each surface.Based on previous simulation studies by Greathouse et al.,11 the Na + ions at 0.5 Å are located above the centers of ditrigonal cavities, and those at 2.0 Å are adsorbed near the Si/Al tetrahedra.At 6.0 Å, where CO 2 begins to fill the pore, more Na + occur above the ditrigonal cavities, and as the pore thickness increases progressively more occur on these sites.At larger pore thicknesses (47.5 Å), the Na + ions located at 2.0 Å from the surfaces are adsorbed only near tetrahedral Al 3+ sites.
As near the basal surfaces bounding the pores, the concentrations of CO 2 and CH 4 on the surfaces of the broken edges of the clay layers are greater than in the control volumes, and CO 2 is preferentially associated with the surface (Fig. 4).The layered structure of the fluid near the broken edge surfaces is also generally similar to that on the basal surfaces, and concentrations greater than in the control regions extending to B20 Å from the surface.The concentrations of CO 2 and CH 4 in the well-defined layers near the broken surfaces are, however, significantly lower than above the basal surfaces (compare values in Fig. 4b with those in Fig. 3), suggesting that both are attracted less strongly to the protonated edge sites than to basal siloxane surfaces.The PDP's show three peaks for CO 2 (1.7 Å, 2.8 Å and 6.3 Å) and two peaks for CH 4 (3.4Å and 6.7 Å).The positions of the two peaks for CH 4 are the same as above the basal surfaces, again illustrating the weak, non-specific interaction of CH 4 with silicate surfaces.In contrast, the positions of the two peaks for CO 2 are significantly closer to the broken edge surfaces.Recent GCMC simulation studies showed a greater preference for CO 2 relative to CH 4 in hydroxylated pores in quartz, which the authors attributed to stronger interaction of CO 2 with the protonates sites. 49The origin of this difference will require substantial additional work, but we note that in the results here the differences between the two broken (010) surfaces with differently substituted octahedral and tetrahedral sites are not significant.Note that in these simulations the structuring of the fluid above the broken surfaces does not depend on pore size, since the broken surfaces face the transition region, for which the thickness is independent of pore size.The PDP's of the fluid species parallel to the broken edges (sum of all molecules within 15.0 Å of the surface) further illustrate the preferential association of CO 2 with the broken edges (Fig. 4a).These plots also demonstrate that the effects of the surface extend only B10 Å into the fluid in the transition region above the opening of the pore.

Conclusions
Computational molecular modeling of the partitioning of supercritical CO 2 and CH 4 between bulk fluid and nano-and meso-pores 4 to 73 Å thick bounded by the basal surfaces of the expandable clay mineral montmorillonite show that both species are preferentially incorporated into the pores with the preference for CO 2 much greater than for CH 4 .This behavior is due to the association of the fluid species with the clay surface and with a much greater preference of the surface for CO 2 relative to CH 4 .Structuring of the fluid extends to B15 Å from the surface, and the CO 2 /CH 4 ratio is greater than in the control region to B20 Å from the surface.The total CO 2 /CH 4 ratio in the pores is greatest at a pore thickness of 6 Å, decreases rapidly to 33 Å, and then gradually approaches the bulk value at larger thickness.Even at a pore thickness of 73 Å, however, it has not reached this value.The behavior of CO 2 and CH 4 on protonated, charge-neutral broken (010) edges of the clay layers is broadly similar, but the concentrations of the fluid species near these surfaces are lower than near the basal surfaces, suggesting that they have weaker interactions with protonated edge sites than with siloxane basal surfaces.
The calculations were done at 323 K and a total fluid pressure of 124 bars using a novel approach (constant reservoir composition molecular dynamics, CRC-MD) that maintains a constant composition in the fluid external to the pore.0][51] This purely MD approach overcomes the difficulties in making stochastic particle insertion-deletion moves in dense fluids encountered in grand canonical Monte Carlo and related hybrid approaches.The simulation cell contains a finite size solid substrate, and the approach can be readily applied to many kinds of porous materials such as zeolites, MOFs and C-based materials.It can also be readily applied to pores bounded by non-porous materials, such as silicates and carbonates.
The results here clearly show that nano-and meso-pores in shales bounded by clay minerals have a preference for CO 2 relative to CH 4 and thus methods of enhanced petroleum production such as CO 2 floods or the use of CO 2 -based fracking fluids should readily displace CH 4 from such pores.The results are in good agreement with previous experimental and computational, high T and P studies that show that there can be an energetic driving force for CO 2 incorporation in the interlayer galleries of expandable clay minerals, whereas CH 4 enters these spaces by a passive, space filling mechanism. 6,31,32,38,39,49The thickness of the surface layer in which the fluid is structured (B15 Å) is similar to that in which H 2 O is structured on oxide and hydroxide surfaces, 72,73 suggesting that computational modeling of geochemically relevant fluids in nano-and mesopores bounded by such materials need only include pores with thicknesses less than a few nm to capture the essential features of their behavior.9][50][51][52] The organic component of shale known as kerogen is commonly thought to contain much of the CH 4 , but recent GCMC studies support the conclusion that although it prefers CO 2 relative to CH 4 , the behavior of these species is strongly dependent on the chemical composition, functionality, moisture content and pore structure of the kerogen. 74,75Overall, the structure, dynamics and energetics of the adsorption environments of CO 2 , CH 4 , other hydrocarbons, and H 2 O in kerogen have not been well explored, and the CRC-MD methods used here can be readily applied to help address these important questions.

Fig. 1
Fig. 1 Schematic representations of the simulation cells used in the constant reservoir composition molecular dynamics, CRC-MD calculations of CO 2 /CH 4 partitioning into pores bounded by montmorillonite basal surfaces.(a) The full simulation cell showing the different regions.Pink arrows represent bias forces.(b) Enlarged image of the silt-like pore and montmorillonite T-O-T layers.Color code: Si -yellow; O -red; Mg -green; Al -pink; H -white; Na -blue; C CO 2 -cyan; C CH 4 -black.Substrate atoms are represented using sticks, and the CO 2 /CH 4 fluids and Na + ions are represented in balls.The directions of the X and Z axis are indicated by the arrows on the lower left.

Fig. 2
Fig. 2 (a) Computed mole fractions of CO 2 (black) and CH 4 (red) molecules and the CO 2 /CH 4 ratio (gold) as functions of pore thickness.Solid and dashed lines correspond to left and right y axis, respectively.(b) Number density of the fluid species within the pore region bounded by the oxygens of the montmorillonite basal surfaces.

Fig. 3
Fig. 3 Computed probability density profiles (PDPs) of Na + (orange), C CO 2 (violet) and C CH 4 (magenta) as functions of distance normal to the pore surfaces (montmorillonite basal surface) with varying pore thicknesses.The vertical brown and blue lines represent the positions of oxygen atoms of the montmorillonite basal surface coordinated to tetrahedral Si sites (O b -blue) and tetrahedral Al sites (O bts -brown).

Fig. 4
Fig. 4 Computed PDP's of C CO 2 (violet) and C CH 4 (magenta) near the broken and protonated (010) surfaces.(a) Average concentration within 15.0 Å of the surfaces and its extension above the pore (box on the top of the figure) plotted parallel to the surface.(b) PDP normal to the broken edges (box to the right of the surface shown).The solid and dashed lines represent the distributions on the left and right sides of the montmorillonite model which expose differently substituted sites.The differences between these surfaces is not significant.Green dotted lines in (a) represent the plane of basal surface oxygen atoms.