Kinetic Monte Carlo Approach To Study Carbonate Dissolution - The

Subscriber access provided by LOCK HAVEN UNIV

Article

A Kinetic Monte Carlo Approach to Study Carbonate Dissolution Inna Kurganskaya, and Andreas Luttge J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.5b10995 • Publication Date (Web): 29 Feb 2016 Downloaded from http://pubs.acs.org on March 5, 2016

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry C 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 33

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

A Kinetic Monte Carlo Approach to Study Carbonate Dissolution Inna Kurganskaya1,* and Andreas Luttge2,3 1

Institue for Geological Sciences, University of Bern, Baltzerstrasse 3, 3012 Bern, Switzerland 2

University of Bremen, MARUM and FB5, PO Box 330 440, 28334 Bremen, Germany 3

Rice University, Houston, TX 77005, USA

ABSTRACT Reactive properties of carbonate minerals and rocks attract a significant attention with regard to modern environmental and industrial problems, including geologic carbon sequestration, toxic waste utilization, cement clinker production, the fate of carbonate shellbearing organisms in acidifying oceans. Despite the ultimate importance of the problem, the number of studies connecting atomistic scale models to experimental observations is limited. In this work we employ Kinetic Monte Carlo (KMC) approach to model carbonate dissolution at the nanometer-micron scale range and to access quantitative relationships between the variety of surface reactive sites and experimentally observed spatio-temporal rate variance. Presented KMC model adequately reproduces experimentally observed dissolution patterns and rates. We examine mechanistic and quantitative relationships between site reactivity, atomic step velocity, etch pit evolution dynamics and macroscopic dissolution rates. The analysis of reactive site statistics shows the leading role of kink sites in the control of the overall rate. The kink site

ACS Paragon Plus Environment

1


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 33

density is closely tied to the controlling types of surface features, e.g., straight and curved steps having different structural orientations, as well as their dynamic interaction. Structurally different kink sites and steps bring different contributions to the total rate. The modeling results indicate that rate values depend on the spatial distribution and type of lattice defects. Greatly inhomogeneous defect distribution commonly observed in natural minerals leads to the orders of magnitude local rate variance. The reason is in the dependence of reactive site density on the local availability of reactive steps sources, e.g. screw dislocations. We conclude that upscaling of the atomistic models of carbonate dissolution and macroscopic rate predictions should be done taking in account the entire variety of reactive sites and surface features, as well as their spatiotemporal distribution.

1. INTRODUCTION Dissolution kinetics of carbonate minerals and their surface reactivity has been a subject of intensive studies with regard to industrial and environmental applications. Carbonates serve as a primary material source for cement clinker production. These minerals are widely used in pharmacy, paint and coating industry. Carbonate rocks contain a great amount of pores and cavities that can be used to trap waste gases and fluids. Geologic sequestration of industrial waste requires an accurate estimate of mineral reactivity changes induced by new geochemical environment. A global response of oceanic system to climate change and subsequent water acidification is closely related to carbonate dissolution/precipitation rates at changing geochemical conditions. Dissolution kinetics of carbonate minerals has been studied over the decades by using a great variety of tools and methodological approaches. The first studies were aimed to measure dissolution rates of fine grain-size material and their dependence on solvent saturation state, pH


2

Page 3 of 33

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


and salinity1–11. In these studies the reacting solid was treated as a uniform phase dissociating at a constant rate depending on solvent composition. The entirely new concept involving surface sites (Surface Speciation Model, SSM12) has been proposed to construct the surface model and explain its reactive properties at different pH conditions. The model implies a uniformly flat ideal solid substrate, while real surfaces of dissolving carbonate minerals have complex geometry outlined by etch pits, curved, straight and macro steps13,14–20. The collective action of the reactive surface features is responsible for the mineral dissolution process, where the undercoordinated sites, e.g. step and kink sites21,22, play a key role in defining reaction kinetics. First Atomic Force Microscope (AFM) observations of calcite surface showed the existence of moving atomic steps14,23. The dissolution model derived from these studies was based on the step movement, considering kink site formation, propagation, and “annihilation”24. This mechanistic model became a conceptual basis for the next generation of studies. Thus, Duckworth and Martin25 combined both structural and SSM approaches into a new reactivity model considering concentrations of reactive step sites. Wolthers et al.26,27 further developed structural SSM approach by incorporating all types of structurally different kink sites. Bracco et al.28,29 and Stack30 inherited the same step-based approach to derive macroscopic rates as a function of solvent saturation state and step density. The atomic scale details regarding mineral-water interface structure, water and solvent ion adsorption became available from Molecular Dynamics calculations27,31–42 and experimental Xray reflectivity studies43–47. This information has an ultimate importance for the understanding of the reaction kinetics at the atomic scale, where a great variety of reactions is possible. Thus, MD calculations showed the substantial difference in dissolution enthalpies for structurally different


3


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 4 of 33

sites33,38 and thus confirmed the assumption about different reactivity directions – for obtuse and acute steps and kinks – initially made by Liang et al.48,49. Studies of the dissolution mechanisms at the larger spatial scales involve kinematics and interaction of bigger surface features, e.g. etch pits and macrosteps. McInnis and Brantley15,50 analyzed etch pit distribution (EPD) at the reacting calcite surface and offered a mechanistic dissolution model based on etch pit growth and coalescence. Duckworth and Martin51 measured microscopic rates based on mass removal from etch pits and compared the results to the macroscopic rates. Although these methods provided an insight into reaction mechanisms at the microscopic scale, they missed a very important rate component – surface normal retreat defined as a rate of a mean height level drop52–54. This process is attributed to material removal by atomic steps at the uppermost surface features. Generally, it may happen without any detectable change of surface geometry, e.g. etch pit deepening and widening. Therefore, while the etch pit volume does not substantially change, the surface may still actively dissolve layer by layer from the top. Luttge and co-workers applied the stepwave model to explain mechanistic relationships between atomic step movement, surface normal retreat, etch pit dynamics, and their effects on macroscopic dissolution rates55–59. They assumed that etch pits work only as step sources, and once a steady-state is achieved, the distribution and density of etch pits are meaningless for the overall rate value. Practically, however, this ideal steady-state situation is problematic to ensure and measure – thus, according to Arvidson et al.52, the rate values obtained from the different laboratories may vary up to two orders of magnitude. Recently, Fischer et al.54 and Luttge et al.60 introduced a hypothesis about an intrinsic variance of mineral dissolution rates that practically cannot be reduced to a single rate constant value. Instead, they suggested using a rate spectra approach, according to which a possible range of


4

Page 5 of 33

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


rates, or spectrum, should be measured or estimated to characterize mineral’s reactivity at given environmental conditions. Fischer et al.54 systematically measured surface roughness and dissolution rate over time for a series of selected areas at calcite single crystal and micrite surfaces. They noticed that depending on the type of the system studied (single crystal vs. micrograin assemblage), rate spectra may show different width and shape, e.g. micrograin system showed wider rate distribution due to the variance in grain size. The observed rate dependence on the system choice shows a major failure of the simple microscopic models to explain kinetic behavior of the macroscopic systems. We suggest that this problem reflects at least in part the limited knowledge concerning the mechanistic relationships between reactivity and scales of length and time. The inter-scale connections, at least at the nm-micron range, can be understood by using Kinetic Monte Carlo (KMC) modeling of dissolution61 and growth62 processes. This approach allows ones to observe dynamics of surface topography change, characterize dissolution rates and their distribution over the surface, and compare simulation results to experimental observations. McCoy and LaFemina63 developed the first well-recognized KMC model of calcite dissolution. This model was parameterized by the experimental AFM data for obtuse and acute step velocities64. One of the important outcomes was the estimate of dissolution rates from different kink and step sites. The results of these simulations were interpreted in terms of a mechanistic model of the straight step movement24. Later, Williford et al.65 updated this model by including reactions of back-precipitation, surface diffusion and diffusion in the boundary layer, in order to test the influence of the fluid flow rate on step velocities. These early models were very successful in explaining experimentally measured step velocities. However, the correct kinetic description of the entire reactive system requires consideration of all kinetically


5


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 6 of 33

important processes and reactive surface features. Jordan and Rammensee66 observed curved steps on carbonate surfaces – these steps move much faster than straight steps and are responsible for most of the material flux from the surface66–68. Curved steps form as a result of the straight step coalescence, that takes place due to the presence of multiple step sources – etch pits. The resulting cumulative etch pit and step action can be measured as normal surface retreat56,58 that gives a precise estimate of the material flux from a given surface. Therefore, the full predictive potential of a KMC model can be achieved if it would be able to reproduce these kinetic events. Although carbonate dissolution had been studied in great details at the microscopic scale, the details about reactive site distributions and their contributions to the overall rate are still not well understood. In the present study we show how a KMC approach can be used to explain mechanistic and quantitative relationships between elementary surface reactions represented by dissolution from structurally different sites, microscopically observed step movement, etch pit interaction, and macroscopically measured dissolution rates. The principal novelty is the use of screw dislocations and 3D crystal model, enabling formation of 3D etch pits and macrosteps. Our goal is to combine all the findings and concepts mentioned above into a one comprehensive approach that potentially should be able to predict dissolution rates for the realistic complex systems studied at both experimental laboratory and field study settings.


6

Page 7 of 33

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. 1. Schematic diagram of basic reactive features at calcite cleavage (104) surface: obtuse (So) and acute (Sa) steps, four types of kink sites, kaa, kao, koa and koo, where the first letter indicates step orientation, the second letter indicates the orientation of a kink within the step. Crystallographic orientations here and thereafter in the paper are given for the hexagonal unit cell. 2. METHODS 2.1. System modeled Carbonate minerals have ionic crystal structure, where [CO3]2- and M2+ ions (M2+=Ca2+, Mg2+, Sr2+, Fe2+, Zn2+, etc.) are the main building blocks. M2+ ions are coordinated by six oxygen atoms belonging to the neighboring carbonate groups. Each carbonate group is coordinated by six M2+ ions in the bulk structure. The bonding network topology between M2+ and carbonate groups is the same as for Kossel crystal21,22 – a generalized crystal model where atomic centers are represented by cubic blocks which centers form simple orthogonal 3D grid. In this crystal type the number of bonds n=6 characterize bulk positions, n=5 terrace, n=4 step and n=3 kink positions. In the case of calcite, each reactive center, either Ca or [CO3]2-, is bonded to the


7


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 8 of 33

surface by n (n=1-5) Ca-O bonds which must be broken in order to release [Ca·(H2O)6]2+ or [CO3]2- into the solvent. Carbonate crystals have rhombohedral cleavage lying in three non-parallel planes forming a rhombohedron. The planes intersect with each other at an inclined angle 101.92⁰, resulting in non-equivalent (obtuse and acute) orientations of exposed cleavage faces. Atomic steps adjacent to these faces, as a result, may attain either obtuse or acute orientation. Depending on step cleavage plane orientation, steps and step sites are classified as obtuse and acute, so and sa. Kink sites formed at each of these steps also can take any of these two orientations in addition to step orientation. Thus, two step sites, sa and so, and four kink sites, kaa, kao, koa and koo are distinguished for carbonate minerals having calcite-type crystal structure (Fig. 1). The orientation of steps and kinks has a substantial influence on their reactivity, reflected in remarkable difference for obtuse and acute step velocities measured in numerous experiments48,49,63,66,67. In present work we modeled dissolution of calcite, although the method can be easily extended for the other carbonates by appropriate parameterization. 2.2. KMC model A KMC algorithm schedules a stochastically defined order of the surface reaction occurrence. An input to a KMC model is thus created as a predefined list of possible reactions with their rates69. The reaction list is constructed for the dissolution events taking place at surface sites which reactivity is defined by the numbers of nearest neighbors. The influence of obtuse and acute orientations can be implemented by assigning difference in bond breaking rates in “fast” (obtuse) and “slow” (acute) directions63,65. In this study we followed the modeling scheme previously designed by McCoy and LaFemina63. We developed an in-house FORTRAN code from scratch where we added three


8

Page 9 of 33

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


new components: (i) direct use of Ca-CO3 bond network instead of CaCO3 block units; (ii) the ability to implement lattice defects – screw dislocations and point defects- into the system; (iii) long-term runs for 3D systems with lateral periodic boundary conditions. The use of 3D instead of common 2D arrays (so-called solid-on-solid model61) enables formation of extra steps at etch pit walls and macrosteps. Visual Molecular Dynamics (VMD 1.9.1) software was used to visualize surface topography obtained from the simulations. We simulated two dislocation types: point defects and screw dislocations. Generally, dislocations can be introduced into a KMC model by defining dislocation parameters, such as core radius and strain energy70,71 affecting dissolution probabilities for selected surface areas. However, Meakin and Rosso71 noted that the same simulation results can be obtained for the models explicitly considering strain field at core regions and those simply incorporating either rapidly dissolving or pre-opened nanotubes. We implemented screw dislocations by pre-opening of the nanotubes having fixed width (2x2 CaCO3 units) and length. The length was taken long enough to ensure the steady-state – 200-300 CaCO3 units. The size of the entire system was 1000x1000x2 CaCO3 units for step velocity measurements and 600x600x300 CaCO3 units for structures involving screw dislocations. In order to simplify the bookkeeping procedures, we used morphological unit cell constructed to form cleavage planes parallel to basic crystallographic planes (see Supporting Information). All simulations were started from the atomically flat cleavage faces lying in (a-b) planes with the hollow cores opened along c axis. We applied periodic boundary conditions (PBC) in lateral (x-y or a-b) directions in order to eliminate the influence of edges and diminish the influence of the system size. The reaction rates for Ca2+ and [CO3]2- units connected via i Ca-O bonds to the surface are defined as follows:


9


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 10 of 33

(1).

Here is dissolution rate for Ca2+ or [CO3]2- attached to the surface by bonds, is the number of reaction attempts, “fundamental frequency” according to Lasaga and Luttge’s terminology58,72,73 or Arrhenius pre-factor, ∆ – is the activation energy of Ca-O-C bond hydrolysis. Anisotropy of bond breaking in obtuse vs acute directions is taken into account by using

activation energy corrections ∆

for obtuse and acute step site orientations (i=4) and four

kink site orientations (i=3). This energy correction is also used as an additional stabilizing term for the terrace sites (i=5). Here j denotes specific orientation for the step site - “a” (acute) or “o” (obtuse) and for the kink site – “aa”, “ao”, “oa” and “oo” (see Fig.1 for details). If i =5 (terrace ! sites), then j is denoted as “T” and ∆ ! stabilization effects (∆

can be used as a term taking in account cleavage face

= 0 if this effect is turned off).

The values of the activation energy parameters for breaking of different bonds in principle can be obtained from Molecular Dynamics or ab initio studies. However, a consistent data set for the reaction rates or activation energies is not available in the moment due to the two primary reasons: (i) a substantial computational cost for estimation of rare event rates74; (ii) dependence of the obtained results from parameterized force fields. Thus, Raiteri et al.39 showed that different parameterization approaches may lead to different results for dissolution free energies. Particularly, they found that the free energy for dissolution of calcium and carbonate ions do not substantially differ from each other, while according to Spagnoli et al.38 free energy of calcium ion dissolution is about 10 times greater comparing to carbonate ion. Despite the problem of exact activation energies estimation, the energetic parameters provided, e.g. dissolution enthalpy and free energy, showed important qualitative information regarding


10

Page 11 of 33

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


site reactivity. Thus, based on the reaction enthalpy, double kink formation and propagation at obtuse edges is more energetically favorable then at acute edges33. Free energy difference calculations for dissolution of ions showed that dissolution from steps is much more favorable than dissolution from terrace sites. We varied modeling parameters in Eq. (1) until the model reproduced experimentally48,49 observed step velocities and vo/va ratios (Table 1) at far-from-equilibrium conditions. The scheme for the step velocity measurements is shown on Fig. 2.

Fig. 2. Measurements of straight step velocities on calcite surface, KMC simulations. A) a monolayer etch pit originated at the cleavage plane center; B) step velocity measurements at different time steps (see Table 1, set C, for simulation parameters). This parameterization technique had been previously used by McCoy and LaFemina63 for calcite. In our first parameterization trial we used their value of fundamental frequency in set A, Table 1. Step separation value L obtained using this set is ~1-5 nm. It is quite small comparing to experimental values75 — 10-50 nm. The value of L can be increased if we increase activation energy of bond breaking and “fundamental frequency” υ. For example, frequency 1012 s-1 corresponds to the water vibration frequency and can be used to approximate frequency of bond hydrolysis reaction attempts76. Thus, parameter set B can produce greater values of step separation—1-13 nm (Table 1).


11


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 12 of 33

An example of obtuse and acute step velocities measured as a function of time is shown on Fig. 2B — these velocities reached constant values 1.2 and 3.1 nm/s for va and vo correspondingly at the step lengths exceeding 100 nm (set B, Table 1). Table 1. Parameter sets used in KMC modeling of calcite dissolution. set

υ, s-1

∆Eb, kT units

∆Esa; ∆Eso, kT units

∆Ek_aa ;∆Ek_ao; , ∆Ek_oa; ∆Ek_oo, nm/s kT units

, nm/s

La, nm

Lo, nm

[Lo], nm

A

5.22·1010

6.35

1.8;0

1.8;1.5;1.5;0

0.9

3.3

0-3

1-5

2.5

B

1012

7.3

1.5;0

1;0.5;0.5;0

1.2

3.1

0-5

1-13

5

C

5.22·1010

6.45

3;0

3;1.5;1.5;0

0.23

2.3

0-1

1-10

4

D

1012

7.3

3.5;0

2.5;1.5;1.5;0

0.2

2.0

0-2

1-11

5

υ-“fundamental frequency”, e.g. frequency of the reaction attempts, ∆Eb– activation energy for Ca-O bond hydrolysis, ∆Esa; ∆Eso –energy corrections for acute and obtuse steps, ∆Ek_aa; ∆Ek_ao; ∆Ek_oa; ∆Ek_oo – energy corrections for four kinks site types – acute and obtuse kinks at acute step; acute and obtuse kinks at obtuse step; -acute step velocity; -obtuse step velocity, Listep separation, [Lo] – mean value for obtuse steps.

Free energy difference calculations for ion dissolution reactions obtained by Raiteri et al.39 give values ~100-300 kJ/mole for dissolution of an ion belonging to the terrace site type. This may roughly correspond to the values ~20-60 kJ/mole (8-24 kT units) for a one bond breaking. Our activation energy values fall into the lowest values in this range. The use of the values > 10 kT units results in either unrealistic values of the fundamental frequency factor (>1015) or too slow step velocities.

3. RESULTS AND DISCUSSION 3.1. Surface topography: general features


12

Page 13 of 33

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


3.1.1. Monolayer pits A monolayer rhomb-shaped pit is an elementary feature frequently observed on the dissolving surface of carbonate minerals16,19,77,78 (Fig. 3A). Monolayer pits can grow around monoatomic holes formed at point defects or as a result of unassisted terrace sites dissolution. The last mechanism does not require any lattice defects to be present. Thus, according to Teng77, a spontaneous formation of monolayer pits can be observed at certain undersaturation state. The same effect can be obtained in KMC simulations, depending on the activation energies for terrace site dissolution. For example, if we use parameter set A (Table 1) without adding any ! stabilizing terms to the terrace sites (e.g., ∆

= 0), the monolayer etch pits will be

spontaneously created (Fig. 3B). As a result, a simple layer-by-layer dissolution mechanism can be observed19.

Fig. 3. Monolayer pits formed on surfaces of carbonate minerals. A. Dolomite, AFM observations19. Reprinted from Geochimica et Cosmochimica Acta, Vol. 80, Urosevich et al. “In situ nanoscale observation of the dissolution of {101$4} dolomite cleavage surfaces”, pp. 1-13, Copyright (2012) with permission from Elsevier; B. Calcite, KMC simulations. 3.1.2. Multilayer pits


13


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 14 of 33

Fig. 4. Morphology of a multilayer etch pit formed on calcite cleavage face, pure water at farfrom-equilibrium conditions. A. AFM data79. Reprinted from Journal of Crystal Growth, Vol. 307 (1), Vinson et al. “Kinetic inhibition of calcite (104) dissolution by aqueous manganese(II)”, pp. 116-125, Copyright (2007) with permission from Elsevier; B,C. KMC simulations, B. parameter set A, C. parameter set B (Table 1). Multilayer pits (Fig. 4A) form around screw dislocation cores that open in the very beginning of the dissolution process. A dislocation hollow core is necessary to start etch pit formation and support localized step and stepwave production. A pit grows as a local region of reactive sites connected via close genetic relationships. It experiences different stages of its existence, as point-bottomed and flat-bottomed pits50. A point-bottomed pit exists as soon as the pit-forming dislocation exists and supplies the retreating pit by the steps. A flat-bottomed pit appears when the dislocation in the lattice closes up and bottom becomes a stable atomically flat plane. Lateral shape of the multilayer rhombohedral pits is defined by terraced straight steps oriented parallel to the cleavage planes. Straight step orientation and overall etch pit’s shape largely depend on activation energies of bond breaking. We observed that the rhombohedral orientations as shown in Fig. 3 and Fig. 4 appear at values ∆ECa-O > 5.5 kT. Pits have anisotropic shape due to the difference in terrace width for obtuse and acute steps - walls made by acute steps are always steeper than walls made by obtuse steps. Terrace width values depend on the parameter set used.


14

Page 15 of 33

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


Thus, narrower terraces and steeper slopes are observed for pits obtained by using parameter set A, Table 1 (Fig. 4B), in comparison to terraces and slopes of the pits obtained by using parameter set B (Fig. 4C). Please, note that regardless the substantial size difference between the pits observed in an experiment (Fig. 4A) and in the simulations (Figs. 4B-C), the overall pit’s morphology is preserved while pit grows bigger in size provided by enough free space and dislocation length to develop. 3.1.3. Pits interaction and surface evolution Dissolution mechanisms involving multilayer pits and screw dislocations are quite complex. Fixed positions of hollow cores create a situation where steps originate at the same points during a prolonged period of time, thus supporting a “steady-state” (Fig. 5). Long-term interaction of the straight steps results in formation of actively dissolving areas covered by the curved steps. These steps are formed as a result of the straight step coalescence at the etch pit’s corners. The reactive corners created in this process progressively grow in length and form rounded steps (Fig. 5A). The curved steps travel from the points of the corner-corner pit intersections to the corners of an outer dominating pit, where they terminate. In this process of the step interaction the bigger pits outlined by the straight steps are formed. The effect of the step bunching – accumulation of the steps due to their slow velocity- creates another interesting effect – formation of the additional surface area and terraces parallel to the other cleavage planes (Fig. 5B). Interaction of the pits possessing this structure results in formation of macro-steps and ridges as well as “staircase” stepped slopes on the mineral’s surface (Fig. 5C). Curved steps also accumulate due to the step bunching process resulting in formation of stepped “ridges” and “capes” (Fig. 5D).


15


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 16 of 33

Fig. 5. Calcite surface topography, KMC simulations (see parameter values in Table 1), crystals with screw dislocations. A. Coalescence of three etch pits and formation of a one shallow pit (set B). A red circle with red arrows show corner-corner interaction of the straight steps, formation of a reactive corner that progressively develop into a curved step; B. step bunching for the system with va/vo ~ 10 and formation of additional terraces at the pit’s walls; C. Formation of macrosteps as a result of step velocity anisotropy and step bunching (set C); D. Bunching of rough and straight (see inset) steps. 3.2. Straight step dynamics Detachment of molecules from kink sites and overall step propagation cause material flux from the surface. At experimental laboratory settings this flux can be detected as increase of the element’s concentration in reactive solvent. In the beginning of the dissolution process, when steps start to form at the dislocation outcrop and move outwards, all kinetic activity is localized at the dislocation centers. A single straight step is an elementary kinetic unit responsible for dissolution at this stage. Total dissolution rate (material flux) shows temporal linear growth (Fig.


16

Page 17 of 33

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


6A) interfered by fluctuations due to the fluctuations of kink site densities along the step (Fig. 6C). This overall rate dynamics perfectly correlates to the stable linear growth of the step site numbers over time (Fig. 6B). Amounts of obtuse and acute step sites are almost equal, because of the pit’s rhombohedral morphology (Fig. 6A). Linear temporal growth of the step site amounts, so and sa, indicates constant step velocity. On contrast, the growth of the kink site amount is much less stable due to quite large fluctuation amplitude (about an order of magnitude). However, the same linear trend, related to etch pit expansion and linear growth of the step length, is observed (Fig. 6C). This trend is quite expectable, since longer steps may bear greater amount of the kink sites. If we normalize kink site number by the step site number for the corresponding kink and step sites (e.g. kao and kaa types normalized by the number of acute sa step sites, and koo and kao by the number of obtuse so sites), we can see that the normalized values fluctuate around some constant value (Fig. 6D). These small scale fluctuations are reflected in the fluctuations of the total rate (Fig. 6A). A quantitative comparison of the kink site concentrations (Fig. 6C-D) clearly shows that koa kink sites are the most abundant. This does not mean that these sites determine the overall dissolution kinetics. Instead, they are preserved in a greater amount as less reactive features, in contrast to koo kink types. Fig. 6E demonstrates contributions of all kink site types into the overall reactive flux: site types kao and koo give 35-40% of the total flux each, site types kaa and kao give 9-15% each. Therefore, regardless the same step lengths for obtuse and acute steps, obtuse steps contribute 70-80% to the total material flux. This result also in a good agreement to general expectations for the contributions of these two types of steps based on vo/va ratio: faster steps remove material much more effectively.


17


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 18 of 33

Fig. 6. Kinetic parameters for the growing monolayer pit shown in Fig. 2A, KMC simulations [abbreviations: a-acute, o-obtuse orientations]. A. Dissolution rate (material flux) from the simulated area vs time. B. % of step sites from total surface sites (step site density). C. % of kink sites from total surface sites (kink site density). D. % of kink sites from the corresponding step sites. E. Contributions of different kink sites to the total reactive flux from the surface. 3.3. Curved step dynamics Jordan and Rammensee66 recognized two kinds of steps on calcite surface: straight “slow” (also referred as “acute” and “obtuse”) and curved “fast”. Curved steps bear quite large amount of kink sites, which make them remarkably reactive features. Figs. 7A-C illustrates the difference in reactivity between straight steps outlining rhomb-shaped pits and curved steps80. While rough steps remove material from the surface, straight step remain relatively immobile (Figs. 7A-C). Jordan et al.67 recognized three types of “fast” steps, depending on the orientations of coalescing straight steps. Thus, coalescence of obtuse steps results in formation of Goo, the


18

Page 19 of 33

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


fastest curved step (Fig. 7D). Two coalescing acute steps give the slowest step in this rough step “family”: Gaa. Step, formed as a result of the acute and obtuse step coalescence, Goa or Gao, has an intermediate velocity. We measured curved step velocities as a distance traveled by these steps divided by the corresponding time interval. The obtained values for Gaa = 7-10 nm/s, Gao=Goa=5-25 nm/s, Goo=25-30 nm/s. Fig. 7E shows a height difference map for the time stepspositive light features there indicate distances traveled by curved steps. The velocity of these steps is not a constant value and depends on the step curvature and orientations of annihilating steps. Obviously, as predicted by Jordan and Rammensee66, v(Gaa)

Kinetic Monte Carlo Approach To Study Carbonate Dissolution - The

Recommend Documents