A computational chemical study of penetration and displacement of water films near mineral surfaces
Geochemical Transactions volume 2, Article number: 35 (2001)
A series of molecular dynamics simulations have been performed on organic–water mixtures near mineral surfaces. These simulations show that, in contrast to apolar compounds, small polar organic compounds such as phenols can penetrate through thin water films to adsorb on these mineral surfaces. Furthermore, additional simulations involving demixing of an organic–water mixture near a surfactant-covered mineral surface demonstrate that even low concentrations of adsorbed polar compounds can induce major changes in mineral surface wettability, allowing sorption of apolar molecules. This strongly supports a two-stage adsorption mechanism for organic solutes, involving initial migration of small polar organic molecules to the mineral surface followed by water film displacement due to co-adsorption of the more apolar organic compounds, thus converting an initial water-wet mineral system to an organic-covered surface. This has profound implications for studies of petroleum reservoir diagenesis and wettability changes.
Adsorption of organic solutes on mineral surfaces is an important process in many natural and engineered environments and is particularly important in controlling the distribution of organic compounds in the subsurface. Molecular transformations related to the chemisorption of organic matter can remove organic contaminants from soils and sediments[1, 2] and have also been associated with diagenetic alterations in biomarker compounds. Physisorption of organic molecules on mineral surfaces controls the wettability of permeable rocks and thus capillary pressure terms in multiphase fluid flow and has been proposed as a factor controlling the composition of migrating petroleum. [5–7] Furthermore, oil wetting of mineral surfaces has been proposed to retard or stop diagenesis and cementation of oil reservoirs. In this paper we report on investigations using molecular dynamics calculations[9, 10] (MD) on a possible chain of events that could change an initially water-wet mineral surface to an oil (or more apolar)-wet surface. We hypothesise that this chain of events involves a penetration of the water film covering the mineral surface by small relatively hydrophilic polar organic molecules which physisorb to reduce mineral surface polarity, followed by an adsorption of apolar organic molecules on the thus preconditioned surface. To evaluate this two-stage adsorption mechanism a series of simulations have been performed studying the stability of water films adsorbed on mineral surfaces as a function of the polarity of an adjacent organic phase. For the organic phase a set of compounds ranging from cyclohexane (highly apolar), carbazole (apolar), phenol (polar) to acetic acid (highly polar) was used. For each organic phase two simulations were performed, one starting with a randomly distributed organic–water mixture and the other with completely demixed water and organic phases, with the water film wetting the mineral surfaces.
The inorganic component of the simulations constituted of a (1,-1,1) calcite and a (1,1,1) α-quartz surface. Both calcite and quartz are important contributors to the sediment mineral matrix and authigenic cements. To facilitate comparison of mineral surface effects on organic–water phase composition a completely hydroxylated α-quartz surface was used in the simulations, instead of replacing part of these hydroxyl groups with oxygen bridges (which is more realistic at neutral pHs). Thus, both the calcite and the α-quartz surface presented a homogeneous surface to the organic/water phases.
Finally, to assess the surfactant behaviour of small polar molecules, a set of simulations was performed in which the calcite surface was selectively and partially covered with phenol molecules. This system was covered with a water film and brought into contact with an apolar cyclohexane phase. Sorption of small polar molecules such as phenols to mineral surfaces has been shown to occur during both simulated and natural petroleum migration[11, 12] through initially water-wet rocks. By comparing the results from these simulations with those from the cyclohexane–water–mineral surface simulations the influence of pre-adsorption of small polar organic molecules on the final wettability of mineral surfaces to apolar components could be monitored.
For the organic molecules charge distributions were calculated using the Electronegativity Equalization Method (EEM). EEM uses atomic electronegativity and hardness to determine a geometry-dependent molecular charge distribution. The EEM parameters for carbon, nitrogen, oxygen and hydrogen were derived by using the fitting procedure described by Njo et al. The EEM parameters for the organic compounds were optimised to reproduce the MKS charges for proteins as reported by Cornell et al. Fig. 1 shows the reproduction of the MKS protein charge distributions using the optimised EEM parameters (Table 1). Fig. 2 shows the EEM charge distributions for carbazole, phenol and acetic acid as used in the simulations.
To allow a proper description of Coulomb interactions between charges qi and qj on neighbouring atoms in the inorganic phases a shielding-term σij was incorporated, adjusting the Coulomb interaction for electron orbital overlap as described in eqn. (1), where qi and qj are the atomic charges (in atomic units), e is the elementary charge, 4πεo is the vacuum permittivity and rij is the interatomic distance (in Å).
For the organic compounds, 1–2 and 1–3 Coulomb interactions were not taken into account. For these compounds shielding parameters of 0.2 were used, effectively removing the shielding influence from the Coulomb potential. Table 2 shows the charge distributions and shielding parameters used for the mineral phases. The α-quartz charge distribution was taken from Demiralp et al., calcite charges were initially taken from Tobias and Klein. However, since Tobias and Klein employ a Lennard-Jones potential to calculate van der Waals interactions while we employ a Morse-potential as described by Demiralp et al. these charges were modified to reproduce bulk densities of calcite and aragonite.
For the α-quartz phase the parameters derived by Demiralp et al. were used, expanded with parameters to describe the surface hydroxyl groups. For the calcite phase force field parameters were determined using the Morse functional form in eqn. (2) as described by Demiralp et al. This Morse-potential was used to describe both the inter-and intramolecular van der Waals interactions as well as the covalent Si-O, C-O and O-H bonds in the mineral phases.
For cyclohexane, a united atom approximation was used in which the CH2-groups are described as a single, charge-neutral, atom. As such, the cyclohexane molecules in our simulations are completely non-polar (more apolar than realistic cyclohexane molecules) and as such provide an extreme case in our partitioning simulations. United CH2-atom parameters were taken from DeBolt and Kollman. For the other organic compounds the AMBER-protein force field was used, alongside EEM-charges based on the MKS-method as described in the previous section. To allow the use of straightforward combination rules to describe inorganic–organic van der Waals interactions the Lennard-Jones potential, used in the AMBER-force field to describe these interactions, was replaced by the Morse-potential used in eqn. (1). Morse parameters for the organic compounds were determined by fitting them to the AMBER Lennard-Jones potential (Table 3). The water molecules were described using the transferable intermolecular potential 3 (TIP3P) of Jorgensen et al. in which the original TIP3P Lennard-Jones parameters were again used to determine suitable Morse parameters (Table 3).
The simulations were performed using the Delphi program. All simulations were performed at a simulated temperature of 298 K using a time-step of 2 fs. The system temperature was held constant using the algorithm described by Berendsen et al. using a temperature damping constant of 2000 fs. Bonds and bond angles were held stationary during the simulation using the approach described by Andersen. Apart from the surface hydroxyl groups, mineral phase atom positions were kept fixed during the simulations.
A cutoff radius of 8.5 Å was used in all simulations. The system energy was corrected for van der Waals interactions beyond this range. To ensure a stable simulation of the Coulomb interactions a 7th order taper function was applied as described by de Vos Burchart.
To allow for a direct comparison between the processes at both the calcite and α-quartz interface a periodic cell (dimensions 24.906 × 40.00 × 19.806Å) was generated containing a slab of both minerals. For α-quartz a hydroxylated slab containing 6 × 1 × 4 orthogonalised unit cells (dimensions: 25.528 × 7 × 19.652Å) was constructed, while for calcite a 4 × 1 × 4 slab was generated, with dimensions 24.284 × 9 × 19.960 Å. Cell parameters of these two mineral phases were averaged in the a- and c-directions to make both fit within the same cell dimensions, and atom coordinates were scaled accordingly. The α-quartz and calcite slabs were separated by 24 Å, thus leaving a volume of 24.906 × 24 × 19.806 Å between them that was filled with the organic–water mixtures. 50 vol.% of this interslab space was filled with organic compounds and the other 50 vol.% with water, both at densities in accordance to their respective bulk densities at standard conditions. Table 4 shows the resulting compositions for the simulated systems. For all these system compositions, both a random (mixed components) and a demixed (separate compound phases) configuration was generated. These configurations were subjected to a short, low-temperature (T = 25 K) MD-run to remove short contacts. Hereafter, MD simulations were performed at T = 298 K until the systems reached steady-state. Molecular trajectories from these steady-state systems, containing 100,000 iterations (equivalent to a time interval of 100,000 × 2 fs = 0.2 ns), were used to determine position-dependent diffusion constants and average molecular y-coordinates. Separate simulations were performed to determine the influence of simulation temperature and system size on the position-dependent diffusion constants in the phenol–water system.
Results and discussion
Demixing simulations on water-wet mineral surfaces
Fig. 3, 4, 5, 6 show the initial (random) and final configurations obtained from the cyclohexane–water, carbazole–water, phenol–water and acetic acid–water MD simulations. A clear influence of the organic phase polarity on the demixing events is observed; in case of the apolar organic phases (cyclohexane and carbazole) a distinct organic phase is formed, in the phenol case the water/organic phase separation is far less distinct whereas in the acetic acid case no clearly distinguishable demixing of acetic acid and water takes place. Furthermore, significant differences in mineral surface wetting behaviour are observed. At the end of the simulation, the cyclohexane phase (Fig. 3) is completely removed from both mineral surfaces. The final carbazole phase (Fig. 4) has a different shape than the cyclohexane phase, penetrating the water layer at both the calcite and the α-quartz surface, although the bulk of the carbazole phase is not in contact with the mineral surfaces. In the phenol and acetic acid cases both mineral surfaces are accessible for sorption of the organic compounds.
Due to the small timescales used in the MD simulations the steady-state configurations obtained from the randomly orientated systems may reflect metastable intermediate configurations rather than the thermodynamic end points. To check this, a second set of simulations was performed in which each organic–water system was initialised with completely water-wet mineral surfaces.
Fig. 7, 8, 9 show the initial and final configurations for these simulations. In the phenol and acetic acid cases end points quite similar to those obtained from the initially random molecular configurations are observed. In the carbazole case, however, we find that, unlike the results from the random configurations, the mineral surfaces remain completely water-wet. Averaged over the last 50,000 iterations (equal to a time interval of 0.1 ns) the system energy of the carbazole–water system started from the random configuration is -52821.7 ± 19.7 kJ mol-1. Over the same interval the initially water-wet carbazole system has a system energy of -52919.2 ± 7.5 kJ mol-1. This difference of over 90 kJ mol-1 in system energy indicates that breaking through the water film is not thermodynamically favoured for the carbazoles and that the end configuration obtained from the random geometry indeed reflects a metastable situation, probably caused by the energy barriers associated with the replacement of an adsorbed carbazole molecule with water.
Position-dependent diffusion coefficient analysis
The analysis presented in the previous section has primarily been qualitative, based on visual observations regarding the relative positions of the organic and water phases with respect to the mineral surfaces. One of the strong features of MD simulations is that by monitoring the development of the molecular position through time it gives the opportunity for a more quantitative analysis. To thus scrutinise the observations made in the previous section we have calculated the molecular diffusion constants for each molecule in the various organic–water–inorganic systems studies. To this end, atom positions were saved at 1000 iteration intervals (equivalent to a time interval of 1000 × 2 fs = 2 ps). From these trajectory files the molecular diffusion constant dm was calculated from the cartesian coordinates of the molecular centres of mass (xm, ym, zm) using eqn. (3).
Fig. 10 shows these molecular diffusion constants for the cyclohexane–water simulation plotted against the average molecular y-coordinate (plotted as the x-axis in Fig. 10, 11, 12) during the entire diffusion analysis time window. The average molecular y-coordinate reflects the molecular distance normal to the mineral surfaces, hence it can be used as a direct indication for the degree of interaction with these surfaces; molecules with average y-coordinates around -12 Å are in direct contact with the calcite surface, molecules roughly between -5 Å and + 5 Å have only indirect surface interactions while molecules placed around y-coordinates of +10 Å can be assumed to be adsorbed on the α-quartz surface. Fig. 10 confirms the qualitative statements made in the previous section: specific regions containing exclusively water molecules or cyclohexane molecules are observed, signifying complete phase separation, and the entire cyclohexane phase is located almost exactly between the mineral surfaces with a water film thickness of about 5–7 Å. Furthermore, Fig. 10 shows greatly reduced diffusion coefficients of the water molecules near the mineral surfaces, indicating that these first layers of the surface water film are tightly adsorbed.
Fig. 11 and 12 confirm the observations that, in contrast to the apolar cyclohexane molecules, small polar compounds such as phenol and acetic acid can migrate to the mineral surface, building up significant concentrations of adsorbed organic species. Fig. 11 also demonstrates that this observation can be made regardless of the starting conditions used in the MD simulations, as the randomly distributed phenol–water and acetic acid–water mixtures yield similar equilibrated concentration profiles compared to the water-wet starting configurations. Molecular diffusion of both phenol and acetic acid is significantly reduced near the mineral surfaces, indicating that these molecules indeed have a physical interaction with mineral surface sites.
Fig. 11 shows that raising the simulated temperature from 298 to 373 K affects the positions and diffusion constants of the phenol molecules in the equilibrated system, as an increased phenol concentration near the calcite surface is observed. Furthermore, the phenols associated to the calcite surface are separated from the other phenol molecules in the system by a water layer of almost 10 Å and phenol molecules non directly adsorbed on a mineral surface show a higher diffusion constant. Apart from these differences the results at 298 and 373 K are not dissimilar; at both temperatures the phenol molecules have access to both mineral surfaces and do not form a distinctly separated phase, distinguishing them in these aspects from the apolar cyclohexane molecules (Fig. 10).
To demonstrate that these results are not biased by the small system sizes, two more phenol–water simulations were performed in which the distance between the α-quartz and calcite slabs was gradually increased and the phenol and water molecule amounts were changed accordingly. The first of these simulations was performed with cell parameters of 24.906 × 72 × 19.806 Å, 79 phenol molecules and 460 water molecules. Fig. 13 shows the final configuration obtained from this simulation, which was commenced with a random phenol–water configuration. This figure, alongside the phenol molecular diffusion constants in Fig. 14, demonstrates that system size does not have a major affect on the simulation outcome. As with the smaller periodic cells we observe that the phenol molecules manage to penetrate the mineral surface water layer and the subsequent distributions of compounds in the equilibrated large phenol–water system are very similar to those observed in the smaller systems. The only difference of note between the larger and smaller periodic cell simulations is an expanded phenol bulk phase (between average y-coordinates of about 5 and 30 Å) in the larger periodic cell. Molecular diffusion constants in this expanded bulk phase are, on average, higher that those observed in the smaller periodic cells. This indicates that, despite not being in direct physical contact with the mineral surfaces, phenol molecules in the bulk-region of the smaller periodic cells (between -5 and + 5 Å) are still affected in their diffusive behaviour by the presence of the mineral phases. A second phenol–water simulation with an even larger system size (cell parameters 24.906 × 119 × 19.806 Å, 147 phenol and 848 water molecules) yielded phenol distributions and molecular diffusion constants very similar to those observed for the 24.906 × 72 × 19.806 Å system (Fig. 14).
Impact of surfactant phenol molecules on cyclohexane-water phase separation
In the previous two sections we have demonstrated that small polar organic compounds like phenol and acetic acid can migrate through thin water films and replace water molecules adsorbed on mineral surfaces. This is in agreement with laboratory studies showing that initially water saturated rocks quickly sorb small polar molecules from hydrocarbon–water systems. Upon adsorption on the mineral surface such organic species could potentially act as surfactants, altering the wetting behaviour of these surfaces. To determine the surfactant influence on phase behaviour two MD simulations were performed in which a random cyclohexane–water phase was demixed in the presence of an uncoated α-quartz surface and a calcite surface coated by, respectively, 6 and 1 phenol molecules. The phenol molecules were kept from migrating away from the calcite surface by introduction of a mild attractive potential between the phenol hydroxyl carbon and an adjacent surface calcium atom. Fig. 15 shows that the final shapes of the water and cyclohexane phases from these phenol-coated simulations are very different from those observed in the uncoated simulations (Fig. 3). In case of the heavily coated calcite surface (6 phenols) almost all water molecules are in the end removed from this surface, ending with a bilayer consisting of, first, the surfactant phenol layer and, thereafter, the entire cyclohexane phase. Even a single adsorbed phenol changes the surface characteristics sufficiently to alter the cyclohexane shape from a band stretching across the periodic cell at maximum distances from both mineral phases (Fig. 3) into a droplet touching the calcite surface at the location of the phenol surfactant (Fig. 15). These observations indicate the major impact of even low concentration of surfactant species on mineral surface wettability. Combined with the earlier observation that small polar organic compounds can migrate across mineral surface water films, this provides support for the hypothesis that initially water-wet mineral surfaces can undergo wettability alterations by a two-stage process, involving initially the establishment of a surfactant layer by small polar organic molecules replacing water at the mineral surface. These surfactant molecules change the surface properties, breaking surface water films and allowing surface access to more apolar compounds, including even high molecular weight components such as asphaltenes, that would have been unable to penetrate an uninterrupted water film on their own.
Comparison between calcite and α-quartz surface
In general terms, the calcite and α-quartz surfaces show comparable behaviour with respect to accessibility to organic compounds. Apolar compounds like cyclohexane do not partition to either surface (Fig. 10) while both phenol and acetic acid do (Fig. 11, 12 and 14). On closer analysis we find that the diffusion constants of acetic acid near the calcite surface are significantly lower than those near the α-quartz surface. This effect, if present, is far less distinct for phenol. This indicates that the highly polar acetic acid molecules get tightly bound to the ionic calcite surface while they associate less strongly to the more covalent α-quartz. This could signify that highly polar mineral surfaces could become oil (or organic) wet in contrast to more apolar surfaces due to their high affinity for polar organic surfactants.
The results in this study have potentially profound implications for several processes central to petroleum geology. Clearly rapid sorption of polar surfactants such as phenols can transform the wetting properties of minerals, thus aiding the development of permanently hydrophobic surfaces following direct sorption of more apolar and/or higher molecular weight compounds from apolar phases. This has implications for the wettability of petroleum reservoir rocks, and thus fluid flow, and also for the preservation of water films necessary in the distribution of silica and other inorganic solutes involved in mineral diagenesis.
A series of MD simulations were performed to study the phase behaviour of organic–water mixtures in the presence of α-quartz and calcite surfaces. Upon changing the polarity of the organic phase by going from a charge-neutral cyclohexane via apolar carbazole and polar phenol to highly polar acetic acid distinct changes in phase behaviour are observed. Cyclohexane and carbazole form discrete organic phases removed from the mineral surfaces by a water film; phenol and acetic acid show a sufficient water solubility to penetrate through these films to compete with the water molecules for mineral surface adsorption sites. Potential biases in these observations, due to simulation conditions, were tested by changing the MD-starting organic–water configurations from completely random to completely phase-separated and by changing the system size and simulation temperature, and neither of these were found to have significant effects on the accessibility of the mineral surfaces for the organic compounds.
Further MD simulations with pre-adsorbed phenol compounds on a calcite surface demonstrate that these small polar molecules can have a profound impact on surface wettability, making the mineral surface accessible to even completely non-polar compounds. This provides strong support for a two-stage process incurring wettability changes in mineral systems, commencing with small polar species migrating to the initially water-wet mineral surface, thus changing its surface characteristics, after which the rest of the water phase gets replaced by more apolar compounds adsorbing on and around these polar surfactants.
Alexander M: Environ Sci Technol. 1995, 29: 2713-10.1021/es00011a003.
Kubicki JD, Itoh MJ, Schroeter LM, Apitz SE: Environ Sci Technol. 1997, 31: 1151-10.1021/es960663+.
Sieskind O, Joly G, Albrecht P: Geochim Cosmochim Acta. 1979, 43: 1675-10.1016/0016-7037(79)90186-8.
Anderson WG: J Pet Technol. 1986, 1125-
Nagy B: Geochim Cosmochim Acta. 1960, 42: 77-95.
Krooss BM, Brothers L, Engel MH: Petroleum migration: special publication of the Geological Society. Edited by: England WA, Fleet AJ. 1991, the Geological Society, London, 59: 149-
Larter SR, Bowler BFJ, Li M, Chen M, Brincat D, Bennett B, Noke K, Donohoe P, Simmons D, Kohnen M, Allan J, Telnaes N, Horstad I: Nature. 1996, 383: 593-10.1038/383593a0.
Worden RH, Oxtoby NH, Smalley PC: Pet Geoscl. 1998, 4: 129-
Allen MP, Tildesley DJ: Computer Simulation of Liquids. 1987, Clarendon Press, Oxford
Frenkel D, Smit B: Understanding molecular simulation. 1996, Academic Press, San Diego
Taylor P, Larter SR, Jones DM, Dale J, Horstad I: Geochim Cosmochim Acta. 1997, 61: 1899-10.1016/S0016-7037(97)00034-3.
Larter SR, Bowler B, Clarke E, Wilson C, Moffat B, Bennett B, Yardley G, Carruthers D: Geochem Trans. 2000, 9:
Mortier WJ, Ghosh SK, Shankar S: J Am Chem Soc. 1986, 108: 4315-10.1021/ja00275a013.
Njo SL, Fan J, van de Graaf B: J Mol Catal A: Chem. 1998, 134: 79-10.1016/S1381-1169(98)00024-7.
Cornell WD, Cieplak P, Bayly C, Gould IR, Merz KM, Ferguson DM, Spellmeyer DC, Fox T, Caldwell JW, Kollman PA: J Am Chem Soc. 1995, 117: 5179-10.1021/ja00124a002.
Demiralp E, Çagin T, Goddard WA: Phys Rev Lett. 1999, 82: 1708-10.1103/PhysRevLett.82.1708.
Tobias DJ, Klein ML: J Phys Chem. 1996, 100: 6637-10.1021/jp951260j.
DeBolt SE, Kollman PA: J Am Chem Soc. 1995, 117: 5316-5340. 10.1021/ja00124a015.
Jorgensen WL, Chandrasekhar J, Madura JD: J Chem Phys. 1983, 79: 926-935. 10.1063/1.445869.
van de Graaf B, Baas JMA: J Comput Chem. 1984, 5: 314-10.1002/jcc.540050406.
Berendsen HJC, Postma JPM, van Gunsteren WF: J Chem Phys. 1984, 81: 3684-10.1063/1.448118.
Andersen HC: J Comput Phys. 1983, 52: 24-10.1016/0021-9991(83)90014-1.
de Vos Burchart E: 1992, PhD Thesis, Delft University of Technology, the Netherlands
Bjørlykke K, Ramm M, Saigal GC: Geol Rundsch. 1989, 78: 243-268. 10.1007/BF01988363.
This work was supported by a Royal Society Research Fellowship and by TMR grant No. ERBFMBICT971871 for ACTvD. We thank the reviewers for their helpful suggestions and remarks.
About this article
Cite this article
van Duin, A.C., Larter, S.R. A computational chemical study of penetration and displacement of water films near mineral surfaces. Geochem Trans 2, 35 (2001). https://doi.org/10.1186/1467-4866-2-35