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.

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.

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

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):

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:

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.

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.

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

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.