Carbonates Dissolution Mechanisms in the Presence of Electrolytes

Subscriber access provided by University at Buffalo Libraries

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Carbonates Dissolution Mechanisms in the Presence of Electrolytes Revealed by Grand Canonical and Kinetic Monte Carlo Modelling Inna Kurganskaya, and Sergey V. Churakov J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b08986 • Publication Date (Web): 26 Nov 2018 Downloaded from http://pubs.acs.org on December 3, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Carbonates Dissolution Mechanisms in the Presence of Electrolytes Revealed by Grand Canonical and Kinetic Monte Carlo Modelling Inna Kurganskaya1,2* and Sergey V. Churakov1,3 1

2

University of Bern, Institute of Geological Sciences, CH-3012, Bern, Switzerland

University of Bremen, MARUM and FB5, Klagenfurter Str. 4, 28334 Bremen, Germany . 3Paul Scherrer Institute, Laboratory for Waste Management, CH-5232, Villigen-PSI, Switzerland

*

correspondence e-mail: [email protected]

ABSTRACT Dissolution of carbonate minerals as is a complex multi-step process, characterized by the particular sequence of steps dependent on pH and background electrolyte concentration. Currently, available dissolution models for carbonates do not consider dependence of the surface speciation on the local surface topography. We have developed a new approach, combining Grand Canonical (GCMC) and Kinetic Monte Carlo (KMC) methods to investigate influence of water pH and electrolyte concentration onto processes of surface charging and dissolution of carbonates. GCMC simulations of the calcite-electrolyte system are used to calculate populations of protonated sites. We consider two basic speciation models characterized by different spatial charge distributions at the surface: “Ionic”, where surface >CO32− sites are represented by “−2” charges at corresponding lattice positions; and “Oxygen”, where surface >CO32− sites are represented by triplets of “−2/3” charges at the positions of oxygen atoms. The speciation of carbonate ion protonation probabilities is found to be controlled by local charge densities and the presence of electrolyte species. In all simulation results, protonation affinity of the surface >CO32− sites followed the trend, kink (most acidic) > step > terrace (least acidic), with the same trend observed with respect to adsorption probabilities of Cl− ions. The influence of protonated site concentrations

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 37

obtained in GCMC simulations were investigated in KMC simulations. The direct comparison of simulated and experimental data showed that the Oxygen model, with an assumption of congruent dissolution, reproduces both the pH dependence of the calcite dissolution rate as well as the morphology of the calcite surface. Based on the considered model, we could identify four key factors that define pH-dependent dissolution mechanisms of calcite: (1) increase of kink site propagation rate at pH < 10; (2) increase of kink site generation frequency at pH 4-7; (3) increase of monolayer pit generation frequency at pH = 2-4; (4) acceleration of kink site propagation and generation at pH 2-4 due to the second protonation step. The combined GCMC + KMC approach shows great potential in resolving surface speciation of carbonates as functions of solvent composition and surface geometry, and their influence on the dissolution mechanisms and rates. Generally, this approach could potentially be applied to any other mineral-fluid system.

1. Introduction Dissolution of carbonate minerals is a ubiquitous process taking place in many natural and manmade systems: coastal marine and oceanic sediments1–4, marine carbonate reefs5, oil and gas reservoirs6, environments involved in the sequestration of greenhouse gases7, toxic and radioactive waste geologic repositories. Kinetic equations relating macroscopic dissolution rate and chemical composition of the fluid phase have systematically showed major failures in reproducibility when different

measurements

techniques

or

sample

preparation

protocols

are

used8.

Dissolution/precipitation reactions of minerals oftentimes take place in the contact with fluids of variable chemical composition in complex heterogeneous environments including a variety of rock textures, pore geometries, fluid flow regimes. The observed rate variance can be attributed to several scale-dependent

phenomena:

unforeseen

variations

in

solvent

chemistry or

hydrodynamics9, or intrinsic variations of reactivity attributed to the crystalline matter10. The intrinsic variance has been shown first for carbonates, where dissolution rates were represented as frequency spectra11. The reactivity of step and kink sites reflect their crystallographic orientation and coordination, and the distribution of these sites is an important source of rate variance at the atomic scale12. The influence of fluid chemical composition on rate variance is a complex problem not fully understood yet, primarily due to the two reasons: (1) adsorption of ions and protons at


2

Page 3 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


carbonate-water interface is difficult to quantify experimentally, and (2) the influence of surface charge on dissolution kinetics is not well-characterized. Our goal is thus to fill this gap by relating processes of carbonate surface protonation, ion adsorption and dissolution. The ionic structure of the calcite-water interface and its influence on calcite dissolution kinetics as a function on pH has been the subject of intense study for many decades13. The first widely used model for surface speciation as a function on pH for carbonates was introduced by Van Cappellen et al14. This model incorporated reactions of protonation, deprotonation, and adsorption of carbonate and calcium ions as major events responsible for surface charge. The dissolution rate was then calculated as a linear combination of rate constants and concentrations of charged surface sites. Duckworth et al.15 extended this approach by introducing corrections for the presence of atomic steps, and related dissolution rate to step velocities. The key terms in their kinetic equations for step velocities are concentrations of protonated carbonate and deprotonated metal sites taken from the previous surface complexation model14. Jordan and Rammensee16 related carbonate dissolution rates with curved steps formed as a result of etch pit interaction, but without effects of protonation and pH influence. Our approach here is to consider all possible kinetic features, e.g. steps moving in different crystallographic directions, etch pits, and the chemistry-dependent speciation of the surface. We employ the probabilistic Kinetic Monte Carlo approach, where all reactive sites explicitly contribute to the total material flux depending on their reactive probabilities17,18. In this study we further extend the previously-developed Kinetic Monte Carlo (KMC) model for dissolution of carbonate minerals in the presence of lattice defects12, incorporating the influence of pH and electrolyte concentration on the dissolution kinetics. Our essential goal is to include molecular scale description of surface speciation into KMC model by taking into account the influence of structural and chemical factors on surface site acidity.


3


Page 4 of 37

Various models of carbonate surface speciation that take site geometry and other details into account have been developed during the last decade. Wolthers et al.19 modified a surface complexation model to fit experimentally measured zeta potentials accommodating different partial charges on terrace, ledge and kink sites. They employed the theoretical CD-MUSIC approach to calculate acidity (pKa) constants for sites with different coordination and orientations (obtuse and acute)20 based on the site partial charges and O-H bond lengths. Andersson and Stipp employed density functional theory based-simulations to directly evaluate pKa values for various site types, taking in account local atomic arrangements of carbonate groups, oxygen, calcium and water molecules21,22. The acidity of surface sites controls their protonation state at different pH conditions. However, the concentrations of ionized sites are strongly affected by the presence, ionic strength, and chemistry of electrolytes, together with the surface density of charged sites23. In order to incorporate all these effects, we developed a new surface speciation model, based on the Grand Canonical Monte Carlo Simulation (GCMC) approach for surface site protonation and ion adsorption developed previously24–26. The dedicated GCMC simulation algorithm25 was used to calculate concentrations of charged sites in the presence of electrolytes and as a function on pH on carbonate mineral surface. The new model fully accommodates the 3D geometry of the dissolving mineral surface, with explicit etch pits, kink and step sites. Interactions in the electrolyte solution and the ion surface interaction are modeled using the Primitive Model (PM) of electrolytes27. In the PM model ionic species are considered as charged, hard spheres, and solvent molecules are replaced with a dielectric continuum. We introduce here a series of models having different levels of granularity and model parameterization of the mineral surface. Initially, we test the role of ion charge localization by either considering the single carbonate groups on the surface as point charge


4

Page 5 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


particles or distributing corresponding point charges on the oxygen sites of >CO32− group explicitly. We refer to these as the Ionic and the Oxygen model, respectively. Next, we investigate the effect of the geometric structure of the surface, the structure of the etch pit and steps. The site speciation statistics was calculated for a series of pH values at electrolyte concentrations representing experimental conditions in the previously published AFM studies28. This information is used to create protonated site populations in the KMC model for carbonates12 updated for taking in account this effect. We perform the KMC simulations with charged site populations obtained from different GCMC models to test their performance by comparison with experimental data available from literature. For this purpose, we review the literature with respect to published data for atomic step velocities, etch pit structures and bulk dissolution rates as a function of pH. Compared to previous studies, for the first time we couple molecular-based GCMC description of surface speciation assuming local equilibrium conditions, with the stochastic KMC modeling of carbonate dissolution process. Use of this innovative approach allows us to elucidate key atomic scale mechanisms responsible for the pH dependence of calcite dissolution rate and to quantify the individual roles of ion adsorption and surface topography. The approach is generic and can be used regardless of the chemistry of the cation. Calcite was chosen as it is by far the best-studied of all carbonate minerals.

2.Methods 2.1. Grand Canonical Monte Carlo model 2.1.1. GCMC approach GCMC simulations are used to obtain surface speciation and ion distributions at the calciteelectrolyte interface, assuming local chemical equilibrium with an external electrolyte solution of a fixed composition and proton activity (e.g. imposed pH conditions). All particles in the system


5


Page 6 of 37

(surface sites) and ions in the electrolyte are considered as charged hard spheres interacting via coulombic potential. The solvent molecules are replaced by a dielectric continuum with dielectric constant of bulk water equal to 78 (dimensionless relative electric permittivity). This is the wellknown primitive model of electrolytes which allows us to take into account the ion-ion correlation effect and steric exclusion. The PM model has been successfully applied to mineral surfaces (silica24, C-S-H25), colloidal systems29–31, proteins, and highly concentrated mixtures of electrolyte solutions. The model is in accurate prediction for ionic activity coefficients32 and experimentally measured zeta-potentials25,33. The equilibrium distribution of ions in the system is obtained by randomly changing an ion’s position and inserting or removing charge-neutral ionic pairs using the standard Metropolis algorithm34. The surface sites are also allowed to protonate and deprotonate according to the following reaction: +

↔

>

+

where

is an arbitrary negatively charged ion involved in the corresponding removal or addition

of a neutral ion pair

(1)

from the system. The protons in the protonation/deprotonation reactions

are treated implicitly via proton activity26. In this approach the Boltzmann factor

− ∆

for

acceptance probability for the de-protonation and protonation steps is calculated according to the equations 2 and 3, respectively: − ∆ − ∆ with

= *+

= =

)

−

− ∆

!

"#$10

+

− ∆

!

"−#$10

− ' ( (2) − ' ( (3),

where * and T are the Boltzmann constant and temperature,

chemical potential of the ion B, , is the number of ions of type B, ∆

!

is the

is the change in

electrostatic potential energy of the simulated system, V is the total volume of the simulation


6

Page 7 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cell, ' = −#-.

' , ' is the intrinsic acidity constant of the corresponding surface

site, assuming a charge neutral system at infinite dilution. Further details of the simulation algorithm and its implementation are described in Churakov et al.25 2.1.2. System setup for GCMC simulations The surface sites are represented by spherical ions with radii of 2 Å placed at fixed lattice positions. 2D Periodic Boundary Conditions are applied in the simulation box in directions parallel to the surface. The simulation box is aperiodic in the direction normal to the surface. For the sake of computational efficiency, the {10.4} rhombohedral calcite plane was mapped to an orthogonal lattice while maintaining the site density (9.8 sites/nm2). The simulation box size used throughout the simulation is 64.25 × 64.25 × 300.00 Å. The long range electrostatic effects are taken into account by charge sheets method35. The calcite (CaCO3) crystal lattice is described as a binary system of positively and negatively charged hard spheres (Fig. 1A-C). We use the commonly accepted Terrace-Ledge-Kink site classification for carbonates, taking in account obtuse and acute step and kink orientation (Fig. 1D). The modeled surfaces follow 3D geometries reproducing various dissolution features occurring on carbonate surfaces: etch pits, hollow cores, atomic steps and kinks (Fig. 2). The simulations were performed for a system in equilibrium with 0.1M bulk electrolyte solution containing NaCl + (NaOH or HCl) similar to the experimental conditions used by Shiraki et al28. A range of pH conditions from pH 2 to pH 12 was modelled by adding corresponding amounts of HCl and NaOH. The concentration of Na+, Cl−, OH−, H3O+ ions in the corresponding bulk fluid is listed in Supporting Information. The chemical potentials of electrolyte species at given pH were calculated separately in NVT ensemble using Widom’s insertion method36,37,25. The calculated ion chemical potentials at every pH point where then used in GCMC simulations of the calcite-fluid


7


Page 8 of 37

interface. The potential presence of dissolved carbonate and Ca2+ ion was neglected to simulate far-from-equilibrium conditions with respect to the calcite dissolution. In the GCMC simulations the surface does not dissolve but only undergoes two-step protonation and deprotonation reactions at anionic (eqs. 4,5) and cationic sites (eqs. 6,7), respectively: ↔ > /012 +

> /01 > /01

2

↔ > /01

> /30

2) 2

> /30

)

↔> /30

) )

+ )

↔> /30 +

(4) (5)

+

)

)

(7)

(6)

The intrinsic acidity constants of surface sites have been taken equal to the acidity constants of aqueous CO32− ion and hydrated Ca2+ ion, pKa(>HCO3−/>CO32−) = 10.3522, pKa(>Ca(OH)2 />CaOHO−) = 13.021. The constants for the second reaction steps are given as pKa(>CO3H20 />CO3H−) = 3.6, pKa(>CaOH+/CaO0) = 16.0. The intrinsic acidity represents the reactivity of the site in a hypothetical system at infinite dilution, assuming an absence of electrostatic interaction with the lattice sites. All surface specific effects on the effective acidity of surface sites are properly represented as the electrostatic interactions of a site with surrounding lattice neighbors and ions in the electrolyte solution. 2.1.3. “Ionic” and “Oxygen” surface speciation model Two types of surface speciation models are considered: a “Ionic” model, where carbonate ion sites are represented by the “−2” charged particles placed at the carbon atom positions (Fig. 1B), and an “Oxygen” models, where carbonate ion sites are explicitly represented by three oxygen atoms carrying “−2/3” charge each and distributed according to the geometry of CO32− group (Fig. 1C). 2.1.4. Surface site crystal lattice interaction


8

Page 9 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Further we have tested the convergence of the system setup (e.g. effective surface protonation and ion sorption) with respect to the electrostatic field created by the bulk lattice substrate. We constructed a set of systems with up to three additional explicit bulk lattice layers placed under the reacting (outermost) surface. We performed these convergence tests using both “Ionic” and “Oxygen” models for the bulk lattice layers (Fig. 1D). The modelling results for the surface speciation at perfect {10.4} cleavage plane of calcite obtained with bulk layers represented by the Ionic and the Oxygen model were found to be identical. The observed convergence of the electrostatic layers with respect to the number of surface-adjacent layers is related to the fact that each layer is charge neutral and non-polar. For the sake of computational efficiency, we therefore represent the bulk layer using the “Ionic” model in all further simulations. The convergence tests using setup with 0,1,2 and 3 “bulk” layers showed that system size convergence is achieved with 1, layers (see SI). Therefore, we used the systems setup with one “bulk” layers in all models

2.1.4. Role of the surface topography We studied the influence of the surface geometry and different arrangements of the reactive features, e.g. etch pits, steps and kinks, on the site protonation and ion adsorption at the interface. For that purpose, we generated eight different configurations representing various topographical environments: flat plane (Fig. 2A), monolayer pit with double kink (Fig. 2B), monolayer pit with single kinks (Fig. 2C), monolayer curved pit with complex rough steps (Fig. 2D), multilayer etch pit with bunching acute steps (Fig. 2E), multilayer pit with curved rough steps (Fig. 2F), two coalescing etch pits with straight outlines and curved steps between them (Fig. 2G), and two coalescing pits with curved rough outlines (Fig. 2H).


9


Page 10 of 37

Fig. 1. Calcite surface models for GCMC titration experiments. Balls represent ions, triangles represent carbonate groups. A. Geometric orientation of steps and kink sites; B. The Ionic model; C. The Oxygen model; D. Side view of the surfaces constructed for (from top to bottom) one bulk layer (L1), three bulk layers (L3), one Ionic bulk layer, one Oxygen bulk layer.

Fig. 2. Surface geometries for GCMC titration models. A. Plane terrace; B-D. Monolayer pits, B. Straight step model 1 (Mono-Str-1); C. Straight step model 2 (Mono-Str-2); D. Curved step model


10

Page 11 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(Mono-Curved); E. Single pit with straight steps (1-Pit-Str), E2 side view; F. Single pit with curved steps (1-Pit-Curved); G. Two coalescing pits with straight steps (2-Pits-Str), G2 side view; H. Two coalescing pits with curved steps (2-Pits-Curved), H2 side view. The positions of surface sites in each simulation are fixed, and no dissolution/precipitation reactions take place in GCMC simulations. The system is in equilibrium with respect to pH and the ionic composition of the external electrolyte solution. Equilibrium distributions of ions in the electrolyte and surface speciation are obtained by statistical averaging of 1000 uncorrelated snapshots from GCMC simulations, using every 1000th configuration. The production runs comprise 106 GCMC steps followed by the same number of pre-equilibration steps. This number was found to be sufficient for these systems to attain equilibrium state suitable for statistical sampling. The ions are considered to be “adsorbed” on a surface site when the radial distance to the surface site is within the 5 Å cutoff radius. This distance corresponds approximately to the mid-point between the first and the second coordination shell in the closest packed lattice (see Supporting Information). The sites are classified according to the common “obtuse-acute” notation (Fig. 1A).

2.2. Kinetic Monte Carlo model 2.2.1 Kinetic Monte Carlo method The Kinetic Monte Carlo method is a widely used approach to simulate temporal evolution of a polyatomic system. We employed the model previously designed for carbonate minerals12. The calcite surface sites are either calcium or carbonate ions. The dissolution rate of each surface group is defined as * 4,6,7

/791

=:∙

?@/?AB CD

E∙

Ca or >CO3

site for protonation state of >CO3 (in the case of >Ca site the protonation states of neighboring >CO3 sites are counted to form j term). Although the de-/protonation reaction both >Ca and >CO3 is considered in the speciation modelling according to the eqs. 4-7, the >Ca site is present as > /30

2) 2

in the entire range of pH conditions studied. No orientation correction is used for kink

sites for reasons of model simplicity and parameter number reduction. From a mechanistical point of view, setting different activation energy barriers for obtuse and acute step site dissolution is enough to create difference in step velocities observed in AFM experiments28,38–41. The reactivity difference for kink sites of different orientations is problematic to obtain from direct experimental observations. If detachment rates at different kink sites will be established experimentally, the difference should be considered in the next generation of KMC models. The model was parameterized to reproduce atomic step velocities in the pH region 8-12 measured by Shiraki et al.28 under far from equilibrium conditions, 0.1M NaCl. Their data exhibit a substantial variance especially in the region of pH < 8, probably due to the uncontrolled conditions with respect to the saturation state of the solvent. These conditions are hard to avoid experimentally for calcite in the acidic region because of the fast dissolution kinetics coupled with the transport of material to and from the surface42. Our goal then was to reproduce general trends for step velocities at pH = 8-12 where dissolution process is dominant and back-precipitation reactions should not take place. The only fitting parameter in our extended KMC model is the in the activation energy for the protonated Ca-HCO3 bond ∆HK

79B

, which is and intrinsic pH


12

Page 13 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


independent variable. We therefore use most reliable dataset (e.g., pH 8-12) for the parameter estimation. The activation energy is a pH independent parameter and the model calibration at pH 8-12 is valid for the entire pH range. 2.2.2 Incongruent and Congruent dissolution models Different approaches were applied for parameterization the KMC model for calcite dissolution. All parameter sets are summarized in Table 1. For the Incongruent model, we assume that the dissolution probability of a given site is influenced by the protonation state of that particular site only. Thus, the protonation state of a carbonate site does not influence the detachment probability of the neighboring calcium sites and vice versa. This Incongruent KMC model was applied to the results of speciation calculations obtained for the Ionic and the Oxygen GCMC models. The essential feature of the Incongruent KMC model is that the dissolution probabilities of carbonate and calcium sites have distinct pH dependence. In the Congruent model the dissolution probability of neighboring Ca and carbonate ions are affected by the protonation state of carbonate sites. The reactivity of two-protonated sites >CO3H20 in all models is calculated by using the eq 7. However, the stability of such sites is a matter of debate, since carbonic acid has a tendency to decompose fast (25 s−1 at T= 25 °C)43 into CO2 and H2O. In order to test how kinetic pathways will change if the doubly-protonated site >CO3H20 instantaneously leaves the surface, we considered two different parameter sets (1 and 2) for the Congruent Oxygen. Set 2 repeated the parameters used in set 1, but the dissolution probability of all doubly-protonated >CO3H20 sites is set to 1. The results for the step velocity parameterization are shown on the plots for both the Ionic and the Oxygen models (incongruent dissolution case) (Fig. 3).


13


Page 14 of 37

Dissolution patterns and material fluxes are calculated for the system sizes 64.25 × 64.25 nm, with one “infinitely long” screw dislocation, and lateral periodic boundary conditions.

Table 1. KMC simulation parameters. 2

∆HI , kT units

∆H , kT units ∆HK 79B , kT ∆HK 7 , units units

0.005

6.35

2.5

−5

0

Ionic congruent

0.005

6.35

2.5

−6

−3

Oxygen congruent

0.01

6.35

2.5

−5

0

Oxygen 0.01 congruent 1,2

6.35

2.5

−3

−3

model

v, 10

Ionic incongruent

L

kT

3. Results 3.1. Grand Canonical Monte Carlo simulations 3.1.1. Ionic vs. Oxygen models The dependence of surface protonation on pH for the Oxygen and Ionic models is shown in Fig. 3. The pH dependence of protonation fractions for the carbonate groups at the terrace sites are nearly identical for both the Oxygen and the Ionic models (Fig. 3A). Less than 1% of these sites undergo the second protonation step. The step sites become protonated at pH=5 and about 15% of these sites become doubly protonated at pH=2 (Fig. 3B). Oxygen and Ionic models predict very different protonation behavior for kink sites. The Ionic model predicts a quasi-linear increase in


14

Page 15 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the fraction of singly protonated kink sites from just 10% at pH=10 to nearly 90% at pH=4. In contrast, the Oxygen model predicts that the fraction of singly protonated sites reaches a maximum at pH=6. At lower pH the fraction of singly protonated sites decreases while the fraction of doubly protonated sites increases. The Oxygen model predicts that the fractions of singly and doubly protonated sites become equal at pH=2.5 (Fig. 3C). The Oxygen model clearly enables the two step protonation mechanism and predicts formation of surface > CO3H20 species, which can be potentially be converted at the surface into H2CO3,aq.

Fig. 3. Protonation fractions for calcite surface sites as functions on pH, a comparison between Ionic and Oxygen models for surface charge distributions. A. Terrace sites; B. Obtuse step sites; C. Kink sites k_oa. 3.1.2. Effect of surface topography on the surface speciation The pH-dependent speciation of surfaces sites with different surface topography (Ionic model) is shown on Fig. 4. Protonation behavior for terrace sites (Fig. 4A-4C) is almost the same for plane and monolayer pit geometries (Fig. 2A-C), while the results for complex geometries show different pH trends. The step site statistics for the monolayer structures is similar (Fig.4B-C), and with the exception of the results obtained for obtuse steps for more complex structures, such as 2-PitsStraight (Fig. 2G) and 2-Pits-Curved (Figs. 2H and 4C). The protonation fractions of kink sites


15


Page 16 of 37

obtained for different surface geometries and site orientations show larger variability (Fig. 4D-F). At the neutral and higher pH the variations can be as large as 4 pK units. This variability becomes much smaller at lower pH. The model predictions for pH below 4 are very similar for all types of kink sites. The results clearly show that second order neighborhood, vicinity of hollow cores, bunched steps, step curvature and other details have influence onto the protonation probabilities.

Fig. 4. Effects of surface geometry onto protonation fraction of calcite surface sites, Ionic model. A. Terrace sites, B. Acute step sites; C. Obtuse step sites; D. Kink sites k_ao, E. Kink sites k_aa, F. Kink sites k_oo.

2.1.3. Influence of site coordination and ions adsorbed at the interface


16

Page 17 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The comparison of protonation fractions for different sites clearly demonstrate dependence of the protonation probability on site coordination (terrace, ledge or kink). Kink sites are the easiest to protonate, and the protonation fraction starts to increase from at pH ≤ 10 (Fig. 5A). The step sites are almost two times less likely to protonate. The pH-dependent speciation curve for steps is shifted with respect to the corresponding curve for kinks by 2-3 pH units towards decreasing pH. The protonation equilibrium curve for terrace sites is shifted by 2 pH units compared to step sites (Fig. 5A). Simulations of a hypothetical system in which the electrostatic interaction between ions and the surface has been switched off show that surface protonation is substantially suppressed, and only the kink sites are able to protonate. The extent of the protonation is reduced by about five times (Fig. 5B). The same trend for the protonation probability, kinks > steps > terrace, is followed, but the increase in the number of protonated sites as pH-decreases is much less steep. This finding demonstrates the pivotal role of ion surface interaction for protonation/deprotonation equilibria at the surface. The same effect has been observed on the surface of amorphous silica by Porus et al.24, who showed that repulsive electrostatic interactions between surface sites is compensated by the ions’ surface interaction at the interface. In addition to this general trend, we also demonstrate the kinetically most important surface sites, e.g. steps and kinks, are less affected by the interaction with the neighboring surface sites. This is likely to be the consequence of typically larger distances between neighboring kink sites. The amount of the sites adsorbing Cl− ions as function of pH is shown on Fig. 5C. The adsorption affinity of Cl- to the surface sites is described by a clear trend kink > step > terrace. The trend is the same as the protonation affinity of the surface sites. The Cl− ions are electrostatically attracted to calcium sites, but reside in the close vicinity of carbonate groups, enhancing their protonation. Comparisons of all three plots (5A-C) show that all pH-


17


Page 18 of 37

dependent protonation is largely due to the presence of adsorbed ions. The pH-dependence of ion adsorption stems from the pH-dependence of carbonate protonation, which creates positive charge attracting Cl− ions. In turn, the adsorption of Cl− ions enhance protonation of the surface sites, attracting even more Cl−, so the process reinforces itself.

Fig. 5. Effects of adsorbed ions on the site protonation probability, statistics is shown separately for each site type, the surface geometry is monolayer curved pit (Fig. 2D). A. Protonation fractions of protonated >CO3H+ sites in the presence of 0.1M NaCl; B. Protonation fractions of sites with ions prevented from approaching the surface; C. Fractions of sites containing at least one Cl− ion within the cutoff radius of 5 Å 3.2. Kinetic Monte Carlo simulations 3.2.1. Surface morphology Etch pit structures and surface morphologies obtained by KMC simulations from the set of parameters listed in Table 1 and surface speciation datasets obtained with the Ionic and the Oxygen model are shown in Fig 6. The predicted etch pit morphology exhibit clear pH-dependence in all models studied, a consequence of the pH dependent change in the surface speciation. The etch pit and atomic step morphology depends on the model and parameter choice, but the overall surface structure change follows a common trend from alkaline pH (8-12), where atomic steps exhibit


18

Page 19 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


straight morphology, to acid pH (2-4), where etch pits exhibit a variety of morphologies depending on model and parameter choice. At pH 12 carbonate groups do not protonate, steps exhibit straight morphologies (Fig. 6C,F,I,L). Etch pits at pH 7 maintain rhombohedral morphology, but the step and kink site density increases (Fig. 6 B,E,H). Pits at pH 2 exhibit curved morphologies with the greatest kink site density (Fig 6A,D,G,J). Incongruent (Fig. 6A-F) and congruent (Fig. 6G-L) models show distinct difference with respect to the density of kink sites along the steps and density of the steps, and the step curvature. The mechanisms generating step curvature at low pH are discussed further below.


19


Page 20 of 37


20

Page 21 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 6. Etch pit morphologies, KMC simulations, shown surface sizes are 60x60 nm. A-C Ionic model, incongruent parameter set 1, D-F: the Oxygen model, incongruent parameter set 2, G-I Ionic model, congruent parameter set 3, J-L the Oxygen model, congruent parameter set 4. The congruent oxygen model with immediate dissolution of >CO3H20 sites results in breaking the mechanism of the step and terrace formation at low pH (Fig. 7). The number of terrace sites with the dissolution probability equal to unity increases with decreasing pH, as a result rough pitted surfaces dominate at this regime. Such a behavior has not been observed in the AFM experiments so far.

Fig. 7. Surface structure and dissolution mechanism at different pH values, KMC simulations, shown surface sizes are 60x60 nm. The oxygen congruent model 2 (instantaneous dissolution of all >CO3H20 sites). 3.2.3. Material fluxes Dissolution rates calculated as “steady-state” material fluxes from the surface are shown on Fig. 8. Experimental observations refer to the flow throw AFM experiments, where rates are calculated by analyzing effluent concentrations. According to our simulation data (all models) the rates in the pH region 5-11 grow linearly spanning a range about two orders of magnitude. Evolution of dissolution rates predicted by Incongruent models exhibit a characteristic plateau at pH 2-4. All fluxes obtained in our simulations refer to the systems with “infinitely” long screw dislocations, except the data shown as gray symbols, which were obtained from the Oxygen


21


Page 22 of 37

Congruent model 1 run without any dislocations. In the latter case dissolution is sustained by spontaneous formation of point defects and monolayer pit mechanism. The rates predicted by using the Oxygen congruent model 2 (instantaneous dissolution of all >CO3H20 sites) span six orders of magnitude.

Fig. 8. Calcite dissolution rates vs. pH. Experimental data and simulations. Rates are calculated for pitted surfaces in the screw dislocation control regime, the special case for the absence of screw dislocations (monolayer pit regime) is shown for the Oxygen congruent model 1, grey squares.

4. Discussion 4.1. Surface speciation and site acidity Acidity constants for calcite surface sites have been a subject of intense study. Several formulations assuming different number of distinct surface sites have been proposed. Early models operated with the two-site models, fitting thermodynamic properties of >Can+ and >CO3m− sites to the experimentally determined surface charge and zeta potentials. The effective site charge n and m are the models specific parameters, which are typically ranging between 1 and 2. Pokrovsky et al. estimated effective pKa = 5.1 for carbonate protonation reaction44 from electrophoretic


22

Page 23 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


measurements, Van Capellen et al. derived pKa = 4.9 from experimental titration data14. Due to the differences in the model formulations, these values should not be compared directly compared to the intrinsic pKa values used in our simulations, which are the same for all sites with identical chemical composition. Instead we can discuss the effective acidity of the different surface sites encountered in our GCMC model due to the specific structural coordination of the surface sites. On the basis of structural environment shown in fig. 2 we can recognize 14 structurally nonequivalent sites. (7 for >Ca2+, 7 for >CO32−). Based on our detailed simulation results, we define an effective pKaeff value for the individual surface species as the value of pH at which the concentration of protonated and deprotonated sites of the same type is equal. In this way, we obtain an effective pKaeff = 4.5 for step and pKaeff = 7 for kink sites. These values are comparable with the estimates of Pokrovsky and Van Capellen assuming a one-site model. In contrast, the effective pKaeff for terrace sites (the most abundant sites) are around 2.5 (Fig. 6). Wolthers et al.20 used bond valence theory (CD-MUSIC) combined with Molecular Dynamics simulations to refine O-H bond lengths for estimation of pKa for sites of different geometry. Their pKa values range from negative to 0, and all experimentally measured zeta potentials of calcite surface were attributed to adsorbed carbonate and calcium ions. Zeta potentials in these studies refer to equilibrated calcite suspensions19, systems which are thus chemically different from our far-from-equilibrium conditions. Andersson and Stipp21,22 calculated pKa values for different sites by using COSMO-RS approach and obtained negative pKa values for carbonate sites on the surface. According to our simulation results, the pKaeff values for carbonate sites in the absence of background electrolyte are negative as well, and protonation is strongly defined by the adsorbed ions. 4.2. Dissolution mechanisms 4.2.1. Etch pit and step morphologies


23


Page 24 of 37

Experimental AFM observations38 of calcite atomic step morphologies show clear dependence on pH in the acidic region 1.7-4.3, (Fig. 9 A-C).

Fig. 9. Surface morphologies and etch pit structure obtained in experiments. A-C 0.1M NaCl, pH adjusted by HCl/NaOH38. Adapted from the original publication: Giudici, G. D. Surface Control vs. Diffusion Control during Calcite Dissolution: Dependence of Step-Edge Velocity upon Solution pH. Am. Mineral. 2002, 87 (10), 1279–1285. Copyright (2002), used with the permission


24

Page 25 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


from the Mineralogical Society of America; D-F water, pH adjusted by HCl/NaOH45. Adapted from the original publication: Atanassova, R.; Cama, J.; Soler, J. M.; Offeddu, F. G.; Queralt, I.; Casanova, I. Calcite Interaction with Acidic Sulphate Solutions: A Vertical Scanning Interferometry and Energy-Dispersive XRF Study. Eur. J. Mineral. 2013, 25 (3), 331–351. Copyright (2013), used with the permission from the Schweizerbart Science Publishers, www.schweizerbart.de/journals/ejm; G-I KMC simulations, The Oxygen congruent model 1. The step roughness increases at low pH as well (Fig. 9). At low pH (1.7-3) normal “Ledge-Kink” step structure is broken, the steps are composed of rough irregular clusters (Fig. 9A), which are curved and bunched (Fig. 9B,D,E). At pH =4.3 and higher the steps exhibit straight regular euhedral morphology (Fig. 9C,F). Comparison of these data with our modelling results (Fig. 9GI) demonstrates that the congruent Oxygen model best reproduces the typical, experimentally observed step morphologies as a function on pH. All other models result in curved steps at pH 47, which have not been observed experimentally. This is a strong argument that the congruent Oxygen model-1 captures the most important molecular mechanism controlling the calcite dissolution in the wide range of pH. We consider here three primary mechanisms that may control step curvature: (1) Straight steps are formed, when probability of kink site dissolution >> probability of step site dissolution; (2) Curved steps are formed when kink sites are less reactive then step sites and accumulate, creating so called step pinning effect46; (3) Rough steps are formed when dissolution rates for kinks and steps are comparable. All these effects may take place on calcite surfaces and contribute to the step morphology. In the Incongruent models (both Ionic and Oxygen) surface sites are substantially less reactive and accumulate on the surface, poisoning reactive sites and creating inhibiting step pinning effect


25


Page 26 of 37

leading to the curved step morphology46,47 (Fig. 10A,C). The formation of the step curvature due to accumulation of >Ca2+ kink sites is the dominant mechanism for this model. This mechanism is even more pronounced for the Oxygen model, where two-protonated sites leave the surface faster and facilitate formation of >Ca2+ kink sites (Fig. 10C). In the Congruent Ionic models, in low pH region (pH=2-4) both kink and step sites are protonated and the ratio between their dissolution probabilities is >>1, so the straight step structure is maintained (Fig. 10B). Mechanism controlling the step morphology in the Congruent Oxygen models is the most complicated: >CaHCO3 stepkink pairs are formed, but Ca2+ step sites are highly reactive and facilitate formation of >CO3H/>CO3H20 kink sites that break down >CaHCO3 pairs (Fig. 10D). As a result, the moderate step curvature is maintained by the dynamic balance between two competing processes shown on Fig. 13B and 13C.


26

Page 27 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 10. Mechanistic block models of calcite dissolution at low pH (CO3H− kink sites, dissolution rate of >CO3H− step sites is about the same as for >Ca2+ kink sites (Rstep(>CO3H−)/ Rkink(>Ca2+) = 0.7), the formation of these sites is highly favorable. B. Ionic Congruent model, straight step regime, step sites are always much more abundant than kink sites, due to the contrastingly different


27


Page 28 of 37

reactivity. C. The Oxygen Incongruent model, pinned curved step regime. The mechanism is almost equivalent to that shown in B, but dominance of >Ca2+ kink sites is even more pronounced, due to the abundance of highly reactive >CO3H20 kink sites, enhanced formation of >Ca2+ kinks in the process of >CO3H20 step site dissolution, dissolution rates of >CO3H− step sites is comparable to dissolution rates of >Ca2+ sites (Rstep(>CO3H−)/Rkink(>Ca2+) = 0.26). D. The Oxygen Congruent model, regime of curved steps with zig-zag pattern, formed by the stable ion pairs >Ca2+ kink site – >CO3H− / >CO3H20 step sites, >Ca2+ step sites are more reactive then dominating >CO3H− step sites, >CO3H20 kink sites are very unstable and do not accumulate. 4.2.2. Rates as material fluxes Experimentally measured rates of calcite dissolution span over two-three orders of magnitude in the pH 2-11 region. Our simulation results predict same order of rate variation, but the absolute values of the dissolution rates obtained in the simulations are about 2-3 orders of magnitude higher (Fig. 8). Please note, that our simulations use an idealized system setup with infinitely long screw dislocations at fixed positions and the rate values are sampled out of the “steady-state” material fluxes representing the highest ranges of rates possible (all etch pits are coalescing and giving the maximal number of steps and kinks possible at these conditions). The simulation for the system without screw dislocations (the Oxygen congruent model) resulted in an order of magnitude lower rates, showing trends and values similar to those observed by Chou et al47. The “plateau” in the pH region 5-12 is commonly observed in a series of experiments8, but in a similar set of experiments the pseudo-linear dependence also present38,47. The Oxygen model predicts linear increase of the dissolution rates in this region. Chou et al.47 noted that both dissolution and precipitation processes are taking place and balance between them control overall material flux for calcite.


28

Page 29 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


We emphasize that the current model does not take into account the potential contribution of reprecipitation process (back reaction). We suggest that the “plateau” or weak pH dependence of the dissolution rate at pH 5-12 could be a result of a competition between protonation-assisted dissolution and fast back precipitation reaction. This balance is then perturbed at lower pH, and dissolution overrides precipitation47. To the best of our knowledge, the plateauing of the dissolution rate in the pH region 2-5 as predicted by Ionic models has not been observed experimentally so far. This can be explained in two ways: (1) current Ionic model does not account for the relevant molecular processes at low pH, e.g. second protonation step as observed in the Oxygen model; (2) published experimental data in pH region 2-5 refer to the condition at which the other processes, e.g. back-precipitation and transport-controlled dissolution may play an important role, and thus not represented in the “pure” dissolution rates which is modelled in our KMC simulations. The potential presence of other reaction mechanisms is highly speculative. They can neither be confirmed nor excluded with certainty on the basis of available observation. We thus conclude that the Oxygen model best predicts the pH-dependent change of the calcite dissolution rate, and the corresponding evolution of the surface morphology as function of pH. 4.2.3. Reaction mechanisms and rate limiting steps pH-dependence mechanism of calcite dissolution is largely defined by the pH dependent speciation of the carbonate groups (e.g. their protonation state) and the corresponding speciation dependent activation energies for dissolution of carbonate and calcium sites. The reactivity of both >CO3H- and >CO3H20 sites strongly affect the model predictions. One would expect that doubly protonated sites >CO3H20 are very unstable and dissociate instantaneously to CO2 and H2O, as happens with carbonic acid in bulk water42,48. In the model where all >CO3H20 sites have


29


Page 30 of 37

dissolution probability equal to 1 (the Congruent Oxygen 2), we cannot reproduce the etch pits observed experimentally at these conditions. One possible explanation is that terrace and step sites are also largely stabilized by lattice coordination, and the transformation rate of carbonic acid molecules into CO2 and water is suppressed. From this study, we can conclude that the dissolution rate is controlled by three pH dependent factors: (1) increase of the step velocity, (2) increase in the step density38 and (3) an increase in monolayer etch pits formation probability as a function of decreasing pH49. The increase of the monolayer pit density with decreasing pH is confirmed experimentally for dolomite49, and it is likely that the same mechanism takes place on the calcite surface at low pH. The key assumption of all the models considered in this study is that protonation of carbonate groups occurs faster than any other process, such as protonation-assisted detachment of step site, either Ca2+ or carbonate. This process controls the densities and spatial distribution of the most reactive sites. These model predictions largely depend on whether congruent or incongruent dissolution process takes place, i.e. whether >Ca2+ sites are affected by the protonation of carbonate sites. In the incongruent models, >Ca2+ site dissolution probability is unaffected by protonation of neighboring carbonate ions. >Ca2+ sites therefore accumulate and limit the dissolution rate (Fig. 10A,C). Therefore, at pH 5-6 and lower, the dissolution process is controlled by the dissolution rate of >Ca2+ sites. The incongruent model, therefore, predicts a plateauing of the dissolution rate at pH 2-6. In the congruent models, probabilities of both >Ca2+ and carbonate sites increase with protonation, and therefore the rate defining steps are kink site generation frequencies, which gradually increase with decreasing pH.


30

Page 31 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The linear increase of the dissolution rate with decreasing pH is observed only for the Congruent Oxygen model 1, because at pH 5 the second protonation step creates new reactive sites at both step and kink positions. This is the only model showing linear pH dependence observed experimentally for carbonate grain systems8, spanning ~2.5 orders of magnitude, compared with the 2-3 orders of magnitude increase observed experimentally. 4.2.4. Role of electrolytes The role of electrolytes on dissolution kinetics is a matter of experimental investigation. According to Ruiz-Agudo et al.50, an increase of ionic strength from 0.001 to 0.1 M NaCl results in five-fold increase of calcite dissolution rate, while etch pit spreading rate increases by factor of two. In contrast, Finneran and Morse51 showed that NaCl retards the dissolution process. Their ionic strength spanned a larger range, from 0 to 6 molal, over which rates decreased by an order of magnitude. Ruiz-Agudo et al.50 suggested that the electrolyte may have a stabilization effect on surface sites by adsorbed ions (inhibitory effect) at ionic strengths < 0.1 M, or by modification of the interfacial water structure at higher ionic strengths (catalytic effect). According to our simulation results, in 0.1 M NaCl solution there is a clear catalytic effect due to the significant increase in density of protonated sites. The adsorption of Na+ ions occurs in such statistically negligible concentrations as to preclude a kinetic effect, whereas Cl− ions adsorbs in significant amounts and may also influence the site dissolution probability, either inhibitory or catalytic. This aspect should be a subject of further investigations.

Conclusions We offer a new approach for the studies of carbonate mineral dissolution kinetics based on a combination of Grand Canonical and Kinetic Monte Carlo methods. The GCMC model captures


31


Page 32 of 37

the effect of electrolyte and local coordination environment of surface sites onto probability of site protonation and deprotonation. We use this information to reflect the populations of surface sites with enhanced reactivities as function of pH in a KMC model. We tested a series of surface and lattice charge models used in GCMC simulations and compared their predictions against experimental observations. The results reveal the strong influence of electrolytes onto the dissolution kinetics through significant enhancement of site protonation probability. Macroscopic dissolution rates calculated as material fluxes span two to three orders of magnitude over the entire pH range, owing to the carbonate protonation mechanism as a function on pH. Among the models for calcite dissolution considered in this study the Oxygen Incongruent model demonstrates the best performance in reproducing the experimentally measured calcite dissolution rate as function of pH and surface morphology. The speciation model used in this study does not rely on fitting experimental data. The outcome of the speciation modelling is determined by the charge distribution of sites on the surface, effective ion size and the intrinsic acidity constant of bicarbonate molecule (an experimental quantity). The current model is developed for dissolution at far from equilibrium conditions and does not take into account the potential contribution of reprecipitation processes. This limitation will be resolved in subsequent work The dissolution kinetics of carbonate minerals as a function of pH has a great relevance to the pH-sensitive carbonate-fluid chemical system involving complex interrelated chemical network between solid and dissolved calcium and carbonate aqueous species, CO2(aq), HCO3-, and atmospheric CO2. According to our findings, presence of electrolytes substantially changes dissolution rates, thus showing the great sensitivity of the carbonate surface to minor changes in pH in the marine and brine-related environments.


32

Page 33 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The developed combined GCMC-KMC modelling approach presented in this study has a great potential in investigation of the role of solvent chemical composition on dissolution mechanisms at the microscopic scale. The examples include but not limited to pH and saturation state dependence, catalytic and inhibitory effects of ions and molecules present in the fluid. This enables the possibility to understand detailed mechanisms governing mineral-fluid interactions in various surface and subsurface environments, i.e., marine and riverine sediments, geological reservoirs and repositories, mine tailings, coral reefs and hydrothermal systems.

ASSOCIATED CONTENT Supporting Information. The Supporting Information is available free of charge. Concentrations of ionized sites on calcite surface for different models, numbers of bulk layers, and surface geometries. Relative concentration of Cl− ions at the calcite-water interface. Electrolyte speciation and concentrations at different pH used in the simulations. Atomic step velocities vs. pH, experimental data and KMC simulations. Effect of the crystal electrostatic field on the surface speciation at different surface topographies. Protonated site fractions from total site populations, GCMC models. Density of protonated sites at calcite surface, KMC model. (PDF) AUTHOR INFORMATION Corresponding Author *e-mail: [email protected]


33


Page 34 of 37

Present Addresses †FB5 and MARUM, University of Bremen, Germany. Klagenfurter Str. 4, 28359 Bremen, Germany. Author Contributions The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Notes The authors declare no competing financial interests. ABBREVIATIONS KMC, Kinetic Monte Carlo; GCMC Grand Canonical Monte Carlo; AFM, Atomic Force Microscopy. REFERENCES (1) (2) (3) (4) (5) (6)

(7)

Arvidson, R. S.; Morse, J. W. Formation and Diagenesis of Carbonate Sediments. In Treatise on Geochemistry (2nd Ed.); Holland, H. D., Turekian, K. K., Eds.; Elsevier: Oxford, 2014; pp 61–101. Morse, J. W.; Arvidson, R. S.; Lüttge, A. Calcium Carbonate Formation and Dissolution. Chem. Rev. 2007, 107 (2), 342–381. Mackenzie, F. T. Carbonate Mineralogy and Geochemistry. In Sedimentology; Encyclopedia of Earth Science; Springer Netherlands, 1978; pp 147–158. Morse, J.W.; Mackenzie, F.T. Geochemistry of Sedimentary Carbonates 1990, in Developments in Sedimentology series, 48. Elsevier Science. Nakamura, T.; Nakamori, T. A. Geochemical Model for Coral Reef Formation. Coral Reefs 2007, 26 (4), 741–755. Huang, S.; Zhang, Y.; Zheng, X.; Zhu, Q.; Shao, G.; Cao, Y.; Chen, X.; Yang, Z.; Bai, X. Types and Characteristics of Carbonate Reservoirs and Their Implication on Hydrocarbon Exploration: A Case Study from the Eastern Tarim Basin, NW China. J. Nat. Gas Geosci. 2017, 2 (1), 73–79. Sanna, A.; Uibu, M.; Caramanna, G.; Kuusik, R.; Maroto-Valer, M. M. A Review of Mineral Carbonation Technologies to Sequester CO2. Chem. Soc. Rev. 2014, 43 (23), 8049– 8080.


34

Page 35 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(8) (9) (10) (11) (12) (13)

(14) (15) (16)

(17) (18) (19)

(20) (21) (22) (23) (24)

Arvidson, R. S.; Ertan, I. E.; Amonette, J. E.; Luttge, A. Variation in Calcite Dissolution Rates:: A Fundamental Problem? Geochim. Cosmochim. Acta 2003, 67 (9), 1623–1634. Colombani, J. The Alkaline Dissolution Rate of Calcite. J. Phys. Chem. Lett. 2016, 7 (13), 2376–2380. Lüttge, A.; Arvidson, R. S.; Fischer, C. A Stochastic Treatment of Crystal Dissolution Kinetics. ELEMENTS 2013, 9 (3), 183–188. Fischer, C.; Arvidson, R. S.; Lüttge, A. How Predictable Are Dissolution Rates of Crystalline Material? Geochim. Cosmochim. Acta 2012, 98, 177–185. Kurganskaya, I.; Luttge, A. Kinetic Monte Carlo Approach To Study Carbonate Dissolution. J. Phys. Chem. C 2016, 120 (12), 6482–6492. Heberling, F.; Bosbach, D.; Eckhardt, J.-D.; Fischer, U.; Glowacky, J.; Haist, M.; Kramar, U.; Loos, S.; Müller, H. S.; Neumann, T.; et al. Reactivity of the Calcite-Water-Interface, from Molecular Scale Processes to Geochemical Engineering. Appl. Geochem. 2014, 45, 158–190. Van Cappellen, P.; Charlet, L.; Stumm, W.; Wersin, P. A Surface Complexation Model of the Carbonate Mineral-Aqueous Solution Interface. Geochim. Cosmochim. Acta 1993, 57 (15), 3505–3518. Duckworth, O. W.; Martin, S. T. Connections between Surface Complexation and Geometric Models of Mineral Dissolution Investigated for Rhodochrosite. Geochim. Cosmochim. Acta 2003, 67 (10), 1787–1801. Jordan, G.; Higgins, S. R.; Eggleston, C. M.; Knauss, K. G.; Schmahl, W. W. Dissolution Kinetics of Magnesite in Acidic Aqueous Solution, a Hydrothermal Atomic Force Microscopy (HAFM) Study: Step Orientation and Kink Dynamics. Geochim. Cosmochim. Acta 2001, 65 (23), 4257–4266. A. Voter. Introduction to the Kinetic Monte Carlo Method. In Radiation Effects in Solids (ed. K.E. Sickafus and E.A. Kotomin); NATO Science; 2007; Vol. 235, pp 1–23. Wehrli, B. Monte Carlo Simulations of Surface Morphologies during Mineral Dissolution. J. Colloid Interface Sci. 1989, 132 (1), 230–242. Wolthers, M.; Charlet, L.; Cappellen, P. V. The Surface Chemistry of Divalent Metal Carbonate Minerals; a Critical Assessment of Surface Charge and Potential Data Using the Charge Distribution Multi-Site Ion Complexation Model. Am. J. Sci. 2008, 308 (8), 905– 941. Wolthers, M.; Di Tommaso, D.; Du, Z.; de Leeuw, N. H. Calcite Surface Structure and Reactivity: Molecular Dynamics Simulations and Macroscopic Surface Modelling of the Calcite-Water Interface. Phys. Chem. Chem. Phys. 2012, 14 (43), 15145–15157. Andersson, M. P.; Stipp, S. L. S. How Acidic Is Water on Calcite? J. Phys. Chem. C 2012, 116 (35), 18779–18787. Andersson, M. P.; Rodriguez-Blanco, J. D.; Stipp, S. L. S. Is Bicarbonate Stable in and on the Calcite Surface? Geochim. Cosmochim. Acta 2016, 176 (Supplement C), 198–205. Labbez, C.; Jönsson, B.; Skarba, M.; Borkovec, M. Ion−Ion Correlation and Charge Reversal at Titrating Solid Interfaces. Langmuir 2009, 25 (13), 7209–7213. Porus, M.; Labbez, C.; Maroni, P.; Borkovec, M. Adsorption of Monovalent and Divalent Cations on Planar Water-Silica Interfaces Studied by Optical Reflectivity and Monte Carlo Simulations. J. Chem. Phys. 2011, 135 (6), 064701.


35


Page 36 of 37

(25) Churakov, S. V.; Labbez, C.; Pegado, L.; Sulpizi, M. Intrinsic Acidity of Surface Sites in Calcium Silicate Hydrates and Its Implication to Their Electrokinetic Properties. J. Phys. Chem. C 2014, 118 (22), 11752–11762. (26) Labbez, C.; Jönsson, B. A New Monte Carlo Method for the Titration of Molecules and Minerals. In Applied Parallel Computing. State of the Art in Scientific Computing; Lecture Notes in Computer Science; Springer, Berlin, Heidelberg, 2006; pp 66–72. (27) Valleau, J. P.; Cohen, L. K. Primitive Model Electrolytes. I. Grand Canonical Monte Carlo Computations. J. Chem. Phys. 1980, 72 (11), 5935–5941. (28) Shiraki, R.; Rock, P. A.; Casey, W. H. Dissolution Kinetics of Calcite in 0.1 M NaCl Solution at Room Temperature: An Atomic Force Microscopic (AFM) Study. Aquat. Geochem. 2000, 6 (1), 87–108. (29) Vibhu, I.; Modak, B.; Patra, C. N.; Ghosh, S. K. Zeta Potential of Colloidal Particle in Solvent Primitive Model Electrolyte Solution: A Density Functional Theory Study. Mol. Phys. 2013, 111 (4), 489–496. (30) Henry, C.; Minier, J.-P.; Pozorski, J.; Lefèvre, G. A New Stochastic Approach for the Simulation of Agglomeration between Colloidal Particles. Langmuir 2013, 29 (45), 13694– 13707. (31) Modak, B.; Patra, C. N.; Ghosh, S. K. Solvent Primitive Model Study of Structure of Colloidal Solution in Highly Charge Asymmetric Electrolytes. AIP Conference Proceedings 2013, 1512 (1), 606–607. (32) Abbas, Z.; Ahlberg, E.; Nordholm, S. Monte Carlo Simulations of Salt Solutions: Exploring the Validity of Primitive Models. J. Phys. Chem. B 2009, 113 (17), 5905–5916. (33) Churakov, S. V.; Labbez, C. Thermodynamics and Molecular Mechanism of Al Incorporation in Calcium Silicate Hydrates. J. Phys. Chem. C 2017, 121 (8), 4412–4419. (34) Understanding Molecular Simulation; Elsevier, 2002. (35) Torrie, G. M.; Valleau, J. P. Electrical Double Layers. I. Monte Carlo Study of a Uniformly Charged Surface. J. Chem. Phys. 1980, 73 (11), 5807–5816. (36) Widom, B. Some Topics in the Theory of Fluids. J. Chem. Phys. 1963, 39 (11), 2808–2812. (37) Svensson, B. R.; Woodward, C. E. Widom’s Method for Uniform and Non-Uniform Electrolyte Solutions. Mol. Phys. 1988, 64, 247–259. (38) Giudici, G. D. Surface Control vs. Diffusion Control during Calcite Dissolution: Dependence of Step-Edge Velocity upon Solution PH. Am. Mineral. 2002, 87 (10), 1279– 1285. (39) Liang, Y.; Baer, D. R. Anisotropic Dissolution at the CaCO3(101M4) — Water Interface. Surf. Sci. 1997, 373 (2–3), 275–287. (40) Liang, Y.; Lea, A. S.; Baer, D. R.; Engelhard, M. H. Structure of the Cleaved CaCO3(101M4) Surface in an Aqueous Environment. Surf. Sci. 1996, 351 (1–3), 172–182. (41) Liang, Y.; Baer, D. R.; McCoy, J. M.; Lafemina, J. P. Interplay between Step Velocity and Morphology during the Dissolution of CaCO3 Surface. J. Vac. Sci. Technol. A 1996, 14 (3), 1368–1375. (42) Chou, L.; Garrels, R. M.; Wollast, R. Comparative Study of the Kinetics and Mechanisms of Dissolution of Carbonate Minerals. Chem. Geol. 1989, 78 (3–4), 269–282. (43) Wang, X.; Conway, W.; Burns, R.; McCann, N.; Maeder, M. Comprehensive Study of the Hydration and Dehydration Reactions of Carbon Dioxide in Aqueous Solution. J. Phys. Chem. A 2010, 114 (4), 1734–1740.


36

Page 37 of 37 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(44) Pokrovsky, O. S.; Schott, J. Surface Chemistry and Dissolution Kinetics of Divalent Metal Carbonates. Environ. Sci. Technol. 2002, 36 (3), 426–432. (45) Atanassova, R.; Cama, J.; Soler, J. M.; Offeddu, F. G.; Queralt, I.; Casanova, I. Calcite Interaction with Acidic Sulphate Solutions: A Vertical Scanning Interferometry and EnergyDispersive XRF Study. Eur. J. Mineral. 2013, 25 (3), 331–351. (46) Arvidson, R. S.; Collier, M.; Davis, K. J.; Vinson, M. D.; Amonette, J. E.; Luttge, A. Magnesium Inhibition of Calcite Dissolution Kinetics. Geochim. Cosmochim. Acta 2006, 70 (3), 583–594. (47) Vinson, M. D.; Arvidson, R. S.; Luttge, A. Kinetic Inhibition of Calcite (1 0 4) Dissolution by Aqueous Manganese(II). J. Cryst. Growth 2007, 307 (1), 116–125. (48) Ghoshal, S.; Hazra, M. K. H2CO3 → CO2 + H2O Decomposition in the Presence of H2O, HCOOH, CH3COOH, H2SO4 and HO2 Radical: Instability of the Gas-Phase H2CO3 Molecule in the Troposphere and Lower Stratosphere. RSC Adv. 2015, 5 (23), 17623–17635. (49) Urosevic, M.; Rodriguez-Navarro, C.; Putnis, C. V.; Cardell, C.; Putnis, A.; Ruiz-Agudo, E. In Situ Nanoscale Observations of the Dissolution of Dolomite Cleavage Surfaces. Geochim. Cosmochim. Acta 2012, 80, 1–13. (50) Ruiz-Agudo, E.; Kowacz, M.; Putnis, C. V.; Putnis, A. The Role of Background Electrolytes on the Kinetics and Mechanism of Calcite Dissolution. Geochim. Cosmochim. Acta 2010, 74 (4), 1256–1267. (51) Finneran, D. W.; Morse, J. W. Calcite Dissolution Kinetics in Saline Waters. Chem. Geol. 2009, 268 (1–2), 137–146.

TOC


37

Carbonates Dissolution Mechanisms in the Presence of Electrolytes

Recommend Documents