Fluid Interface

Article pubs.acs.org/JPCC

Insights into Microscopic Diffusion Processes at a Solid/Fluid Interface under Supercritical Conditions: A Study of the Aqueous Calcite (101̅4) Surface Chun-Yaung Lu,*,† Danny Perez,† Donald D. Hickmott,‡ and Arthur F. Voter† †

Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, United States Earth and Environmental Sciences Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, United States

‡

ABSTRACT: The diffusion of adsorbed molecules at a solid/ fluid interface is a phenomenon of great fundamental interest and technological importance. We use molecular dynamics and kinetic Monte Carlo simulations to investigate the diffusion processes at an aqueous calcite (101̅4) interface under various geologically relevant supercritical conditions [101.3 MPa (1000 atm), 300−800 K]. Between 600 and 700 K, the adsorption configuration changes from an outer- to inner-sphere surface complex, and the desorption free energy barrier, ΔGd, increases greatly. The ratio ΔGd/kBT takes a minimum value at 600 K, making it the temperature with the greatest desorption tendency. Consequently, the temperature dependence of the mean lateral diffusivity of a finite bulk system reverses at 600 K. More interestingly, the system exhibits typical characteristics of bulk fluidmediated surface diffusion. For example, the mean-squared displacement shows different scaling properties with changes in temperature or bulk thickness. The optimized bulk thickness which maximizes the effective surface diffusivity was also found to vary with temperature.

I. INTRODUCTION Phenomena at mineral interfaces in supercritical fluids are fundamental in a wide range of Earth systems. They occur during the formation of ore deposits,1 within ocean-floor hydrothermal systems, during metamorphism within the continental lithosphere,2 during geologic carbon sequestration,3,4 and during geothermal energy production.5 Kinetic processes and mechanisms in such conditions are ill-constrained because experimental studies of solid/liquid surface dynamics at relevant elevated pressures (P) are difficult. In particular, geothermal energy production, including in enhanced geothermal systems,5 can occur at temperatures (T) and Ps above the critical point of water. Recent developments in extraction techniques aim to improve the economics of geothermal energy.6 In the case of the Iceland Deep Drilling Project at the Krafla caldera in Iceland, it has been speculated that supercritical conditions (Tc = 647.1 K, Pc = 22.064 MPa) for pure water would be reached at about 3.5 km depth and could be deeper for geothermal fluid.7 Given the fact that the Earth’s P gradient is roughly 27 MPa/km in the continental crust,8 the pressure can reach 100 MPa at 4 km depth, and the temperature at a similar depth can be as high as 800 K. Understanding the reactions occurring at mineral interfaces in geothermal fluids under such high P−T conditions is thus important for geothermal energy production. It is also of great fundamental interest to understand how extreme P−T environments affect interfacial phenomena. Since high P−T measurements require special apparatus, relevant experiments are challenging and limited by technical difficulties. Considering the feasibility and the cost in time and money to conduct © 2012 American Chemical Society

experiments, computer simulations are a promising tool to study systems under extreme conditions.4,9 It is known that adsorbed particles (ions or molecules) on a solid surface can undergo two types of surface translational motion. One is in-surface self-diffusion, and the other is bulkmediated surface diffusion (BMSD, bulk here referring to the fluid).10−12 Compared to bulk diffusion, in-surface diffusion is usually much slower due to the rough energy landscape; i.e., it takes time for the adsorbed particle to overcome the energy barriers hampering lateral mobility. Because the energy landscape is determined by interactions with the surface, it is not isotropic in space and also changes with time. BMSD occurs when fluid is present above the solid surface. The adsorbed particles can then desorb from the surface and readsorb elsewhere; i.e., the excursions into the bulk fluid lead to effective surface displacements. It has been shown that where the characteristic readsorption time is much less than the retention time (the time for the adsorbent to become permanently lost to the bulk), the effective surface diffusion becomes non-Fickian and has anomalous scaling properties. More precisely, the particles perform “Lévy walks” on the surface.10−12 In general, reactant molecules need to reach an active site to initiate a reaction. Therefore, the kinetics of the reaction are influenced by the transport properties of the reactants. For mineral systems, the active site can be a kink site or a vacancy on the surface, and dissolved mineral ions or Received: October 19, 2012 Revised: November 12, 2012 Published: November 29, 2012 25934

dx.doi.org/10.1021/jp310391d | J. Phys. Chem. C 2012, 116, 25934−25942

The Journal of Physical Chemistry C

Article

in the literature.21,25−30 We adopted the force field recently developed by Raiteri et al.,29 which has been shown to provide accurate prediction of the solvation free energy of the component ions under ambient conditions. While the force field used in this work was not specifically designed to model water−calcite interaction at high T, we believe that the phenomena we discuss below are generic and should be applicable to a wide range of mineral/fluid systems. The organization of this paper is as follows. In section II we introduce the methodologies and tools used in this work. In section III.A we discuss the results of desorption free energy calculations and the system structure under different P−T conditions. In section III.B1 we discuss additional observations from the simulations and the details of the KMC modeling. In section III.B2 we present the results of mean diffusivities calculated for finite bulk systems. In section III.B3 we discuss the BMSD properties of the system. Finally, in section IV we offer some conclusions.

impurity species may combine with the vacancy, causing surface reconstruction. Since bulk fluid diffusion is much faster than insurface diffusion, the process can be enhanced when bulk diffusion is involved. In high temperature regimes, both the insurface and the BMSD mechanisms become faster and compete with each other. As will be shown below, the overall diffusion behavior is sensitive to the surface and fluid properties, which are T and P dependent. Theoretical understanding of the physics of diffusion can be complex due to the interactions between the adsorbent particle, surface, and solvent molecules. Large-scale fully atomistic modeling, while avoiding the use of many simplifying assumptions, is usually computationally expensive, making long-time simulations impracticable. One common strategy to extend the simulation time scale is to adopt a continuum approach, which can effectively reduce the computational costs by ignoring some microscopic details while still considering the important physics. Continuum methods are typically based on the solution of partial differential equations such as the Fokker−Planck equation.13 In the overdamped limit, the time evolution of the probability distribution function for diffusive dynamics can be described by the Smoluchowski equation. The first passage time distribution for given initial and final states (for example, from the adsorbed state to the desorbed state) can also be generated.13−16 For a metastable multistate system, the transitions between stable configurations usually occur on a long time scale compared to thermal vibrations. In this case, it can be more convenient to describe the time evolution of the system using discrete state-to-state dynamics. One example of a coarsegrained stochastic simulation method is the well-known kinetic Monte Carlo (KMC) algorithm.17−20 In each KMC step, a new state and the duration (passage) time in the old state are randomly drawn from appropriate distributions. For a surface desorption process, several milestone states can be defined along the reaction coordinate (e.g., the vertical distance between the adsorbent and the surface). The transition time from one state to another satisfies a first passage time distribution which can be obtained from solutions of the Smoluchowski equation. In the present work, the aqueous calcite (101̅4) surface was chosen as a model system for studying diffusion at a solid/fluid interface under supercritical aqueous conditions. Calcite (CaCO3) is an abundant mineral species and is a primary carbon reservoir in the continental crust. Its surface interactions with fluids are important in oil and gas reservoirs, in geologic carbon sequestration environments, and in geothermal systems. The (101̅4) surface largely dominates the calcite morphology and therefore is important for crystal growth and dissolution.21 In aqueous environments, CO32− is in equilibrium with HCO3− and H2CO3. The relative concentrations depend on the pH value, T, and P. More accurate modeling, involving chemical reactions such as the protonation and deprotonation of carbonate/water species and the effect of the ionic strength of the solution, can be obtained through first principles methods;22−24 however, this is only achieved through an enormous increase in simulation complexity and computational costs. Since our main purpose is not to very accurately reproduce the behavior of a specific system but to gain general insights into diffusion at mineral/water interfaces in supercritical conditions, we chose to rely on classical atomistic simulations with an empirical potential. Many empirical potentials for the aqueous calcite system have been proposed

II. THEORETICAL METHODS A. Molecular Dynamics Simulations. In this work, molecular dynamics (MD) simulations were performed in the canonical (constant-NVT) ensemble using a Nosé−Hoover chain31 with a relaxation time of 0.1 ps. The velocity-Verlet algorithm32 with a 2 fs time step was used to integrate the equations of motion. We verified that the simulation results are not sensitive to the settings of relaxation time (in the range 0.1−1.0 ps) and time step size (1−2 fs). We used the force field developed by Raiteri et al.29 to simulate the interactions in the aqueous calcite system. The force field contains two-body, three-body, four-body, and long-range electrostatic interactions. Both water and carbonate structures are considered to be flexible. The SPC/Fw model was used for water, which has two harmonic O−H bonds and one harmonic H−O−H angle.33 Carbonate is described by three harmonic C−O bonds, three harmonic O−C−O angles, and one out-of-plane anharmonic C−O3 bond. All MD simulations were performed with the LAMMPS package (version: March 2012).34 The long-range Ewald sums for Coulombic interactions were computed using the “Ewald/n” style with the precision set to 10−6. B. System Details and Geometry. The simulation cell with periodic boundary conditions consisted of 516 SPC/Fw water molecules above a triclinic 5 × 6 × 7 (CaCO3 units) substrate originally created from the (101̅4) cleavage of calcite. The normal to the (1014̅ ) surface was aligned with the z-axis. The carbon atoms in the bottom two layers were held frozen at their equilibrium bulk position while the top surface was exposed to water. The critical temperature and pressure of the SPC/Fw water model are35 Tc = 634.5 K (expt 647.1 K) and Pc = 17.25 MPa (expt 22.064 MPa), respectively. The optimized lattice constants of bulk calcite for given target T and P were determined using a gradient-based free energy minimization method described in ref 36. In this framework, the free energy derivative with respect to the element of the edge vector matrix h is −

⎛ δij ⎞ ∂G = V {[(p − P)h−1]ij + [(p − P)h−1]ji }⎜1 − ⎟ 2⎠ ∂hij ⎝ (1)

where G is the free energy, P is the target pressure, p is the internal pressure tensor, which can be estimated from a constant-h MD simulation (all atoms are allowed to move), and 25935

dx.doi.org/10.1021/jp310391d | J. Phys. Chem. C 2012, 116, 25934−25942

The Journal of Physical Chemistry C

Article

δij is a Kronecker delta function. The system (with water) was again optimized solely in the z-direction using the same method until the desired hydrostatic pressure was reached. Whenever a new CaCO3 unit was introduced to the system, the cell geometry was reoptimized in the z-direction. Once the cell was optimized, it was held constant during the production simulations. C. Desorption Free Energy. Let the perpendicular distance to the surface, z, be the reaction coordinate. The free energy difference between two planes, za and zb, is ⎛Q ⎞ ΔG = G(zb) − G(za) = −kBT ln⎜⎜ b ⎟⎟ ⎝Qa ⎠

trajectories. The diffusion processes in the perpendicular (to the surface) and lateral directions are treated independently. The perpendicular diffusion is described by eq 3 while the lateral diffusion is considered as free 2-D Brownian motion with an effective diffusivity appropriate for that state. To illustrate our approach, consider the model system shown in Figure 1, which consists of an energy well and a flat free

(2)

where G denotes the free energy, kB is the Boltzmann constant, T is the absolute temperature, and Qa and Qb are the corresponding partition functions at za and zb. Since the partition function is proportional to the population, the ratio Qa/Qb is equal to the population ratio ρa/ρb, which can be estimated from a MD trajectory equilibrated in the region of interest. If ΔG is high compared to the thermal energy kBT, sampling the free energy surface with direct MD simulations becomes inefficient and more advanced techniques (such as umbrella sampling37) would be needed. In this work, the desorption barriers are generally less than 6kBT, and thus sampling with direct MD simulations is adequate. D. Numerical Solution of the Smoluchowski Equation. The desorption dynamics can be characterized by a onedimensional diffusion process on a free energy landscape, G(z), where z denotes the reaction coordinate, which is the perpendicular distance of the adsorbent measured from the interface. The time evolution of the probability density for the position of the adsorbent, 7(z , t ), can be described by the Smoluchowski equation:38 ∂7(z , t ) ∂ ∂ = D e−βG(z) [e βG(z)7(z , t )] ∂z ∂z ∂t

Figure 1. One-dimensional model system of the adsorption/ desorption processes at T = 500 K with diffusivity D = 0.596 Å2 ps−1. (a) The potential energy (red): V(z) = (z4 − 2z2 + 1)/5 for z < 0 and V(z) = 1/5 for z ≥ 0. z1, z2, and z3 are the coordinates of the energy minimum (z = −1), the well boundary (z = 0), and a point (z = 1) in the flat-energy region, respectively. (b−d) Passage time distributions for the transitions defined in (a) (also see the details in the text).

energy region representing the adsorbed and desorbed state, respectively. The passage time, τ, is defined as the elapsed time until the system, starting from a given state, reaches a neighboring state for the first time. The distribution of τ in the well can be estimated by solving eq 3 with initial condition 7(z , 0) = δ(z − z1). The probability density of the passage time from z1 to z2, Pr1→2(τ), is obtained from the (normalized) flux J+(z2,τ) passing through the absorbing boundary at z2. Similarly, by setting the initial condition 7(z , 0) = δ(z − z 2), the distributions Pr2→1(τ) and Pr2→3(τ) can be obtained from the flux J−(z1,τ) and J+(z3,τ), respectively. The corresponding branching ratios B2→1 and B2→3 are defined as the probability of reaching either state 1 or 3 first, starting from 2. The last piece of the kinetic framework is the distribution Pr3→2(τ), which is mathematically equivalent to the first passage time distribution of a 1-D Brownian motion that starts at z3 and ends at z2. The solution39 is

(3)

where β−1 = kBT and D is a coordinate-independent diffusion coefficient. The Smoluchowski equation can be solved numerically by discretizing the coordinate space (with mesh size d = zn+1 − zn) and recasting eq 3 into a master equation:38 ∂7(zn , t ) = k(n|n + 1)7(zn + 1 , t ) ∂t + k(n|n − 1)7(zn − 1 , t ) − [k(n + 1|n) + k(n − 1|n)]7(zn , t )

(4)

where k(m = n ± 1|n) denotes the rate of probability flow from grid point n to m: k(m|n) =

⎧ −β[G(zm) − G(zn)] ⎫ D ⎬ exp⎨ 2 ⎩ ⎭ 2 d

Pr

(5)

The quantity k(m|n)7(zn , t ) ≡ J ±(zn , t ) is the probability flux from grid point n to m at time t, where the superscript ± indicates whether the flux is forward (+, zm > zn) or backward (−, zm < zn). If zn is an absorbing boundary, k(m|n), the rate associated with the flux leaving zn, is set to zero. E. Desorption Dynamics. A complex surface desorption process can have several intermediate steps. We utilized KMC with a set of predefined states and their corresponding firstpassage time distributions to generate long state-to-state

1 ∂ (τ ) = 2 ∂τ

3→2

∞

∑

( −1)n fn (τ )

n =−∞

(6)

where fn (τ ) = −erf[α(2nL − z′)] + erf{α[(2n + 1)L − z′]} + erf[α(2nL + z′)] − erf{α[(2n + 1)L + z′]} (7)

erf(x) = 2π−1/2∫ x0dx exp(−x2) represents the error function, α = (4Dτ)−1/2, z′ = z3 − z2, and L = zmax − z2, and zmax is a reflecting boundary representing the upper bound of the 25936

dx.doi.org/10.1021/jp310391d | J. Phys. Chem. C 2012, 116, 25934−25942

The Journal of Physical Chemistry C

Article

desorbed state (i.e., the surface of the fluid). In the limit zmax → ∞, eq 6 reduces to Pr 3 → 2(τ ) =

⎛ z′ 2 ⎞ z′ exp⎜ − ⎟ τ 4πDτ ⎝ 4Dτ ⎠

simulations, no evidence of dissociation was found for T = 300−800 K on time scales of tens of nanoseconds. Thus, in this work CaCO3 was considered as an inseparable and nonreactive unit. As shown in Figure 2a, the free energy profile changes only slightly from 300 to 600 K. However, for T > 600 K, a dramatic increase of the desorption barrier is observed. A closer inspection reveals important features of the interface structure. The mass densities of water oxygen as a function of vertical distance at various temperatures are summarized in Figure 2b. The three density peaks at 2.1, 3.4, and 4.8 Å correspond to the first three layers of surface water. Similar estimations of the peak positions by either experiments or simulations have been reported previously.27,30 The three free energy minima at 2.5, 3.4, and 5.0 Å roughly coincide with the water density peaks, indicating that the stable adsorbed CaCO3 configurations are associated with surface water molecules. As T increases from 300 to 800 K, the bulk water density (measured at >10 Å above the calcite surface) decreases from 1.030 to 0.477 g/cm3 (1.037 and 0.482 g/cm3, respectively, for bulk water densities computed by using the IAPWS formula with the same P and T conditions41). The relative height of the second and the third density peaks decrease with T, indicating a decrease in ordering of the water at the interface. The decrease in density at high T also reduces water’s solvation power. As a result, dissolved CaCO3 becomes less favorable and the equilibrium position corresponding to the free energy minimum moves closer to the surface as T increases. The free energy minimum is originally at ≈5 Å above the surface for T ≤ 600 K but gradually shifts to ≈3.4 Å at 800 K. The two different adsorbing configurations of CaCO3 at T = 500 and 800 K are illustrated in Figures 3a and 3b, respectively. In

(8) 1→2

2→1

With the passage time distributions [Pr (τ), Pr (τ), Pr2→3(τ), Pr3→2(τ)] and the branching ratios (B2→1, B2→3), we are able to generate a state-to-state trajectory by randomly drawing a new state and a first-passage time from the corresponding distributions. In practice, we numerically integrate Pr(τ) and use inverse transform sampling40 to draw an appropriate time. The kinetic framework we just described can be extended to any number of intermediate states, as required. Note that our approach differs from standard KMC as the state-to-state dynamics are not assumed to be Markovian. Therefore, the choice of the next visited state and the residence time are correlated.

III. RESULTS AND DISCUSSION A. Desorption Free Energy. The free energy profiles of a CaCO3 molecule desorbed from the calcite (101̅4) surface were computed for T = 300−800 K at 101.3 MPa. The results are summarized in Figure 2a. The reaction coordinate is the vertical distance between the adsorbent carbon atom and the calcite surface (defined as the average z position of the top-layer Ca atoms). Although CaCO3 can dissociate into Ca2+ and CO32−, previous studies indicated that the ion pair (CaCO3) is stable in the aqueous environment at room temperature.30 In our

Figure 3. Snapshots (cropped) of MD simulations at 101.3 MPa with temperatures (a) 500 and (b) 800 K. Atoms are colored by the element type: calcium (yellow), carbon (black), carbonate oxygen (red), and water oxygen (cyan). The carbonate oxygens in the calcite slab and the water hydrogens are not shown for clarity.

Figure 3a, water molecules are present between the bound CaCO3 and the surface. This type of configuration is usually referred to as the “outer-sphere surface complex” in contrast with the “inner-sphere surface complex” where the CaCO3 binds directly to the surface as shown in Figure 3b.42 One important consequence of the changes in the free energy profile is to modify the desorption rate dependence on temperature. In Figure 2c, we compare the ratio of the desorption barrier (ΔGd) to thermal energy (kBT) for T = 300−800 K. The ratio decreases from 300 to 600 K, indicating that the desorption barrier grows slower than the increment of thermal energy. In this region, the tendency for desorption therefore increases with T. However, the trend reverses at 600 K, reflecting the sudden increase of the desorption barrier. As discussed above, this is related to changes in the solvent state.

Figure 2. (a) Free energy curves for a CaCO3 molecule desorbed from the aqueous calcite (101̅4) surface at various T and 101.3 MPa. The curves of 400 and 500 K are not shown in the plot for visual clarity (similar to 300 and 600 K). The desorption barrier, ΔGd, is defined as the energy difference between the minimum and the plateau level. The three local minima are numbered and marked by short arrows. (b) Mass density of water oxygen as a function of vertical distance at various T and 101.3 MPa. All density values were transformed to the corresponding water density [multiplied by (18.01/16.0)]. The positions of the first three peaks are numbered and marked by short arrows. (c) The ratio ΔGd/kBT as a function of T. 25937

dx.doi.org/10.1021/jp310391d | J. Phys. Chem. C 2012, 116, 25934−25942

The Journal of Physical Chemistry C

Article

definitions). There is a notable transition between 600 and 700 K as expected from the changes in the free energy profile, which predicts the dominance of the inner-type adsorption above 600 K. Since at 700 and 800 K the adsorbent is directly bound to the surface, the motion of the CaCO3 molecule is strongly affected by the surface structure and the pattern of the underlying lattice can be easily seen in the figures. A closer inspection indicates that the CaCO3 molecule in the inner layer diffuses preferentially along the (42̅ 6̅ 1) direction (cf. Figure 4, green). The molecule can also makes excursions through the outer layer (Figure 4, red), resulting in more isotropic displacements. At lower temperatures, the molecule is weakly bound and the surface diffusion is consequently more isotropic. In order to study the diffusion processes over longer time scales, we adopted a continuum description of the system and modeled the kinetics using the KMC algorithm previously described. The kinetic framework for the aqueous calcite system is illustrated in Figure 5, where z1 − z5 represents the z positions of the predefined states used in the modeling (most of them are either local maxima or minima). The space above the calcite surface is partitioned into three sections: inner layer, outer layer, and bulk fluid. The locations of the inner/outer and outer/bulk boundaries were determined from the free energy profiles as illustrated in Figure 5. The passage time distributions, Pri→j(τ), were obtained by solving the Smoluchowski equation (eq 3) or by using the Green’s function (eq 6 or eq 8). The lateral diffusion was assumed to be independent of the diffusion along the z direction. Since a CaCO3 molecule in the inner or outer layer interacts with the surface, its lateral motion is no longer isotropic. The diagonal elements of the diagonalized diffusivity tensor calculated for the inner and outer layers using molecular dynamics are summarized in Figure 6a,b.

This effect also plays an important role in the diffusion process at the interface which will be discussed in the following section. B. Diffusion of a CaCO3 Molecule at the Interface. 1. Observations and Modeling Details. As discussed above, at 101.3 MPa, the free energy profile of a CaCO3 molecule desorbed from the calcite (101̅4) surface varies significantly with T. To understand how this change affects the dynamics at the interface, in Figure 4 we compare trajectories generated by

Figure 4. Trajectories of a carbonate carbon atom diffusing in a confined volume above the calcite/water interface at various temperatures and 101.3 MPa. The (101̅4) direction is perpendicular to the plane, and two in-plane directions, (011̅0) and (4̅2̅61), are shown in the plot. In the simulation, trajectories were reflected back toward the surface by a wall located at the local free energy maximum around 6.3 Å (≈ 6.8 Å for 700 and 800 K; see Figure 2a). The total simulation time is 10 ns with the interval between two consecutive points 2 ps. The boundary of the inner (green) and outer state (red) is the maximum between minimum 2 and 3 shown in Figures 2a and 5.

the MD simulations at different Ts. During the simulations, all trajectories were constrained within the adsorbed state to prevent desorption from the surface. The area visited by a trajectory is roughly proportional to the mean diffusivity. The trajectories are color-coded according to whether they are in the inner (green) or outer (red) layer (also see Figure 5 for

Figure 6. Diffusivity calculated for T = 300−800 K and 101.3 MPa. (a) and (b) summarize the diagonalized diffusion coefficients for the inner and outer layer in two orthogonal directions x′ (x″) and y′ (y″), respectively. The prime notations are used to differentiate the new (diagonalized) directions from the original simulation coordinates. (c) 1-D diffusion coefficient determined from the bulk water under the same temperature and pressure conditions.

The 1-D diffusivities of the bulk water under the same conditions are shown in Figure 6c. The inner-layer diffusivity is roughly 5−11 times smaller than the outer-layer one and 25− 80 times smaller than the bulk diffusivity at the same T and P. Once a passage time τ is drawn from the distribution, a lateral displacement, r (in the x and y directions), is generated using the appropriate 1-D Green’s function for free diffusion:

Figure 5. Illustration of the kinetic framework used to describe a CaCO3 molecule desorbed from the aqueous calcite (101̅4) surface. The solid black curve is the desorption free energy profile for T = 300 K and P = 101.3 MPa (same as Figure 2a). z1−z5 denote the coordinates of the predefined states. The regions z < z2 and z2 ≤ z < z4 are defined as the inner (green) and outer layer (red), respectively. The fluid above z4 is treated as bulk fluid (blue). τi→j denotes the passage time from state i to j. 25938

dx.doi.org/10.1021/jp310391d | J. Phys. Chem. C 2012, 116, 25934−25942

The Journal of Physical Chemistry C .(r , t ) =

⎛ r2 ⎞ 1 exp⎜ − ⎟ 4πD*t ⎝ 4D*t ⎠

Article

atoms) placed in the water layer above the calcite surface. The agreement is very good, which demonstrates that our model captures the essential ingredients controlling the diffusion process. 3. Bulk Mediated Surface Diffusion. Another way to describe the diffusion at a solid/fluid interface is to focus on the motions in the adsorbed state and consider the desorption/ adsorption processes as excursions into the bulk fluid which also generate surface displacements. The effective surface diffusion will scale differently from the regular Fickian diffusion law due to the bulk excursions. The studies of so-called bulk mediated surface diffusion are concerned with the kinetics of desorption/adsorption processes and their effect on the diffusion on the surface. As discussed previously, the desorption tendency is mainly controlled by the free energy barrier, but other factors, such as contact fluid thickness and in-surface, bulk diffusivities also affect the desorption/adsorption kinetics. For the aqueous calcite system we are interested in, the desorption free energy barrier exhibits an unusual trend (see Figure 2) which makes it an excellent model system to study various diffusive regimes of BMSD. The qth moment of the effective surface displacement r is defined as

(9)

where . denotes the Green’s function and D* is the diffusivity in the specified direction. Given the time t = τ, the displacement r can be drawn from the distribution .(r , τ ). 2. Mean Lateral Diffusivities of Finite Bulk Systems. The state-to-state trajectories generated from the KMC algorithm were used to study the diffusion of a CaCO3 molecule at the interface. All trajectories were initialized in the adsorbed state at t = 0. The mean lateral diffusivities, D̅ , were estimated from the slope of the mean-squared displacement (MSD) curves using the relation MSD = ⟨rx2 + ry2⟩ = 2D̅ t, where r denotes the displacements in the specified direction. The results for T = 300 K, P = 101.3 MPa with a 12 Å contact water thickness are summarized in Figure 7a. The MSD curve (red) is compared

⟨r q(t )⟩ =

1 N0

∞

∫−∞ r qn(r , t ) dr

(10)

where N0 is the total number of adsorbents initially on the surface and n(r,t) is the number of the adsorbents on the surface with surface displacement r at time t given n(0,0) = N0. ⟨r2(t)⟩ is usually referred to as the effective mean-squared displacement along the surface (the same abbreviation, MSD, is used below). The quantity

Figure 7. (a) MSD of a CaCO3 molecule diffusing in the water layer above the calcite (101̅4) surface at 300 K and 101.3 MPa (red solid surve). Two black dashed lines represent the same quantities for bulk and in-surface diffusion. The blue dotted line (with slope 0.16 Å2 ps−1) was linearly extrapolated from the long-time data of the red curve. The green line (slope = 0.17) was generated based on the bulk and surface diffusion coefficients weighted by the corresponding probabilities. (b) The values of diffusivities were obtained from the mean-square displacement plots [as in (a)] for various temperatures and water thicknesses. The black dots with error bars were estimated from direct MD simulations with contact water thickness 12 Å.

Γ(t ) =

1 N0

∞

∫−∞ n(r , t ) dr

(11)

measures the survival probability in the adsorbed state at time t. The survival probabilities calculated for a CaCO3 molecule diffusing on the aqueous calcite (101̅4) surface are summarized in Figures 8a and 8b for the systems of infinite and 19 Å contact water thickness, respectively. In the case of infinite contact water, Γ(t) eventually decays to zero due to entropic effects. If the water thickness is finite, the probability instead approaches a nonzero equilibrium value, and the order of survival probabilities reflects the tendency of being in the adsorbed state.

with the predicted MSD lines generated using the bulk and insurface (weighted sum of the inner- and outer-layer contributions) diffusivities. Since the system is bounded in the z direction, the bulk to surface probability ratio estimated from a dynamical trajectory will converge to its equilibrium value at long time. The green line was generated by interpolating the bulk and surface lines using their equilibrium probabilities [estimated from ΔGd(z)]. The slope of the red curve (≈0.16 Å2 ps−1) estimated at long time agrees with the slope of the green line (0.17), indicating that the long-time limit of the mean lateral diffusivity is also an equilibrium quantity. The values of D̅ estimated from the MSD slopes at long time for various temperatures and contact water thicknesses are summarized in Figure 7b. As seen in the figure, D̅ increases with contact water thickness due to increased contributions from the bulk diffusion. D̅ also increases with temperature and reaches its maximum value at 600 K. When T > 600 K, the system switches to inner-sphere type adsorption and D̅ therefore decreases. In order to validate the continuum model, we compare with direct MD simulations for a thickness of 12 Å (cf. black dots in Figure 7b). In the MD simulations, the water thickness is controlled by a reflecting boundary (which affects only carbonate carbon

Figure 8. Survival probability in the adsorbed state as a function of time at various temperatures and 101.3 MPa with (a) infinite and (b) 19.0 Å contact water thickness. 25939

dx.doi.org/10.1021/jp310391d | J. Phys. Chem. C 2012, 116, 25934−25942

The Journal of Physical Chemistry C

Article

It is known that BMSD exhibits non-Fickian diffusion behaviors with anomalous scaling properties. At finite times, the MSD curves can be fitted by a power law, f(t) = Atξ, where the exponent, ξ, varies with the fitting range.43 The case ξ = 1 is referred to as normal (Fickian) diffusion, while the cases ξ < 1 and ξ > 1 are usually called sub- and superdiffusion, respectively. For t → ∞, the survival probability for an infinite bulk system decays as Γ(t) ∼ t−1/2, which leads to a subdiffusive surface diffusion with MSD ∼ t1/2. The long-time ξ fitted from the MSD and Γ(t) are summarized in Figure 9a. The results

Figure 10. Short-time effective MSD (red solid lines) as a function of time for various temperatures. The calcite (101̅4) surface is in contact with infinite water. Black solid lines represent the same quantities for in-surface diffusion generated using the corresponding surface diffusivities.

Figure 9. Power exponents, ξ, obtained from the least-squres fits for MSD and Γ(t) in the (a) long-time (90−100 ns) and (b) short-time (0−0.5 ns) regions at various temperatures. The fitting function has the form f(t) = c + Atξ, where c was held at zero in the short-time fits.

agree with the predicted limits, MSD ∼ t1/2 and Γ(t) ∼ t−1/2. The 300 and 800 K exponents for Γ(t) are slightly higher than −0.5, indicating the curves have not yet reached the long-time limit. The deviations are consistent with the slow decay observed for the curves at the same temperatures in Figure 8a. Unlike the long-time limits of ξ for MSD and Γ(t), which are independent of the desorption rates, the short-time scaling reveals more detailed information about the desorption kinetics. It has been shown that for strongly adsorbing systems, on time scales much shorter than the retention time, th, the qth moment of surface displacement ⟨rq(t)⟩ scales as ∼t(q+1)/2, ∼tq, and ∼t ln(th/t) for the cases of q > 1, q < 1, and q = 1 respectively.10,12 For an aqueous calcite system with infinite contact water, the short-time (∼1 ns) effective MSD as a function of time at various Ts are summarized in Figure 10. The exponents ξ fitted from the MSD curves in the range 0−0.5 ns are shown in Figure 9b. It is clear from the ξ values and by comparing with the in-surface diffusion lines (black) that the MSD curves at 300, 700, and 800 K exhibit a superdiffusive trend (ξ > 1) while others are subdiffusive (ξ < 1). On a much shorter time scale, say 1 (data not shown). However, on such short time scales for the present (three-state) system, the scaling is more complicated than what was predicted from a simple two-state surface/bulk model. From Figure 10 we know the desorption/adsorption kinetics affects the scaling properties of effective surface diffusion. It has been reported that the diffusion trend can also be modified by changing the contact fluid thickness.44 We investigate this effect by varying the contact water thickness in the calcite system. The MSD curves calculated from the systems of different bulk thicknesses and temperatures are summarized in Figure 11. For a given temperature, 500 K for example, the effective surface diffusion can be faster or slower than the in-surface diffusion depending on the water thickness. It is interesting that the

Figure 11. Effective MSD as a function of time for various temperatures and contact water thicknesses. The order of MSD curves and the corresponding contact water thicknesses are summarized in the legends. MSD lines estimated from surface diffusivities are also included in the plots (black solid lines).

optimized water thickness that maximizes the effective surface diffusivity also varies with T. For T = 400, 500, and 600 K, the adsorbent tends to be lost to the bulk, and therefore the effective surface diffusivity is maximized with thinner contact water. In contrast, for T = 300, 700, and 800 K, the adsorbent tends to stick on the surface. Thus, thicker contact water can provide more space for bulk excursions to take place and makes the surface mass transport more efficient. 25940

dx.doi.org/10.1021/jp310391d | J. Phys. Chem. C 2012, 116, 25934−25942

The Journal of Physical Chemistry C

Article

(3) Cole, D. R.; Chialvo, A. A.; Rother, G.; Vlcek, L.; Cummings, P. T. Philos. Mag. 2010, 90, 2339−2363. (4) Kerisit, S.; Weare, J. H.; Felmy, A. R. Geochim. Cosmochim. Acta 2012, 84, 137−151. (5) Tester, J. W.; Anderson, B. J.; Batchelor, A. S.; Blackwell, D. D.; DiPippo, R.; Drake, E. M.; Garnish, J.; Livesay, B.; Moore, M. C.; Nichols, K.; Petty, S.; Toksöz, M. N.; Veatch, R. W. The Future of Geothermal Energy: Impact of Enhanced Geothermal Systems (EGS) on the United States in the 21st Century: Interdisciplinary panel report; US Department of Energy, 2006. (6) Aines, R. D. et al. Geothermal Tomorrow 2008; US Department of Energy, 2008. (7) Fridleifsson, G. O.; Elders, W. A. Sci. Drill. 2007, 4, 26−29. (8) Best, M. G. Igneous and Metamorphic Petrology, 2nd ed.; John Wiley & Sons: New York, 2002. (9) Cygan, R. T.; Kubicki, J. D. Molecular Modeling Theory: Applications in the Geosciences, 1st ed.; The Mineralogical Society of America, 2001. (10) Bychuk, O. V.; O’Shaughnessy, B. J. Phys. II 1994, 4, 1135− 1156. (11) Bychuk, O. V.; O’Shaughnessy, B. J. Colloid Interface Sci. 1994, 167, 193−203. (12) Bychuk, O. V.; O’Shaughnessy, B. Phys. Rev. Lett. 1995, 74, 1795−1798. (13) Kampen, N. V. Stochastic Processes in Physics and Chemistry, 3rd ed.; Elsevier B.V.: Oxford, UK, 2007. (14) Szabo, A.; Schulten, K.; Schulten, Z. J. Chem. Phys. 1980, 72, 4350−4357. (15) Zwanzig, R. Proc. Natl. Acad. Sci. U. S. A. 1988, 85, 2029−2030. (16) Banushkina, P.; Meuwly, M. J. Chem. Phys. 2007, 127, 1351011−135101-7. (17) Bortz, A. B.; Kalos, M. H.; Lebowitz, J. L. J. Comput. Phys. 1975, 17, 10−18. (18) Gillespie, D. T. J. Comput. Phys. 1976, 22, 403−434. (19) Gillespie, D. T. J. Comput. Phys. 1978, 28, 395−407. (20) Gilmer, G. H. Science 1980, 208, 355−363. (21) de Leeuw, N. H.; Parker, S. C. Faraday Trans. 1997, 93, 467− 475. (22) Lardge, J. S.; Duffy, D. M.; Gillan, M. J. J. Phys. Chem. C 2009, 113, 7207−7212. (23) Lardge, J. S.; Duffy, D. M.; Gillan, M. J.; Watkins, M. J. Phys. Chem. C 2010, 114, 2664−2668. (24) Tommaso, D. D.; de Leeuw, N. H. Cryst. Growth Des. 2010, 10, 4292−4302. (25) de Leeuw, N. H.; Parker, S. C. J. Phys. Chem. B 1998, 102, 2914−2922. (26) Kerisit, S.; Parker, S. C. J. Am. Chem. Soc. 2004, 126, 10152− 10161. (27) Perry, T. D., IV; Cygan, R. T.; Mitchell, R. Geochim. Cosmochim. Acta 2007, 71, 5876−5887. (28) Bruneval, F.; Donadio, D.; Parrinello, M. J. Phys. Chem. B 2007, 111, 12219−12227. (29) Raiteri, P.; Gale, J. D. J. Am. Chem. Soc. 2010, 132, 17623− 17634. (30) Raiteri, P.; Gale, J. D.; Quigley, D.; Rodger, P. M. J. Phys. Chem. C 2010, 114, 5997−6010. (31) Hoover, W. G. Phys. Rev. A 1985, 31, 1695−1697. (32) Verlet, L. Phys. Rev. 1967, 31, 98−103. (33) Wu, Y.; Tepper, H. L.; Voth, G. A. J. Chem. Phys. 2006, 124, 204503−204503. (34) Plimpton, S. J. Comput. Phys. 1995, 117, 1−19. (35) Raabe, G.; Sadus, R. J. J. Chem. Phys. 2007, 126, 044701-1− 044701-7. (36) Martoňaḱ , R.; Laio, A.; Parrinello, M. Phys. Rev. Lett. 2003, 90, 075503−075503. (37) Torrie, G. M.; Valleau, J. P. J. Comput Phys 1977, 23, 187−199. (38) Bicout, D. J.; Szabo, A. J. Chem. Phys. 1998, 109, 2325−2338. (39) Usenko, A. S.; Zagorodny, A. G. Mol. Phys. 1987, 61, 1213− 1246.

IV. CONCLUSION In this paper, we studied how the diffusion dynamics at the calcite/water interface change under different supercritical conditions (101.3 MPa, 300−800 K). From MD simulations, we found that the water conditions greatly affect the desorption free energy profile of a CaCO3 molecule on the surface. Above 600 K, the adsorption configuration makes a transition from an outer- to inner-sphere surface complex and the desorption barrier increases dramatically. As a result, the tendency for desorption no longer increases with temperature but exhibits a reversed trend above 600 K (cf. Figure 2c). This is a good illustration of the tunable solvation power of supercritical fluids, and of its effect on the diffusion dynamics of absorbed species. The use of a KMC algorithm enabled us to reach much longer time scales than would have been possible using MD simulations alone, while maintaining good accuracy. We showed that, for T > 600 K, the adsorbent molecule stays closer to the surface (a region of much smaller diffusivity), and the mean diffusivity decreases with further increase in temperature, leading to a counterintuitive reverse trend (Figure 7b). More interestingly, the diffusion on the calcite surface was found to have the typical characteristics of BMSD. We clearly showed that not only the local diffusivities but the desorption tendency and the bulk size also play important roles in BMSD. Depending on the ambient conditions, the effective surface diffusion can be superdiffusive or subdiffusive in the range of time we are interested in. We also showed that the optimized bulk thickness which maximizes the effective surface diffusivity is T dependent. We believe this factor will affect the efficiency of local surface mass transport. Overall, under supercritical condition, the solvent properties can change dramatically by adjusting T and P. Interactions among the adsorbents, surface and solvent molecules decide the dynamics at the interface. Even the simple aqueous calcite system discussed above exhibits complicated dynamical behavior, and our work provides a new strategy to study the fundamental issues of those interfacial phenomena under conditions where direct experimental interrogation is difficult. Such phenomena must be understood to evaluate reaction kinetics and fluid transport, including the relative importance of surface diffusion versus diffusion in bulk fluids,2 in a range of crustal geologic and energy systems.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

■

REFERENCES

This work was supported by the Laboratory Directed Research and Development at Los Alamos National Laboratory (LANL) under Project 20110363ER. LANL is operated by Los Alamos National Security, LLC, for the National Nuclear Security Administration of the U.S. DOE under Contract DE-AC5206NA25396.

(1) Hedenquist, J. W.; Lowenstern, J. B. Nature 1994, 370, 519−527. (2) Walther, J. V.; Wood, B. J. Contrib. Mineral. Petrol 1984, 88, 246− 259. 25941

dx.doi.org/10.1021/jp310391d | J. Phys. Chem. C 2012, 116, 25934−25942

The Journal of Physical Chemistry C

Article

(40) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes 3rd ed.: The Art of Scientific Computing, 3rd ed.; Cambridge University Press: New York, 2007. (41) Cooper, J. R. Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam, 2007. (42) Sposito, G. The Chemistry of Soils, 1st ed.; Oxford University Press: New York, 1989. (43) Revelli, J. A.; Budde, C. E.; Prato, D.; Wio, H. S. Eur. Phys. J. B 2003, 36, 245−251. (44) Revelli, J. A.; Budde, C. E.; Prato, D.; Wio, H. S. Eur. Phys. J. B 2004, 37, 205−212.

25942

dx.doi.org/10.1021/jp310391d | J. Phys. Chem. C 2012, 116, 25934−25942