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

Bond-breaking activation energy EB is nowadays a key parameter for understanding and modeling crystal dissolution processes. However, a methodology to estimate EB 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 1 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: where k is the rate constant at a given temperature [s −1 ], ν is the pre-exponential or frequency factor [s −1 ], 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 2 , 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 3 suggested that the ratelimiting 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. 4 ).
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 5 ., halides e.g 6 ., carbonates e.g 7 ., 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 15 . 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 16 . 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.
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 32 , 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 1 (e.g. 37 ) 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 23 . 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 2 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 setup: a cube of different sizes (from 50 × 50 × 50 to 1000 × 1000 × 1000 sites), a frequency factor of 1.0 × 10 12 s −1 , 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).
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%.
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).
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 i P i and shown in Fig. 3. Interestingly, the respective contribution of each n i P i value does not depend on the crystal size, at least for the sizes analyzed in this work.
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 29 , 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 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), where n 0 S is the total relative number of surface sites. Numerical simulations with different E B values showed that n 0 S ' n 0 3 þ n 0 4 þ n 0 5 and Eq. (13) can be rewritten as where p B is the bond-breaking probability (see Table 2) and 0 < κ; < n 0 s . κ 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 κ=n 0 3 % 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 Due to the property of κ, we have 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 16 , 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 ν was set to 10 12 s −1 for all simulations, the minimum value to 10 11 s −1 , and the maximum to 10 13 s −1 . 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′ 6 .
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 E B h i at 100 different times, and the corresponding standard deviation σ EB . The ratio σ EB = E B h i provides a first order evaluation of the estimation accuracy. The same type of outputs is provided for the frequency factor ν.
The optimization procedure provides a quite accurate estimation of the bond-breaking activation energy and the frequency factor. The variation coefficient σ EB = E B h i 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)

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 3 ) represents an ideal target in that respect, because (i) its reactivity can be described following a Kossel geometry 29 , (ii) bond-breaking activation energies have been previously estimated following ab initio methods (Table 1)  The experiments from the following studies were selected: Smith et al. 40 , corresponding to batch experiments where both solute concentration in the reactor and VSI data were used for estimating the dissolution rates; Cubillas et al. 42 and Xu et al. 43 , corresponding to flow through reactor experiments, where the dissolution rates were estimated using the outlet solute concentration; Bouissonnié et al. 41 , corresponding to flow-through reactor experiments, where dissolution rates were estimated using VSI.
Some experimental conditions are summarized in Table 4. 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. 40 for their experiments (CDE2 and CDE3). Only one experiment (named X1) from Xu et al. 43 was selected due to the lack of information on the specific surface areas. The experiments from Bouissonnié et al. 41 are named AB1 to AB5. Smith et al. 40  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 40 or 41 ): 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.   Grain size (mm) V (cm 3 ) T (°C) S (m 2 /g) Cubillas et al. 42 1.0-1. Bouissonnié et al. 41 3.0-6.0 50 22 -Smith et al. 40 10.0 40 20 - 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 flowthrough experiments when dissolution rates are estimated using solute concentrations in the fluid. 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. 41 and batch reactors for Smith et al. 40 .
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. 44 and Noiriel et al. 45 .
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 log 10 v ð Þ h i± σ log 10 ðvÞ where log 10 v ð Þ h iis the average value of the decimal logarithm of the frequency factor and σ log 10 v ð Þ 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 12 s −1 . 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. 43 and VSI are very close (between 32.2 and 33.6 kJ/mol).
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. 7 -(see Table 1), close to the lower bound of 20 kJ/mol, in agreement with Kurganskaya and Luttge 20 who consider that the lowest values of the range given by Raiteri et al. 7 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.
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 −7 . The surface site densities listed in Table 7 are used to compute the ratio κ=n 0 3 (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.

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.

Stochastic simulation of crystal dissolution
We assume that a crystal can be described as a Kossel-type crystal 46 also called Terrace Ledge Kink system (TLK-47 ) 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 29 .
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 1 and 48 , the Transition State Theory (TST) provides a simple expression for the rate constant, known as Eyring (or Eyring-Polanyi) equation 2 : where K ij is the rate constant [s −1 ], ΔE ij the activation barrier of the process [J/mol], T the temperature [K] and R the gas constant [J/mol/K]. The term ν is called the pre-exponential factor or frequency factor [s −1 ]. 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 where k B is the Boltzmann constant, h the Planck constant and T the  49,50 . For T = 295 K, the corresponding value is v = 6.15 × 10 12 s −1 .
The value usually adopted for the characteristic frequency (or pre-exponential frequency) is 10 12 s −1 , the intermolecular vibrational frequency of bulk water 51 . Other values can be found in the literature such as 5.22 10 10 s −1 52 or values ranging from 1 to 8 × 10 10 s −1 estimated by model calibration on AFM data obtained during calcite dissolution experiments 29 . Dissolution is described by a solid-on-solid approach and the dissolution rate is defined by 53 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 54 : The computation of a time step τ i for the ith iteration requires the generation of a random exponential distribution with an average of Δt h i ¼ 1=k ij . This is obtained by 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 55 .

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 16 .
The crystal mass depends on the total number of sites inside the crystal and at its surface: where ρ is the crystal density [M/L 3 ], ℓ the edge length of a single site [L], n T is the total number of crystal sites, n s ¼ P 5 i¼1 n i is the number of sites in contact with the fluid, and n v = n 6 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 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 i B the probability of breaking the i bonds. n 0 i ¼ 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 n 0 i ; i ¼ 1::6 À Á , E B , and the frequency factor ν. This optimization problem is thus defined by The first equation of the system is the objective function for matching r M , the computed, andr 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 56 . The algorithm starts with numerous initial solutions (in our case, an initial solution is a set of n′ i , E B and ν 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 ofr M À r M ð Þ 2 for this particle, called local best) and the best value ofr M À r M ð Þ 2 for all particles (called the global best). Moving a particle consists in changing the values of n′ i , E B and ν following: where ψ j k 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], ψ j k;pbest is the local best value and ψ j gbest 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 57 . We chose this one for its simplicity. The iterative algorithm is stopped either whenr M À r M ð Þ 2 <τ, 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 τ ¼ 10 À8r 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 timeindependent 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. 42 provided flow rates and solute concentration, which allows computingr M the measured dissolution rates following: where q is the flow rate [L 3 /T] and C the solute concentration [M/L 3 ]. We P. Ackerer et al.
used the initial mass to compute r M , the simulated dissolution rate (see Eq. (13)).
The flow-through experiments reported in Xu et al. 43 and Smith et al. 40 provided dissolution rates expressed in mol/m 2 /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 writê 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

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.