A Kinetic Monte Carlo Approach to Model Barite Dissolution: The Role of Reactive Site Geometry

Abstract

Barite (Ba[SO₄]) is one of the promising candidates for sequestration of radioactive waste. Barite can incorporate radium (Ra) and form ideal solid solutions, i.e., (Ba,Ra)[SO₄]. Together with isostructural celestite (Sr[SO₄]), ternary solid solutions, (Ba,Sr,Ra)[SO₄], may exist in natural conditions. Our fundamental understanding of the dissolution kinetics of isostructural sulfates is critically important for a better risk assessment of nuclear waste repositories utilizing this mineral for sequestration. So far, the barite-water interface has been studied with experimental methods and atomistic computer simulations. The direct connection between the molecular scale details of the interface structure and experimental observations at the microscopic scale is not yet well understood. Here, we began to investigate this connection by using a kinetic Monte Carlo approach to simulate the barite dissolution process. We constructed a microkinetic model for the dissolution process and identified the reactive sites. Identification of these sites is important for an improved understanding of the dissolution, adsorption, and crystal growth mechanisms at the barite–water interface. We parameterized the molecular detachment rates by using the experimentally observed etch pit morphologies and atomic step velocities. Our parameterization attempts demonstrated that local lattice coordination is not sufficient to differentiate between the kinetically important sites and estimate their detachment rates. We suggest that the water structure and dynamics at identified sites should substantially influence the detachment rates. However, it will require more work to improve the parameterization of the model by means of Molecular Dynamics and ab initio calculations.

1. Introduction

Barite (BaSO₄) plays an important role in modern industry. It is used in plastic, paints and oil production [1,2], paper making, chemical manufacturing, and offshore oil extraction [3]. Due to the high neutron absorption of barium, barite became also one of the components of neutron shielding materials [4]. Radium is one of the late ²³⁸U decay products, which can become the dominant source of radioactivity in nuclear waste repositories [5]. The ability of uptaking radium makes barite a perspective material that can potentially be used to significantly improve the safety of nuclear waste repositories. Investigation of the barite dissolution process is, therefore, highly important for a better understanding of marine sedimentation [6,7], as well as for the construction of safer nuclear waste repositories.

Previously, barite dissolution was studied primarily by means of experimental [8,9,10,11,12,13,14,15,16,17] techniques as well as by atomistic models [18,19,20,21]. Experimental studies provided an in-depth insight into catalytic effects of salts [8,17] and chelating agents [16,22] onto morphological changes of reactive surface topography and dissolution rates. However, dissolution in pure water is exceedingly harder to investigate with conventional microscopic techniques, such as atomic force microscopy (AFM). Kuwahara (2011,2012) provided a detailed description of barite dissolution in pure water at different temperatures [14,15]. Moreover, some data are available from other studies [17]. Stack et al. [19,20,21] conducted a detailed investigation of the barite-water interface and reactive sites at the molecular scale. In particular, Stack et al. [19] developed a Molecular Dynamics approach to calculate reaction rates of surface site detachment and tested their method successfully on a Ba kink site. The details of the barite-water interface structure were also studied by using X-ray reflectivity techniques [18,23]. These studies showed that water adsorbed at the interface has a four-layered organized structure and coordinates surface ions. Another important observation is a significant relaxation of surface SO₄ groups in the presence of water [23]. These studies indicate that the mineral–water interface structure, especially in the vicinity of reactive sites, e.g., step and kink sites, should play an important role in mineral dissolution and growth kinetics.

Despite these detailed studies of the dissolution mechanism(s) at the microscopic and atomic scales, the direct connection between the scales is not yet well understood. Reactions at the atomic scale and the interface structure are not well-linked with the overall mechanistic picture. For example, in order to calculate step retreat velocity, one would need to obtain at least two, commonly three, detachment rates from different surface sites (Figure 1A): (1) detachment from a step site and formation of a “negative” double kink which would require breaking of four bonds; and (2) detachment from a kink site, left or right (they may be not necessarily equivalent, depending on the lattice structure) which would require breaking of three bonds. The step advancement velocity taking place at oversaturated conditions during crystal growth, would require two–three rates (Figure 1B): (1) formation of a “positive” double kink via formation of two new bonds; and (2) attachment at either left or right kink site. Barite surfaces have at least three types of atomic steps oriented along different crystallographic directions [10,14,15,16,17,22], two chemical types of lattice sites, Ba and SO₄. Therefore, the minimum number of required kinetic parameters for a microkinetic model of dissolution is six, the same for growth. At least nine parameters are required for dissolution or growth at near equilibrium conditions if attachment and detachment rates can be coupled via chemical potential differences [24,25,26]. At near equilibrium conditions both attachment and detachment processes must be considered, as well as parameters for surface diffusion if relevant [26,27].

Although the elementary reactive events in the Kossel type of crystal (Figure 1) seem to be easily identifiable, it may be not easy to recognize geometry, location, and bonding topology of process-controlling sites for an arbitrary mineral structure. The aim of this study is to construct a multi-scale model of barite dissolution by relating process-controlling elementary molecular reactions and microscopic scale processes, such as etch pit formation and atomic step propagation. We achieve this goal by using a kinetic Monte Carlo (KMC) heuristic approach that we adopted for two decades to understand the dissolution of various minerals, such as carbonates, quartz, feldspars, and phyllosilicates [28].

Figure 1. A schematic drawing for elementary processes on crystalline surface, Kossel crystal model, first introduced in [29,30]. (A): Atomic step retreat during crystal dissolution; (B) Atomic step advancement during crystal growth. The letters “R” and “L” denote right and left kinks, that may have non-equivalent propagation rates in non-Kossel structures.

Figure 1. A schematic drawing for elementary processes on crystalline surface, Kossel crystal model, first introduced in [29,30]. (A): Atomic step retreat during crystal dissolution; (B) Atomic step advancement during crystal growth. The letters “R” and “L” denote right and left kinks, that may have non-equivalent propagation rates in non-Kossel structures.

In principle, a KMC modelling approach is a powerful instrument to study the dissolution process at the microscopic scale and to test working hypotheses regarding kinetically important reaction controls. The method allows us to reveal the evolution of the reactive system from one surface configuration to another as a function of time and to determine the influence of kinetic parameters on the surface morphology. One of the most important fundamental challenges related to the KMC approach is the delicate interplay between the model’s complexity, computational feasibility, and robustness with regard to modelling parameters [31,32,33,34,35]. The KMC model complexity is reflected in the number of factors influencing site reactivity. The computational feasibility requires reasonable assumptions for classifying surface sites into distinct reaction groups. The robustness of the model means that quantitative relationships between the model parameters should have reasonable grounds from the physical chemistry of the process. In other words, the KMC model should not be taken as a black-box simulator mimicking some surface features, but rather serve as a heuristic tool for construction of the overall mechanistic picture of the process.

In this paper, we studied mechanisms of barite dissolution by using the KMC method. As a main result, we identified the minimum set of kinetically relevant reactive sites necessary to construct a microkinetic model for dissolution. We have also found that common KMC approaches to parameterize dissolution rates by using site lattice coordination strikingly fail to reproduce experimentally observed step velocities and etch pit morphologies altogether. We suggest that the water structure and dynamics near the reactive sites along different crystallographic directions should be substantially different and influence on the molecular detachment rates. These molecular details should be investigated further with complementary atomistic methods such as Molecular Dynamics and ab initio calculations. Such studies could provide a better understanding of the dissolution, the ion adsorption, and the crystal growth processes on barite surfaces.

2. Materials and Methods

2.1. Simulated System

The simulated system is constructed by translations of the barite unit cell (a = 8.884 Å, b = 5.458 Å, c = 7.153 Å, α = 90°, β = 90°, γ = 90°, Pnma space group [36]) along three crystallographic axes (Figure 2) and making a cut at the (001) face. Atomic coordinates for translation in the crystallographic basis set were taken as output of the XtalDraw 1.0 program [37]. The simulated system size 600 × 600 × 10 unit cells, which means 2.9 × 10⁷ lattice sites in total. This system size is sufficient to reproduce experimentally observed etch pit morphologies and calculate step velocities independent of a system size (see Supplementary Information for details).

2.2. The Kinetic Monte Carlo Model

The kinetic Monte Carlo approach relates rates of molecular reactions to their probabilities in a specified ensemble of possible reactions and reactive species. This coarse-grain method propagates a reactive system between states before and after a specific reaction happens and, thus, simulates a sequence of reactive events. The decision to perform a specific reaction is made based on a random number generator and reaction probabilities. The method essentially lacks any molecular details describing the reaction mechanism at the molecular scale. The reliability of the model strongly depends on a suitable definition of a reactive ensemble and the correct estimation of reaction probabilities.

2.2.1. Reaction Rates and Probabilities

In the present study the rates of dissolution reactions were defined as functions on the numbers of Ba-O-S (i) and Ba-O-Ba (j) bonds for Ba²⁺ ions, and S-O-Ba bonds for SO₄²⁻ ions:

$$k_i^{SO4}=v·W_T·W_s·exp\left(-i\frac{∆E\left(S-O-Ba\right)}{kT}\right)$$

(1)

$$k_{ij}^{Ba}=v·W_T·W_s·exp\left(-i\frac{∆E\left(Ba-O-S\right)}{kT}+j\frac{∆E\left(Ba-O-Ba\right)}{kT}\right)$$

(2)

Here, $v$ is the reaction attempt frequency in Hz (10¹² Hz for this model as an order of water vibration frequency, this approximation was introduced by Pelmenschikov et al. [38]), $W_T$ is the factor used for the terrace sites to reduce their reactivity in studies of single etch pit morphologies and set 1 to all other sites, $W_s$ is a site-specific correction factor for steric factors, T is the temperature, k is Boltzmann constant, and $∆E$ represents activation energy factors for bond breaking. The corrections for the steric factors were calculated as “activation” functions on parameters corresponding to specific ions (either Ba²⁺ or $S{O}_4^{2-}$ ion in the uppercase index) and their coordination (lowercase index):

$$W_s=exp\left(-\frac{w_{coordination}^{ion}}{kT}\right)$$

(3)

The default $w_{coordination}^{ion}$ value for all sites is 0, with exception for sites with a special coordination environments: (1) SO₄ kink sites with coordination 4 and 5 (see Figure 4 in the Results section); and (2) Ba step sites (i = 4) with 6 and 7 s order (Ba-O-Ba) neighbors (j = 6 or 7) and with two terrace neighbors in the uppermost atomic layer (see Figure 5 in the Results section). The reason for implementation of this parameter is twofold: (1) to enhance dissolution rates for SO₄ kink sites, otherwise the triangular shape of the pits cannot be reproduced; and (2) to recognize Ba step sites along the [120] direction apart from Ba step sites at [010] fast and slow steps. The sites along [120] direction have differently coordinated SO₄ neighbors, and, as a result, different steric environments (see Figure 5 in the Results section). This site distinction is necessary to reproduce straight steps along the [120] direction.

A similar formula was used to reduce probabilities for terrace sites removal ($w_T=0$ for all non-terrace sites):

$$W_T=exp\left(-\frac{w_T}{kT}\right)$$

(4)

This parameter can be omitted for general purposes; here, it was used only to suppress formation of additional pits that interfere with step velocity calculations.

The rates of atomic attachment are excluded from the present simulation approach because here we are focused on the dissolution process only at far-from-equilibrium conditions. In these conditions, dissolution products are assumed to be removed immediately from the surface by the flowing fluid, so we didn’t include surface diffusion processes.

2.2.2. The Algorithm

An in-house KMC program was used to implement the model according to the BKL adaptive time step algorithm [39], also known as the “divide-and-conquer” algorithm [40] involving running sums [41]. The BKL algorithm is essentially rejection-free adaptive time step algorithm [39,40,41]; so at each iteration step, a one dissolution reaction type is randomly chosen according to assigned event probabilities. This step is followed by random choice of a surface site of a certain type. The time is propagated by using a standard formula [41]:

$$∆t=-\frac{1}{Q}\ln \left(r\right)$$

(5)

where r is a uniformly distributed random number, Q is the sum of all reaction rates. Probabilities are essentially normalized as follows:

$$P_i=\frac{N_i{k}_i}{{{\displaystyle \sum }}_{i=1}^{N_{steps}}{N}_i{k}_i}=\frac{N_i{k}_i}{Q}$$

(6)

Details of code organization are shown on the workflow chart (Figure 3). System geometry data, rate parameters, and output details are read first from an input file. Then the program generates a cleavage surface and opens dislocation hollow cores in the middle of the surface. The size of the hollow core for a point defect is essentially one Ba atom, the size of the hollow core for screw dislocation is 3 × 3 × 5 unit cells. There is no strain energy considered in this model, the hollow cores were opened in the beginning of the simulation. This approach was suggested by Meakin and Rosso [40], who arrived at identical results for etch pit growth mechanisms in cases where hollow cores opened first due to excess strain energy and in cases when hollow cores were pre-opened before the simulation start. The simulation run is performed iteratively, where at each iteration step reaction, probabilities are calculated, a random number is generated, and a decision regarding a reaction type is made. A surface site belonging to a chosen reaction type is randomly selected and is then removed from the reactive site list. Neighborhood is updated accordingly, and output data are tabulated at a specified number of iteration steps. Periodic boundary conditions were used to avoid effects of system size limitations.

2.2.3. Parameterization

Calculations of reaction rates at the molecular scale cab be done by using either ab initio or Molecular Dynamics calculations. The KMC method cannot provide molecular detachment or attachment rates and requires them as input parameters [41]. A “classical” approach to scale detachment rates by the number of bonds to be broken [42] may not work for some systems, e.g., the second coordination sphere [43,44,45] or reactive site geometry [46,47,48] has to be considered as important factors influencing the dissolution kinetics. The values of kinetic parameters (e.g., activation energies and reaction attempt frequencies) can be obtained by fitting them to experimentally measured step velocities (e.g., as in [46,47,48]). However, it is important to remember that fitted rates do not necessarily provide a correct estimation of the molecular detachment rates, since the number of fitted parameters is often larger than the number of data for fitting (for example, the rates of obtuse and acute kink site detachments, as well as formation of a double kink constitute three parameters vs. one value of measured step velocity [46,47,48]). The straight step velocity can be calculated as function on kink site propagation and double kink generation rates [49] (Figure 1), which makes these parameters inseparable if no other additional data are available. The major factor responsible for a correct estimation of kink site density along a step is the ratio between rates of kink propagation and a double kink generation [50]. Although this parameter would provide helpful information for KMC model parameterization, estimation of kink site densities from experimental data is typically not feasible, with a few exceptions [51].

The major purpose of a KMC model fitted by experimental data is its utilization as a heuristic approach to identify primary mechanisms and kinetic factors driving dissolution, growth and adsorption processes [28], as well as surface-induced catalysis [31,32]. The values of the KMC parameters should not be taken as “true” molecular scale rates, which should be calculated by using atomistic simulations. The aim of our present parameterization approach is to reproduce the atomic step velocities and experimentally observed etch pit morphologies altogether [14,17], while the ultimate research goal is to reveal bonding topology of most important reactive sites.

First, parameterization attempts revealed the ultimate importance of the second coordination sphere for Ba²⁺ ions (Ba-O-Ba bonds) for dissolution rates from Ba²⁺ surface sites. The use of the parameter set I (

Table 1. Parameter sets for the Kinetic Monte Carlo model of barite dissolution; all parameters are provided in kT units, except $v$, the reaction attempt frequency factor, which is given in Hz.

Parameter	Set I	Set II	Set III	Set IV 1
$v$1	1012	1012	1012	1012
$w_T$	30	30	30	30
$∆E\left(S-O-Ba\right)$	7	7	7	7.7
$∆E\left(Ba-O-S\right)$	7	7	7	7.7
$∆E\left(Ba-O-Ba\right)$	0	1	1	1.2
$w_{i=5}^{SO4 }$	0	0	−5	−5
$w_{i=4}^{SO4}$	0	0	−5	−7
$w_{4,7-2T}^{Ba}$	0	0	−3	−6
$w_{4,6-2T}^{Ba}$	0	0	−5	−7

¹ Parameters are fitted to roughly reproduce experimental data in pure water (see Table 2).

), which incorporated only Ba-O-S and S-O-Ba bond breaking energies, resulted in the formation of an oval-shaped pit (Figure 4A). Parameter set II (Table 1) included activation energy for Ba-O-Ba bond breaking. This parameter set enabled the formation of trapezoidal pits (Figure 4B). Parameter set III included non-zero parameters for steric factors by taking into account specific geometric arrangements of surface sites (see Results for detailed description). The use of these steric factors resulted in formation of triangular pits (Figure 4C). The second coordination sphere for SO₄ sites (neighbors at S-O-Ba-O-S links) was not included in the parameterization routine, because otherwise it would stabilize SO₄ sites in comparison to Ba sites, which is an unlikely process according to thermodynamic calculations [20]. This coordination sphere was used to recognize terrace sites from the rest of the reactive sites and assign them the $w_T$ factor reducing removal probability.

Figure 4. Morphologies of simulated monolayer pits on barite (001) face by using different parameter sets and a schematic explanation for step velocities calculations. (A): Parameter set I excluding Ba-O-Ba coordination sphere, (B): Parameter set II including Ba-O-Ba coordination sphere, (C): Parameter set III including steric factors $w_s^i$. White lines here are used as auxiliary lines perpendicular to the [010], $\left[\overline{2}10\right],$ and $\left[\overline{2}10\right]$-oriented atomic steps. Intersection points of those lines with step edges at two consecutive times steps were used to calculate step propagation increments $∆l_i$.

Figure 4. Morphologies of simulated monolayer pits on barite (001) face by using different parameter sets and a schematic explanation for step velocities calculations. (A): Parameter set I excluding Ba-O-Ba coordination sphere, (B): Parameter set II including Ba-O-Ba coordination sphere, (C): Parameter set III including steric factors $w_s^i$. White lines here are used as auxiliary lines perpendicular to the [010], $\left[\overline{2}10\right],$ and $\left[\overline{2}10\right]$-oriented atomic steps. Intersection points of those lines with step edges at two consecutive times steps were used to calculate step propagation increments $∆l_i$.

The step velocities were calculated by using a simple formula:

$$v_{step, i}=\frac{∆l_i}{∆t}$$

(7)

where $∆t$ is the time step between two measurements, and $∆l_i$ is local step propagation increment calculated as a distance between two points located at the intersection of step edges and auxiliary lines perpendicular to the [010], $\left[\overline{2}10\right],$ and $\left[\overline{2}10\right]$ crystallographic directions (Figure 4C). The results for step velocity calculations are shown in the

Table 2. Step velocities on barite surface, KMC simulations and experimental measurements at 30 °C [14], 10⁻¹ nm/s.

	[010] Slow	[010] Fast	[120]
Pure water, experiments [14]	0.09 ± 0.01	-	0.18 0.01
Pure water, experiments [17]	0.3 ± 0.1	0.8 ± 0.2	-
Pure water, simulations 1	0.16 ± 0.04 2	0.29 ± 0.06	0.30 ± 0.05 3

¹ Parameter set IV, see Table 1; ² Each value represents an average over 10 points, ±standard deviation; ³ The value for step in a monolayer pit.

. The step velocities produced in the simulations and triangular step morphologies are difficult to reproduce altogether due to unknown true values of dissolution rates from different sites. Our values fall into the ranges reported in two separate experimental studies for dissolution of barite in pure water (Table 2). Our [010] slow step velocity is slightly larger than the experimentally measured value by Kuwahara [14] (0.16 vs. 0.09 × 10⁻¹ nm/s) and twice lower then reported by Risthaus et al. [17] (0.16 vs. 0.3 × 10⁻¹ nm/s). Our step velocities for [010] fast and [120] are almost identical because their motion is measured within a monolayer pit. We obtained a stable ratio 2 between the [010] “fast”/[120] step velocity and the slow [010] step velocity. The same ratio is observed by Kuwahara [14] and a similar ratio is observed by Risthaus et al. [17].

2.2.4. Data Visualization

Surface topographies are visualized as height colormaps and 3D ball-and-stick models by using VMD (Visual Molecular Dynamics) software, version 1.9.3 [52]. XYZ data for every Ba, S and O atoms were calculated to generate input files for VMD. 3D XYZ files were turned into height colormaps by switching color scheme based on atom types to the color scheme based on coordinate change (Z-direction).

3. Results

3.1. Morhoplogy of Monolayer Pits

Monolayer pits of triangular and trapezoidal shapes are characteristic dissolution features of the barite (001) face [14,15,17,22]. Triangular pits (Figure 5A) form during dissolution in pure water [14,15,17], corresponding to the simulation conditions in the present study. The pit’s morphology is defined by the [010], [120], and $\left[\overline{2}10\right]$ crystallographic directions (Figure 5A,B). The molecular structure of atomic steps forming along these directions can be understood by making cuts along those directions (Figure 5C). The (001) face can be described as a combination of rectangular Ba and SO₄ sublattices. SO₄ terrace sites form alternating vertical rows of sites with coordination 5 and 6. The left [010] (slow and long) step forms along rows of 6-coordinated sites, while the right (fast and short) [010] step forms along 5-coordinated sites. Sites of higher coordination have lower dissolution probabilities and thus stabilize the step. Rows of Ba ions coordinated by surface oxygen atoms belonging to SO₄ groups form the frontlines of both [010] steps. The [hk0] steps form along directions with alternating Ba and SO₄ sites.

Figure 5. Morphology of monolayer pits on barite surface. (A): Experimental data [14] (Kuwahara, Y, In situ hot-stage AFM study of the dissolution of the barite (001) surface in water at 30–55 °C, (2012), Am. Min., 97, 1564–1573, Copyright (2012), MSA, reproduced with the permission of the Mineralogical Society of America); (B): KMC simulations (this work, the parameter Set IV in the Table 1); (C): Primary crystallographic directions on the barite (001) face, white atoms in SO₄ groups correspond to 6-coordinated sites, while cyan atoms correspond to the 5-coordinated sites.

Figure 5. Morphology of monolayer pits on barite surface. (A): Experimental data [14] (Kuwahara, Y, In situ hot-stage AFM study of the dissolution of the barite (001) surface in water at 30–55 °C, (2012), Am. Min., 97, 1564–1573, Copyright (2012), MSA, reproduced with the permission of the Mineralogical Society of America); (B): KMC simulations (this work, the parameter Set IV in the Table 1); (C): Primary crystallographic directions on the barite (001) face, white atoms in SO₄ groups correspond to 6-coordinated sites, while cyan atoms correspond to the 5-coordinated sites.

The triangular morphology of pits was reproduced only after assigning steric factors to specific surface sites, otherwise only trapezoidal pits formed (Figure 4B). The necessity and role of the steric factors can be understood after a detailed analysis of the step and kink site’s local coordination. Further, we present the results of this analysis and discuss step propagation mechanisms.

3.2. Reactive Sites Coordination

3.2.1. Step Sites

Coordination of SO₄ step sites depends on crystallographic direction of an atomic step (Figure 6). A remarkable result for this system is that site coordination cannot be directly related to a structural status, “terrace”, “ledge”, or “kink”. Ba terrace sites have coordination 5, while step sites may have coordination either 4 or 5 (Figure 7). SO₄ terrace and step sites may have coordination either 6 or 5 as well (Figure 6). The second coordination sphere helps to distinguish terrace sites for Ba as (5,6) (Ba-O-Ba neighbors) as well as (5,15) and (6,17) for SO₄ (S-O-Ba-O-S) sites. As we mentioned above, the second coordination sphere for SO₄ ions was not included into parameterization due to thermodynamic considerations [20], but some preliminary trials to set a parameter (1 kT unit per neighbor) for this sphere resulted in rotation of the pit by 180° (See Supplementary Information file). This result cannot be verified by experimental data, because the same rotation of pits takes place in consecutive half-layers parallel to the (001) plane due to 2₁ stacking rotation axis. We assume that this scenario is not realistic because the nature of bonding in barite is primarily ionic [20].

The SO₄ step sites at the left (slow) [010] step have coordination 6 (Figure 6A). This coordination makes them quite unreactive in comparison to Ba sites along the same step, which have coordination 4. The step dissolution process is initiated by a removal of a Ba step site and generation of SO₄ double kink with coordination 5 (Figure 6A). The SO₄ step sites at the right (fast) [010] step have coordination 5 (Figure 6B). Removal of Ba ion from this step generates SO₄ kinks with coordination 4. The [hk0] steps have quite irregular morphology, although they are in general oriented along the $\left[\overline{2}10\right]$ and [120] directions. There are two types of segments forming along this direction: a straight “classic” [120] cut with alternating Ba and SO₄ sites along the step (Figure 6C). Since this is a diagonal cut through the terrace made up by 5 and 6-coordinated SO₄ sites, the coordination of step sites is alternating from 4 to 5. Formation of a corner segment (Figure 6C) disrupts this structure due to preferential removal of 4-coordinated sites and results in formation of segments dominated by corners (Figure 6D).

Ba step sites in most cases have coordination 4 (Figure 7), thus, typically becoming the first sites to dissolve from the step, in comparison to SO₄ sites. Ba sites at the left (slow) [010] step have coordination (4,7), which is slightly higher than coordination (4,6) at the other steps (Figure 7B–D). The (4,7) site binds to the two SO₄ sites at the upper terrace. The same binding geometry is observed for (4,6) sites at the right (fast) [010] step (Figure 7B) and a corner segment along the [120] step (Figure 7C). This type of “2T” binding geometry makes Ba to be sterically less hindered for the surrounding water molecules than the “3T” sites forming along straight [hk0] segments (Figure 7D). We suggested that this “2T” coordination may increase probability for the Ba-O-S bond hydrolysis, so we assigned a weighting coefficient to those sites (see Table 1 and Methods).

3.2.2. Kink Sites

Removal of Ba ions from step sites initiates the formation of SO₄ double kinks (Figure 7A). Propagation of 5-coordinated SO₄ kink along the left [010] face substantially influences the step retreat rate (Figure 8A). Although 4-coordinated kinks at the neighboring SO₄ row occasionally form, they are in general less stable. SO₄ kink sites at right [010] step have coordination 4 (Figure 8B). Kinks of a similar structure and the same coordination appear at [hk0] steps (Figure 8C). SO₄ sites with coordination 4 also appear at a variety of geometric positions at [hk0] steps, such as corners, overhangs, and terrace sites. In general, there is no specific coordination for SO₄ kink sites. Strictly speaking, a “terrace-ledge-kink” structure is applicable only to the left and right [010] steps. [hk0] steps dissolve as small stable segments and clusters randomly forming along [120] directions. Ba kink sites typically have coordination 3 (Figure 8D) in contrast to SO₄ kink sites of a higher coordination. Depending on a sensitive parameter balance, either Ba or SO₄ kink sites along [010] appear as stable sites (compare configurations on Figure 8A,B,D,E). The parameter set IV in the Table 1 results in formation of stable Ba kink sites along the [010] steps. Ba kink sites are in general not stable on [hk0] faces, but still can occur to some extent as overhangs (Figure 8F). The balance between stable Ba and SO₄ kink sites in general largely depends on the parameter choice and cannot be verified by only knowing etch pit morphology and step velocities.

3.3. Morphology of Mutlilayer Pits

Multilayer pits form around screw dislocation hollow cores. The morphology of the pits is defined by a geometric superposition of triangular pits rotated by 180° in each consecutive layer. As a result, pits have hexagonal structure (Figure 9A). We reproduced the shape of these pits in KMC simulations, although their morphology is not precisely euhedral (Figure 9B). Probably this discrepancy is caused by the scale difference and parameterization issues. The structure of multilayered pits clearly indicates the slow steps in the uppermost atomic layer limit motion of the fast steps at deeper layers. The step density is quite high due to a limited terrace width (about a few unit cells). As a result, the pits have steep walls. Since fast and slow steps superimpose due to stacking rotation, they commonly form sequences of bunched steps with 2-layer macrosteps (Figure 9). Kuwahara (2011, 2012) [14,15] described in detail the formation of these bunched structures. He found that etch pit spreading rates are about an order of magnitude slower than the spreading of unrestricted steps. This result matches to our observations of etch pit kinematics where the slowest step determines the etch pit spreading rate.

4. Discussion

4.1. Rate Limiting Step

Studies of mineral dissolution commonly include discussions about a rate limiting step or a rate-limiting molecular reaction which dissolution rate can be extrapolated to macroscopic systems. Our present study reveals that the barite surface contains quite a diverse number of reactive sites. A large variety of possible geometric kink and step site arrangements and sensitivity of their removal rates to the local neighborhood indicate a fundamental issue in identification of a rate limiting reaction. The simulation outcome in terms of etch pit morphology and step velocities is highly sensitive to the parameter values as well as site identification and differentiation procedure. Similar etch pit structures, for example, can be reproduced if we make Ba or SO₄ kink sites to be a little more stable relative to each other. Clearly, dissolution of Ba from the [010] step sites has a large importance for defining etch pit morphology, although it is not quite clear whether formation of double kinks on [hk0] faces starts from Ba or SO₄ step site removal. Step morphology and propagation are defined by a delicate interplay between dissolution rates at a large tapestry of surface sites. These rates apparently are highly sensitive, not only to the local geometric neighborhood, but also to local water structure at these sites. Our results indicate a necessity for more detailed investigation of dissolution mechanisms from identified rate-determining sites ensemble by means of Molecular Dynamics and Quantum Mechanical methods.

4.2. Recommendations for Molecular Dynamics Calculations

Molecular Dynamics and ab initio calculations may potentially shed light onto molecular controls of site reactivity and role of local site geometry and water configuration. Some successful attempts to understand the dissolution of barite were conducted by Stack et al. [19], who developed an efficient way to calculate reaction rates of surface site detachment. They used a combination of Molecular Dynamics approach powered by a well-developed force field [53] and rare event theory to simulate dissolution of Ba 3-coordinated site via step-wise bond breaking. A similar approach can be very beneficial for further investigations of dissolution mechanisms from the other kinetically important sites. The tentative list includes but is not limited to:

(1)

Dissolution of Ba step sites from the left [010] step where Ba has coordination (4,7) and from the right site where Ba has coordination (4,6);

(2)

Dissolution of Ba step site from the straight [120] step;

(3)

Dissolution of SO₄ kink sites of coordination 5 at the left [010] step and coordination 4 at the right [010] step;

(4)

Extra calculations for checking preferential removal of SO₄ and Ba at step sites. These calculations should ensure whether Ba or SO₄ kink sites along each kinetically relevant atomic step influence onto step propagation rates.

5. Summary and Conclusions

Despite the importance of the barite–celestite system for the nuclear waste community, the number of studies on the dissolution mechanisms is quite limited. Here, we presented a first kinetic Monte Carlo simulation approach of barite dissolution that can be extrapolated to isostructural sulfates (e.g., celestite and anglesite). Our study revealed that common approaches to parameterize a KMC model by classifying reactive sites based on their lattice coordination is not sufficient for reproducing experimental data. We suggest that local structure and dynamics of water in the vicinity of reactive sites should substantially affect molecular detachment rates. Moreover, the influence of the barite–water interface structure and dynamics is likely to overweight the influence of lattice coordination. We identified a minimum set of reactive sites that might control the overall reaction mechanism of dissolution. In principle, this set can be expanded further on for more realistic modelling of the dissolution process. Molecular Dynamics simulations and/or quantum mechanical calculations for the reaction mechanisms and rates at these different energetic sites will be ultimately required for a comprehensive understanding of this important mineral system.