Article pubs.acs.org/Langmuir

Adsorption Properties of Triethylene Glycol on a Hydrated {101̅4} Calcite Surface and Its Effect on Adsorbed Water Richard Olsen,* Kim N. Leirvik, Bjørn Kvamme, and Tatiana Kuznetsova Department of Physics and Technology, University of Bergen, Allégaten 55, 5007 Bergen, Norway S Supporting Information *

ABSTRACT: Molecular dynamics (MD) and Born−Oppenheimer MD (BOMD) simulations were employed to investigate adsorption of aqueous triethylene glycol (TEG) on a hydrated {101̅4} calcite surface at 298 K. We analyzed the orientation of TEG adsorbed on calcite, as well as the impact of TEG on the water density and adsorption free energy. The adsorption energies of TEG, free energy profiles for TEG, details of hydrogen bonding between water and adsorbed TEG, and dihedral angle distribution of adsorbed TEG were estimated. We found that while the first layer of water was mostly unaffected by the presence of adsorbed TEG, the density of the second water layer was decreased by 71% at 75% surface coverage of TEG. TEG primarily attached to the calcite surface via two adjacent adsorption sites. Hydrogen bonds between water and adsorbed TEG in the second layer almost exclusively involved the hydroxyl oxygen of TEG. The adsorption energy of TEG on calcite in a vacuum environment calculated by classical MD amounted to 217 kJ/mol, which agreed very well with estimates found by using BOMD. Adsorption on hydrated calcite yielded a drastically lower value of 33 kJ/mol, with the corresponding adsorption free energy of 55.3 kJ/mol, giving an entropy increase of 22.3 kJ/mol due to adsorption. We found that the presence of TEG resulted in a decreased magnitude of the adsorption free energy of water, thus decreasing the calcite wettability. This effect can have a profound effect on oil and gas reservoir properties and must be carefully considered when evaluating the risk of hydrate nucleation.

■

INTRODUCTION Calcite is a mineral that has been studied extensively in the literature through experiments,1 as well as through various theoretical methods such as molecular dynamics (MD) simulations.2 It is extremely abundant in the Earth’s crust, and is used within a large number of industries, including CO2 storage; carbonate rock makes a significant portion of all oil and gas reservoirs. Early studies using atomistic simulations (MD in particular) were focused solely on bulk calcium carbonate and were primarily concerned with correct reproduction of its structural, elastic, and vibrational properties.3 With increasingly reliable MD models being developed for the calcium phase, the focus has been turned toward the calcite−water interface. This effort has provided a microscopic insight into the highly structured layers that water forms at the {101̅4} calcite surface.4 There has also been an increasing number of MD studies of other molecular species near the calcite surface, such as investigations into alcohol−calcite interactions,5,6 investigations into adsorption of metal ions on calcite,7 and simulations of peptides near the calcite surface.8,9 Vaterite (a polymorph of calcium carbonate) is metastable, and below 40 °C it transforms into calcite.10 However, above 40 °C vaterite transforms into aragonite.11 It has been found © XXXX American Chemical Society

that use of both monoethylene glycol and tetraethylene glycol stabilizes vaterite.12 Thus, pointing toward the affinity of these diols toward vaterite, it can be assumed that an affinity toward calcite is likely. Furthermore, studies have shown that also other alcohols, such as 1-propanol, ethanol, 2-propanol, and methanol, have an affinity toward calcite.13 The increasing use of water flooding for enhanced oil recovery frequently involves carbonate reservoirs (e.g., Ekofisk offshore Norway, which has been in production since 1970 and still has a high production rate after 25 years of water flooding). Like many offshore hydrocarbon systems, the thermodynamic conditions in the upper few hundred meters are inside a hydrate formation region.14 To prevent formation of hydrate, one may add salts into the reservoir,15 but salts cause corrosion. As an alternative, thermodynamic hydrate inhibitors, such as triethylene glycol (TEG), can be injected into critical parts of a reservoir.16 Barriers for heterogeneous hydrate nucleation on a surface are dependent on wettability. Furthermore, hydrogen bonds toward TEG are the main mechanism preventing hydrate formation.15 Thus, it is of importance to study how adsorbed Received: January 15, 2015

A

DOI: 10.1021/acs.langmuir.5b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir TEG on calcite affects the adsorption free energy of water (a high adsorption free energy magnitude implies a high degree of wetting). It is also of interest to know which adsorption sites are used by adsorbed water versus which are used by TEG to determine possible nucleation mechanisms. The TEG orientation on adsorbed calcite, as well as the chemical potential difference between TEG on the surface and TEG in the bulk, will determine how many hydroxyl groups and ether oxygens are available for hydrogen bonding to water (and therefore how effective TEG is as a hydrate inhibitor close to calcite). Earlier research has found that wettability defines the oil behavior in reservoirs.17 It has been suggested that if forces between water and calcite are strong enough, oil might not be able to displace adsorbed water.18 Thus, within the field of enhanced oil recovery, it is of great importance to know how the wetting properties are altered in a reservoir that has been affected by TEG. Experimental studies have shown that thermodynamic inhibitors such as monoethylene glycol (MEG) and TEG increase calcium carbonate scale formation due to a lowered solubility of calcium carbonate in water.19,20 Knowledge of how adsorbed glycols orient on calcite surfaces can help to obtain a deeper understanding of this effect. Besides being used as a thermodynamic hydrate inhibitor, TEG is also frequently used as a dehydration agent in various processes, for example, in the drying of natural gas,21 which will lead to trace amounts of TEG following industrial pipelines. In our work, we studied how adsorbed TEG on hydrated {101̅4} calcite surfaces (at 298 K and 1 g/cm3 density) affects the surrounding water. We estimated the free energy changes of adsorbed water, studied the orientations of adsorbed TEG, estimated the TEG adsorption energies, and finally estimated the adsorption free energies of TEG. This was achieved through MD simulations by employing a combination of MD force fields reported in the literature (see the Supporting Information and the Computational Details). The adsorption energies of TEG and orientations of adsorbed TEG found using MD were compared to results obtained using ab initio Born− Oppenheimer molecular dynamics (BOMD) simulations. The second section of this paper contains computational details, while the third contains results from postprocessing of simulation data, as well as a discussion of these results. Our findings are summarized in the last section.

■

K ir(ri − r0, i)2

∑

E bond =

i ∈ bonds

and

∑

Eang =

K iθ(θi − θ0, i)2

i ∈ angles

respectively, where r0,i and θ0,i are the equilibrium bond lengths and equilibrium covalent angles, while Kri and Kθi are the respective force constants. The dihedral angles of the molecules were assigned the energy Edih =

1 2

4

∑ ∑ Vn, i[1 + (−1)n+ 1 cos(nϕi)] i ∈ dih. n = 1

where ϕi is dihedral angle i of the system. Figure 1 shows the structure of TEG and the OPLS atom types that have been assigned. Lennard-Jones parameters εi and σi and associated

Figure 1. Structure of the TEG model, where HO, HC1, and HC2 (hydrogen), OH and OS (oxygen), and CT1 and CT2 (carbon) are the atom types.

partial charges qi of HO, OH, OS, and CT1, as well as all bond stretching parameters, covalent angle parameters, and some dihedral angle parameters, were taken from ref 47, while the remaining parameters were found in refs 44−52 or in papers cited within these references. For TEG self-interactions as well as interactions toward H2O and calcite, Lennard-Jones 12−6 short-range potentials combined with Coulomb interactions (parametrized through partial charges, qi, of the atoms) were used. Thus E LJ + C

1 = 2

∑ ∑ (, ) ∀ i ∈ (

⎡ 2 qq ∑ ⎢⎢ e i j 4π ϵ0 rij ∀ j∈) ⎢ ⎣

⎧⎛ ⎞12 ⎛ ⎞6 ⎫⎤ σij ⎪⎥ ⎪ σij + 4ϵij⎨⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎬⎥fij ⎪⎝ rij ⎠ ⎝ rij ⎠ ⎪ ⎩ ⎭⎥⎦

(1)

where ( and ) are sets of atomic indices representing two molecules, respectively, and qi is in units of the electron charge |e|. The distance between atoms i and j is represented by rij = [(ri − rj)2]1/2, where ri and rj are the respective positions of two atoms. Furthermore, f ij = 0.5 if ( = ) (intramolecular interactions) with i ∈ ( and j ∈ ) separated by three bonds, while f ij = 0 if ( = ) with i ∈ ( and j ∈ ) separated by less than three bonds; otherwise, f ij = 1. For different atoms, the geometric mixing rules, σij = (σiσj)1/2 and εij = (εiεj)1/2, were applied. Water Model: Flexible SPC. Water was simulated using the flexible simple point charge (fSPC) model.55 In this model, shortrange interactions and Coulomb interactions are represented by the Lennard-Jones and Coulomb potential of eq 1. Bond stretching for a single water molecule in this model is described by

COMPUTATIONAL DETAILS

Software Packages. MD simulations of this work were performed using the Aug 12, 2013, and Sept 5, 2014, versions of LAMMPS,22,23 while the Nov 21, 2014, version of the collective variables module, Colvars,24−26 was used for constrained dynamics and calculation of free energy profiles through the adaptive biasing force (ABF) method. VMD (Visual Molecular Dynamics) version 1.9.127 was used to represent and study the molecular structures, as well as the molecular motions of the systems. We used the PackMol software package28 to create initial configurations of TEG in the water phase and the GDIS software package29 to build necessary crystal structures. BOMD simulations were performed using CP2K version 2.5.1.30−43 TEG Model. A TEG molecule has the structural formula HO− (CH2)2−O−(CH2)2−O−(CH2)2−OH. To model TEG, the all-atom (AA) OPLS force field44−54 was utilized. The energies of the bonds and covalent angles are therefore represented by the sums of the harmonic potentials

E bond = Dm{[1 − e−ρm Δr1]2 + [1 − e−ρm Δr2]2 } + c{Δr1 + Δr2}Δr3 + b(Δr3)2 + dΔr1Δr2

(2)

where Δr1 and Δr2 are the two H−O bond length perturbations, while perturbations of H−H atomic separations are described by Δr3. These perturbations have the form Δri = ri − ri,0 B

DOI: 10.1021/acs.langmuir.5b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir where ri,0 describe the two H−O equilibrium bond lengths for i = 1 and 2 and the equilibrium H−H separation for i = 3. We used an unmodified version of LAMMPS, which does not allow for higher order cross-terms Δr1Δr3, Δr2Δr3, and Δr1Δr2. These terms were therefore neglected. Thus, c → 0 and d → 0 in eq 2. For interactions toward TEG, we used the Lennard-Jones potential of eq 1 and fSPC model together with geometric mixing rules. Calcite−water interactions were modeled using Buckingham potentials EB =

1 2

⎧ Cij ⎫ ⎪ ⎪ ⎨Aij e−rij / ρij − 6 ⎬ ⎪ rij ⎪ ⎭ ∀ i∈) ⎩

∑ ∑ ∑ (, ) ∀ i ∈ (

(3)

where we used geometric mixing rules Aij = (AiAj) , ρij = (ρiρj)1/2, and Cij = (CiCj)1/2. Water oxygen parameters were taken from de Leeuw et al.,56,57 and water hydrogen parameters were back-calculated using geometric mixing rules.58 Calcite Model. Calcite has the chemical formula CaCO3. Experimental evidence suggests that surface carbonate ions on a hydrated calcite surface only relax approximately 4° away from the bulk positions.59,60 Exact density profiles of water were not of primary concern in our work. The focus was on the overall behavior of the water density profiles, and how they were affected by the presence of adsorbed TEG. Therefore, to reduce the computational requirements, we assumed an idealized calcite structure taken from crystallographic data without further optimization.61 We modeled this by excluding calcium carbonate atoms from MD integrations to fix the crystal lattice in place. Furthermore, short-range interactions between calcium carbonate atoms were turned off to avoid unnecessary excess strain. Short-range interactions between water and calcite were modeled using the Buckingham potential of eq 3 with geometric mixing rules for the parameters. We used Buckingham parameters from Hwang et al.,62 while the partial charges used were previously reparametrized using density functional theory (DFT).63 To estimate the TEG−calcite interaction parameters, we employed the method described by Freeman et al.,64 which is based on existing Buckingham parameters of the pure calcite model. For example, to find the Lennard-Jones parameters for calcite Ca−TEG O interactions, we first fit a Lennard-Jones curve to the calcite Ca−calcite Ca Buckingham curve by solving 1/2

Figure 2. Fits of the Lennard-Jones potentials (dashed lines) to the Buckingham potentials (solid lines) of the calcite Ca−Ca (a, b), O−O (c, d), and C−C (e, f) interactions. (b), (d), and (f) are zoomed in to capture the local minima of (a), (c), and (e), respectively.

slab lengths in the z-direction to achieve a target density, ρtarget, of approximately 1 g/cm3. Thus, the length of the water slabs was set to

∂E BCa − Ca |r = r0 = 0 ∂r with respect to r0 (choosing the local minimum solution) and then requiring that both Ca − Ca ∂E LJ

∂r

Δz =

m TEGNTEG + m H2ONH2O (39.7 Å)(40.3 Å)ρtarget

where m and N are the mass and number of molecules, respectively. We ran six different MD simulations in the NVT ensemble at 298 K, with 0, 1, 16, 32, 64, and 128 TEG molecules present in the bulk slab. Temperature control was achieved using a Nosé− Hoover thermostat40 with a 100 fs thermostat relaxation time. Once the systems were stabile (achieved through short runs using Langevin dynamics combined with a cutoff on the velocity of the atoms) and equilibrated, we ran the system with no TEG for 4 ns, the system with one TEG molecule for 20 ns, and the remaining systems for 40 ns. All systems were run with a 1 fs time step and a 2 Å added distance beyond the force cutoff, while the cutoff was set to 12 Å for Coulomb potentials, Lennard-Jones potentials, and Buckingham potentials. We employed the method of Ewald summation to handle electrostatic contributions across periodic boundary conditions (PBCs). After 1.7 ns the TEG molecule of the system containing a single TEG (with the initial center of mass position 20 Å from the surface calcium atoms) adsorbed on the calcite surface, as shown in Figure 3, which also shows the orientation of water molecules near the calcite surface. A visual inspection shows a

|r = r0 = 0

and ECa−Ca (r0) = ECa−Ca (r0). Thus, we aimed to reproduce the local B LJ minimum of the Buckingham potential. Subsequently, the calcite Ca Lennard-Jones parameters were used with geometric mixing rules toward TEG O. The resulting Lennard-Jones fits to the Buckingham potentials are shown in Figure 2. The above procedure resulted in εCa = 2.5137 kJ/mol, σCa = 4.0525 Å, εO = 0.40057 kJ/mol, σO = 3.0329 Å, εC = 0.39856 kJ/mol, and σC = 3.4704 Å.

■

RESULTS AND DISCUSSION Simulation Details: MD. Systems consisting of a 39.7 Å × 40.3 Å × 28.9 Å rigid wall of calcite (800 CaCO3 units), followed by a water phase of 3820 water molecules and differing amounts of TEG molecules (between 1 and 128 TEG molecules) were constructed and extended using periodic boundary conditions. Our calcite slab was built using GDIS29 from crystallographic information found in ref 61. The initial configurations of water and TEG were set up using PackMol28 to evenly distribute the TEG molecules inside the half of bulk water closest to the {101̅4} surface. Depending on how many TEG molecules were added to the system, we varied the bulk C

DOI: 10.1021/acs.langmuir.5b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

were treated as fully flexible. Furthermore, all systems were assigned a PBC box of 14.97 Å × 24.795 Å × 32.0 Å. For carbon, oxygen, and hydrogen, we used the TZV2PGTH basis set (i.e., triple-ζ valence and two polarizations, also known as 6-311(2df,2pd)), while for calcium we used the DZVP-GTH-PADE-q10 basis set (i.e., double-ζ valence and single polarization, also known as 6-31G*). These were used with GTH-PADE-q4, GTH-PADE-q6, POTENTIAL GTHPADE-q1, and GTH-PADE-q10 Goedecker−Teter−Hutter pseudopotentials for carbon, oxygen, hydrogen, and calcium, respectively.32,34,37 DFT calculations were done with five multigrids and a cutoff of 300 Ry, while the target accuracy for the SCF convergence was set to 10−6. Simulations of calcite with a single adsorbed TEG, calcite in a vacuum, and a single TEG in a vacuum were run for 10.6, 12, and 7.4 ps, respectively, with a time step of 0.5 fs. In addition, we ran the same simulations again, but including a Grimme D2 dispersion correction using a BLYP reference function with a potential cutoff of 20 Å.33 These simulations ran for 9.7, 12.6, and 9.9 ps, respectively. TEG Impact on Water. Water Density Profile. Using trajectories from time periods where TEG molecules were adsorbed on the surface, we constructed water density profiles (with a bin size of 0.15 Å). Figure 4a shows the density profiles for various amounts of TEG in the system, obtained by

Figure 3. Distribution plots (also known as heat maps) of water oxygen (red), water hydrogen (blue), TEG hydroxyl (orange), and the remainder of TEG (black) taken over the last 1000 frames (400 ps) of simulation containing a single TEG molecule adsorbed to the {101̅4} calcite surface.

high oxygen affinity toward the calcite surface, a feature that has also been found by other groups.65,66 In the last 4 ns of each simulation run, 1, 9, 19, 23.75, and 29.75 TEG molecules had adsorbed (on average) to the lefthand surface (at z = 0 Å) within systems of 1, 16, 32, 64, and 128 TEG molecules, respectively. Furthermore, 6, 12, 18.5, and 19.25 TEG molecules had adsorbed (on average) to the righthand surface within systems of 16, 32, 64, and 128 TEG molecules, respectively. Estimation of the adsorption energies of TEG on hydrated calcite and calcite in a vacuum environment also involved performing a 4 ns MD simulation of a TEG molecule in a vacuum using the same simulation parameters as in the case of a single TEG on hydrated calcite. Similarly, the adsorption energies of TEG on a calcite surface surrounded by a vacuum were evaluated from an additional 4 ns MD simulation of TEG adsorbed on calcite (without water) using the same simulation parameters. The total energy of TEG at infinite dilution in water was determined by comparing the total energy of the aqueous system with that of a pure water slab with the same dimensions and the same number of water molecules. Since our calcite structure was fixed, the energy of calcite in a vacuum was available, did not require any time averaging, and only consisted of electrostatic contributions. Simulation Details: BOMD. BOMD simulations were used to study TEG orientations on the {1014̅ } calcite surface, as well as TEG adsorption energies on the same surface, in a vacuum environment at fixed NVT, where T = 298 K. We employed a Nosé−Hoover thermostat40,41 with a 50 fs thermostat relaxation time, a chain length of 3, a third-order Yoshida integrator, and two multiple-time steps. We constructed three separate systems aimed at estimating adsorption energies: calcite with a single adsorbed TEG, calcite in a vacuum, and a single TEG in a vacuum. The calcite slab dimensions were set to 14.97 Å × 24.5 Å × 7.6 Å (54 CaCO3 units), and the slab was constructed using GDIS29 from crystallographic information found in ref 61. All atoms of calcite

Figure 4. (a) Water density profiles with and without TEG adsorbed on the surface. The Ca atoms of the calcite surface are placed at z = 0 Å, while the shaded area shows approximately where the TEG molecule is located in the z-direction. (b) Density of the second and third layers of water as a function of the density of adsorbed TEG on the calcite surface at 298 K. Filled circles and triangles are the amplitudes of the density peaks at 4 Å (see (a)). Dashed and dotted lines are from the analytical expression of eq 7. A surface density of 62.4 ng/cm2 corresponds to a surface saturated by TEG. (c) Number of H-bond donors (per frame) of water distributed along the z-axis into bins of 0.05 Å. Inset: Hw−Ow H-bond profile for one and nine adsorbed TEG molecules. (d) Number of H-bonds (per frame) distributed into energy bins of 2 kJ/mol. Hw and Ow are hydrogen and oxygen of water, and Ooh and Oet are hydroxyl and ether oxygen of TEG, respectively. D

DOI: 10.1021/acs.langmuir.5b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

We solve for ( and ) by requiring that the water density reaches a minimum, ρHmin2O, at some maximum amount of adsorbed TEG, ρmax ads , and that its first derivative with respect to ρads vanishes at this point to ensure continuity. The maximum of adsorbed TEG can be estimated roughly using the radius of gyration of TEG at the calcite surface

averaging over the last 4 ns of the respective simulations. As can be seen from the figure, the distance to the first layer of water is approximately 1.88 Å. Furthermore, the density due to water oxygens oriented toward the surface (i.e., the first peak) and the density of the hydrogens of water (i.e., the second peak) in the first water layer were not significantly altered by the presence of TEG on the surface. This was most likely because the first layer of water was adsorbed to calcium atoms while TEG was adsorbed to oxygens and therefore was not affected by the adsorbed TEG (see Figure 3). However, the water density in the second layer of water (the peak at approximately 3.75 Å) was affected by adsorbed TEG molecules, where the density was decreased by approximately 4%, 28%, 58%, 69%, and 71% due to 1, 9, 19, 24, and 30 adsorbed TEG molecules (on average), respectively. The density of the third layer of water was also decreased due to the presence of adsorbed TEG. For the system with no TEG, we found that one water molecule adsorbed per Ca atom. We also estimated the average number of adsorbed water molecules per Ca atom of calcite in the first two layers of water to be nw =

A 10 NCam H2O 24

∫0

5

N

rg =

2

where A, NCa = 80, mH2O, and ρH0 2O(z) are the area of the calcite surface in square angstroms, number of calcium atoms at the surface, mass of water in grams, and density profile of pure water (where z = 0 is located at the position of the surface calcium atoms), respectively. This result is consistent with reported experiments of water films on calcite.18 Our estimate of the water density profile has features similar to those obtained using ReaxFF,67,68 including a small second peak due to the strong orientation of oxygens toward the surface.4 Note that the large height of peaks in the water density can be attributed to the small bin size used. Semiempirical Expression for Water Density. In the case in which adsorbed TEG molecules are far apart (and z is at a distance from the surface where the water density is affected), we may assume that

max max ρads = γNads /A

The minimum mass of water that can reside in any neighborhood of z can be roughly approximated by m Hmin = (Vvac − VTEG)ρHe O 2O 2

where Vvac = 2Arg is the volume surrounding the TEG molecules, VTEG = (4/3)πrg3Nmax ads is the volume of adsorbed TEG molecules, and ρHe 2O ≈ 0.997 g/cm3 is the density of water. On the basis of this assumption, the minimum amount of water at maximum adsorption of TEG is ρHminO 2

(4)

for the corresponding densities, ρ, of water, where ρads is the surface density of adsorbed TEG and α̃ is a constant related to α. However, close to already adsorbed TEG molecules, there will be fewer water molecules available to push away; i.e., α will be lower. Thus, we modify eq 4 to ρH O (z) = ρH0 O (z) − α(̃ ρads )ρads 2

2

(6)

with M = 14mH + 4mO + 6mC, where mH, mO, and mC are the atomic masses of hydrogen, oxygen, and carbon, given in grams per mole (NA is Avogadro’s number). We denote the conversion factor by γ. Thus

where NH2O(z) and NH0 2O(z) are the number of water molecules located in the neighborhood of z, with and without any adsorbed TEG molecules, respectively. NTEG is the number of adsorbed TEG molecules, and α is the number of water molecules that are removed per TEG adsorption. Equivalently 2

= 2.64 ± 0.06 Å

9 x TEG molecules x 10 M /NA ng = A A 10−16 cm 2 Å2 ng x = 2492.08 A cm 2

NH2O(z) = NH0 2O(z) − αNTEG

2

N

∑i =m1 ωi

where Nm = 24 is the number of atoms in TEG, rc is the weighted center of the molecule, and ri is the position of the ith atom of TEG. We set weight factors, ωi, to the appropriate atomic masses. If we assume that each adsorbed TEG requires approximately an area ΔA = 4rg2 on the surface, where each TEG molecule is organized in a grid of squares, we find that N = A/ΔA ≈ 57.3 is the maximum number of TEG molecules that can fit on the surface of area A. However, there are only 80 CaCO3 groups on our calcite surface, where we observed from our MD simulations that each TEG adsorbed onto two adsorption sites, leaving a maximum of 40 adsorbed TEG molecules. Thus, we assumed that Nmax ads = 40. Conversion from x TEG molecules per area, A (Å2), was achieved by

dz ρH0 O (z) ≈ 2.7 per Ca

ρH O (z) = ρH0 O (z) − αρ̃ ads

∑i =m1 ωi(rc − ri)2

=

m Hmin 2O Vvac

max ⎞ ⎛ 2πrg 2Nads ⎜ ⎟ρ e = ⎜1 − ⎟ H 2O 3A ⎝ ⎠

After some algebra and evaluating the coefficients ( and ) , we find the complete expression ⎧ ⎞ ⎪ A ⎛ Aρads ρH O (z , ρads ) = ⎨1 + max ⎜ max − 2⎟ 2 ⎪ γNads ⎝ γNads ⎠ ⎩ max ⎤ ⎛ ⎞ ⎫ ⎡ ρe 2πrg 2Nads ⎪ ⎥ H2O ⎟ρ ⎬ρ 0 (z) × ⎜⎜1 − ⎢1 − ads H 2O 0 ⎟ ⎢⎣ ⎥⎦ ρH O (z) ⎠ ⎪ 3A ⎝ ⎭ 2

(5)

where we for simplicity ignore the z dependence of α̃ (ρads) and choose a simple linear form:

(7)

where only the radius of gyration (rg) and water density peak at z (which was found to be ρH0 2O ≈ 2.849 g/cm3 for the second

α(̃ z , ρads ) = (ρads + ) E

DOI: 10.1021/acs.langmuir.5b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir layer of water and ρH0 2O ≈ 1.178 g/cm3 for the third layer of water from our simulations) need to be determined. Figure 4b shows the density of the second and third layers of water as a function of the density of adsorbed TEG as calculated using eq 7 together with values obtained from our MD simulations. We found a good agreement between the two (the standard deviation was found to be 0.058 g/cm3 for the second layer and 0.085 g/cm3 for the third layer). Impact on the Free Energy Profile of Water. The free energy difference between water M at z and water M at z0 can be estimated by (see the Supporting Information) ⎡ P(z , ρ ) ⎤ ads ⎥ ΔA H2O(z , ρads ) = −RT ln⎢ ⎢⎣ P(z 0 , ρads ) ⎥⎦

where P(z,ρads) is the probability of water M being present at z with adsorbed TEG density ρads. Since this probability is also proportional to the density ρH2O(z,ρads), we can estimate the free energy difference, ΔAH2O, between water at z and water in the bulk by69 ⎡ ρ (z , ρ ) ⎤ HO ads ⎥ ΔA H2O(z , ρads ) = −RT ln⎢ 2 bulk ⎢ ⎥ ρ ⎣ ⎦ H 2O

Figure 5. Perturbation of the free energy profile of water, ΔApert, close to the calcite surface (where the water density is affected by TEG) due to adsorbed TEG at 298 K, where ρH0 2O(z) is the density profile of water with no adsorbed TEG on the calcite surface and ρads is the surface density of TEG. A surface density of 62.4 ng/cm2 corresponds to a surface saturated by TEG.

(8)

where ρHbulk is the water density in the bulk phase and 2O ρH2O(z,ρads) is described in eq 7. It is important to realize that eq 8 is not valid when calculating a free energy difference between two systems with different functional forms of the potential and kinetic energy terms entering into partition functions used to calculate P. This would result in probabilities, P, with different proportionality factors to the density, ρ, which would not cancel in a probability fraction. Equation 8 can be rewritten as

energy of water, at maximum adsorption of TEG, has a positive value (above 0.63 g/cm3 at ρads = 62.4 ng/cm2, the perturbation is positive; see Figure 5) and lowers the magnitude of the adsorption free energy of water. Hence, water wetting decreases. Hydrogen Bonding between Water and TEG. We investigated how water oriented itself close to adsorbed TEG, in an MD simulation with one adsorbed TEG, by looking at the hydrogen bond (H-bond) statistics (based on the last 6 ns of simulation). Oxygens of water (Ow) acted as donors, and hydroxyl oxygens (Ooh), as well as ether oxygens (Oet), of adsorbed TEG acted as acceptors. We employed purely geometric H-bond criteria,70,71 where the presence of a Hbond was assumed if (1) the distance between the donor (d) and acceptor (a), da···d, was less than dca···d, (2) the distance between the donor hydrogen (h) and acceptor, da···h, was less than dca···h, and (3) the angle, αa···h−d, was greater than 150°. The cutoff distances dca···d and dca···h were set to the first minima of the donor−acceptor and donor−hydrogen−acceptor radial distribution functions (RDFs), respectively. We then found dcOoh···Ow ≈ 3.5 Å, dcOet···Ow ≈ 3.5 Å, dcOoh···Hw ≈ 2.5 Å, dcOet···Hw ≈ 2.2 Å, and dOw···Hoh ≈ 4 Å, where Hw and Hoh denote the hydrogen of water and hydroxyl oxygen of TEG, respectively. No H-bonds were found between the hydroxyl hydrogens of TEG and water. Thus, the oxygen of water did not orient strongly toward the hydroxyl hydrogens of TEG. Figure 4c shows how water formed H-bonds to both hydroxyl oxygens (Ooh) and ether oxygens (Oet) of TEG in layers outside the calcite surface (in the z-direction). By comparison with Figure 4a, it can be seen that the hydrogens of the first layer of water bind (by between 0.2 and 0.3 H-bond per time frame) to both Ooh and Oet and are therefore likely to orient toward them. However, in the second layer of water, the hydrogens of water only have a significant orientation toward the hydroxyl oxygens of TEG. The distribution in Figure 4d shows the strength of these H-bonds to be on the order of −130 to −230 kJ/mol.

⎡ ρ (z , ρ ) ⎤ HO ads ⎥ ΔA H2O(z , ρads ) = −RT ln⎢ 2 0 ⎢⎣ ρH O (z) ⎥⎦ 2 ⎡ ρ 0 (z ) ⎤ HO − RT ln⎢ 2bulk ⎥ ⎢ ρ ⎥ ⎣ H 2O ⎦

Thus ΔA H2O(z , ρads ) = ΔA H0 2O(z) + ΔA pert (z , ρads )

(9)

where ⎡ ρ (z , ρ ) ⎤ HO ads ⎥ ΔA pert (z , ρads ) ≡ −RT ln⎢ 2 0 ⎢⎣ ρH O (z) ⎥⎦ 2

is a perturbation of the water free energy profile due to adsorbed TEG and ΔAH0 2O(z) is the free energy profile of water when no TEG is adsorbed. Note that eq 9 is only valid in the range of z where the water density is affected by TEG adsorption. Outside this range, ΔAH0 2O(z) applies. Figure 5 shows a plot of ΔApert(ρH0 2O(z),ρads) at 298 K which implies that the free energy changes of water due to adsorbed TEG should be in the range of ±4 kJ/mol. The majority of points in the water density profile (with no TEG adsorbed to the surface) were found to be greater than 0.63 g/cm3. Therefore, our estimated perturbation to the free F

DOI: 10.1021/acs.langmuir.5b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

Figure 6. TEG orientation toward the {101̅4} calcite surface for (a) an MD simulation of completely adsorbed TEG, (b) an MD simulation of partially adsorbed TEG, (c) a BOMD simulation of completely adsorbed TEG, and (d) a similar BOMD simulation with dispersion correction. The parameters di and αi are shown in Table 1. The nomenclature for distances and angles in (b)−(d) is equivalent to that in (a). (e) Dihedral angle distributions for six dihedrals of TEG. Distribution plots are normalized such that the area under the curves yields unity. Gray solid lines are due to a single TEG molecule in bulk water, while black solid lines are due to a TEG molecule adsorbed on the calcite surface.

Table 1. Values of v from Figure 6a−c with Standard Deviations, σ, Based on the Last 400 ps of MD Simulations for (I) TEG Where Both Hydroxyls Are Adsorbed on Hydrated Calcite, (II) TEG Where Only One of the Hydroxyls Has Adsorbed to the Hydrated Surface, and (III) TEG Adsorbed on Calcite in a Vacuum, (IV) the Last 6 ps of the BOMD Simulation Where TEG Is Adsorbed onto Calcite in a Vacuum, and (V) a Similar BOMD Simulation with Dispersion Correction I

II

III

IV

V

param

v

σ

v

σ

v

σ

v

σ

v

σ

d1 (Å) α1 (deg) d2 (Å) α2 (deg) d3 (Å) d4 (Å) d5 (Å) d6 (Å)

2.5 150.0 1.49 169.0 3.4 4.1 4.6 5.4

0.1 6.7 0.09 5.7 0.1 0.2 0.3 0.5

2.4 151.4 1.41 171.0 3.3 4.0

0.1 5.7 0.08 4.9 0.1 0.2

2.5 131.4 1.6 151.0 3.2 3.8 3.3 4.6

0.1 7.4 0.1 8.8 0.1 0.2 0.1 0.3

2.6 141.0 2.0 150.7 3.4 2.6 2.7 4.6

0.2 8.2 0.2 12.1 0.1 0.1 0.1 0.3

3.1 139.6 2.0 128.8 3.9 2.5 2.6 5.4

0.3 14.8 0.3 14.3 0.2 0.1 0.1 0.2

As water−water H-bond criteria, we have used the first minima of the RDFs, which we found to be 3.3 Å between the oxygens of water and 2.4 Å between the hydrogens and oxygens of water (excluding hydrogens on the same water molecule). The inset of Figure 4c shows water−water H-bonding as a function of the position, z, of the donor atoms, both with one TEG molecule and with nine TEG molecules adsorbed to the surface (where the statistics were based on the last 800 ps of the simulations). It can be seen that H-bonds between water molecules were reduced in the presence of adsorbed TEG. TEG Structuring. Orientation Relative to the Surface. Snapshots of TEG after adsorption to calcite in MD simulations and in BOMD simulations are shown in Figure 6. When only a single TEG was present (in parts a, c, and d), we observed a complete adsorption of both TEG hydroxyl groups. However, when enough TEG had adsorbed onto the surface, TEG molecules were allowed to adsorb using only one of the two hydroxyl groups (see part b). As alluded to in the Introduction, tetraethylene glycol has an affinity toward vaterite (a

polymorph of calcium carbonate). In addition, alcohols tend to have an affinity toward calcite. Thus, it is not unreasonable that TEG adsorbed to calcite in our simulations. Figure 3 contains density plots of fully adsorbed TEG, water, and TEG hydroxyl groups to illustrate molecular motion during the last 400 ps of MD simulation when the snapshot for Figure 6a was taken. It can be seen that TEG hydroxyls only deviated slightly from the positions shown in Figure 6. This is further quantified in Table 1, which tabulates the time-averaged distances, di, and angles, αi, between the TEG atoms and atoms of calcite together with the corresponding standard deviations. The comparisons of the distances and angles (di and αi) given by BOMD simulations (IV, V) with those obtained via MD simulations of TEG and calcite in a vacuum (III) indicated that only the distances between TEG oxygen and calcite calcium, as well as the distances between TEG hydroxyl hydrogen and carbonate oxygen, differed significantly (d4 by 46%, d5 by 22%, and d2 by 20%). Furthermore, in the case of BOMD simulation without dispersion correction, both TEG G

DOI: 10.1021/acs.langmuir.5b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

used, where we in addition included pressure regulation using a Nosé−Hoover barostat74 with a 5 ps relaxation time. From this simulation we estimated the dihedral distributions of six dihedrals of TEG in bulk water. These distributions, S(θ), were divided into 720 bins of width Δθ = 0.25° and normalized such that

hydroxyls attached to a single adsorption site, in direct contrast to our MD simulations where each hydroxyl group of fully adsorbed TEG coordinated to its separate adsorption site (as checked by counting the carbonate oxygens within 2.6 Å of the hydroxyl hydrogens of TEG and subsequently verifying that both hydrogens adsorbed to separate oxygens). However, when dispersion correction was applied in BOMD simulations, each TEG hydroxyl hydrogen adsorbed onto a separate carbonate oxygen (see Figure 6d and Table 1 (simulation V)). For both MD simulations of TEG adsorption on calcite in a vacuum and BOMD without dispersion correction, adsorption took place on the oxygens of carbonate, and (due to the flexibility of the TEG backbone) in both cases each adsorbed TEG occupied approximately the same space on the surface (i.e., two carbonate groups). The radius of gyration for MD vacuum simulation of TEG on calcite was found to be 2.57 ± 0.07 Å, using eq 6, while for the BOMD simulation it was found to be 2.45 ± 0.05 Å; these values are in very good agreement. However, for BOMD simulation with dispersion correction, the radius of gyration was found to be larger (3.2 ± 0.05 Å) due to the adsorption sites not being located on nearest neighbor carbonate groups. From Table 1 we see that for MD simulations of TEG on hydrated calcite surfaces (I and II) the ether oxygens of TEG are further away from the calcite calcium (see d5 and d6) than for MD simulation of TEG on calcite in a vacuum (III). This effect can be explained by water molecules adsorbing onto calcite between the adsorption sites of TEG (see Figure 3). It is also interesting to compare the hydroxyl placements of TEG with those of an ethanol hydroxyl group. Simulations done by other groups6 show a distance of 1.4 Å between the hydroxyl hydrogen of ethanol and the carbonate oxygen and a distance of 2.4 Å between the hydroxyl oxygen and closest calcite calcium (corresponds to d2 and d4 of Figure 6, respectively). This corresponds well with the distance of 1.45 Å we found between the hydroxyl hydrogen of TEG and the closest carbonate oxygen. However, we found that the distance between the hydroxyl oxygen of TEG and the closest calcite calcium atom (d4 ≈ 3.8) was much longer than the corresponding distance for ethanol. Thus, the hydroxyl groups of TEG will orient themselves differently toward the calcite surface than the hydroxyl group of ethanol, but will be located at approximately the same distance from the surface. It has been reported that the solubility of calcium carbonate (calcite) in water decreases if the water contains glycol.19,20,72 A possible explanation could be that the hydroxyl groups, responsible for the hydrophilic behavior, point toward calcite. Thus, the “hydrophobic part” is left exposed to water, yielding a lower solubility of calcite. From Figure 6b we see that, under certain conditions, TEG molecules may arrange themselves with one hydroxyl pointing toward calcite and the other toward the water phase. Thus, under conditions where enough TEG is present and is packed tightly enough on the surface, instead the “hydrophilic part of TEG” would be exposed to water, similar to a Langmuir− Blodgett film, resulting in a different effect on the surrounding water than what we see in Figure 4a,b for moderate amounts of TEG at 298 K and 1 g/cm3 density. Dihedral Distributions. We performed a 37 ns NPT simulation at 298 K and 1 atm (corresponds to approximately 1 g/cm3 density) with 1 TEG molecule and 1500 water molecules (but using MDynaMix73 version 5.1). The simulation parameters of the system with TEG on hydrated calcite were

∫0

π

d θ S(θ ) = 1

Thus, in the discrete limit, we calculated 720

C=

∑ SjΔθ j=0

and subsequently let Sj → Sj/C. The resulting distributions are shown in Figure 6e (as gray solid lines) and are compared with the dihedral distributions of TEG adsorbed on the calcite surface (which are shown as black solid lines). It can be seen that the largest difference between the dihedral distributions of TEG adsorbed on the surface and TEG in the bulk appears in the OCCO dihedrals where the 180° peak vanishes for the adsorbed molecule. A second effect that can be observed is the additional peak at 150° in the HOCC dihedral. This could be due to additional vibrations on the hydroxyl groups because of their attachment to the calcite surface. The remaining dihedral distributions for TEG adsorbed on the calcite surface did not differ significantly compared to those for TEG in bulk water. TEG−Surface Energies. Adsorption Energy. Adsorption energies were approximated both at hydrated calcite and in a vacuum environment. For hydrated calcite, we calculated the time-averaged total energy (including long-range interactions with a 12 Å cutoff and corresponding corrections from Ewald summation) of the simulation with hydrated calcite, Ecal+H2O, TEG in water, EH2 O+TEG, pure water, EH2O, and hydrated calcite including a single adsorbed TEG, Ecal+H2O+TEG. Thus, the adsorption energy was calculated as (h) Eads = (Ecal + H 2O + [E H 2O + TEG − E H 2O]) − Ecal + H 2O + TEG

(10)

In a vacuum environment, using nomenclature equivalent to eq 10 (v) Eads = (Ecal + E TEG) − Ecal + TEG

(11)

For comparison we also estimated the time-averaged nonbonded energy contributions between TEG and the atoms of calcite in simulations containing a single TEG molecule: UTEG − calcite =

∑ i ∈ calcite

⎧ ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σ σ ⎪ ∑ ⎨4εij⎢⎢⎜⎜ ij ⎟⎟ − ⎜⎜ ij ⎟⎟ ⎥⎥ r ⎝ rij ⎠ ⎦ j ∈ TEG ⎪ ⎩ ⎣⎝ ij ⎠

⎫ e 2 qiqj ⎪ ⎬f + E EWALD + 4πε0 rij ⎪ ij ⎭

(12)

where EEWALD is the contribution that originates from periodic boundaries through the method of Ewald summation. Similarly, we estimated the interactions between TEG and water, UTEG−water. H

DOI: 10.1021/acs.langmuir.5b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

constant boundaries were both set to k = 50.208 kJ/(mol·m2). The three simulations were run with z ∈ [−35, −25] Å, z ∈ [−25, −15] Å, and z ∈ [−15, 0] Å, respectively. This resulted in gradient and sample information in the three ranges that were used as a starting point for a 60 ns ABF simulation with z ∈ [−35, 0] Å. The resulting free energy profile is shown in Figure 7. We see no evidence of any significant local energy

Table 2 shows the estimated adsorption energies. Averages were taken over the last 4−19 ns of the MD simulations and Table 2. Adsorption Energies of TEG Found from MD Simulations in a Vacuum Environment (v), MD Simulations of Hydrated Calcite (h), BOMD Simulations in a Vacuum, and Dispersion-Corrected BOMD Simulations in a Vacuum (d)a E(h) ads (kJ/mol) MD (v) MD (h) BOMD BOMD (d) a

E(v) ads (kJ/mol) 217

33

UTEG−calcite (kJ/mol) −262 −184

UTEG−water (kJ/mol) −303

119 280

Energies are described in eq 10, eq 11, and eq 12.

over the last 6−12 ps of the BOMD simulations (depending on the length of the simulations and equilibration period). It can be seen that the adsorption energy of TEG found from MD simulations in a vacuum is 82% higher than that found from BOMD simulations in a vacuum and 22.5% lower than that found from BOMD simulations including dispersion corrections, indicating good agreement. Hydrating the calcite surface resulted in a 184 kJ/mol decrease in the adsorption energy of TEG, thus making adsorption of TEG more favorable when in a vacuum environment. This can also be seen from TEG−calcite pair interaction energies, which have a lower magnitude at hydrated calcite (see Table 2). As seen from Table 2, interactions between TEG and water are strong, indicating that adsorbed water partly screens TEG from the surface. Our adsorption energies can also be compared to the 106 kJ/ mol adsorption energy for 1-propanol, estimated using DFT by Okhrimenko et al.,13 which is also of the same order of magnitude. Free Energy Profile. TEG was pulled perpendicularly off the {101̅4} surface by running the simulation of adsorbed TEG on hydrated calcite again, but with a harmonic potential 1 /harm = k (z − z 0)2 2 added to the system Hamiltonian to restrain the TEG molecule. The variable z is defined as the projected distance along the zaxis (perpendicular to the {101̅4} surface of calcite) between the geometric center of TEG and the geometric center of all water in the system. The restraint was allowed to move by gradually allowing z0 to change from z0 = −34.5 (corresponding to ∼4.3 Å relative to the surface calciums of calcite) to z0 = 0 (i.e., ∼38.8 Å relative to the surface calciums of calcite) throughout a 2 ns simulation period. The force constant was set to k = 50.208 kJ/(mol·m2). Constrained dynamics was achieved by using the Colvars24 software package. Accumulated work at 40 Å was found to be approximately 400 kJ/mol, which serves as a lower limit (but not the lowest) of energy required to pull TEG off the surface. Information from the trajectory resulting from pulling TEG off the surface was used as input in subsequent ABF simulations. Next, we employed the Colvars module to run three separate 20 ns ABF simulations25,26 of hydrated calcite with adsorbed TEG where the collective variable, z, was defined with a grid spacing of 0.1 Å and fullSamples was set to 500 (allowing for scaling of the biasing force). The lower and upper force

Figure 7. Free energy profile obtained using adaptive biasing force (ABF) dynamics (z = 0 is the location of the surface Ca atoms).

maxima of the free energy near the surface. Averaging the graph in Figure 7 from 11.2 Å, we found an adsorption free energy of ΔA = 55.3 ± 0.5 kJ/mol

for TEG on the hydrated calcite surface. Using the adsorption energy for TEG on hydrated calcite (see Table 2), we find an entropy decrease of ΔU − ΔA ≈ −0.07 kJ/(mol ·K) T Thus, the entropy increases for a system where TEG adsorbs to calcite (ΔU = 33 kJ/mol and T = 298 K). The overall increase in the system entropy may be explained by an increase in the water and TEG degrees of freedom (by removal of TEG from the bulk), which is larger than the corresponding decrease of the degrees of freedom due to constraining TEG on the surface. ΔS =

■

SUMMARY AND CONCLUSIONS Performing MD simulations of hydrated calcite with differing amounts of TEG on a {1014̅ } calcite surface, at 298 K and 1 g/ cm3 density, we found that the overall water density profile showed a 71% decrease in the second water layer with 75% of the available adsorption sites occupied. Furthermore, it was found that the first layer of water was unaffected by the presence of adsorbed TEG (while both the second and third layers were attenuated). This was attributed to TEG and water having different sites of adsorption on the calcite surface. At low TEG mole fractions, each hydroxyl group of TEG adsorbed to adjacent adsorption sites, resulting in the hydrophobic part of TEG pointing toward the bulk phase of water, thus offering an explanation to the reported reduction in solubility of calcium carbonate in water solutions with TEG. Furthermore, such an orientation of TEG would decrease the number of hydroxyl groups available for hydrate inhibition. It was also found that, at higher mole fractions of TEG, adsorption could occur by only one of the TEG hydroxyl groups adsorbing, leaving the hydrophilic part pointing toward bulk water. Thus, it may be that a Langmuir−Blodgett-like film could form on the surface under the right conditions. The I

DOI: 10.1021/acs.langmuir.5b02228 Langmuir XXXX, XXX, XXX−XXX

Langmuir orientations of TEG adsorbed on calcite found by MD simulations were found to correspond well with those found by ab initio simulations performed on similar systems. We found perturbations of the free energy profile of water, due to adsorbed TEG, to be positive but below 4 kJ/mol (in the range where the water density was affected), thus resulting in a lowered magnitude of the adsorption free energy of water and therefore decreased water wetting, which is important information both to determine the risk of hydrate nucleation and to predict oil behavior in reservoirs consisting of hydrated calcite with adsorbed TEG. It was found that water molecules of both the first and second water layers could orient their hydrogens toward the adsorbed TEG hydroxyl oxygens and form H-bonds (with nonbonded energies of approximately −230 kJ/mol). However, H-bonds of water to adsorbed TEG ether oxygens were almost exclusively by waters in the first layer (with nonbonded energies of approximately −133 kJ/mol). H-bonds to TEG from the first layer may be attributed to the specific orientation of water adsorbed to calcite. Dihedral distribution peaks did not differ significantly between TEG in bulk water and TEG adsorbed on hydrated calcite. The most significant difference was the vanishing of the 180° peak of the OCCO dihedrals, which was due to bending of the TEG molecule resulting from both hydroxyls attaching to the calcite surface. The adsorption energy of TEG adsorbed to hydrated calcite was found to be 33 kJ/mol using MD simulations. MD simulations of TEG on calcite in a vacuum environment yielded an adsorption energy equal to 217 kJ/mol, which was only 22.5% lower than the corresponding energy found from ab initio simulations. We estimated the adsorption free energy of TEG to be approximately 55.3 kJ/mol, thus yielding an entropy increase, TΔS, of approximately 22.3 kJ/mol by adsorption to calcite at 298 K. In this work we applied only one model for TEG (the OPLS model), one model for water (fSPC), and one model for calcite. It is not given that the chosen combination of models is the most accurate at predicting physical reality. Thus, further work is required to verify the predictability of the chosen models, and to possibly perform necessary adjustments. However, our work still yields valuable insight into the adsorption behavior of TEG onto calcite and may pave the way to construct improved models for similar systems.

■

■

ACKNOWLEDGMENTS

■

REFERENCES

We acknowledge financial support from Statoil ASA, the academy agreement between Statoil and the University of Bergen, and the Research Council of Norway through the project CLIMIT “Safe long term sealing of CO2 in hydrate” (Project Number 224857), as well as computational resources provided by the Bergen Center for Computational Science (BCCS).

(1) Heberling, F.; et al. Reactivity of the calcitewater-interface, from molecular scale processes to geochemical engineering. Appl. Geochem. 2014, 45, 158−190. (2) Hamm, L. M.; Bourg, I. C.; Wallace, A. F.; Rotenberg, B. Molecular Simulation of CO2- and CO3-Brine-Mineral Systems. Rev. Mineral. Geochem. 2013, 77, 189−228. (3) Pavese, A.; Catti, M.; Price, G.; Jackson, R. Interatomic Potentials for CaCO3 Polymorphs (Calcite and Aragonite), Fitted to Elastic and Vibrational Data. Phys. Chem. Miner. 1992, 19, 80−87. (4) Fenter, P.; Kerisit, S.; Raiteri, P.; Gale, J. D. Is the Calcite-Water Interface Understood? Direct Comparisons of Molecular Dynamics Simulations with Specular X-ray Reflectivity Data. J. Phys. Chem. C 2013, 117, 5028−5042. (5) Sand, K. K.; Yang, M.; Makovicky, E.; Cooke, D. J.; Hassenkam, T.; Bechgaard, K.; Stipp, S. L. S. Binding of Ethanol on Calcite: The Role of the OH Bond and Its Relevance to Biomineralization. Langmuir 2010, 26, 15239−15247. (6) Cooke, D. J.; Gray, R. J.; Sand, K. K.; Stipp, S. L. S.; Elliott, J. A. Interaction of Ethanol and Water with the {101¯4} Surface of Calcite. Langmuir 2010, 26, 14520−14529. (7) Kerisit, S.; Parker, S. C. Free Energy of Adsorption of Water and Metal Ions on the {101¯4} Calcite Surface. J. Am. Chem. Soc. 2004, 126, 10152−10161. (8) Yang, M.; Stipp, S. L. S.; Harding, J. Biological Control on Calcite Crystallization by Polysaccharides. Cryst. Growth Des. 2008, 8, 4066− 4074. (9) Yang, M.; Rodger, P. M.; Harding, J. H.; Stipp, S. Molecular dynamics simulations of peptides on calcite surface. Mol. Simul. 2009, 35, 547−553. (10) Rodriguez-Blanco, J. D.; Shaw, S.; Benning, L. G. The kinetics and mechanisms of amorphous calcium carbonate (ACC) crystallization to calcite, via vaterite. Nanoscale 2011, 3, 265−271. (11) Ogino, T.; Suzuki, T.; Sawada, K. The formation and transformation mechanism of calcium carbonate in water. Geochim. Cosmochim. Acta 1987, 51, 2757−2767. (12) S̆kapin, S. D.; Sondi, I. Synthesis and characterization of calcite and aragonite in polyol liquids: Control over structure and morphology. J. Colloid Interface Sci. 2010, 347, 221−226. (13) Okhrimenko, D. V.; Nissenbaum, J.; Andersson, M. P.; Olsson, M. H. M.; Stipp, S. L. S. Energies of the Adsorption of Functional Groups to Calcium Carbonate Polymorphs: The Importance of -OH and -COOH Groups. Langmuir 2013, 29, 11062−11073. (14) Hermansen, H.; Landa, G.; Sylte, J.; Thomas, L. Experiences after 10 years of waterflooding the Ekofisk Field, Norway. J. Pet. Sci. Eng. 2000, 26, 11−18. (15) Sloan, E. D., Koh, C. A. Clathrate Hydrates of Natural Gases; CRC Press: Boca Raton, FL, 2008. (16) Mohammadi, A. H.; Eslamimanesh, A.; Yazdizadeh, M.; Richon, D. Glycol Loss in a Gaseous System: Thermodynamic Assessment Test of Experimental Solubility Data. J. Chem. Eng. Data 2011, 56, 4012−4016. (17) Hassenkam, T.; Skovbjerg, L. L.; Stipp, S. L. S. Probing the intrinsically oil-wet surfaces of pores in North Sea chalk at subpore resolution. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 6071−6076. (18) Bohr, J.; Wogelius, R.; Morris, P.; Stipp, S. Thickness and structure of the water film deposited from vapour on calcite surfaces. Geochim. Cosmochim. Acta 2010, 74, 5985−5999.

ASSOCIATED CONTENT

S Supporting Information *

Force field parameters used in our work, discussion of nonbonded energy distributions, derivation of the free energy expression, and a note on the effects of reducing the calcite slab size. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.langmuir.5b02228.

■

Article

AUTHOR INFORMATION

Corresponding Author

*Phone: +47 55583310. E-mail: [email protected]. Notes

The authors declare no competing financial interest. J

DOI: 10.1021/acs.langmuir.5b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir (19) Tomson, M. B.; Kan, A. T.; Fu, G.; Al-Thubaiti, M.; Shen, D.; Shipley, H. J. Scale Formation and Prevention in the Presence of Hydrate Inhibitors. SPE J. (Soc. Pet. Eng.) 2006, 11, 248−258. (20) Kaasa, B.; Sandengen, K.; Østvold, T. Thermodynamic Predictions of Scale Potential, pH, and Gas Solubility in GlycolContaining Systems. SPE International Symposium on Oilfield Scale, Aberdeen, U.K., May 11−12, 2005; SPE International: Richardson, TX, 2005; DOI: 10.2523/95075-MS. (21) Elsarrag, E.; Magzoub, E. E. M.; Jain, S. Mass-Transfer Correlations for Dehumidification of Air by Triethylene Glycol in a Structured Packed Column. Ind. Eng. Chem. Res. 2004, 43, 7676−7681. (22) Plimpton, S. Fast Parallel Algorithmsfor Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (23) Plimpton, S.; Thompson, A.; Crozier, P. LAMMPS Molecular Dynamics Simulator, 2014. http://lammps.sandia.gov; accessed Sept 5, 2014. (24) Fiorin, G.; Klein, M. L.; Hénin, J. Using collective variables to drive molecular dynamics simulations. Mol. Phys. 2013, 111, 3345− 3362. (25) Darve, E.; Rodríguez-Gómez, D.; Pohorille, A. Adaptive biasing force method for scalar and vector free energy calculations. J. Chem. Phys. 2008, 128, 144120. (26) Hénin, J.; Fiorin, G.; Chipot, C.; Klein, M. L. Exploring Multidimensional Free Energy Landscapes Using Time-Dependent Biases on Collective Variables. J. Chem. Theory Comput. 2010, 6, 35− 47. (27) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (28) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A package for building initial configurations for molecular dynamics simulations. J. Comput. Chem. 2009, 30, 2157− 2164. (29) Fleming, S.; Rohl, A. GDIS: a visualization program for molecular and periodic systems. Z. Kristallogr. 2005, 220, 580−584. (30) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (31) Frigo, M.; Johnson, S. The Design and Implementation of FFTW3. Proc. IEEE 2005, 93, 216−231. (32) Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 1703−1710. (33) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787−1799. (34) Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic separable dual-space Gaussian pseudopotentials from H to Rn. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 3641−3662. (35) Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J. cp2k: atomistic simulations of condensed matter systems. Comput. Mol. Sci. 2014, 4, 15−25. (36) Kolafa, J. Time-Reversible Always Stable PredictorCorrectorMethod for Molecular Dynamics of Polarizable Molecules. J. Comput. Chem. 2004, 25, 335−342. (37) Krack, M. Pseudopotentials for H to Kr optimized for gradientcorrected exchange-correlation functionals. Theor. Chem. Acc. 2005, 114, 145−152. (38) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (39) Lippert, G.; Hutter, J.; Parrinello, M. A hybrid Gaussian and plane wave density functional scheme. Mol. Phys. 1997, 92, 477−488. (40) Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 1984, 52, 255−268. (41) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511−519. (42) VandeVondele, J.; Hutter, J. An efficient orbital transformation method for electronic structure calculations. J. Chem. Phys. 2003, 118, 4365−4369.

(43) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. QUICKSTEP: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Comput. Phys. Commun. 2005, 167, 103−128. (44) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, 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, 11225−11236. (45) Jorgensen, W. L.; McDonald, N. A. Development of an all-atom force field for heterocycles. Properties of liquid pyridine and diazenes. J. Mol. Struct.: THEOCHEM 1998, 424, 145−155. (46) McDonald, N. A.; Jorgensen, W. L. Development of an AllAtom Force Field for Heterocycles. Properties of Liquid Pyrrole, Furan, Diazoles, and Oxazoles. J. Phys. Chem. B 1998, 102, 8049− 8059. (47) Damm, W.; Frontera, A.; Tirado-Rives, J.; Jorgensen, W. L. OPLS All-Atom Force Field for Carbohydrates. J. Comput. Chem. 1997, 18, 1955−1970. (48) Rizzo, R. C.; Jorgensen, W. L. OPLS All-Atom Model for Amines: Resolution of the Amine Hydration Problem. J. Am. Chem. Soc. 1999, 121, 4827−4836. (49) Watkins, E. K.; Jorgensen, W. L. Perfluoroalkanes: Conformational Analysis and Liquid-State Properties from ab Initio and Monte Carlo Calculations. J. Phys. Chem. A 2001, 105, 4118−4125. (50) Price, M. L. P.; Ostrovsky, D.; Jorgensen, W. L. Gas-Phase and Liquid-State Properties of Esters, Nitriles, and Nitro Compounds with the OPLS-AA Force Field. J. Comput. Chem. 2001, 22, 1340−1352. (51) Jorgensen, W. L.; Ulmschneider, J. P.; Tirado-Rives, J. Free Energies of Hydration from a Generalized Born Model and an AllAtom Force Field. J. Phys. Chem. B 2004, 108, 16264−16270. (52) Jensen, K. P.; Jorgensen, W. L. Halide, Ammonium, and Alkali Metal Ion Parameters for Modeling Aqueous Solutions. J. Chem. Theory Comput. 2006, 2, 1499−1509. (53) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Optimized intermolecular potential functions for liquid hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638−6646. (54) Jorgensen, W. L.; Tirado-Rives, J. The OPLS Potential Functions for Proteins. Energy Minimizations for Crystals of Cyclic Peptides and Crambin. J. Am. Chem. Soc. 1988, 110, 1657−1666. (55) Toukan, K.; Rahman, A. Molecular-dynamics study of atomic motions in water. Phys. Rev. B: Condens. Matter Mater. Phys. 1985, 31, 2643−2648. (56) de Leeuw, N. H.; Parker, S. C. Molecular-dynamics simulation of MgO surfaces in liquid water using a shell-model potential for water. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 13901−13908. (57) de Leeuw, N. H.; Cooper, T. G. Surface simulation studies of the hydration of white rust Fe(OH)2, goethite α-FeO(OH) and hematite α-Fe2O3. Geochim. Cosmochim. Acta 2007, 71, 1655−1673. (58) Kvamme, B.; Kuznetsova, T.; Kivelæ, P.-H. Adsorption of water and carbon dioxide on hematite and consequences for possible hydrate formation. Phys. Chem. Chem. Phys. 2012, 14, 4410−4424. (59) Heberling, F.; Trainor, T. P.; Lützenkirchen, J.; Eng, P.; Denecke, M. A.; Bosbach, D. Structure and reactivity of the calcitewater interface. J. Colloid Interface Sci. 2011, 354, 843. (60) Ukrainczyk, M.; Greiner, M.; Elts, E.; Briesen, H. Simulating preferential sorption of tartrate on prismatic calcite surfaces. CrystEngComm 2015, 17, 149−159. (61) Markgraf, S. A.; Reeder, R. J. High-temperature structure refinements of calcite and magnesite. Am. Mineral. 1985, 70, 590−600. (62) Hwang, S.; Blanco, M.; Goddard, W. A., III Atomistic emulations of Corrotion Inhibitors Adsorbed on Calcite Surfaces I. Force field Parameters for Calcite. J. Phys. Chem. B 2001, 105, 10746− 10752. (63) Cuong, P. V.; Kvamme, B.; Kuznetsova, T.; Jensen, B. Adsorption of water and CO2 on calcite and clathrate hydrate: the effect of short-range forces and temperature. Int. J. Energy Environ. 2012, 6, 301−309. (64) Freeman, C. L.; Harding, J. H.; Cooke, D. J.; Elliott, J. A.; Lardge, J. S.; Duffy, D. M. New Forcefields for Modeling K

DOI: 10.1021/acs.langmuir.5b02228 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir Biomineralization Processes. J. Phys. Chem. C 2007, 111, 11943− 11951. (65) Raiteri, P.; Gale, J. D.; Quigley, D.; Rodger, P. M. Derivation of an Accurate Force-Field for Simulating the Growth of Calcium Carbonate from Aqueous Solution: A New Model for the CalciteWater Interface. J. Phys. Chem. C 2010, 114, 5997−6010. (66) Gale, J. D.; Raiteri, P.; van Duin, A. C. T. A reactive force field for aqueous-calcium carbonate systems. Phys. Chem. Chem. Phys. 2011, 13, 16666−16679. (67) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A., III ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396−9409. (68) Rahaman, O.; van Duin, A. C. T.; Bryantsev, V. S.; Mueller, J. E.; Solares, S. D.; Goddard, W. A., III; Doren, D. J. Development of a ReaxFF Reactive Force Field for Aqueous Chloride and Copper Chloride. J. Phys. Chem. A 2010, 114, 3556−3568. (69) Marrink, S.-J.; Berendsen, H. J. C. Simulation of Water Transport through a Lipid Membrane. J. Phys. Chem. 1994, 98, 4155− 4168. (70) Saiz, L.; Padró, J. A.; Guàrdia, E. Structure of liquid ethylene glycol: A molecular dynamics simulation study with different force fields. J. Chem. Phys. 2001, 114, 3187−3199. (71) Padró, J.; Saiz, L.; Guàrdia, E. Hydrogen bonding in liquid alcohols: a computer simulation study. J. Mol. Struct. 1997, 416, 243− 248. (72) Flaten, E. M.; Seiersten, M.; Andreassen, J.-P. Polymorphism and morphology of calcium carbonate precipitated in mixed solvents of ethylene glycol and water. J. Cryst. Growth 2009, 311, 3533−3538. (73) Lyubartsev, A. P.; Laaksonen, A. M.DynaMix - a scalable portable parallel MD simulation package for arbitrary molecular mixtures. Comput. Phys. Commun. 2000, 128, 565−589. (74) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant pressure molecular dynamics algorithms. J. Chem. Phys. 1994, 101, 4177−4189.

L

DOI: 10.1021/acs.langmuir.5b02228 Langmuir XXXX, XXX, XXX−XXX