Estimating the activation energy of bond hydrolysis by time-resolved weighing of dissolving crystals

Abstract

Bond-breaking activation energy E_B is nowadays a key parameter for understanding and modeling crystal dissolution processes. However, a methodology to estimate E_B based on classical dissolution experiments still does not exist. We developed a new method based on the calibration of a Kossel type dissolution model on measured dissolution rates obtained by mass (or volume) variations over time. The dissolution model does not depend on the geometry of the crystal surface but only on the density of the different types of sites (kink, step, terrace, bulk). The calibration method was applied to different experimental setups (flow through and batch) with different ways of estimating the dissolution rates (solute concentration in the fluid, surface topography) for calcite crystals. Despite the variety of experimental conditions, the estimated bond-breaking activation energies were very close to each other (between 31 and 35 kJ/mol) and in good agreement with ab initio calculations.

Introduction

It has long been known that chemical reactions are thermally activated. In his pioneering work, Arrhenius¹ formalized his experimental observations to provide one of the most fundamental empirical relations in chemical kinetics, well known under the term ‘Arrhenius equation’. The corresponding formula, which describes the temperature dependence of reaction rates, introduced the concept of ‘activation energy’ (E_a), which may be viewed as the minimum amount of energy required to trigger a given chemical reaction:

$$k = \nu \exp ( - E_{\rm{a}}/RT)$$

(1)

where k is the rate constant at a given temperature [s⁻¹], $\nu$ is the pre-exponential or frequency factor [s⁻¹], which is specific to the considered chemical reaction, R is the gas constant [J/mol/K] and T is the absolute temperature [K].

It was only about 50 years later that the Arrhenius equation received a theoretical support for elementary reactions. Following from the transition state theory developed by Eyring², the activation energy was then shown to correspond to the enthalpy of the formation of the activated complex (or, in other words, to the difference between the binding energies involved in the formation of the activated complex and the binding energies of the reactants to be broken). As a consequence, transposing this relation to dissolution reactions, Lasaga³ suggested that the rate-limiting steps of mineral dissolution may be determined by comparing the temperature dependence of a dissolution reaction rate with that of elementary step reactions. The activation energy of dissolution of a wide range of minerals including oxides, halides, sulfates, carbonates or silicates has been determined experimentally since then, showing values ranging from a few kJ/mol to up to ~100 kJ/mol (see compilation in e.g.⁴).

In parallel, the advent of ab initio methods and their application to the field of Earth sciences/mineralogy from the early 1990s paved the way to the calculation of the activation energy of bond hydrolysis (or “bond-breaking activation energy”) from first principles. The energetics of bond hydrolysis (E_B) have been investigated for a wide range of minerals and structures, including oxides e.g⁵., halides e.g⁶., carbonates e.g⁷., silicates and in particular, tectosilicates^8,9,10,11, and orthosilicates^12,13. A synthetic overview of E_B values is provided in Table 1. Among the outcomes of these studies, some of them have tried to link the activation energies derived from ab initio calculations to the activation energies determined experimentally, concluding for instance that the rate-limiting step of all silicates was the cleavage of the last Si–O–Si bond linking an Si atom to the surface^14,9. The use of activation energies determined experimentally to infer reaction mechanisms has however been questioned by several studies, pointing out that the overall rate dependence on temperature actually stems from the contribution of two terms, i.e., the enthalpy of formation of the rate-controlling surface sites, and the activation energy of the activated complex itself¹⁵. As a consequence, the measured activation energy is often referred to as ‘apparent’ or ‘bulk activation energy’, resulting from the sum of the enthalpy of proton adsorption and that of bond hydrolysis¹⁶. Whether or not direct insights into the energetics of bond hydrolysis can directly be retrieved from the measurement of ‘apparent activation energies’ has therefore remained an open question.

Table 1 Theoretical bond-breaking activation energies of bonds encountered in common minerals.

The determination of bond-breaking activation energies from first principles has received renewed interest in recent years due to the development of a generation of stochastic dissolution models at the atomic scale, where the bond-breaking probability is scaled to the bond-breaking activation energy (e.g.^{17,18,19,20,21,22,23,24,25,26,27,28,29,30,31}). Over the last decade, this generation of models has gradually superseded the conventional approach describing mineral dissolution kinetics using closed-form rate equations inherited from the early work of Aagaard and Helgeson³², for which the theoretical basis has been regularly questioned (see e.g^33,34). The main interests and successes of atomic-scale stochastic dissolution models are manifold and include: (i) a better account for the intrinsic reactivity of minerals, including rate variability^17,21,22 and dissolution anisotropy^35,26; (ii) a fine description of topographical features such as etch pit morphology^36,19,27 and evolution of grain morphology^35,28, (iii) the clarification of reaction mechanisms, such as the formation of surface layers on dissolved silicates (e.g.³⁷) and the impact of the saturation state on etch pit nucleation, step-wave propagation and reaction rates (e.g.^35,36). Ultimately, stochastic dissolution models might play an important role for the upscaling of reaction rates and the development of a next generation of reactive transport codes with a stronger mechanistic basis²³. This is of interest for a wide range of fields of application dealing with materials degradation, ranging from metal corrosion to chemical weathering or geological storage of CO₂ and nuclear waste.

In this paper, we present an original method to estimate E_B values from the monitoring of the dissolution of crystals. This method can also be applied to check the quality of dissolution parameters resulting from ab initio simulations. The developed methodology assumes that the crystals are of Kossel type, and consists in calibrating the corresponding dissolution model parameters which includes E_B to fit the output of dissolution experiments. These calibrations are applied to a wide range of E_B using dissolution rates computed by a Kinetic Monte Carlo (KMC) method to evaluate the methodology reliability. This method for evaluating E_B is then applied to numerous calcite dissolution experiments where dissolution rates were estimated at fixed saturation indices through the monitoring of cation release in the aqueous fluid or by nanoscale topography measurements at the crystal surface.

Results and discussion

Main observations drawn from the simulations

Several numerical experiments were performed to grasp the main characteristics associated to the dissolution of Kossel crystals. The numerical experiments were performed using the following set-up: a cube of different sizes (from 50 × 50 × 50 to 1000 × 1000 × 1000 sites), a frequency factor of 1.0 × 10¹² s⁻¹, and a temperature of 295 K. The bond-breaking activation energies vary from 5.0 to 100.0 kJ/mol (see Table 1). This range of activation energies leads to a very wide range of bond-breaking probabilities (Table 2).

Table 2 Activation energies and corresponding bond-breaking probabilities used in the simulations.

In the literature, variables like atom release or solute concentrations are presented versus time or versus reaction advancement (or dissolution progress) defined by the cumulated number of removed sites divided by the total number of initial sites. When possible, even if reaction advancements are very convenient for comparing simulations, we will avoid using dissolution progress because results interpretation may be quite difficult and eventually misleading. Actually, the relation between the dissolution progress and time is not always linear and depends on the activation energy (Fig. 1). In the provided example, the first 10% of the dissolution requires about 5% of the time needed to dissolve the entire crystal for E_B= 5 kJ/mol, and 23% for E_B greater than 15 kJ/mol. For E_B greater than 15 kJ/mol, the linear relationship between time and dissolution progress is valid only for a dissolution progress ranging between 20% and 80%.

Fig. 1: Time evolution as a function of the dissolution progress.

N_r is the cumulated number of dissolved sites, N₀ the initial number of sites and t_end, the simulated time necessary to dissolve the entire crystal.

The differences in E_B have also a significant impact on the evolution of the crystal geometry during dissolution. The geometry remains unchanged in average for low bond-breaking activation energies (lower than 15 kJ/mol) whereas it significantly changes for higher bond-breaking activation energies (Fig. 2).

Fig. 2: Evolution of the crystal geometry during dissolution.

Images (a) and (c) for E_B = 5 kJ/mol and (b) and (d) for E_B = 15 kJ/mol. Images (a) and (b) correspond to a reaction advancement of 10%, and images (c) and (d) to 50%. Only 1/8 of the crystal is represented. The white solid line is the initial size of such a one eight of the crystal.

The evolution of the geometry shows that, for a crystal of this size, the contributions of the terrace, step and kink sites to dissolution are quite similar, which explains the more or less unchanged geometry for E_B = 5 kJ/mol. For E_B ≥ 15 kJ/mol, the geometry of the entire crystal evolves from a cubic to an octahedral shape. Kink sites are the predominant sites observed at the surface.

In order to evaluate whether deriving bond-breaking activation energy from a series of (blind) numerical simulations (or from dissolution experiments) is doable, we performed a first order sensitivity analysis by studying the effects of activation energy on the surface topography. The results underlined by the morphology evolution (i.e., the observation of two distinct dissolution regimes depending on whether E_B is smaller or greater than 10–15 kJ/mol) are confirmed by a more detailed analysis of the contribution of each type of sites located at the crystal surface. This contribution is quantified by n_iPⁱ and shown in Fig. 3. Interestingly, the respective contribution of each n_iPⁱ value does not depend on the crystal size, at least for the sizes analyzed in this work.

Fig. 3: Contribution of the different types of sites to the crystal dissolution rate.

N_n is the number of nearest neighbors: N_n = 3 for kinks, N_n = 4 for steps, N_n = 5 for terraces. Random fluctuations are due to the stochastic nature of the KMC method.

The sites with less than 3 bonds do not contribute appreciably to the dissolution because they are too few in number. Consistently with previous studies (e.g.^17,24 and references therein²⁹, for nanosized grains), kink sites are by far the main contributors. They are the only contributors for E_B ≥ 15 kJ/mol. Step sites have a significant contribution for low E_B values (around 20% for E_B = 5) and a small contribution for E_B = 10. Terrace sites contribute for less than 5% of the total dissolution rate, and for low E_B values only. Interestingly, the different contributions reach steady-state quite rapidly for E_B = 5 and 15 kJ/mol (in less than 10% of the time required to dissolve 90% of the crystal), meaning that the relative number of sites of different types remains constant over most of the dissolution process. For intermediate E_B values (10 kJ/mol), the evolution of the number of sites requires significantly more time to reach this steady state. Moreover, significant oscillations occur during the first stages of the dissolution until equilibrium between the different types of sites is reached.

Simulations also showed that the contribution of sites with 1 or 2 neighbors to crystal dissolution can be neglected. These results can be used to define the relationship between E_B and the bulk or apparent activation energy, E_A. We call the apparent activation energy the energy required to have a dissolution rate r_M without taking care of the various sites (kink, step, terrace) at the surface (which is often the case for classical dissolution experiments),

$$r_{M} = m\;\nu n_S^\prime {\rm{e}}^{\left( { - E_{\rm{A}}/RT} \right)}$$

(2)

where $n_S^\prime$ is the total relative number of surface sites.

Numerical simulations with different E_B values showed that $n_S^\prime\simeq n_3^\prime + n_4^\prime + n_{^5}^\prime$ and Eq. (13) can be rewritten as

$$r_{M} = m\;\nu \kappa {\rm{e}}^{\left( { - 3E_{\rm{B}}/RT} \right)}$$

(3)

with $\kappa = n_3^\prime + n_4^\prime p_{\rm{B}} + n_5^\prime p_{\rm{B}}^2$ where p_B is the bond-breaking probability (see Table 2) and $0\, <\, \kappa, < \,n_s^\prime$. $\kappa$ is the sum of the different site contributions to dissolution. Depending on the E_B value and the shape of the crystal surface, we can have $\kappa /n_3^\prime \approx 1.0$ meaning that only kink sites contribute to dissolution. Considering the value of bond-breaking probabilities (Table 2), kink sites are the main contributor for high E_B values.

Combining Eqs. (2) and (3) leads to

$$E_{A} = 3E_{B} - RT\;{Ln}\left( {\frac{\kappa }{{n_S^\prime }}} \right)$$

(4)

Due to the property of $\kappa$, we have

$$E_{A} \ge 3E_{B}$$

(5)

We recall here that from physical and practical standpoints, the bulk activation energy measured experimentally is linked to the macroscopic activation energy (E_A) through the additional contribution of, e.g., the enthalpy of proton adsorption in the acidic pH range¹⁶, which is neglected here. If the enthalpy of adsorption of the reactive species is known, the apparent activation energy derived experimentally may be used to extract E_A and ultimately, to define an upper limit for the bond-breaking activation energy estimation for a Kossel type dissolution model.

Estimation of bond-breaking activation energies and frequency factor parameters

Dissolution rates from KMC simulations screening E_B from 5 to 100 kJ/mol are used as measured values. The robustness and accuracy of the method are evaluated by comparing the estimated parameters (bond-breaking activation energy and frequency factor) with the parameters used for the corresponding simulations.

The boundary values for bond-breaking activation energies are set to [1,30] kJ/mol for expected values lower than 20 kJ/mol and to [10,130] kJ/mol for expected values greater than 10 kJ/mol. These constraints are set to avoid potential impact of a priori knowledge on the optimization procedure. Of course, if a priori knowledge exists, the constraints can be refined by reducing the min-max interval. The frequency factor $\nu$ was set to 10¹² s⁻¹ for all simulations, the minimum value to 10¹¹ s⁻¹, and the maximum to 10¹³ s⁻¹. The sites with one and two bonds are not considered because their contribution to the dissolution is negligible (see Fig. 3), which reduces the number of estimated parameters to 6 (frequency factor, E_B, and the number of kink, step, terrace and bulk sites). Numerical simulations also show that the concentration of bulk sites is greater than 0.95, even 0.99 for high activation energies. A lower limit of 0.90 was prescribed for n′₆.

The bond-breaking activation energies can be estimated whatever the time step or the dissolution progress. 100 time steps regularly distributed over a reaction progress ranging from 20% to 90% are used to estimate the parameters. Table 3 indicates the E_B values used to simulate the crystal volume evolution over time, the average estimated value $\left\langle {E_{\rm{B}}} \right\rangle$ at 100 different times, and the corresponding standard deviation $\sigma _{E_{\rm{B}}}$. The ratio $\sigma _{E_{\rm{B}}}/\left\langle {E_{\rm{B}}} \right\rangle$ provides a first order evaluation of the estimation accuracy. The same type of outputs is provided for the frequency factor ν.

Table 3 Parameter estimation results for activation energies and frequency factors.

The optimization procedure provides a quite accurate estimation of the bond-breaking activation energy and the frequency factor. The variation coefficient $\sigma _{E_{{B}}}/\left\langle {E_{{B}}} \right\rangle$ indicates that the accuracy for E_B = 5 kJ/mol is smaller than the others. A detailed analysis of the numerical data shows that sites with one or two neighbors contribute slightly to dissolution for that particular case and these sites are not taken into account in the optimization process.

The accuracy of the estimated decimal logarithm of the frequency factor is independent of the experimental conditions (E_B values). In average, the estimation of the frequency factor is less accurate than the estimation of E_B values (Table 3). A more detailed analysis is depicted in Figs. 4, 5. There is no trend in the estimated values of E_B and ν for the 100 different time steps (samples). The quality of both estimations depends on neither the E_B value nor the dissolution stage.

Fig. 4: Estimated E_B values at different steps of the dissolution process.

Values indicated in the legend are the exact values.

Fig. 5: Estimated log10 frequency factor values at different steps of the dissolution process.

Colors represent the different E_B values. The value used for the simulations is ν = 10¹² s⁻¹.

Bond-breaking activation energy estimation based on laboratory experiments

Laboratory experiments dedicated to the estimation of dissolution rates are based on flow through or batch reactors associated with the analyses of the chemical composition of the water and/or the topography of the sample surface by, for example, atomic force microscopy (AFM) or vertical scanning interferometry (VSI).

The feasibility and robustness of our approach were tested for these kinds of laboratory experiments studies dedicated to the dissolution of calcite. Calcite (CaCO₃) represents an ideal target in that respect, because (i) its reactivity can be described following a Kossel geometry²⁹, (ii) bond-breaking activation energies have been previously estimated following ab initio methods (Table 1) and (iii) the dissolution kinetics of calcite have been extensively studied in the literature.

Amongst the numerous published experiments, we selected experimental conditions consistent with the assumptions made in the dissolution model, i.e.:

• neutral or basic conditions where pH variations do not affect the dissolution rate³⁸;

• saturation index greater than 0.40, where the formation of etch pits at the crystal surface is hindered³⁹ and where the dissolution rate can be assumed to vary linearly with the saturation index^40,41.

The experiments from the following studies were selected:

• Smith et al.⁴⁰, corresponding to batch experiments where both solute concentration in the reactor and VSI data were used for estimating the dissolution rates;

• Cubillas et al.⁴² and Xu et al.⁴³, corresponding to flow through reactor experiments, where the dissolution rates were estimated using the outlet solute concentration;

• Bouissonnié et al.⁴¹, corresponding to flow-through reactor experiments, where dissolution rates were estimated using VSI.

Some experimental conditions are summarized in Table 4.

Table 4 General experimental conditions (V: reactor volume, T: fluid temperature, S: specific surface area).

The experiments were performed by changing parameters such as the flow rate and/or the saturation index. Due to model assumptions, we selected only experiments with provided saturation index higher than 0.40. We kept the naming from Smith et al.⁴⁰ for their experiments (CDE2 and CDE3). Only one experiment (named X1) from Xu et al.⁴³ was selected due to the lack of information on the specific surface areas. The experiments from Bouissonnié et al.⁴¹ are named AB1 to AB5. Smith et al.⁴⁰ analyzed the dissolution of two planes ({104} and {001}) of the crystal. Only data related to the {104} plane were used (experiments D, E and F). The {001} face was not taken into account as calcite can be represented by a Kossel crystal (for which faces represent flat/terrace surfaces) only if the crystal habit is shaped with {104} faces.

In order to take into account the impact of the saturation index of the solution, the dissolution rate constant of calcite (k) was calculated based on the assumption that the rate obeys the transition state theory, which is reasonable in this range of saturation states (see⁴⁰ or⁴¹):

$$k = r/\left( {1 - {\Omega}} \right)$$

(6)

with Ω the saturation index.

Data related to the estimation of dissolution rates based on the calcium concentration at the outflow are summarized in Table 5. These experiments differ in the initial mass of calcite and the residence time in the reactor defined by the ratio between the volume of the fluid in the reactor and the injection rate. When the fluid volume was not provided, we assume that it is equal to the volume of the reactor. The kinetic rate constants k are very close for the same type of experiments but vary by a factor close to 5 between the highest and smallest values, considering all data. Batch experiments provide significantly lower values than flow-through experiments when dissolution rates are estimated using solute concentrations in the fluid.

Table 5 Experimental conditions with m the initial mass of calcite, q the flow rate and Res the residence time of the solution in the reactor.

Data related to the estimation of dissolution rates based on surface retreat measured by VSI are given in Table 6. The rate constants are quite similar (except for experiment F) despite different experimental setups: flow-through reactor for Bouissonnié et al.⁴¹ and batch reactors for Smith et al.⁴⁰.

Table 6 Experimental conditions with dh/dt the retreat velocity, q the flow rate and Res the residence time of the solution in the reactor.

On average, the dissolution rates based on the solute concentration are higher than dissolution rates based on the surface retreat. For calcite, this is in agreement with studies of Arvidson et al.⁴⁴ and Noiriel et al.⁴⁵.

The first optimization was performed considering that the frequency factor is the same for all experiments. The estimated value is 11.96 ± 0.50 for $\left\langle {{{{\mathrm{log}}}}_{10}\left( v \right)} \right\rangle \pm \sigma _{{{{\mathrm{log}}}}_{{{{\mathrm{10}}}}}(v)}$ where $\left\langle {{{{\mathrm{log}}}}_{10}\left( v \right)} \right\rangle$ is the average value of the decimal logarithm of the frequency factor and $\sigma _{{{{\mathrm{log}}}}_{10}\left( v \right)}$ is the corresponding standard deviation. This value is estimated with good accuracy and is consistent with the value computed by Eq. (8), v = 6.15 × 10¹²s⁻¹. We therefore prescribed in the following this value for the frequency factor to limit its effects on the other estimated parameter values.

Bond-breaking activation energies were estimated for each experiment. More than 99% of the optimizations reached convergence, i.e., the relative difference between measured dissolution rate and estimated dissolution rate was smaller than 0.01%. The estimated bond-breaking activation energy values lie between 30.0 and 34.3 kJ/mol (Fig. 6). Data from Xu et al.⁴³ and VSI are very close (between 32.2 and 33.6 kJ/mol).

Fig. 6: Estimated bond-breaking activation energy for the 16 different experimental dissolution rate values.

Average values are computed from the 200 optimizations (symbol) as the corresponding standard deviation (bars). Experiment numbers (x-axis) are given in Table 5, 6.

The estimated bond-breaking activation energies are similar for the different experiments despite differences in experimental conditions (flow rate, grain size and crystal pre-treatment), in the methodology (monitoring of the chemical composition of the water or topography analyses of the crystal surface by VSI) and the advancement of the dissolution (crystal geometry). Similar values were obtained by Raiteri et al.⁷—(see Table 1), close to the lower bound of 20 kJ/mol, in agreement with Kurganskaya and Luttge²⁰ who consider that the lowest values of the range given by Raiteri et al.⁷ are the most probable.

Considering the consistent values of the bond-breaking activation energy, we also analyzed the surface geometry obtained for each optimization (relative number of kink, step, terrace and bulk sites, see Table 7). The distributions of the different sites on the calcite surface are very similar for all experiments, except for experiments 6–8, where the dissolution process was monitored through solute concentration in batch experiments. The uncertainty in the estimation of the concentration of each surface site is quite large (about 50%) due to the number of degrees of freedom in the optimization procedure.

Table 7 Average and related standard deviation for the solid concentration of the different sites.

For an average value of E_B equals 34 kJ/mol and a fluid temperature of 20 °C, the bond-breaking probability is 8.77 × 10⁻⁷. The surface site densities listed in Table 7 are used to compute the ratio $\kappa /n_3^\prime$ (see Eq. (3)) which is very close to 1.0 for all experiments. The kink sites can therefore be considered as the sole contributors to dissolution.

Dissolution rates estimated through VSI data are determined on a very small crystal surface compared to dissolution rates estimated by solute concentration in the reactor. One could consider that VSI data are not appropriate, because of the very small crystal surface area investigated, which may not be representative of the total crystal surface. Our results demonstrate the opposite: the estimated bond-breaking activation energy is consistent with the other estimated values (see Fig. 6 and Table 7). Assuming that VSI data are representative of the dissolution rate, the differences in the estimated rates is due to the crystal surface geometry of the surface explored by the VSI and the surface of the crystal in contact with the fluid, as underlined by the bulk concentration differences in the different sites (see Table 7).

Results summary

A methodology for estimating the bond-breaking activation energy E_B and frequency factor of crystal dissolution has been developed and tested for Kossel type crystals. It consists of an optimization procedure that aims at minimizing the differences between measured and modeled dissolution rates. The measured dissolution rates are obtained following mass (or volume) variations of the crystal over time, which are common data in experimental studies. The modeled dissolution rate is based on Kossel assumptions and requires 8 parameters, i.e., the frequency factor, E_B and the number of different types of sites describing the crystal. The dissolution model does not depend on the geometry of the crystal surface but only on the density of the different types of sites. It allowed also providing a link between the bulk activation energy and E_B.

A first set of numerical simulations were performed over a wide range of E_B (from 5 to 100 kJ/mol) to explore possible simplifications of the dissolution model. This sensitivity analysis of the dissolution rates to E_B showed that:

• 1. The evolution of the crystal geometry and the number of the various sites (especially kink sites) during dissolution observe two different regimes depending on a threshold value for the bond-breaking activation energy of about 15 kJ/mol.

• 2. The kink sites are the main contributors to the dissolution in all cases, and almost the only contributors to dissolution when the bond-breaking activation energy is greater than 15 kJ/mol.

• 3. Whatever E_B, the contribution of sites with less than 3 bonds can be neglected, which reduces the number of parameters to estimate to 6.

The methodology was applied to different experimental setups (flow-through and batch experiments) with different methods for estimating the dissolution rates (solute concentration in the fluid, surface topography) of calcite. Despite the differences in experimental conditions and methods used to estimate the dissolution rates, the following results were obtained:

• 1. The frequency factor defined by its formulation based on Boltzmann and Planck constants and the temperature can be used to analyze experiments dedicated to mineral dissolution kinetics.

• 2. Estimated bond-breaking activation energies lie in a very narrow range (between 31 and 35 kJ/mol), whatever the considered study. These values are in good agreement with ab initio calculations.

• 3. The densities of the different sites (kink, step, terrace, bulk) were also very similar for all saturation indices higher than 0.4, and the kink sites are the main contributors to dissolution.

This method should also be considered as a way to evaluate the reliability of ab initio calculation or as an alternative way of estimating E_B values. It is also easier to use and does not account of complicated physics involved in ab initio calculation. Moreover, the estimation is performed for the actual crystal geometry, on the contrary to quite numerous ab initio studies.

Methods

Stochastic simulation of crystal dissolution

We assume that a crystal can be described as a Kossel-type crystal⁴⁶ also called Terrace Ledge Kink system (TLK-⁴⁷) i.e. a compact stack of sites representing a reactive unit (atom, molecules,…) linked to their neighbors by chemical bonds. This very simple geometry allows handling topographical details such as flat surfaces or terraces (sites with 5 bonds), ledges or steps (sites with 4 bonds) and kinks (sites with 3 bonds). Sites with 6 bonds are not in contact with the fluid and named bulk sites²⁹.

The mathematical model describing the dissolution kinetic is developed in a stochastic framework and the numerical solution is obtained using a Kinetic Monte Carlo (KMC) method. Based on initial ideas from¹ and⁴⁸, the Transition State Theory (TST) provides a simple expression for the rate constant, known as Eyring (or Eyring–Polanyi) equation²:

$$k_{ij} = \nu \exp \left( { - \frac{{{\Delta}E_{ij}}}{{RT}}} \right)$$

(7)

where K_ij is the rate constant [s⁻¹], ΔE_ij the activation barrier of the process [J/mol], T the temperature [K] and R the gas constant [J/mol/K]. The term $\nu$ is called the pre-exponential factor or frequency factor [s⁻¹]. This frequency factor is expected to be a reasonable approximation for the number of reaction attempts in one unit of time during the chemical reactions at the mesoscale. It can be defined by

$$\nu = \frac{{k_{\rm{B}}}}{\hbar }T$$

(8)

where k_B is the Boltzmann constant, $\hbar$ the Planck constant and T the temperature^49,50. For T = 295 K, the corresponding value is v = 6.15 × 10¹² s⁻¹. The value usually adopted for the characteristic frequency (or pre-exponential frequency) is 10¹² s⁻¹, the intermolecular vibrational frequency of bulk water⁵¹. Other values can be found in the literature such as 5.22 10¹⁰ s^{−1 52} or values ranging from 1 to 8 × 10¹⁰ s⁻¹ estimated by model calibration on AFM data obtained during calcite dissolution experiments²⁹.

Dissolution is described by a solid-on-solid approach and the dissolution rate is defined by⁵³

$$r_n = \nu \mathop {\prod}\limits_{i = 1}^n {p_{B}} = \nu p_{{B}}^n\;{{{\mathrm{with}}}}\;p_{{B}} = \exp \left( { - \frac{{E_{{B}}}}{{RT}}} \right)$$

(9)

where r_n is the dissolution rate of sites with n bonds, n is the number of bonds (or the number of nearest first neighbors also called coordination number) of the site, p_B is the bond-breaking probability. In first order, it is assumed that the probability of breaking one bond is independent of the number of bonds.

KMC requires a rigorous estimation of the duration of a single ‘jump’ since this stochastic approach is also an upscaling of time. Time events at the atomic scale are not modeled explicitly and the ‘macroscopic’ time step has to be defined properly. It has been shown that the macroscopic time step follows an exponential distribution (assuming that the time step is ‘sufficiently’ small—even if sufficiently is not well defined) with average⁵⁴:

$$\left\langle {{\Delta}t} \right\rangle = 1/\mathop {\sum}\limits_{j = 1}^5 {\nu n_jp_{\rm{B}}^j}$$

(10)

The computation of a time step τ_i for the ith iteration requires the generation of a random exponential distribution with an average of $\left\langle {{\Delta}t} \right\rangle = 1/k_{ij}$. This is obtained by

$$\tau _i = - {\rm{L}}n(Z)/\mathop {\sum}\limits_{j = 1}^5 {\nu n_jp_{{B}}^j}$$

(11)

with Z a random number with a uniform distribution between [0,1]. However, for sufficiently long simulation time, the generation of random time step can be avoided and the average time step can be used⁵⁵.

Bond-breaking activation energy estimation

In this section, we describe the method we followed for estimating the bond-breaking activation energy based on knowledge of the dissolution rate of a given single mineral determined experimentally. Theoretical developments and E_B estimation methodology are based on the following assumptions (Kossel dissolution model):

• the crystal can be described by an ensemble of identical sites of cubical shape;

• the sites are linked through their common face (maximum of 6 bonds);

• only the sites in contact with the fluid (with less than 6 bonds) can be dissolved;

• a site is removed from the crystal when the 6 bonds are broken;

• for a given site, the probability of dissolving a bond is independent of the site, of its location and of the number of bonds that links this site to its neighbors.

These assumptions concerning the solid are completed by assumptions concerning the fluid: The effect of the enthalpy of proton adsorption is assumed to be negligible, corresponding to experimental conditions close to the Point of Zero Net Proton Charge¹⁶.

The crystal mass depends on the total number of sites inside the crystal and at its surface:

$$m = \rho V = \rho n_T\ell ^3 = \rho (n_s + n_v)\ell ^3$$

(12)

where ρ is the crystal density [M/L³], $\ell$ the edge length of a single site [L], n_T is the total number of crystal sites, $n_s = \mathop {\sum}\nolimits_{i = 1}^5 {n_i}$ is the number of sites in contact with the fluid, and n_v = n₆ the number of bulk sites.

The dissolution rate of the crystal r_M [M/T] depends on the number of sites located at the crystal surface and their dissolution rates and defined by

$$r_M = - \frac{{{\rm{d}}m}}{{{\rm{d}}t}} = \nu \rho \ell ^3\mathop {\sum}\limits_i^5 {n_ip_{{B}}^i} = \nu \rho \left( {n_T\ell ^3} \right)\mathop {\sum}\limits_i^5 {\frac{{n_i}}{{n_T}}p_{{B}}^i} = m\;\nu \mathop {\sum}\limits_i^5 {n_i^\prime p_{{B}}^i}$$

(13)

with n_i the number of sites with i bonds (i lower than 6 because the site is located at the crystal surface), n_T is the total number of sites and $p_{{B}}^i$ the probability of breaking the i bonds. $n_i^\prime\, = \,n_i/n_T$ is the relative number (or solid concentration) of sites having i neighbors.

The estimation of E_B is considered as a non-linear optimization problem with constraints. It consists in estimating the relative number of different sites $\left( {n_i^\prime ,\;i = 1..6} \right)$, E_B, and the frequency factor $\nu$. This optimization problem is thus defined by

$$\left\{ {\begin{array}{l} {{{min}}\left[ {\left( {\hat r_M - m\;\nu \mathop {\sum}\limits_i^5 {n_i^\prime p_{{B}}^i} } \right)^2} \right]} \\ {\mathop {\sum}\limits_{i = 1}^6 {n_i^\prime } = 1.0} \\ {0 \le n_i^\prime \le 1} \\ {E_{{{B}},\min } \le E_{{B}} \le E_{{{B}},\max }} \\ {\nu _{\min } \le \nu \le \nu _{\max }} \end{array}} \right.$$

(14)

The first equation of the system is the objective function for matching r_M, the computed, and $\hat r_M$, the measured, dissolution rates. The other equations of the system are the associated constraints. The sum of the solid site concentrations has to be equal to 1 and each concentration is smaller than one by definition. The E_B,min and E_B,max (respectively v_min and v_max) parameters are a priori lower and upper limits of the bond-breaking activation energy (resp. the frequency factor). The solution of this minimization problem is obtained using a standard Particle Swarm Optimization (PSO) method⁵⁶. The algorithm starts with numerous initial solutions (in our case, an initial solution is a set of n′_i, E_B and $\nu$ values chosen randomly within the ranges as defined by Eq. (14)). Each solution is called a particle. Particles are moved randomly while taking care of their best solution (the smallest value of $\left( {\hat r_M - r_M} \right)^2$ for this particle, called local best) and the best value of $\left( {\hat r_M - r_M} \right)^2$ for all particles (called the global best). Moving a particle consists in changing the values of n′_i, E_B and $\nu$ following:

$$\left\{ {\begin{array}{*{20}{l}} {\psi _{k + 1}^j = \psi _k^j + {\Delta}\psi _k^j} \hfill \\ {{\Delta}\psi _k^j = \omega _1{\Delta}\psi _k^{j - 1} + {\rm Z}_1\omega _2(\psi _{k,{{pbest}}}^j - \psi _k^j) + {\rm Z}_2\omega _3(\psi _{{{gbest}}}^j - \psi _k^j)} \hfill \end{array}} \right.$$

(15)

where $\psi _k^j$ is the jth parameter of particle k, ω_i (i = 1...3) are user’s defined constants, Z_i (i = 1, 2) are random numbers uniformly distributed over [0,1], $\psi _{k,{{pbest}}}^j$ is the local best value and $\psi _{{{gbest}}}^j$ is the global best value. ω_i is a weight that allows for balancing the influence of the previous particle move (i = 1), the local best solution (i = 2) and the global best solution (i = 3). Numerous alternative algorithms exist⁵⁷. We chose this one for its simplicity. The iterative algorithm is stopped either when $\left( {\hat r_M - r_M} \right)^2 < \tau$, where τ is a user-defined tolerance, or when a user-defined maximum number of iterations is reached. 200,000 particles are used to solve the optimization problem and the tolerance is fixed to $\tau = 10^{ - 8}\hat r$ to take into account the wide range of reactions rates. The maximum number of iterations is set to 1000.

Many parameters have to be estimated (bond-breaking activation energy, frequency factor, relative number of sites of coordination 3 to 6) and the solution of the optimization problem will not be unique. Particle Swarm Optimization allows to explore the space of optimal parameters and estimating the parameter uncertainties by running the algorithm several times for the same measured data, due to its random components (see Eq. (15)). The optimization procedure was solved 200 times leading to 200 different parameter sets for each dissolution rate. This number of optimizations allows reaching stable average and standard deviation of estimated parameters values, i.e., these statistical values do not change with additional optimization runs.

Notice also that the estimated dissolution rate is proportional to the frequency factor for a given number of sites and a given bond-breaking energy. This proportionality may make the estimation of the frequency factor difficult, but possible because the frequency factor and E_B are time-independent for a given experiment.

Application of the method to a real case study: calcite dissolution

The mass loss over time, required for the optimization (see Eq. (14)), is very seldom provided in publications dedicated to measurements of mineral dissolution kinetics. Cubillas et al.⁴² provided flow rates and solute concentration, which allows computing $\hat r_M$ the measured dissolution rates following:

$$\hat r_M = qC$$

(16)

where q is the flow rate [L³/T] and C the solute concentration [M/L³]. We used the initial mass to compute r_M, the simulated dissolution rate (see Eq. (13)).

The flow-through experiments reported in Xu et al.⁴³ and Smith et al.⁴⁰ provided dissolution rates expressed in mol/m²/s. These dissolution rates were multiplied by the surface area given in their paper and by the molar mass of calcite to compute r_M.

VSI observations provide the surface retreat over time and we can write

$$\left\{ {\begin{array}{*{20}{l}} {\hat r_M = - \rho S\frac{{{\rm{d}}h}}{{{\rm{d}}t}}} \hfill \\ {r_M = \rho Sh\;\nu \mathop {\sum}\limits_i^5 {n_i^\prime p_{{B}}^i} } \hfill \end{array}} \right.$$

(17)

where S is a reference surface area and h the surface retreat. The minimization of the quadratic difference between measured and computed dissolution rates does not depend on S and ρ, and the first equation in (14) can be rewritten as

$${\rm{min}}\left[ {\left( { - \frac{{{\rm{d}}h}}{{{\rm{d}}t}} - h\;\nu \mathop {\sum}\limits_i^5 {n_i^\prime p_{{B}}^i} } \right)^2} \right]$$

(18)