Probability distributions of mineral dissolution rates: the role of lattice defects

Kurganskaya, Inna; Luttge, Andreas

doi:10.3389/frwa.2023.1225837

ORIGINAL RESEARCH article

Front. Water, 28 September 2023

Sec. Water and Critical Zone

Volume 5 - 2023 | https://doi.org/10.3389/frwa.2023.1225837

This article is part of the Research Topic Reactivity at Pore Scale: from Molecular Mechanisms to Reaction Rates View all 3 articles

Probability distributions of mineral dissolution rates: the role of lattice defects

$\r\nInna Kurganskaya,$ Inna Kurganskaya^1,2^*

Andreas Luttge^1,2,3

¹Geoscience Department (FB5), University of Bremen, Bremen, Germany
²MAPEX, Center for Materials and Processes, University of Bremen, Bremen, Germany
³MARUM-Center for Marine Environmental Sciences, Bremen, Germany

The correct quantification of mineral dissolution rates is a critical task for macroscopic reactive transport modeling. Previous studies showed a substantial rate variability of about two orders of magnitude, which cannot be explained by variance of external environmental parameters alone. If the rate cannot be predicted as a constant parameter, then the critical question is whether it can be predicted as a stable reproducible probability distribution. Although a large variety of factors may contribute to the overall variance across the scales, the effect of defect density and defect spatial distribution can be considered as one of the key variance sources. Here, we tested the reproducibility of probability distributions for Kossel crystals with a different amount and spatial configurations of lattice dislocations. We ran several tests on systems with the same configurations and calculated the probabilities of material flux. Surprisingly, we discovered that the density of dislocations has minimal impact on the probability distributions. However, the spatial location of dislocations has a substantial influence on the rate distributions reproducibility. In cases where multiple etch pits operate simultaneously, reproducible rate distributions are found regardless of the number of dislocations. In cases where dislocations formed clusters, one large etch pit controlled the entire surface, and sets of reproducible probability distributions were detected. Then, more complex statistical behavior is expected, since the result is path-dependent. These results have serious consequences for the implementation of rate distributions in reactive transport models. Further studies, however, are needed to provide clear guidance on relating surface morphologies, dislocation distributions, and dissolution rate variance. The role of material-specific properties, such as crystallographic structure and bonding, in rate distributions, should be additionally addressed. The role of grain boundaries, crystal size and crystal habit, including nanoparticulate forms, in rate variance, also should be addressed for practical applications.

1. Introduction

The growing demand of modern societies for reliable and predictable numerical models of critical geochemical processes can now be fulfilled with computational power. The ideas to solve important environmental problems, e.g., CO₂ sequestration (Balashov et al., 2013; Daval et al., 2013; Hellmann et al., 2013; Jun et al., 2013, 2017; Smit, 2016; Arif et al., 2017; Daval, 2018; Wild et al., 2019; Loganathan et al., 2020; Deng et al., 2022; Shabani et al., 2022; Urych et al., 2022), water cleaning (Bhattacharyya and Gupta, 2007; Barry et al., 2021), soil enrichment and remediation from pollution (Aredes et al., 2012; Björneholm et al., 2016; Niazi et al., 2023), prediction of oil and gas reservoir behavior (Browning and Murphy, 2003; Steefel et al., 2015), or the safety of nuclear waste repositories (Payne et al., 2013; Kalinichev et al., 2017; Leal et al., 2017; Vinograd et al., 2018; Androniuk and Kalinichev, 2020; Wieland et al., 2020; Cygan et al., 2021; Claret et al., 2022; Liu et al., 2022), geothermal resource modeling (Wilson et al., 2001; Yapparova et al., 2014, 2019, 2023; Whitaker and Frazer, 2018; Lamy-Chappuis et al., 2022), in silico via computer simulations seem to be both promising and economically efficient. Understanding the process-controlling reactions at mineral surfaces, i.e., dissolution, adsorption, nucleation, and crystal growth, constitutes a milestone in geochemical process modeling (Brantley et al., 2008; Lichtner et al., 2018; Xiao et al., 2018; Poonoosamy et al., 2021; Prill et al., 2021; Deng et al., 2022). The combination of such models with sophisticated experimental techniques and complementary field studies propels us into a new era of multidisciplinary science. Communication between the communities of geochemists, physicists, chemists, environmental engineers, and applied mathematicians, to name only a few, becomes especially important in this perspective.

The ideas of pioneers such as Vernadsky and Goldschmidt (Müller, 2014) presented all natural processes on Earth as perpetual movement of chemical elements between geochemical reservoirs: the Earth mantle and crust, soils, rivers, watersheds, swamps, lakes, and oceans. Each reservoir is characterized by their own predominant thermodynamic and kinetic mechanisms of material production (nucleation and growth) and destruction (mechanical decomposition and dissolution). Each location within the respective reservoir is a complex chemical microcosm powered by fundamental physical forces operating at different scales, i.e., gradients in temperature and pressure, gravity, chemical potentials, and electrostatic interactions. The fundamental laws of interactions, such as Newtonian equations of kinetic motion, or Coulomb interactions, alone cannot always capture the complex chemical machinery of chemical element conversion and transport from immobile (solids) to mobile (ions in fluids) forms. However, we note that classical chemical kinetics providing reaction rates for very simple molecular processes in many practical applications cannot be considered as a reliable choice for macroscopic modeling of geochemical processes (Lüttge et al., 2013; Luttge et al., 2019). A vast variety of molecular reactions taking place at solid-fluid interfaces of complex solid materials may happen in such systems (Brantley et al., 2008). Material transport is another factor that adds complexity to reaction mechanisms at different scales. The fundamental problem is, thus, how to model geochemical processes across the scales taking into account the complexity of geochemical reactions and material transport (Zhu, 2009).

It is necessary to obtain data about a reactive solid-fluid system at a variety of scales, i.e., from the molecular scale via molecular modeling (Kubicki, 2016) and spectroscopic techniques (Hawthorne, 2018), to the reservoir scale via modeling porosity networks and making massive CT scans (Lichtner et al., 2018). This comprehensive approach provides knowledge about geochemical reaction mechanisms, although the process of data collection is tedious and expensive. The critical issue related to such studies is transferability of knowledge obtained for each process of interest at a specific location to other similar processes at other locations. A remarkable example of geochemical process complexity is an intrinsic variance of mineral dissolution rates, which cannot be attributed to the variance of external environmental conditions (Arvidson et al., 2003; Fischer et al., 2012, 2014; Lüttge et al., 2013; Luttge et al., 2019). The intrinsic rate variance is reported to be as large as two orders of magnitude and attributed to variance of reactive surface site density and reactivity of different surface features (Arvidson et al., 2003; Fischer et al., 2012; Levenson and Emmanuel, 2013; Lüttge et al., 2013; Noiriel et al., 2018). The full description of variance sources is a non-trivial task: thus, molecular scale controls, such as water exchange frequency at different sites (Wolthers et al., 2013), different crystallographic face orientations (Godinho et al., 2012, 2014), different crystal habits (Bouissonnié et al., 2021), or etch pit morphology evolution (Brantley et al., 1986; Bandstra and Brantley, 2008; Pollet-Villard et al., 2016).

The underlying work hypothesis is the intrinsic stochasticity of the dissolution rates. A fundamental problem of practical importance is how to predict dissolution rates for macroscopic scale modeling (Li et al., 2006; Balashov et al., 2013, 2015; Lugo-Méndez et al., 2015; Steefel et al., 2015; Liu et al., 2017; Lichtner et al., 2018; Erfani et al., 2019; Agrawal et al., 2021; Poonoosamy et al., 2021; Prill et al., 2021). The fundamental question is thus if we cannot model geochemical processes in a deterministic way, is a probabilistic approach the alternative?

A probabilistic approach that assumes an ensemble of possible events with related probabilities of their occurrences is widely used in physics, physical chemistry, and geosciences. There are two well-known probabilistic algorithms, the Metropolis Monte Carlo (Frenkel and Smith, 2002; Binder and Heermann, 2010) and Kinetic Monte Carlo algorithms (Blum and Lasaga, 1987; Jansen, 2012; Andersen et al., 2019). The Metropolis Monte Carlo algorithm is a standard tool of statistical physics that allows direct statistical sampling over energetic states which have different occurrence probabilities. In this way important macroscopic parameters, such as systems' energy, can be calculated. The Kinetic Monte Carlo (KMC) algorithm relates reaction rates in a reactive ensemble with probabilities of their occurrences (Gillespie, 1977; Gilmer, 1980; Cheng, 1993; Voter, 2007; Hess and Over, 2017). Estimation of mineral resources employ geostatistics, which relies on random sampling approaches and implies a careful separation of a meaningful signal from stochastic noise. The microscopic complexity of any system is an inherent property that cannot be ignored, but instead can be used to quantify and predict the macroscopic behavior if a proper stochastic model of a process and related statistical sampling approach are established.

The first question of practical interest is whether there is a reproducible probability distribution of dissolution rates which can be directly implemented into macroscopic scale models, e.g., reactive transport at the pore scale. Probabilities in many cases can be approximated by counting and normalizing event frequencies. In this manuscript we deliver the statistical results of numerical KMC experiments on a simple Kossel cubic crystal with screw dislocations, where we counted frequencies of dissolution rate occurrences. We compared rate distributions and dissolution surface features, i.e., etch pits, which generate these distributions. We then arrive at the conclusion regarding statistical behavior of dissolving surfaces depending on spatial location of dislocations.

2. Methods

We used the standard Kinetic Monte Carlo simulation algorithm as it was developed originally for Kossel crystals by Gilmer and Bennema (1972), Gilmer (1976, 1980). A Kossel crystal has simple cubic lattice structure and atomic units are represented by cubic blocks (Kossel, 1927; Stranski, 1928). Kossel crystals are commonly used as a generic Terrace-Ledge-Kink representation of crystalline solid surfaces (Mutaftschiev, 2001) (Figure 1A). If a bulk crystal (coordination number, CN = 6) is cut in two halves along a, b, or c crystallographic axes, atomically flat terrace sites with five neighbors (number of neighbors N = 5) are formed. Atomic rows parallel to the other two axes form ledge sites (N = 4) which are adjacent to the terrace sites and form atomic steps. Removal of atoms from ledge sites generate kink sites (N = 3) which are self-reproducing units if removed sequentially from the same atomic row. In a numerical procedure, filled positions can take value 1 and unfilled value 0, so the entire crystal can be represented as a 3-D array of 0/1 values, or, alternatively, as a 2-D array of surface heights for solid-on-solid models (SOS) where overhangs are excluded. We used the SOS model in the present study.

FIGURE 1

Figure 1. A basic framework for crystal dissolution simulations by the Kinetic Monte Carlo method. (A) Surface of Kossel crystals (bulk coordination number CN = 6) with surface sites: terrace (N = 5), ledge (N = 4), and kink (N = 3) sites. Each surface site has an associated dissolution probability depending on its dissolution rate k_i. (B) Construction of the running sum for the divide-and-conquer (Meakin and Rosso, 2008) or BKL (Bortz et al., 1975) algorithm. The sum of rates for all sites is divided into intervals which lengths corresponds to dissolution rates. Sites with higher rates are selected more frequently than the sites with lower rates.

The KMC algorithm normalizes probabilities for dissolution rates for surface sites with different dissolution rates according to a simple relation (Bortz et al., 1975):

\begin{array}{l} P_{i} = \frac{k_{i}}{\sum_{i = 1}^{M} k_{i}} = \frac{N_{i} \exp (- i β φ)}{\sum_{i = 1}^{M} N_{i} \exp (- i β φ)}, & (1) \end{array}

where M is the total number of possible reactive events, k_i is the reaction rate for dissolution event of an atom with i neighbors (Figure 1A), N_i is the number of sites with i neighbors, φ is a bond breaking activation energy, β = kT is the Boltzmann parameter, where k is the Boltzmann constant, T is the temperature. We used β = 6 in this study, as a parameter providing straight steps which emerge in the majority of mineral-fluid systems. This parameter can be used to model dissolution of calcium carbonate (Kurganskaya and Luttge, 2016), one of the best studied minerals with regard to its dissolution kinetics (Morse et al., 2007).

The rejection-free algorithm, known as divide-and-conquer (Meakin and Rosso, 2008) or BKL (Bortz-Kalos-Lebowitz) algorithm (Bortz et al., 1975) is used in this study as it is indicated in the right hand side of the Equation (1). According to this algorithm, a type of reaction is chosen first, then a site where reaction is performed is chosen at random. The running sums (Voter, 2007; Meakin and Rosso, 2008) are constructed as sequential sums of dissolution rates k_i, and a random number falling into an interval (k_i; k_i+ k_i+1) is used to decide upon a reaction type i (Figure 1B).

The time step between reactive events is calculated as follows (Voter, 2007):

\begin{array}{l} t = - \frac{\ln (r)}{\sum_{i = 1}^{M} N_{i} \exp (- i β φ)}, & (2) \end{array}

where r is a random number. Surface dissolution rates for rate distribution analysis (n data points per distribution) were calculated as an inverted time T_n required to dissolve axb surface sites, where a is the length and b is the width of the simulated system. In this study a = b = 200 atomic units:

\begin{array}{l} R (T_{n}) = \frac{1}{T_{n}}, & (3) \end{array}

\begin{array}{l} T_{n} = \sum_{j = 1}^{a • b} t_{j} . & (4) \end{array}

Screw dislocations were placed into the crystal model as opened hollow cores according to the procedure discussed in previous publications (Meakin and Rosso, 2008; Fischer et al., 2014; Kurganskaya and Luttge, 2016). The location of screw dislocations were chosen at random and then simulations were run at the same dislocation configuration five times in independent runs. We run simulations on systems with two, three, five, ten, twenty, thirty, and forty dislocations. An independent test study for dislocations placed on a regular 5 × 5 grid was performed to assess the role of dislocation positioning. Surface height maps in this study are shown as colormaps where each color represents an atomic layer at a constant height.

Bin size was chosen in the same way in each simulation run, where the representative rate interval was chosen between (2ab)/k₃ and (0.5ab)/k₃ and the total number of bins is 50, so the bin size equals to (0.03ab)/k₃. Rate frequencies falling into these bins were calculated for simulation trajectories equals to 2•10⁷ atomic dissolution events. This binning approach provides ~500 data points for each rate distribution curve.

3. Results and discussion

3.1. Surface configurations

3.1.1. Interaction of two stepwave sources

Presence of screw dislocations at mineral surfaces normally result in opening of etch pits if a system of interest is not close to equilibrium (Lasaga and Blum, 1986; Lasaga, 1998; Lasaga and Luttge, 2001; Lüttge, 2006). Lasaga and Luttge developed a model of crystal dissolution based on etch-pit opening process and formation of concentric step trains, so-called stepwaves (Lasaga and Luttge, 2001; Lasaga and Lüttge, 2003). The stepwaves are constantly generated at the outcrops of screw dislocation hollow cores. Each stepwave travels across the surface and dissolves the crystal layer-by-layer. The stepwaves can collide into each other and form curved steps which propagate significantly faster than the straight rectangular stepswaves (Jordan et al., 2001; Vinson and Lüttge, 2005). Kinetic Monte Carlo simulations of Kossel crystals with multiple screw dislocations demonstrated that the peak-and-valley surface morphology stems from interaction between the etch pits and the stepwaves they emit (Meakin and Rosso, 2008). Curved steps are in general more reactive in comparison to the straight steps. The contribution of curved steps into the overall dissolution rate, as well as variance in their reactive properties were usually not discussed. In general, two diagonally opposite curved steps form when two straight steps emanating from different sources d1 and d2 collide at any point on a d1-d2 line if we assume square-shaped stepvawes (Figure 2). If the stepwaves are rectangular, which may accidentally happen, they may collide at any point within a square which diagonal is d1-d2 line.

FIGURE 2

Figure 2. Possible etch pit configurations arising in the same simulated Kossel crystal with fixed positions of two screw dislocations, system size: 200 × 200 atomic units. Labels d1 and d2 indicate locations of the hollow cores. (A) Stepawaves from two dislocations coalesce at some point around a line connecting dislocation centers. Two etch pits define dynamics of the dissolution process; (B) Stepawaves from two dislocations coalesce close to the d1 hollow core. The etch pit around d2 starts to dominate process dynamics; (C) Etch pit around d1 completely disappears because its stepwaves become eradicated by the stepwaves coming from the d2 source. The pit around d2 is the only one pit present on the surface; (D) The same dynamics as in (C), but the only one present pit is formed around d1 instead of d1 following the same mechanism.

The location at which curved steps form is never precisely defined because each stepwave velocity varies stochastically locally. As a result, each new couple of circular steps may form at a new location. Depending upon this location, reactivity of these steps varies due to the different densities of kink sites they bear. This variance contributes to the overall dissolution rate variance. If two stepwaves meet closer to d1 dislocation, the upper etch pit dominates the surface (Figure 2A), while if they meet closer to d2 dislocation, the lower pit will dominate the surface (Figure 2B). In the cases when the stepwaves from one of the sources slow down, they can be overridden by the stepwaves from another source. As a result, one of the sources becomes inactive and only one of the etch pits is functioning (Figures 2C, D). If screw dislocations are long enough, these two situations can randomly reverse to I and II. Overall, the system fluctuates between all possible configurations I–IV thus giving temporal rate variance. The system shown on Figure 2 evolves in time from the configuration I to the configuration IV via configurations II and III.

3.1.2. Interaction of scattered stepwave sources

The interaction between multiple stepwave sources is analogs to the interaction of the two stepwave sources presented above. The only difference is that the possible number of etch pits and possible surface morphologies vary to a greater extent. Thus, for a system with ten sources (Figure 3) any number of etch pits between one and ten may form: five major pits (Figure 3A), three major pits (Figure 3E), two major pits (Figures 3B, F), and one major pit (Figures 3C, D). The system gradually evolves from configuration I to configuration VI, but other sequences may randomly occur.

FIGURE 3

Figure 3. Different etch pit configurations arising in the same simulated Kossel crystal with fixed positions of ten randomly placed screw dislocations, system size: 200 × 200 atomic units. Labels di indicate locations of the hollow cores. (A) Five etch pits control the dissolution process; (B) Two major etch pits control the process; (C) One major etch pit with a few suppresed pits controls the process; (D) The cluster of d1-d4-d5 dislocations form one dominating pit; (E) Three etch pits merging into one larger pit control the dissolution process; (F) Two interacting etch pits control the process.

The mechanism for different configurations emergence is the same as for the system with two dislocations competition between the arrival times for stepwaves coming from different sources at the same location. In general, the number of pits is usually smaller than the number of dislocations. The number of etch pits forming depends on the geometrical arrangements of the dislocation sources. Thus, closely spaced dislocations have more chances to form a dominating pit because the arrival time for the stepwaves for their mutual collision is shorter than the arrival times for the more distance sources. The game seems to be defined in the very beginning (Figures 3A–D): three major reactivity zones—(1) d1-d4-d5, (2) d2-d3-d6-d7, and (3) d8-d9-d10—are initiated at the very beginning based on relative distances between the dislocations (Figure 3A). Then the third zone (d8-d9-d10) run out of the game because the stepwaves emanating from it get wiped out by the major stepwaves coming from the other two zones (Figure 3B). Finally, the first zone (d1-d4-d5) wins the game because it produces faster traveling major stepwaves which wipe out the second zone (d2-d3-d6-d7) (Figure 3C). One dominating pit occupies the surface (Figure 3D). This situation seems to be final, but it spontaneously reverses to formation of a new zone (d1-d2-d3-d4-d5) which competes with a zone formed by d8 solely (Figure 3E). In the next step these two zones completely disappear with the two new competing zones to be established: d1-d4-d5 and d8-d9-d10. As we can see from this example, coalescence of etch pits into a one major etch pit is in principle reversible into a multiple etch pit regime, but it is not possible to predict exactly which of the screw dislocations will become new major stepwave sources.

3.1.3. Interaction of clustered stepwave sources

The systems with ten or more screw dislocations showed a distinctly different behavior (Figure 4): the dissolution process started as normal from formation of multiple pits around screw dislocation hollow cores (Figure 4A). The major pit on the right started to take over the entire surface due to larger local density of dislocations (Figures 4A, B).

FIGURE 4

Figure 4. Different etch pit configurations arising in the same simulated Kossel crystal with fixed positions of thirty randomly placed screw dislocations, system size: 200 × 200 atomic units. The configurations are shown as time sequence corresponding to n*105 iteration steps; (A) n = 2; (B) n = 4; (C) n = 10; (D) n = 12; (E) n = 106; (F) n = 188.

The process of major right pit dominance continued, where other sources formed short-living tiny pits on terraces emanating from the major pit's source (Figure 4C). The system is stacked in this one major pit configuration (Figures 4D, F). The width and the density of the straight step trains varied over time thus contributing to the overall rate variance.

3.2. Rate distributions

Dissolution rate frequencies showed quite stable reproducible distributions between the runs for systems with the number of dislocations ranging from two to five (Figures 5A, B). Systems with ten dislocations typically exhibited stable distributions (Figure 5C), but some systems can occasionally exhibit a deviant distribution in a new run (Figure 5D). The systems with twenty dislocations (Figures 5E, F) showed multiple sets of stable distributions depending on a dislocation configuration: the configuration set 1 produced three distinct distributions in five simulation runs (Figure 5E), while the configuration set 2 produced two distinct distributions in five runs (Figure 5F). The same behavior is observed for the systems with thirty dislocations (Figure 5G), while the system with forty dislocations exhibited stable distributions again (Figure 5H).

FIGURE 5

Figure 5. Frequencies of dissolution rates calculated as inverted time steps required to dissolve one atomic layer (1/Δt) in Kinetic Monte Carlo simulations of Kossel crystal dissolution. Each simulation was run five times for the same configuration of screw dislocations. System size in each simulation is 200 × 200 atomic units. (A) 2 dislocations; (B) 5 dislocations; (C) 10 dislocations, configuration set 1; (D) 10 dislocations, configuration set 2; (E) 20 dislocations, configuration set 1; (F) 20 dislocations, configuration set 1; (G) 30 dislocations; (H) 40 dislocations.

The comparison of rate variance amongst systems with different dislocation configurations (Figure 6) showed that dislocation density has a very non-trivial influence on rate distributions. Roughly, there are two types of dislocation configurations which result in different behavior of rate distributions: (1) “low” density where multiple etch pits operate simultaneously on the surface. Rate distributions are stable, reproducible, and there is no or very weak dependence of rate values from the number/density of dislocations. This result is counterintuitive because it may seem that higher dislocation density should result in higher dissolution rates due to greater extent of pitting. We observed this regime for the systems with two, three, five, and ten (set 1) dislocations; (2) “high” density (from 20 to 40 dislocations) where etch pit clustering effect results in formation of a one dominating kinetic feature and the rest of dislocations do not develop into independent pits. In this case rate distributions may significantly vary from run to run, thus introducing an overall variance about half of order of magnitude. As we can see from this result, the enhanced dissolution at larger density of dislocations does not stem from larger number of pits, as it may seem intuitively. Instead, a surface dominating pit from a dislocation cluster defines rate distributions which show a fluctuating behavior. The overall fluctuations in rate distributions and difference in average or modal rate values largely depend on relative distances between the dislocations as etch pit sources, instead of dislocation density as a number of defects per unit area.

FIGURE 6

Figure 6. Compilation of synthetic rate probability distributions for systems with different numbers of screw dislocations, KMC simulations of Kossel crystals, the system size is 200 × 200 atomic units.

3.3. The influence of surface feature organization on rate distribution

As we can see from the data on Figures 5, 6, the efficiency of material removal by a couple of pits is almost identical to the efficiency of tens of pits operating simultaneously. Similar efficiency, however, is valid only for “steady-state” systems where all stepwaves from different sources came into contact with each other. The reason of similar efficiency is simple: the cumulative amount of kink sites is quantitatively similar for a large number of “small” curved steps and a few “large” curved steps formed as a result of straight step interaction (see Figures 2, 3) operating simultaneously on the surface. This number of sources, however, largely controls local rate variance. The global rate variance is controlled by the amount of material flux coming from the entire surface. Therefore, smaller number of dislocations is overcompensated by the larger reactive features. This situation changes when one large structure, such as etch pit, controls the entire surface. In each new run the etch pit sets itself at one of possible probability distributions. A possible explanation is the existence of long-range correlations between propagation of the atomic step trains. The observed convergence of rate distributions at substantially large dislocation density as it is shown on Figure 5H indicates that this density is capable of breaking these long-range correlations. As a result, the system returns to a unimodal distribution regime.

To make sure that this result is not a random artifact, we ran an independent test on a system with 25 dislocations placed on a regular grid and compared it with a system with 20 randomly placed dislocations (Figure 7). The system with regularly spaced dislocations generated reproducible rate distributions as systems with smaller number of dislocations shown above. A large variety of possible etch pit configurations were formed (Figure 7A1), but the probability distributions were stable between the runs (Figure 7A2). On contrast, the system with 20 randomly placed dislocations formed a process-dominating pit (Figure 7B1) with three distinct probability distributions appeared in five independent runs (Figure 7B2). This result indicates ultimate importance of screw dislocation positioning: clustered dislocations have great chances to become a dominant etch pit formation center, while the other dislocations will be not active as step sources.

FIGURE 7

Figure 7. The influence of surface organization on rate distributions, KMC simulations of Kossel crystals with screw dislocations, the system size is 200 × 200 atomic units. Regular grid of 25 dislocations, (A1) surface map, (A2) rate distributions for five independent runs. Randomly spaced 20 dislocations, (B1) surface map, (B2) rate distributions for five independent runs.

3.4. Comparison to experimental data

Direct comparison between data sets shown in the present study and experimentally measured rate distributions is difficult because our system is simple and idealistic. In general, crystal lattice, molecular structure of mineral-water interface, chemical composition of reacting fluid, and other details, may influence onto the results. Most experimental studies for rate distributions were done for calcite (calcium carbonate) crystals, which have anisotropy in step velocities affecting etch pit morphology and statistics for atomic step interactions. The rate spectra collected for single crystals in many cases have mixed sources of variance: dislocations and grain edges. Grain edges dissolve crystal according to a different mechanism: they generate rough steps with enhanced reactivity. Reactivity of these steps is similar to reactivity of circular rough steps discussed in Section 3.1. The source is different, and their statistical behavior also can be different. The contribution to the overall rate distributions from this source is significant because these features are large and very reactive. Another issue is the absence of reproducibility control for systems with identical locations of dislocations. Despite these issues, some aspects relevant to this study can be discussed.

The first aspect is stability of measured rate distributions over time for single crystals where dissolution is controlled by etch pits (Brand et al., 2017; Bibi et al., 2018; Noiriel et al., 2018; Guren et al., 2023). On contrast, calcite with micritic or poly-grain structures showed distributions changing in time (Bollermann and Fischer, 2020; Guren et al., 2023). Limited time of experiments and scale/resolution limitations do not allow us to ensure the eternal stability of rate distributions for those systems. Important information can be obtained by comparison of modal values for rate distributions obtained in different laboratories using different samples and analytical techniques at different scales for single crystals with etch pits. Thus, the values 0.1–0.2 × 10⁻⁶ mol m⁻²s⁻¹ were obtained by Brand et al. (2017) and Bibi et al. (2018) for pH 7, and values 5–8 × 10⁻⁶ mol m⁻²s⁻¹ by Noiriel et al. (2018), Noiriel et al. (2020), and Guren et al. (2023) for pH 4. The values for pH 4 are expected to be higher because carbonate dissolve faster in acidic water. The samples they used also contained macrofeatures, e.g., macrosteps and crevices, obtained from mechanical cleaving of single crystals. Their probability distributions for pitted surfaces appear to be quite reproducible as well, taking into account possible differences in microfeature reactivity.

The second aspect is occurrence of surface morphologies with process-controlling pits formed around dislocation clusters and with competing etch pits of similar sizes (as two contrasting cases on the Figure 7). Both cases are common on dissolving surfaces especially for calcite which bonding topology is the same as for Kossel crystal. Formation of large pits around pit's clusters is commonly observed, e.g., in studies by Callagon et al. (2014), Noiriel et al. (2020), and Guren et al. (2023), as well as etch pit morphologies controlled by many pits, e.g., in Arvidson et al. (2003), Fischer et al. (2012), and Brand et al. (2017). If these two contrasting cases indeed result in different, but reproducible rate distributions in one case, and a finite set of reproducible distributions in the other, the issue of defect density and location becomes indeed solvable in the statistical sense.

4. Conclusions

The reproducibility of dissolution rates is the fundamental problem of practical importance for applications involving reactive transport modeling. Intrinsic rate variance issues in macroscopic models can be addressed by using probabilistic approaches in cases when statistical behavior of rates is well-known. Reproducibility and robustness of rate probability distributions in this scope becomes a fundamentally important question, where the influence of dislocation density need to be clarified. As we established in this study, at least for simple Kossel crystals, the numbers or density of dislocations has little relevance to probability distributions. Instead, their spatial locations, or geometric arrangements, have critical influence on overall reactive system dynamics and rate reproducibility. In this study we addressed rate variance generated by different surface morphologies stemming from various etch pit geometric superpositions. We discovered that in cases when multiple etch pits controls surface reactivity, a stable probability distribution can be expected. In cases when dislocations form clusters, formation of one giant macro-pit is expected. In these cases several possible stable distributions can be observed due to the possible existence of long-range correlations in movement of atomic step trains. This result has a remarkable implication for understanding rate variance issues and for quantitative predictions of dissolution rates for modeling of geochemical processes. In particular, long-term behavior of etch pit clusters of various geometric arrangements should be better understood, as well as possible rate distributions they generate for time trajectories of different length. The good news is that despite inherent randomness, the surface dynamics follows some clear trends predictable to some extent in the statistical sense. The issues of local rate variability vs. global temporal variance of material flux should be investigated in more detail. The role of mineral crystallographic structure and bonding in determining etch pit interactions and corresponding rate distributions should be carefully evaluated, as well as the influence of “extrinsic” parameters saturation state, pH, transport, leached layers on silicates, etc., on the kinematics of atomic steps and etch pits. Grain boundaries, crystal size, and crystal habit, including nanoparticulate forms, constitute another important source of rate variance relevant to practical applications. The dominant sources of rate variance should be determined before practical implementations of rate distributions.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

IK: conceptualization, programming, data collection and analysis, and manuscript writing. AL: conceptualization, manuscript editing, and proofreading. All authors contributed to the article and approved the submitted version.

Funding

We would like to acknowledge the Bundesministerium für Bildung und Forschung (BMBF) under Grant No. 03G0900B within the joint project ResKin_Move.

Acknowledgments

We thank Ricarda D. Rohlfs for her enthusiastic scientific discussions, as well as Cornelius Fischer and Rolf S. Arvidson for their discussions and inspiring previous publications on the rate variance issue.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Agrawal, P., Bollermann, T., Raoof, A., Iliev, O., Fischer, C., and Wolthers, M. (2021). The contribution of hydrodynamic processes to calcite dissolution rates and rate spectra. Geochim. Cosmochim. Acta 307, 338–350. doi: 10.1016/j.gca.2021.05.003