Arnaud Bouissonnié
email: arnaud.bouissonnie@etu.unistra.fr

Damien Daval

Philippe Ackerer

The Dissolution Anisotropy of Pyroxenes: A Surrogate Model for Steady-State Enstatite Dissolution Resulting From Stochastic Simulations of the Hydrolysis Process

Over the past decades, the field of mineral dissolution kinetics has undergone a spectacular evolution, with an increasingly detailed description of the atomic scale mechanisms of fluidsolid interactions. The development of probabilistic dissolution models has played a prominent role in this evolution, as they allow for bridging the outputs of ab initio development of surrogate models such as those reported here represents one of the possible 48 strategies for upscaling dissolution processes from the atomic scale to the micron scale.

calculations to macroscopic observables such as dissolution rates and nanoscale surface features. It is however admitted that these models cannot be easily adapted to simulate natural systems at large space and time scales due to the restricted dimensions and durations that can be simulated numerically. In the present study, we demonstrate that the steady-state outputs of the face-specific stochastic treatment of enstatite dissolution, which was experimentally validated in a previous paper, can be boiled down to a single analytical expression under the form: stands for the contribution of the dislocations to the overall dissolution rate.

The derivation of simple analytical expressions to get steady-state rate data that are similar to those obtained using stochastic dissolution models, may contribute to parametrize efficiently the bond-breaking probability of various atoms for pyroxene solid solutions, and raises the question of the extension of such surrogate expressions to other silicate structures. Finally, the

Introduction

Many complementary approaches have been developed to predict the fate of waterrock interactions, at various space and time scales. Historically, the conceptual framework underpinning the corresponding models has been closely linked to the disciplinary field in which they were developed.

When dealing with the Earth system, the first attempts to develop weathering rate laws for large space (km²) and time (Myr) scales were undoubtedly phenomenological, and strongly relied on empirical relations. The pioneering work of Berner and coworkers, aimed at reconstructing the evolution of the CO2 concentration in the atmosphere during the Phanerozoic through the development of the BLAG [START_REF] Berner | The Carbonate-Silicate Cycle and Its Effect on Atmospheric Carbon Dioxide over the Past 100 Millions Years[END_REF] and GEOCARB [START_REF] Berner | A Model for Atmospheric Co2 over Phanerozoic Time[END_REF][START_REF] Berner | 3GEOCARB-II -a Revised Model of Atmospheric Co2 over Phanerozoic Time[END_REF][START_REF] Berner | GEOCARB III: A Revised Model of Atmospheric Co2 over Phanerozoic Time[END_REF] codes, was based on model weathering reactions, whose rates were largely parametrized according to empirical rate-runoff, rate-pCO2 and rate-temperature relations inferred from field measurements. Since the early 90s, this approach has been gradually superseded by reactive transport models, where the source-terms of the classical reaction-dispersion-advection equation are described following kinetic rate laws derived from dissolution experiments conducted on powdered single-crystals. Although the corresponding rate equations, which relate the dissolution rate to fundamental parameters such as temperature, pH, surface area or solution saturation state, were originally claimed to result from the transition state theory (TST) [START_REF] Lasaga | Transition State Theory[END_REF][START_REF] Aagaard | Thermodynamic and Kinetic Constraints on Reaction-Rates among Minerals and Aqueous-Solutions .1. Theoretical Considerations[END_REF][START_REF] Oelkers | The Effect of Aluminum, Ph, and Chemical Affinity on the Rates of Aluminosilicate Dissolution Reactions[END_REF] , several studies subsequently questioned the theoretical validity of the extension of TST to silicate dissolution [START_REF] Lasaga | Fundamental Approaches in Describing Mineral Dissolution and Precipitation Rates[END_REF][START_REF] Gin | Theoretical Consideration on the Application of the Aagaard-Helgeson Rate Law to the Dissolution of Silicate Minerals and Glasses[END_REF] , the experimental validation of which being described as nothing more than a fitting exercise with little physical basis by some authors [START_REF] Luttge | Crystal Dissolution Kinetics and Gibbs Free Energy[END_REF] . More generally, a growing number of experimental studies have underlined the limitations of TST-based relations, which fail to accurately relate dissolution rates to solution saturation state [START_REF] Burch | Free Energydependence of Albite Dissolution Kinetics at 80°C and Ph 8.8[END_REF][START_REF] Taylor | The Dependence of Labradorite Dissolution and Sr Isotope Release Rates on Solution Saturation State[END_REF][START_REF] Beig | Albite Dissolution Kinetics as a Function of Distance from Equilibrium: Implications for Natural Feldspar Weathering[END_REF][START_REF] Hellmann | Dissolution Kinetics as a Function of the Gibbs Free Energy of Reaction: An Experimental Study Based on Albite Feldspar[END_REF][START_REF] Dixit | Effect of Solution Saturation State and Temperature on Diopside Dissolution[END_REF][START_REF] Arvidson | Mineral Dissolution Kinetics as a Function of Distance from Equilibrium -New Experimental Results[END_REF][START_REF] Daval | Dissolution Kinetics of Diopside as a Function of Solution Saturation State: Macroscopic Measurements and Implications for Modeling of Geological Storage of Co2[END_REF][START_REF] Nicoleau | The Di-and Tricalcium Silicate Dissolutions[END_REF][START_REF] Pollet-Villard | Does Crystallographic Anisotropy Prevent the Conventional Treatment of Aqueous Mineral Reactivity? A Case Study Based on K-Feldspar Dissolution Kinetics[END_REF] , and are unable to account for the variability (heterogeneity and anisotropy) of crystal dissolution rates [START_REF] Fischer | How Predictable Are Dissolution Rates of Crystalline Material?[END_REF][START_REF] Godinho | Effect of Surface Orientation on Dissolution Rates and Topography of Caf2[END_REF][START_REF] Daval | Linking Nm-Scale Measurements of the Anisotropy of Silicate Surface Reactivity to Macroscopic Dissolution Rate Laws: New Insights Based on Diopside[END_REF][START_REF] Smith | Effects of Crystal Orientation on the Dissolution of Calcite by Chemical and Microscopic Analysis[END_REF][START_REF] Wild | Ph-Dependent Control of Feldspar Dissolution Rate by Altered Surface Layers[END_REF][START_REF] Kurganskaya | Kinetic Monte Carlo Approach to Study Carbonate Dissolution[END_REF][START_REF] Daval | Dynamics of Altered Surface Layer Formation on Dissolving Silicates[END_REF][START_REF] Fischer | Pulsating Dissolution of Crystalline Matter[END_REF][START_REF] Zhang | Aluminosilicate Dissolution Kinetics: A General Stochastic Model[END_REF] . With respect to the relation between dissolution rate and fluid saturation state, the deviation from a TST-behavior has often been attributed to the nucleation and opening of etch pits where dislocations outcrop at the mineral's surface [START_REF] Burch | Free Energydependence of Albite Dissolution Kinetics at 80°C and Ph 8.8[END_REF][START_REF] Pollet-Villard | Does Crystallographic Anisotropy Prevent the Conventional Treatment of Aqueous Mineral Reactivity? A Case Study Based on K-Feldspar Dissolution Kinetics[END_REF][START_REF] Smith | Effects of Crystal Orientation on the Dissolution of Calcite by Chemical and Microscopic Analysis[END_REF][START_REF] Brantley | Dissolution at Dislocation Etch Pits in Quartz[END_REF][START_REF] Lasaga | Surface Chemistry, Etch Pits and Mineral-Water Reactions[END_REF][START_REF] Lasaga | Variation of Crystal Dissolution Rate Based on a Dissolution Stepwave Model[END_REF][START_REF] Bouissonnié | From Mixed Flow Reactor to Column Experiments and Modeling: Upscaling of Calcite Dissolution Rate[END_REF] . Etch pits have been extensively

reported for a wide variety of minerals [START_REF] Daval | Linking Nm-Scale Measurements of the Anisotropy of Silicate Surface Reactivity to Macroscopic Dissolution Rate Laws: New Insights Based on Diopside[END_REF][START_REF] Smith | Effects of Crystal Orientation on the Dissolution of Calcite by Chemical and Microscopic Analysis[END_REF][START_REF] Bouissonnié | From Mixed Flow Reactor to Column Experiments and Modeling: Upscaling of Calcite Dissolution Rate[END_REF][START_REF] Amelinckx | The Direct Observation of Dislocations[END_REF][START_REF] Teng | Controls by Saturation State on Etch Pit Formation During Calcite Dissolution[END_REF][START_REF] Kurganskaya | Does the Stepwave Model Predict Mica Dissolution Kinetics?[END_REF][START_REF] Kurganskaya | A Comprehensive Stochastic Model of Phyllosilicate Dissolution: Structure and Kinematics of Etch Pits Formed on Muscovite Basal Face[END_REF][START_REF] Pollet-Villard | Influence of Etch Pit Development on the Surface Area and Dissolution Kinetics of the Orthoclase (001) Surface[END_REF][START_REF] Perez | Comparing the Reactivity of Glasses with Their Crystalline Equivalents: The Case Study of Plagioclase Feldspar[END_REF] , when the fluid under saturation is large enough and does not exceed a critical value of the Gibbs free energy. The efficiency of these sites to enhance mineral dissolution rate has been showed both experimentally and numerically, and is accounted for by the stepwave model [START_REF] Lasaga | Variation of Crystal Dissolution Rate Based on a Dissolution Stepwave Model[END_REF] . At the atomic-scale, the opening of etch pits results in an increase of the density of kink and step sites. This process exposes atoms with a lower coordination to the fluid which are, as a result, more rapidly dissolved.

More generally, the reactivity of minerals is highly heterogeneous, with hotspots of reactivity such as kinks, steps (and ultimately crystal edges and corners) and can be accurately captured by the "rate spectra" concept [START_REF] Fischer | How Predictable Are Dissolution Rates of Crystalline Material?[END_REF][START_REF] Luttge | Kinetic Concepts for Quantitative Prediction of Fluid-Solid Interactions[END_REF] . Finally, several studies have pointed out that mineral dissolution is anisotropic 19, 22-23, 38, 40-41 , resulting from the anisotropic distribution of atomic positions in the crystal lattice. These observations have emphasized the need for deeper experimental and theoretical studies of fluid-mineral interactions from the atomic-to the nmscale to unravel the dissolution mechanisms.

Precisely at the other end of the spatiotemporal spectrum, the advent of microscopic and spectroscopic techniques of characterization of the fluid-mineral interface offered new avenues to understand mineral reactivity [START_REF] Daval | Dynamics of Altered Surface Layer Formation on Dissolving Silicates[END_REF][START_REF] Schott | X-Ray Photoelectron Studies of the Mechanism of Iron Silicate Dissolution During Weathering[END_REF][START_REF] Eggleston | Sample Preparation and Aging Effects on the Dissolution Rate and Surface-Composition of Diopside[END_REF][START_REF] Gratz | Atomic Force Microscopy of Atomic-Scale Ledges and Etch Pits Formed During Dissolution of Quartz[END_REF][START_REF] Jordan | Dissolution Rates and Activation Energy for Dissolution of Brucite (001) : A New Method Based on the Microtopography of Crystal Surfaces[END_REF][START_REF] Lüttge | An Interferometric Study of the Dissolution Kinetics of Anorthite; the Role of Reactive Surface Area[END_REF][START_REF] Wirth | Focused Ion Beam (Fib): A Novel Technology for Advanced Application of Microand Nanoanalysis in Geosciences and Applied Mineralogy[END_REF][START_REF] Cailleteau | Insight into Silicate-Glass Corrosion Mechanisms[END_REF][START_REF] Noiriel | Resolving Time-Dependent Evolution of Pore-Scale Structure, Permeability and Reactivity Using X-Ray Microtomography[END_REF][START_REF] Hellmann | Nanometre-Scale Evidence for Interfacial Dissolution-Reprecipitation Control of Silicate Glass Corrosion[END_REF][START_REF] Leonard | Exploring Dynamic Surface Processes During Silicate Mineral (Wollastonite) Dissolution with Liquid Cell Tem[END_REF][START_REF] Ruiz-Agudo | Control of Silicate Weathering by Interface-Coupled Dissolution-Precipitation Processes at the Mineral-Solution Interface[END_REF][START_REF] Geisler | Real-Time in Situ Observations of Reaction and Transport Phenomena During Silicate Glass Corrosion by Fluid-Cell Raman Spectroscopy[END_REF] . Studies dedicated to the molecular-scale description of the silicate dissolution process following ab initio quantum mechanical calculations and molecular dynamics (MD) models emerged some 20 years ago [START_REF] Xiao | Ab Initio Quantum Mechanical Studies of the Kinetics and Mechanisms of Silicate Dissolution: H+(H3o+) Catalysis[END_REF][START_REF] Pelmenschikov | Lattice Resistance to Hydrolysis of Si-O-Si Bonds of Silicate Minerals: Ab Initio Calculations of a Single Water Attack onto the (001) and (111) Β-Cristobalite Surfaces[END_REF][START_REF] Criscenti | Silicate Glass and Mineral Dissolution: Calculated Reaction Paths and Activation Energies for Hydrolysis of a Q(3) Si by H3o+ Using Ab Initio Methods[END_REF][START_REF] Morrow | Ab Initio Investigation of Dissolution Mechanisms in Aluminosilicate Minerals[END_REF][START_REF] Morrow | Description of Mg2+ Release from Forsterite Using Ab Initio Methods[END_REF][START_REF] Morrow | Quantum Mechanical Modeling of Hydrolysis and H2o-Exchange in Mg-, Ca-, and Ni-Silicate Clusters: Implications for Dissolution Mechanisms of Olivine Minerals[END_REF] . These studies provided essential theoretical insights into the dissolution process. Admittedly, the corresponding studies have been of limited direct interest (and were not intended) to model chemical weathering in natural settings, due to the restricted durations and dimensions of the system that can be simulated numerically.

A first intermediate step in the upscaling process of dissolution reaction rates has been successfully reached with the emergence of stochastic models of silicate and carbonate dissolution [START_REF] Kurganskaya | Kinetic Monte Carlo Approach to Study Carbonate Dissolution[END_REF][START_REF] Zhang | Aluminosilicate Dissolution Kinetics: A General Stochastic Model[END_REF][START_REF] Kurganskaya | A Comprehensive Stochastic Model of Phyllosilicate Dissolution: Structure and Kinematics of Etch Pits Formed on Muscovite Basal Face[END_REF][START_REF] Cailleteau | Insight into Silicate-Glass Corrosion Mechanisms[END_REF][START_REF] Geisler | Real-Time in Situ Observations of Reaction and Transport Phenomena During Silicate Glass Corrosion by Fluid-Cell Raman Spectroscopy[END_REF][START_REF] Morrow | Ab Initio Investigation of Dissolution Mechanisms in Aluminosilicate Minerals[END_REF][START_REF] Lasaga | Mineralogical Approaches to Fundamental Crystal Dissolution Kinetics -Dissolution of an a(3)B Structure[END_REF][START_REF] Devreux | Leaching of Borosilicate Glasses. Ii. Model and Monte-Carlo Simulations[END_REF][START_REF] Zhang | Si Order in Albite and Its Effect on Albite Dissolution Processes: A Monte Carlo Study[END_REF][START_REF] Meakin | Simple Kinetic Monte Carlo Models for Dissolution Pitting Induced by Crystal Defects[END_REF][START_REF] Kurganskaya | Kinetic Monte Carlo Simulations of Silicate Dissolution: Model Complexity and Parametrization[END_REF][START_REF] Briese | The Effect of Crystal Size Variation on the Rate of Dissolution -a Kinetic Monte Carlo Study[END_REF][START_REF] De Assis | Dissolution of Minerals with Rough Surfaces[END_REF][START_REF] Lasaga | Mineralogical Approaches to Fundamental Crystal Dissolution Kinetics[END_REF][START_REF] Zhang | Morphological Evolution of Dissolving Feldspar Particles with Anisotropic Surface Kinetics and Implications for Dissolution Rate Normalization and Grain Size Dependence: A Kinetic Modeling Study[END_REF][START_REF] Suarez | Short-and Long-Term Weathering Rates of a Feldspar Fraction Isolated from an Arid Zone Soil[END_REF] , whose parametrization may strongly rely on the results of ab initio calculations mentioned above, with considerable simplifications of the reaction details.

Such an approach dramatically increased the size of the studied systems from clusters of a few atoms to nanocrystals of several tens of millions of atoms, simulating time scales of up to several weeks [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF][START_REF] Lüttge | A Stochastic Treatment of Crystal Dissolution Kinetics[END_REF] . However, implementing such models into reactive transport codes is not an easy task since they usually do not provide the rate laws that are required by these codes to predict the rates of water-rock interactions at large space and time scales [START_REF] Knauss | Evaluation of the Impact of Co2, Co-Contaminant Gas, Aqueous Fluid and Reservoir Rock Interactions on the Geologic Sequestration of Co2[END_REF][START_REF] Godderis | Modelling Weathering Processes at the Catchment Scale: The Witch Numerical Model[END_REF] . Promising strategies combining Monte Carlo simulations and Voronoi methods have been recently developed to circumvent this issue [START_REF] Luttge | Kinetic Concepts for Quantitative Prediction of Fluid-Solid Interactions[END_REF][START_REF] Rohlfs | Crystal Dissolution Kinetics Studied by a Combination of Monte Carlo and Voronoi Methods[END_REF] , offering the perspective to reach the next step in the upscaling process, i.e., the continuum scale, where reactive transport can be applied. Such recent breakthroughs emphasize the crucial importance of maintaining efforts to investigate into more details the various theoretical and empirical emerging relations of stochastic dissolution models, which represents one of the goals of the present paper.

In a previous study [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] , we developed a probabilistic model of enstatite dissolution which successfully reproduced the measured anisotropy of enstatite dissolution rates and associated surface features. In the present study, hundreds of simulations were conducted with this model by varying the bond-breaking probabilities of the crystal to explore the analytical relation that may link together the detachment of individual atoms from the enstatite surface to the overall face-specific dissolution rates. This work was in part motivated by the fact that mineral dissolution kinetics is often treated by mixing macroscopic parameters with microscopic ones, resulting in rate laws for which the theoretical basis remains questionable (see discussion in e.g. [START_REF] Gin | Theoretical Consideration on the Application of the Aagaard-Helgeson Rate Law to the Dissolution of Silicate Minerals and Glasses[END_REF] ). Here we illustrate how overall rate constants may be derived from atomic-scale parameters, considering especially the bond-breaking probability of hydrolysis processes. In particular, we show that, at steady-state conditions, an analytical relation can be derived to express the overall dissolution rate constant, whose parameters can be adjusted following the outputs of numerical simulations performed at the atomic scale. We then show that the resulting analytical relation can be used to predict the steady-state enstatite facespecific dissolution rates over a wide range of bond-breaking probabilities. This relation is then further extended from defect-free surfaces to defective surfaces impacted by screw dislocations. Finally, we discuss the implications of such surrogate expressions, both from the perspective of the mechanisms of pyroxene dissolution and for the upscaling of silicate dissolution rates.

Methods

Model description and algorithms

A complete description of the model is given in Bouissonnié et al. [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] . In short, the positions of Mg, Si and O atoms provided by Hugh-Jones and Angel [START_REF]A Compressional Study of Mgsio3 Orthoenstatite up to 8.5 Gpa[END_REF] are used to create the enstatite cell thanks to the symmetry elements of the Pbca space group. Mg and Si atoms are then linked to the 6 and 4 nearest O atoms respectively. Since probabilistic dissolution models generally consider the bond-breaking probabilities of M-O-M' bonds [START_REF] Lasaga | Surface Chemistry, Etch Pits and Mineral-Water Reactions[END_REF][START_REF] Kurganskaya | A Comprehensive Stochastic Model of Phyllosilicate Dissolution: Structure and Kinematics of Etch Pits Formed on Muscovite Basal Face[END_REF][START_REF] Suarez | Short-and Long-Term Weathering Rates of a Feldspar Fraction Isolated from an Arid Zone Soil[END_REF] , each Mg and Si atom is connected to its first coordination sphere [START_REF] Kurganskaya | A Comprehensive Stochastic Model of Phyllosilicate Dissolution: Structure and Kinematics of Etch Pits Formed on Muscovite Basal Face[END_REF] .

Usually, the probability attributed with one event (i.e. the bond-breaking probability) is written as follows [START_REF] Lasaga | Surface Chemistry, Etch Pits and Mineral-Water Reactions[END_REF] :

a B E k T P e   (1)
where 𝑃 stands for the bond-breaking probability, 𝐸 for the activation energy of the bond hydrolysis (J), 𝑘 for the Boltzmann constant (J/K) and 𝑇 for the temperature (K). In the present study, the most likely bond-breaking probability (i.e., the bond hydrolysis associated to with the lowest activation energy, Mg-O-Mg) was arbitrarily set to 0.99 and the two other probabilities were scaled according to the differences between their corresponding activation energies following:

( ) A B a a B E E k T A B P e P    ( 2 )
where A and B represent two different bonds.

In agreement with previous studies [START_REF] Cailleteau | Insight into Silicate-Glass Corrosion Mechanisms[END_REF][START_REF] Devreux | Leaching of Borosilicate Glasses. Ii. Model and Monte-Carlo Simulations[END_REF][START_REF] Zhang | Si Order in Albite and Its Effect on Albite Dissolution Processes: A Monte Carlo Study[END_REF][START_REF] Meakin | Simple Kinetic Monte Carlo Models for Dissolution Pitting Induced by Crystal Defects[END_REF][START_REF] Kurganskaya | Kinetic Monte Carlo Simulations of Silicate Dissolution: Model Complexity and Parametrization[END_REF][START_REF] Briese | The Effect of Crystal Size Variation on the Rate of Dissolution -a Kinetic Monte Carlo Study[END_REF][START_REF] De Assis | Dissolution of Minerals with Rough Surfaces[END_REF] , an "all-or-none" approach has been used.

This approach considers the hydrolysis of atoms instead of individual bonds. In other words, for a given simulated step, an atom is released only when all the remaining bonds that connect the atom to the surface of the crystal are broken simultaneously. Otherwise, the corresponding atom remains connected to the surface. The probability of an atom M to be released in the fluid depends on the number and type of atomic bonds present in its first coordination sphere.

In the case of the enstatite, this probability can be written as follows:

M O Mg M O Si B B M O Mg M O Si E E n m k n m T kT M P P e P e             (3)
where 𝑛 and 𝑚 stand for the number of bonds that the considered atom shares with Mg and Si atoms respectively. The dissolution algorithm was described in details in Bouissonnié et al. [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] and is schematically represented in Fig. 1. Such algorithms have previously proven successful to simulate the dissolution of silicate materials (e.g. [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF][START_REF] Cailleteau | Insight into Silicate-Glass Corrosion Mechanisms[END_REF][START_REF] Devreux | Leaching of Borosilicate Glasses. Ii. Model and Monte-Carlo Simulations[END_REF] ). In short, a random number Z uniformly distributed between 0 and 1 is drawn for each atom at the surface defined by an incomplete coordination sphere (i.e. an atom at the crystal surface). If Z < PM, all atomic bounds are broken and the atom is considered as dissolved (removed from the surface); if not, the atom stays at the surface. The first coordination spheres of the remaining atoms are then updated if necessary (i.e., when a removed atom was belonging to the coordination spheres of one of the neighboring remaining atoms). Because coordination spheres are modified, the probability of dissolution PM of the corresponding atoms increases. In such kinds of algorithms, each iteration step corresponds to a given simulated duration, and the time-step is therefore constant for all simulations. The proportionality factor to convert the number of iterations into time is not known a priori, and can be determined by comparing the outputs of a simulation with experimental data (see [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] for details). This aspect contrasts with kinetic Monte Carlo algorithms based on a "divide and conquer" approach [START_REF] Kurganskaya | A Comprehensive Stochastic Model of Phyllosilicate Dissolution: Structure and Kinematics of Etch Pits Formed on Muscovite Basal Face[END_REF][START_REF] Meakin | Simple Kinetic Monte Carlo Models for Dissolution Pitting Induced by Crystal Defects[END_REF] , which use an adaptive time step that varies at each iteration as a function of the random number that is drawn, corresponding to a specific surface site that is dissolved during the iteration. The main strengths of the "divide and conquer" approach are that they allow for a direct link between time and the number of iterations, while precluding the existence of "dead" iterations, reducing accordingly the computing time required to simulate a given process. Conversely, "dead" iterations may sometimes occur over the course of the algorithm that we used (essentially, at the beginning of the simulations), but as opposed to the "divide and conquer" approach, it allows for simultaneous events to occur over a given iteration (and time step).

Although both methods should ultimately provide identical results (see theoretical derivation in [START_REF] Bortz | A New Algorithm for Monte Carlo Simulation of Ising Spin Systems[END_REF] ), to the best of our knowledge, a comparison between both types of algorithms has never been proposed for mineral dissolution, and is out of the scope of the present study.

The impact of dislocation density has also been studied by running additional simulations for the (100) face. Dislocations were introduced as lines running parallel to the [100] axis, in agreement with previous observations [START_REF] Skrotzki | Defect Structure and Deformation Mechanisms in Naturally Deformed Augite and Enstatite[END_REF] . To simulate the faster dissolution at dislocation sites (placed at the center of the reactive surface), the bond-breaking probabilities associated to with the atoms belonging to the dislocation line were set to 1 following Kurganskaya and LuttgeKurganskaya and Luttge [START_REF] Kurganskaya | A Comprehensive Stochastic Model of Phyllosilicate Dissolution: Structure and Kinematics of Etch Pits Formed on Muscovite Basal Face[END_REF] .

One of the goals of this study was to understand how the dissolution rate evolves as a function of the probabilities used as input parameters. To perform this sensitivity analysis, each of the three fundamental probabilities (PMg-O-Mg, PSi-O-Mg, or PSi-O-Si) was sampled in a given range while the other two remained constant.

Finally, the boundary conditions (BC) were set as follows: atoms that are part of the sides of the simulated volume cannot be dissolved. Since it can induce a dependence of the dissolution rate on the dimensions of the simulated surface, we verified that the simulated surface area was large enough to avoid any impact of the BC on the dissolution rate. In order to determine the dissolution rate, the number of released Mg and Si atoms is stored at each step (iteration). In this study, the dissolution rate calculation is based on the release of Si atoms (equivalent to the one based on Mg atoms departure and mean surface height at steady-state conditions [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] ):

Outputs of the model

𝑟 ( ) = 𝑉 𝑖𝑡𝑒𝑟 × 𝑁 , × 𝑆 𝑁 , (4)
where 𝑟 ( ) stands for the dissolution rate (Å/it) of a given face (here, either (100), ( 010) or (001)), 𝑁 , for the number of Si atoms released at each iteration i,

, Si cell N
for the number of Si atoms in enstatite cell (16), iter for the total number of iterations, and 𝑆 for the geometric surface area of the crystal face (Ų). The second output of the model is related to the chemical environment of the atoms when they are released from the surface. In particular, knowledge of the first coordination sphere of each atom that has been dissolved allows for the determination of the number of Mg and Si atoms that were connected to the released M atoms. As opposed to other similar numerical studies (e.g. [START_REF] Fischer | Beyond the Conventional Understanding of Water-Rock Reactivity[END_REF][START_REF] Fischer | Inherited Control of Crystal Surface Reactivity[END_REF] ), our study was not aimed to capture the intrinsic heterogeneous distribution of reaction rates at the simulated crystal surface. Instead, our primary goal was to provide an analytical expression for the overall facespecific steady state dissolution rate as a function of bond-breaking probabilities, which has to be achieved at the scale of a given face, echoing the face-specific dissolution rates frequently reported from dissolution experiments (see [START_REF] Luttge | Kinetic Concepts for Quantitative Prediction of Fluid-Solid Interactions[END_REF][START_REF] Noiriel | Pore-Scale Geochemical Reactivity Associated with Co2 Storage: New Frontiers at the Fluid-Solid Interface[END_REF] for an overview). 1. Legends correspond to the code indicated for each group of simulations. These codes are used in Table 1 to indicate which simulation corresponds to which groups of data. 𝑟 = 𝑓(𝑃 ) cases respectively). These different groups and the full dataset are listed in the Table 1.

Results

Dissolution rate evolution as a function of bond-breaking probabilities

13 Mg O Mg P   Mg O Si P   Si O Si P   ∆E (Mg-O-Mg-Mg-O-Si) (kJ/mol) ∆E (Mg-O-Si-Si-O-Si) (kJ/mol) ∆E (Mg-O-Mg-Si-O-Si) (kJ/
P   Mg O Si P   Si O Si P   ∆E (Mg-O-Mg-Mg-O-Si) (kJ/mol) ∆E (Mg-O-Si-Si-O-Si) (kJ/mol) ∆E (Mg-O-Mg-Si-O-Si) (kJ/
A linear relationship between the logarithm of the dissolution rate and the logarithm of the bond-breaking probabilities can be observed for the different simulations (see Fig. 2). This linear model matches results from numerical simulations quite satisfactorily. Moreover, the slope remains constant when changing the activation energy for the hydrolysis of one bond, whatever the difference in the activation energy of the two other bonds (i.e. the slopes are the same for the results shown in each Fig. 2A, 2B and 2C). Slopes, intercepts and R² of these linear regressions are listed in Table . 1.

Mean first coordination sphere of released atoms

The attainment of a mean surface configuration when steady-state dissolution rates are reached represents a common trait of all simulations. The corresponding mean first coordination spheres of the M atoms at the iteration step of their detachment from the surface are depicted in Fig.

Dislocation density

The results of simulations conducted with dislocations are reported in Fig. 4. The dislocation density has been investigated by varying the surface area (i.e., all simulations were conducted using a single dislocation line, while varying the simulated surface area).

Discussion

Bulk dissolution rate evolution as a function of bond-breaking probabilities

From the numerical experiments described in Section 3.1, the steady state dissolution rates of enstatite (defined as congruent and constant dissolution rates with time) are proportional to each individual bond-breaking probability raised at a given power inferred from the linear regressions depicted in a log-log diagram (Fig. 2). Moreover, since the slope of the linear regression between Log(Pi) (i stands for Mg-O-Mg or Mg-O-Si or Si-O-Si) and Log(r) is constant whatever the difference between the two other involved probabilities; a more general function of the dissolution rate can be written as:

  Mg O Mg Mg O Si Si O S hkl bulk i k P P P r           ( 5
)
where k is a constant and , ,    are the corresponding slopes listed in Table 2 for the various regressions. The parameter k is then obtained using equation 5 and is constant for each face.

The different parameter values are summarized in If the parameter α can be considered as constant for (100) and ( 010) faces (Table 1), the slope that corresponds to a variation of PMg-O-Mg decreases from 7.46 to 5.05 when the difference between EMg-O-Si and ESi-O-Si decreases, i.e., when EMg-O-Si gets closer to ESi-O-Si. For these cases, the assumptions used in the concept of the surrogate model may not be fulfilled (e.g., k is constant and does not depend on the different probabilities). This assumption will be discussed in the next section. Despite this approximation and the apparent trend for the (001)

face, the proposed model described by Eq. 5 is capable of satisfactorily estimating the simulated dissolution rates (Fig. 5) for the three different faces. Note that Eq. ( 5) has been derived from simulations where only the first coordination sphere of surface atoms is considered for expressing the dissolution probability of individual atoms. Several studies [START_REF] Kurganskaya | A Comprehensive Stochastic Model of Phyllosilicate Dissolution: Structure and Kinematics of Etch Pits Formed on Muscovite Basal Face[END_REF][START_REF] Kurganskaya | Kinetic Monte Carlo Simulations of Silicate Dissolution: Model Complexity and Parametrization[END_REF] have shown that second coordination sphere may also play a role in the dissolution, and may impact the shape of the developed etch pits. However, we showed that considering only the first coordination sphere is enough to reproduce satisfactorily the face-specific shape of etch pits observed on enstatite [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] , explaining why we stick to this model in the present study.

Considering the impact of the second coordination sphere would make Eq. 5 more sophisticated, since it would probably incorporate specific terms related to the second coordination spheres.

Theoretical interpretation of the surrogate model

The development of the surrogate expression given in Eq. 5 is derived from a statistical analysis of the outputs of the stochastic simulations of enstatite dissolution, with no preconception of the mathematical form that should be used to relate the steady-state dissolution rates to the individual bond-breaking probabilities describing enstatite hydrolysis.

In the present section, we provide a possible theoretical explanation of Eq. 5.

The dissolution rate can be defined as the derivative of the number of atoms that are dissolved with time:

, M d dN r dt  (6)
where NM,d stands for the amount of dissolved atoms belonging to the M species during the time interval dt. Arguably, NM,d should depend on two parameters: the amount of atoms located at the mineral surface (NM,S) and the intrinsic detachment rate of these atoms.

Importantly, classical theories of dissolution kinetics suggest that at steady-state conditions, NM,S is constant and proportional to the considered surface area [START_REF] Aagaard | Thermodynamic and Kinetic Constraints on Reaction-Rates among Minerals and Aqueous-Solutions .1. Theoretical Considerations[END_REF] . This aspect was numerically verified in our previous study [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] where we showed that the amount of surface atoms level to a plateau as soon as the dissolution becomes congruent and the dissolution rate is constant. At steady state, the number of dissolved atoms of Mg or Si is given by:

ˆMg Mg Si Si N N P N P   ( 7 )
where N is the number of dissolved atoms of Mg or Si, ˆA N is the number of atoms A (Mg or Si) at the mineral surface, A P is the probability that atom A is dissolved. Following the previous equation, N can be rewritten as:

ˆMg Mg Si Si N N P N P  (8)
, Comparison between the dissolution rates estimated with the surrogate model (r theory ; Eq. 5) and dissolution rates provided by the numerical simulations (r mod ) for faces (100) (A), (010) (B) and (001) (C). The comparison is made for all the different simulations listed in Table 1. For instance, (A) presents all the simulations shown in Fig. 2, and the rate variability results from the corresponding bond-breaking probabilities that are reported in Table 1.

i i A i A A A B P P P      (9)
where i  and i  are the number of bonds connecting the considered A atom to another A or to a B atom (Si or Mg accordingly), respectively. At the solid surface, we consider that there are NC different coordination spheres and i n atoms with a coordination sphere i. Therefore, the probability to dissolve the atom A over the mineral face is:

1 1 ˆC i i N A i A A A B i A P nP P N        and 1 ˆC N A i i N n    (10)
Or

1 1 ˆC i A i A A A A A N A i A A A B A A A B A A A A B i A P nP P P P k P P N                           (11) with 1 1 ˆC i
A i A N A iA A A B i A k nP P N            (12)
and the number of dissolved atoms is:

ˆˆˆM g Mg Si Si Mg Si Mg MgOMg MgOSi Si SiOSi MgOSi Mg Mg Si Si MgOMg MgOSi SiOSi N N N k P P k P P N k N k P P P          (13)
with , ,

Mg M g S i S i            .
The statistical analysis of the numerical experiments provides some values of the parameters , , , k    (see Table 2). Surprisingly, the parameter k does not depend on the probabilities k can be rewritten as:

1, 1 ˆˆC i j i j N j A i A A A B i i j A A n k n P P N N              (14)
Therefore, A k can be a constant value if 1,

C i j i j N j iAA A B i i j n nP P             which is possible when i j    and/or i j  
 , for having positive exponents for the probabilities A A P  and/or A B P  .

Given that the exponents are positive, considering k as a constant is a reasonable assumption if the difference between A A P  and A B P  is significant. The validity of this condition is difficult to appreciate if A A P  and A B P  are close, as it is the case for the simulations run for (001) face (see section 4.1 and Table 1). Finally, this demonstrates that the dissolution rate is therefore compatible with Eq. 5:

Mg O Mg Mg O Si Si O Si kP P P N r W W dt dt            with 𝑊 = × × (15)
with cmx N standing for the number of "complexes" (see 4.4) in the enstatite cell.

Accounting for the impact of dislocation density on the dissolution rates

Dislocations have a measurable impact on the dissolution rate, whatever the set of probabilities that was tested (Fig 6C). This result is consistent with numerical and experimental studies that showed that the presence of dislocations outcropping at mineral surfaces globally increases the dissolution rate when the nucleation of etch pits is thermodynamically favorable (i.e., at far-from-equilibrium conditions [START_REF] Burch | Free Energydependence of Albite Dissolution Kinetics at 80°C and Ph 8.8[END_REF][START_REF] Hellmann | Dissolution Kinetics as a Function of the Gibbs Free Energy of Reaction: An Experimental Study Based on Albite Feldspar[END_REF][START_REF] Arvidson | Mineral Dissolution Kinetics as a Function of Distance from Equilibrium -New Experimental Results[END_REF][START_REF] Smith | Effects of Crystal Orientation on the Dissolution of Calcite by Chemical and Microscopic Analysis[END_REF][START_REF] Lasaga | Variation of Crystal Dissolution Rate Based on a Dissolution Stepwave Model[END_REF][START_REF] Pollet-Villard | Influence of Etch Pit Development on the Surface Area and Dissolution Kinetics of the Orthoclase (001) Surface[END_REF] ).

However, when dealing with natural samples, the dislocation density remains a parameter impossible or difficult to control a priori, which complicates the prediction of its impact on dissolution rate. Furthermore, several studies showed that the mineral dissolution rate is not a strictly increasing function of the dislocation density [START_REF] Lasaga | Variation of Crystal Dissolution Rate Based on a Dissolution Stepwave Model[END_REF][START_REF] Pollet-Villard | Influence of Etch Pit Development on the Surface Area and Dissolution Kinetics of the Orthoclase (001) Surface[END_REF][START_REF] Schott | Dissolution Kinetics of Strained Calcite[END_REF] since above a given threshold, the dissolution rate no longer varies with the dislocation density.

To qualitatively comply with the above-mentioned studies, the results for the (100) face (Fig. 4) were fitted using the following empirical relation between dislocation density and the additional dissolution rate resulting from etch pit opening at dislocation outcrop (see Section 3.3.), which verifies that the dissolution rate levels to a plateau value when the dislocation density tends towards infinity:

(100)

1 d dislocation app r k e            (16)
where app k and ω (≈ 5.51 x 10 10 ) are empirical parameters that were obtained by calibration on the outputs of the simulations, and 𝜌 is the dislocation density. Whereas ω does not seem to vary with the probabilities used as input parameters, it clearly appears that app k (Fig. 4) depends on the individual bond-breaking probabilities. From a physical standpoint, this observation may result from the fact the first coordination spheres of atoms in the vicinity of a dislocation line differ from those of atoms considered for defect-free surface.

Considering this apparent relationship between kapp and the probabilities, a first tentative to link the global rate to the probabilities and the dislocation density was performed using the same exponent for each probability as the "bulk" case. Since the global rate is equivalent to 𝑟 ( ) = 𝑟 ( ) + 𝑟 ( ) , it is possible to assume the following relation:

  1 d Mg O Mg Mg O Si Si O Si hkl dis r k k e P P P                              (17)
A comparison between this relation and the experimental results is given in Fig. 6A and the value of kdis is given in Table 3. Although the relationship between the surrogate model and the results of the simulations was proven efficient to simulate the dissolution of defect-free surfaces (R² > 0.9 in all cases), it is not the case when using Eq. 17, for which R² = 0.76.

In order to improve the agreement between the surrogate model of the global rate and the results of the simulations, numerical experiments were run for face (100) with dislocations similarly to Section 4.1. The dissolution rates were fitted using the following function:

  1 dis dis dis d Mg O Mg Mg hkl disloca O Si Si O S tion dis i k P P P e r k                     (18)
The values of , , dis dis dis    are lower than those obtained in the "bulk" (defect-free) case (Table 3), which makes sense since the coordination spheres of atoms are arguably lower in the vicinity of the dislocation line, resulting in greater measured dissolution rates. , resulting in the following relation:

k dis dis  dis  dis  γ
  1 dis dis d dis hkl di Mg O Mg Mg O Si Si O Si Mg O Mg Mg O Si Si O Si s P P r P k P P P e k                               (19)
Therefore, the agreement between the results of the simulations and the surrogate model is improved and its correlation coefficient increases from 0.72 to 0.86 (Fig. 6B).

If this type of relations can be extended to other minerals and other types of defects (e.g. vacancies), they may represent an interesting method to link together dissolution rates derived from considerations at the atomic scale and those derived from observations at the mineral scale. This type of relations can also represent a significant progress for reactive transport models that simulate steady-state dissolution reactions. Indeed, even though the present work represents a preliminary step to link relations at the atomic scale to macroscopic dissolution rates, it shows promise as a mean to express face-specific steady-state dissolution rates as a function of the bond-breaking probabilities (and therefore, of the activation energies of hydrolysis) of the different types of bonds that exist for a given mineral. While some studies have questioned the possibility to use the outputs of Monte Carlo simulations in reactive-transport codes because of the limited space-and time-scales investigated following such stochastic treatments of the dissolution process, here we show that steady-state dissolution rates may be satisfactorily calculated from considerations at the atomic scale. As a consequence, rate equations similar to Eq. 17 or Eq. 19 may be simply implemented in reactive transport codes as the rate-constant of the source term of the classical reactionadvection-diffusion equation. This may be of interest for modeling processes for which the fluid composition can be considered as unchanged over the considered simulated time, such as chemical weathering in rivers or in some aquifers with a constant fluid velocity, or fluid circulation in deep fractured geothermal reservoirs, where the fluid composition is approximately constant and poorly affected by geothermal power plant functioning (e.g. [START_REF] Sanjuan | Main Characteristics of the Deep Geothermal Brine (5km) at Soultz-Sous-Forêts (France) Determined Using Geochemical and Tracer Test Data[END_REF] ).

More work remains however needed to assess whether this result is limited to enstatite or the Agreement using (A) Eq. 17 and (B) Eq. 19. (C) Comparison between dislocation-driven dissolution and bulk (dislocation-free) dissolution. The presence of dislocations increases dissolution rates for all the simulations. The data points represent the simulations conducted using sets of probabilities and dislocation densities given in Table 4.

pyroxene group (there is a priori no reason to think so), as well as expanding the relation to the transient regime, where the dissolution rate is neither stoichiometric nor constant. As long as this latter task remains not fulfilled, alternate approaches such as the use of Voronoi distance maps [START_REF] Luttge | Kinetic Concepts for Quantitative Prediction of Fluid-Solid Interactions[END_REF][START_REF] Rohlfs | Crystal Dissolution Kinetics Studied by a Combination of Monte Carlo and Voronoi Methods[END_REF] probably represents the most promising compromise to upscale kinetic Monte Carlo simulations. This strategy nonetheless remains far more complicated than just using a surrogate expression as proposed above.

Relation between the coordination of atoms leaving the surface and the surrogate expression

Interestingly, when atoms are released from the mineral surface, they seem to have a specific average coordination. This assertion is supported by the analysis of the environment of Mg and Si when they leave the surface, since the mean values of Mg-O-Mg, Mg-O-Si+Mg-O-Si and Si-O-Si bonds remain constant and unaffected by the set of probabilities used to run the simulations for all faces (see Section 3.2). However, the numerical values of the parameters used in the surrogate expression of the dissolution rate are face-specific. This apparent paradox can be easily explained by recalling that the histograms reported in Fig. 3 actually take into account all atoms leaving the surface, including those whose departure is not rate-limiting of the dissolution process.

By analogy with Eq. 3, the surrogate expression given by Eq. 17 or Eq. 19 may reflect the average coordination of atoms that control the dissolution process, as kink sites would do.

The fact that the numerical values of α, β and γ are face-specific indicates that this average rate-controlling configuration is not unique for a given mineral. Most likely, this result reflects the fact that kink sites, which are much more difficult to define using a real crystal lattice than simplified isotropic cubic Kossel crystals, differ from one face to another when dealing with anisotropic structures. This explanation would in turn be consistent with the observed anisotropic reactivity of the pyroxene structure.

Conclusion

In this study, we presented the results of hundreds of probabilistic simulations of enstatite dissolution to link the overall dissolution rate to the bond-breaking probabilities used as input parameters. By varying independently each probability, we showed that it is possible to build a surrogate model that links the different probabilities to the dissolution rate following a power law. This result contributes to the general effort of upscaling of mineral dissolution kinetics, since this surrogate expression is based on a mechanistic approach developed from considerations at the atomic-scale, from which the resulting dissolution rate constants can be used as a source term in reactive transport simulations. However, the relation remains valid at steady-state conditions only, and the transient regime must be treated following other upscaling approaches, such as the use of Voronoï distance maps [START_REF] Luttge | Kinetic Concepts for Quantitative Prediction of Fluid-Solid Interactions[END_REF][START_REF] Rohlfs | Crystal Dissolution Kinetics Studied by a Combination of Monte Carlo and Voronoi Methods[END_REF] .

The various simulations conducted with dislocations have shown that it is possible to link the dissolution rate to the dislocation density by introducing an exponential factor to the global 𝑟 = 𝑓 ∏ 𝑃 relation. This relation further extends the interest of probabilistic simulations of mineral dissolution to account for the impact of some important parameters that are hardly controlled in experimental studies. In addition, in case of a bond hydrolysis with significant higher activation energy than the others, the dissolution rate is strongly correlated to this activation energy, and the surrogate model may be used to estimate the value of this parameter.

However, the surrogate model also presents some limitations, particularly when M-O-Si and Si-O-Si hydrolysis probabilities are getting close to each other (where M is a divalent cation, i.e., Mg in the case of enstatite). In this specific case, the results may be out of the limits of the theoretical framework that supports the development of the surrogate expression.

Of note, such cases are however those for which the input parameters are physically unrealistic, as the activation energy of Si-O-Si hydrolysis is admitted to be much lower than any other M-O-Si bond.

)

  is the steady-state dissolution rate of a defect-free (hkl) face [Å/iter], PM-O-M' stands for the bond-breaking probability of the M-O-M' bond, and k, α, β and γ are fitting parameters adjusted following the outputs of the stochastic simulations. When dislocations outcrop at the surface of a given (hkl) face of enstatite, the relation then becomes:

Figure 1

 1 Figure 1 Scheme of the numerical algorithm of the probabilistic model.

Fig. 2 .

 2 Fig.2. Results of the different simulations for the (100) face. Each graph represents the evolution of the steadystate dissolution rate as a function of the logarithm of: (A) P Mg-O-Mg , (B) P Mg-O-Si and (C) P Si-O-Si . Each color represents a different ∆Ea between the two other probabilities which remains fixed. Slopes and intercepts are listed in Table1. Legends correspond to the code indicated for each group of simulations. These codes are used in

  varies). The slope (a), intercept (b) and correlation coefficient (R²) of the linear regression of 𝐥𝐨𝐠(𝐫) = 𝐟(𝐥𝐨𝐠(𝐏)) are given in the last three columns (see text for details). When P Mg-O-Mg was left constant, its value was set to 0.99. This value allows for minimizing the simulation time required while having the possibility to associate a corresponding physical ∆E a value. When the relation between dissolution rate and P Mg-O-Mg was studied, P Mg-O-Si was arbitrarily set to 0.3. Finally, the chosen probabilities were constrained by P Mg-O-Mg > P Mg-O-Si > P Si-O-Si . In total, 345 simulations have been conducted (115 per face) to explore a wide range of activation energies. Over these 115 simulations, 30 were dedicated to unravelling the dependence of the dissolution rate evolution on PMg-O-Mg , 25 on PMg-O-Si , and 60 on PSi-O-Si .

Fig. 3 .

 3 Fig.3. Distribution of the different coordination spheres of an atom right before its dissolution. The first two lines depict Mg atoms, and the following two lines, Si. The proportion is defined by the number of atoms at the surface having m or n bounds divided by the total number of atoms, so that the sum of all modes is 1. Blue and red bars correspond to a set of low (P Mg-O-Mg = 0.5, P Mg-O-Si = 0.3 and P Si-O-Si = 0.0109) and high (P Mg-O-Mg = 0.99, P Mg-O-Si = 0.3 and P Si-O-Si = 0.1547) probability values, respectively. Magenta, red and green rectangle stand for the simulations on (100), (010) and (001) faces respectively.

3 .

 3 The coordination spheres will be referred as M-(O-Mg)n or M-(O-Si)m in the following (n and m standing for the number of Mg and Si neighbors, respectively). Whereas the coordination of Si-(O-Si)m is dominated by a single mode for which m = 1, others are more dispersed around one or more principal modes. Regarding Mg-(O-Mg)n, n is always an even number, regardless of the set of probabilities, consistent with the structure of enstatite. The value of the principal modes of Mg-(O-Mg)n, Mg-(O-Si)m, and Si-(O-Mg)m depends on the probability and differs from one face to another, as illustrated in Fig. 3, where the blue and red bars corresponds to two distinct sets of probabilities ([PMg-O-Mg = 0.99, PMg-O-Si = 0.3, PSi-O-Si = 0.0109] and [PMg-O-Mg = 0.7, PMg-O-Si = 0.3 and PSi-O-Si = 0.1547]) selected from the slowest and fastest range of dissolution rates respectively. However, these distributions also reveal that the mean coordination spheres remain constant, regardless of the set of probabilities used to run the simulation or the face that is considered (4.02 ± 0.01; 6.64 ± 0.02 and 1.002 ± 0.002 for Mg-O-Mg; Mg-O-Si+Mg-O-Si and Si-O-Si bonds, respectively).

Fig 4 .

 4 Fig 4. Contribution of dissolution rate specific to dislocations (𝒓 𝒅𝒊𝒔𝒍𝒐𝒄𝒂𝒕𝒊𝒐𝒏 (𝟏𝟎𝟎)

  Simulations containing dislocations have been conducted on the (100) face. The presented result show the relation between the dissolution rate and the dislocation density after subtraction of the "bulk" contribution to the overall dissolution rate (i.e., The results highlight an asymptotic relationship between dissolution rate and dislocation density.

Fig 5 .

 5 Fig 5. Comparison between the dissolution rates estimated with the surrogate model (r theory ; Eq. 5) and

P

   and A B P  in the range of the numerical values used for the experiments. An appropriate choice of   ,

Figure 6 .

 6 Figure 6. Comparison of experimental dissolution rates (r mod ) and dissolution rates determined with the surrogate model (r theory ) with the presence of one dislocation (dislocation density is changed by varying the surface area).

Fig. 5 .

 5 Bolded lines represent the input parameters used in the simulations shown in Fig.3.

Table 1 .

 1 Different ∆𝐄 𝐚 values used in this study (the "-" symbol indicates that ∆𝐄 𝐚

	a	b	R²
	mol)

Table 2 ,

 2 the , ,    values are obtained by taking the average values of the corresponding slopes and the dissolution rates calculated with the surrogate model (Eq. 5) are compared with the simulated one in Fig.5.

	Face	k	α	β	γ
	(100)	8.87	2.46	4.62	0.94
	(010)	7.80	2.22	4.65	0.75
	(001)	4.72	6.36	4.19	0.14

Table 2 .

 2 Values of the different face-specific parameters of the surrogate model.

Table 3 .

 3 Values of the parameters of the model described by Eq. 17 and Eq. 18 and the correlation coefficient R 2 . The bold characters represent the fitted parameters.

	R 2

Table 4 .

 4 Dislocation densities and probabilities used in the simulations used for

	Dislocation density	MgOMg	MgOSi	SiOSi
	1.99E+11	0.99	0.3	0.0109
	1.99E+11	0.99	0.25	0.0096
	1.99E+11	0.99	0.9	0.0096
	1.99E+11	0.99	0.25	0.0091
	1.99E+11	0.99	0.25	0.0669
	4.88E+10	0.99	0.3	0.0109
	4.88E+10	0.99	0.25	0.0096
	4.88E+10	0.99	0.9	0.0096
	4.88E+10	0.99	0.25	0.0091
	4.88E+10	0.99	0.25	0.0669
	2.21E+10	0.99	0.3	0.0109
	2.21E+10	0.99	0.25	0.0096
	2.21E+10	0.99	0.9	0.0096
	2.21E+10	0.99	0.25	0.0091
	2.21E+10	0.99	0.25	0.0669
	1.21E+10	0.99	0.3	0.0109
	1.21E+10	0.99	0.25	0.0096
	1.21E+10	0.99	0.9	0.0096
	1.21E+10	0.99	0.25	0.0091
	1.21E+10	0.99	0.25	0.0669
	3.00E+09	0.99	0.3	0.0109
	3.00E+09	0.99	0.25	0.0096
	3.00E+09	0.99	0.9	0.0096
	3.00E+09	0.99	0.25	0.0091
	3.00E+09	0.99	0.25	0.0669
	1.9899E+11	0.99	0.5	0.0096
	1.9899E+11	0.8	0.3	0.0109
	1.9899E+11	0.7	0.3	0.0109
	4.8773E+10	0.99	0.5	0.0096
	4.8773E+10	0.8	0.3	0.0109
	4.8773E+10	0.7	0.3	0.0109
	2.211E+10	0.99	0.5	0.0096
	2.211E+10	0.8	0.3	0.0109
	2.211E+10	0.7	0.3	0.0109

Finally, this study represents a first step towards the development of surrogate expressions of mineral dissolution rates. Its application to natural environment, its extension to other groups of minerals as well as its extension to transient state still have to be explored.

Acknowledgements

A.B. thanks the University of Strasbourg for funding his PhD grant. We are also grateful for the careful reviews and detailed suggestions made by two anonymous reviewers, which improved an earlier version of the present paper.

TOC Graphic