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

Damien Daval

François Guyot

Philippe Ackerer

The dissolution anisotropy of pyroxenes: Experimental validation of a stochastic dissolution model based on enstatite weathering

Keywords: stochastic model, pyroxene dissolution anisotropy, dissolution kinetics

teaching and research institutions in France or abroad, or from public or private research centers.   The dissolution anisotropy of pyroxenes: Experimental validation of a stochastic dissolution model based on enstatite weathering

Introduction

The interactions between aqueous fluids and silicates play a central role in many major geological, environmental and industrial processes. Silicate dissolution and precipitation reactions control on various space and time scales many aspects of the continental surfaces such as providing nutrients to living organisms [START_REF] Huntington | Calcium Depletion in a Southeastern United States Forest Ecosystem[END_REF] , participating to the surface water chemistry, mobilizing [START_REF] Frugier | Son68 Nuclear Glass Dissolution Kinetics: Current State of Knowledge and Basis of the New Graal Model[END_REF] or sequestrating natural and anthropic contaminants [START_REF] Johnson | Acid Mine Drainage Remediation Options: A Review[END_REF] , shaping landscapes [START_REF] Brantley | Twelve Testable Hypotheses on the Geobiology of Weathering[END_REF] and ultimately regulating the long-term carbon cycle and climate [START_REF] Hir | Scenario for the Evolution of Atmospheric Pco(2) During a Snowball Earth[END_REF][START_REF] Beaulieu | High Sensitivity of the Continental-Weathering Carbon Dioxide Sink to Future Climate Change[END_REF] . On an equal footing, predicting the reactivity of silicate surfaces is further important for industrial purposes, such as cement hydration [START_REF] Nicoleau | Analytical Model for the Alite (C3s) Dissolution Topography[END_REF][START_REF] Robin | Experimental Study and Numerical Simulation of the Dissolution Anisotropy of Tricalcium Silicate[END_REF] and kinetics of catalyzed reactions [START_REF] Kurganskaya | A Comprehensive Stochastic Model of Phyllosilicate Dissolution: Structure and Kinematics of Etch Pits Formed on Muscovite Basal Face[END_REF] . As a consequence, the dissolution rates of silicates have now been investigated for more than half a century, giving rise to a dense body of literature (e.g. Brantley et al. [START_REF] Brantley | 7.3 -Reaction Kinetics of Primary Rock-Forming Minerals under Ambient Conditions A2 -Holland, Heinrich D[END_REF] for a review).

From a mechanistic standpoint, the conventional treatment of mineral reactivity has long been based on the fundamental assumption that silicate dissolution may be described within the framework of the transition state theory (TST), early developed by Eyring 11 for elementary and homogeneous reactions. The extension of TST by Aagaard and Helgeson [START_REF] Aagaard | Thermodynamic and Kinetic Constraints on Reaction-Rates among Minerals and Aqueous-Solutions .1. Theoretical Considerations[END_REF] to heterogeneous fluid-solid reactions supposes that silicate dissolution reactions are controlled by a unique rate-limiting step, implicitly requiring that the surface reactivity is homogeneous, isotropic and immutable for a given fluid composition [START_REF] Noiriel | Pore-Scale Geochemical Reactivity Associated with Co2 Storage: New Frontiers at the Fluid-Solid Interface[END_REF] . Although the resulting 'kinetic rate laws' have been implemented into reactive transport codes for several decades, an increasing number of studies has consistently challenged this conceptual framework: First, it has been theoretically and experimentally demonstrated that mineral reactivity is highly heterogeneous, resulting from the distribution of reactive sites at the fluid-solid interface [START_REF] Luttge | Kinetic Concepts for Quantitative Prediction of Fluid-Solid Interactions[END_REF] . The 'rate spectra' concept [START_REF] Fischer | How Predictable Are Dissolution Rates of Crystalline Material?[END_REF] , which provides a measurement of the rate distribution at the mineral surface in the spatial domain, was developed to indirectly capture this energetic landscape.

Second, the advent of nanoscale imaging of surface topography has resulted in the demonstration that crystal dissolution is a highly anisotropic process [START_REF] Robin | Experimental Study and Numerical Simulation of the Dissolution Anisotropy of Tricalcium Silicate[END_REF][START_REF] Gratz | Dissolution of Quartz in Aqueous Basic Solution, 106-236°C: Surface Kinetics of "Perfect" Crystallographic Faces[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] 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] 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] Perez | Comparing the Reactivity of Glasses with Their Crystalline Equivalents: The Case Study of Plagioclase Feldspar[END_REF] . These studies showed that the anisotropic reactivity of minerals results from the anisotropic distribution of atomic positions in the crystal lattice, and may be described in the framework of the periodic bond chain (PBC) theory [START_REF] Hartman | On the Relations between Structure and Morphology of Crystals[END_REF] . This theory suggests that crystal faces may be sorted following the number of chains of energetically strong bonds (PBCs) they contain, with the faces containing the largest numbers of PBCs dissolving slower, and vice versa. Pollet-Villard et al. [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] further proposed that the recognition of the anisotropy of crystal reactivity challenged the idea of a unique relation between the overall dissolution rate of a crystal and the Gibbs free energy of reaction (∆G), which is constitutive of the TST formulation, while questioning the existence of a unique intrinsic dissolution rate constant for crystals.

Finally, silicate dissolution is accompanied with gradual physicochemical and structural modifications of their surface, resulting in the formation of nm-to µm-thick silica-rich surface layers [START_REF] Daval | Mechanism of Wollastonite Carbonation Deduced from Micro-to Nanometer Length Scale Observations[END_REF][START_REF] Kwak | Metal Carbonation of Forsterite in Supercritical Co2 and H2o Using Solid State 29si, 13c Nmr Spectroscopy[END_REF][START_REF] Gin | Nuclear Glass Durability: New Insight into Alteration Layer Properties[END_REF][START_REF] Gin | Origin and Consequences of Silicate Glass Passivation by Surface Layers[END_REF][START_REF] Hellmann | Nanometre-Scale Evidence for Interfacial Dissolution-Reprecipitation Control of Silicate Glass Corrosion[END_REF][START_REF] Daval | Time-Dependent Feldspar Dissolution Rates Resulting from Surface Passivation: Experimental Evidence and Geochemical Implications[END_REF] . Accounting for the transport properties of such layers, especially when developed under mildly aggressive conditions (circum-neutral pH, low temperature), was also proposed to be essential to develop reliable models of silicate dissolution [START_REF] Daval | Dynamics of Altered Surface Layer Formation on Dissolving Silicates[END_REF][START_REF] Rébiscoul | Water Dynamics in Nanoporous Alteration Layer Coming from Glass Alteration: An Experimental Approach[END_REF] .

In a nutshell, all these studies converge towards the same conclusion: providing a better theoretical description of mineral dissolution crucially requires gaining knowledge on the relation between the atomic configuration of (and therefore, the energetic distribution at) the dissolving mineral surface and the overall dissolution rate of a given crystal.

Driven by this assertion, studies dedicated to the molecular-scale description of the silicate dissolution process following ab initio quantum mechanical calculations 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] . Such studies shed light on essential theoretical information regarding the dissolution process such as suggesting reaction paths, determining transition state configurations (i.e., bond lengths and angles) and calculating activation energies of the hydrolysis reactions in various pH domains. They also independently confirmed that the dissolution process cannot be boiled down to the kinetics of a single surface reaction.

Actually, the large number of possible configurations of dangling atoms at the silicate surface over the course of the dissolution translates into as many possible different hydrolysis reactions, justifying the need for specific numerical simulations of the dissolution process at a scale exceeding clusters of a few atoms. This was precisely one of the motivations for the development of kinetic models using a Monte Carlo approach (often abbreviated as kMC) [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 Approach to Study Carbonate Dissolution[END_REF][START_REF] Lasaga | Mineralogical Approaches to Fundamental Crystal Dissolution Kinetics[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] Zhang | Aluminosilicate Dissolution Kinetics: A General Stochastic Model[END_REF][START_REF] Meakin | Simple Kinetic Monte Carlo Models for Dissolution Pitting Induced by Crystal Defects[END_REF][START_REF] Cailleteau | Insight into Silicate-Glass Corrosion Mechanisms[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] , whose parametrization may strongly rely on the results of ab initio calculations mentioned above. However, with the noticeable exception of a few recent 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 Approach to Study Carbonate Dissolution[END_REF][START_REF] Kurganskaya | Kinetic Monte Carlo Simulations of Silicate Dissolution: Model Complexity and Parametrization[END_REF] , to date, most of these modelling studies were conducted with simple Kossel crystals [START_REF] Lüttge | A Stochastic Treatment of Crystal Dissolution Kinetics[END_REF] , making difficult the direct comparison of the outputs of the simulations with experimental data.

In the present study, we combined experimental measurements of the anisotropy of silicate dissolution (face-specific topographic features and dissolution rates) with stochastic simulations to critically assess the validity of a simple dissolution model. Because of its very anisotropic crystallographic structure made of chains of silica tetrahedrons running parallel to one of the main crystallographic axes, pyroxene was chosen as a case study. It was previously shown that the anisotropic structure of this family of silicates (also known as inosilicates) resulted in a strong anisotropy of dissolution [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] Daval | Dynamics of Altered Surface Layer Formation on Dissolving Silicates[END_REF][START_REF] Lorenzo | The Carbonation of Wollastonite: A Model Reaction to Test Natural and Biomimetic Catalysts for Enhanced Co2 Sequestration[END_REF] , intensifying the differences of reactivity between the various faces that define their usual crystal habit. Owing to its simple chemical composition (MgSiO3) and structure (orthorhombic), enstatite was chosen as a model within the pyroxene group.

After a brief description of the methods used and the ensuing information, we will discuss the experimental results (i.e., the anisotropy of enstatite dissolution rates and surface features observed using electron microscopy) in regard of pre-existing studies based on pyroxene dissolution. These results will then be compared to the outputs of the simulations, which will lead to a validation of the dissolution model and a discussion of the wider implications and conclusions for pyroxene dissolution and dissolution mechanisms at the fluid-silicate interface.

Materials and methods

Starting materials and preparation

The enstatite samples used in this study come from Kilosa (Tanzania), and were supplied by Mawingu Gems Company. Negligible amounts of clinoenstatite (around a few wt. %) were detected using X-ray diffraction (Fig. S1). Monolithic enstatite samples (up to several tens of mm on a side) were cut with a diamond blade saw and polished down to the nanometer scale through a multi-step abrasive sequence. The crystallographic orientation of each sample was determined using an electron backscatter diffraction (EBSD) detector attached to a scanning electron microscope (SEM) (TESCAN® Vega 2). The initial roughness of each surface was measured at different scales with vertical scanning interferometer (VSI; ZYGO® NewView 7300). For each sample, the initial average arithmetic roughness (Ra), defined as the arithmetic average of the absolute values of the roughness profile, ranged between 5 nm and 15 nm, and Rq (the roughness calculated with the root-mean-square of the vertical distance between the surface elevation and the mean surface line) ranged between 19 nm and 70 nm.

These initial roughness parameters were measured on 270 × 360 µm 2 images obtained using VSI.

Dissolution experiments and analytical procedures

Oriented samples were mounted in a passivated titanium support as designed by Daval et al. [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] , where a Viton® disk held with a screw creates a non-wetted reference area. Partially masked samples were dissolved in a synthetic acidic solution in 120 ml-PFA reactors at 90 °C for durations ranging from 6 hours to 3 weeks to check for the linearity of the surface retreat as a function of time.

The aqueous solutions were made with ultrapure deionized water (18.2 MΩ.cm) and analytical grade nitric acid. The pH was adjusted to 0.0 ± 0.2. Nitric acid was chosen to avoid the dissolution of the titanium support that was observed in preliminary tests conducted with hydrochloric acid under similar conditions. Solution sampling was regularly carried out for the analyses of Si, Mg and Fe using ICP-AES (Thermo ICAP 6000). The CHESS code [START_REF] Van Der Lee | Chess Tutorial and Cookbook[END_REF] was used to determine the in situ pH and saturation indices with respect to enstatite and secondary phases. It was verified that all solutions were undersaturated with respect to any secondary phase.

For each orientation, at least two solid samples were analyzed for topography measurements. While the first one was periodically recovered to perform the analyses and put back in reactor, the second sample was left in the reactor all over the duration of the experiment. This protocol was followed to make sure that the reactivity of samples that experience wetting/drying cycles to get time-resolved topography measurements remains similar to the reactivity of samples continuously immersed in solution. The topography measurements were performed using VSI in stitching mode (magnification: ×5) to estimate the average changes in height between the unreacted reference surface and the reacted mineral surface. For each oriented sample, dissolution rates were calculated according to the following equation (see e.g. Lüttge et al. [START_REF] Lüttge | An Interferometric Study of the Dissolution Kinetics of Anorthite; the Role of Reactive Surface Area[END_REF] ):

𝑟 ( ) = ∆ℎ ∆𝑡 𝑉 (1)
where 𝑟 ( ) (mol.m -2 .s -1 ) is the dissolution rate of the (hkl) face, Δh is the surface retreat (m), Δt is the alteration duration (s) and 𝑉 (m 3 .mol -1 ) is the molar volume of enstatite. As detailed below, the minute lateral extension of etch pits (~ a few tens of nm) prevented the investigation of individual pits by VSI even at high magnification (×100; lateral resolution ~ e.g. rate spectra.

After completion of the surface topography measurements, the formation of nm-thick, amorphous Si-rich layers (ASSL) was verified to make sure that the layer thickness was negligible compared to the corresponding surface retreat measured by VSI. The samples were either platinum-or carbon-coated and thin sections perpendicular to the investigated faces were prepared by focused ion beam (FIB) milling using the FEI Helios 600 Nanolab dualbeam operating at CP2M (Marseille, France). The FIB Ga ion milling was carried out at an ion beam voltage of 30 kV and with beam currents from 3 nA to 10 pA for the final steps.

Micrometer-thick sections were lifted out in situ using an Omniprobe 200 micromanipulator and transferred to a half copper grid for final ion milling to electron transparency (final thickness of about 100 nm). High-resolution transmission electron microscopy (HRTEM) and scanning transmission electron microscopy (STEM) observations were performed on FIB foils using a JEOL 2100F microscope operating at 200 kV, equipped with a field emission gun. Energy dispersive X-ray (EDX) spectra were acquired in STEM mode to probe the chemical compositions, with a focused electron beam (1 nm). When the FIB thin section intersected etch pits (i.e., for the (100) and (010) samples), TEM was further used to characterize their morphologies.

Modeling strategy

The stochastic approach for mineral dissolution studies has emerged in the mid-80's [START_REF] Lasaga | Surface Chemistry, Etch Pits and Mineral-Water Reactions[END_REF] .

Since then, numerous studies have used stochastic simulations to describe mineral dissolution at the atomic scale, principally using simple Kossel crystals [START_REF] Meakin | Simple Kinetic Monte Carlo Models for Dissolution Pitting Induced by Crystal Defects[END_REF][START_REF] De Assis | Dissolution of Minerals with Rough Surfaces[END_REF][START_REF] Lüttge | A Stochastic Treatment of Crystal Dissolution Kinetics[END_REF] (simple cubic crystal lattice; 6-fold coordinated atoms) but some have also applied them on more complex mineral 9, 22, 41 and glass [START_REF] Devreux | Leaching of Borosilicate Glasses. Ii. Model and Monte-Carlo Simulations[END_REF][START_REF] Cailleteau | Insight into Silicate-Glass Corrosion Mechanisms[END_REF][START_REF] Cailleteau | Why Do Certain Glasses with a High Dissolution Rate Undergo a Low Degree of Corrosion?[END_REF] structures. These studies have shown the robustness of these models through their ability to reproduce etch pit shapes [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 Approach to Study Carbonate Dissolution[END_REF] and sometimes, by comparing modeled and measured dissolution rates [START_REF] Kurganskaya | Kinetic Monte Carlo Approach to Study Carbonate Dissolution[END_REF] .

These models are based on a simple relation that links the bond-breaking probability to the activation energy of the hydrolysis of this bond [START_REF] Lasaga | Surface Chemistry, Etch Pits and Mineral-Water Reactions[END_REF] :

𝑃 = 𝑒 ∆ ( 2
)
where P is the bond-breaking probability, ∆Ea, the activation energy (J), kB, the Boltzmann constant (J/K) and T the temperature (K).

Modeling enstatite crystal structure

Enstatite is an inosilicate made of chains of silicate tetrahedrons running parallel to the caxis. The crystal lattice of orthoenstatite (used in the present study) is orthorhombic and belongs to the Pbca space group. Hugh-Jones et al. [START_REF]A Compressional Study of Mgsio3 Orthoenstatite up to 8.5 Gpa[END_REF] have provided the coordinates of the atomic positions (Table 1) that allow forming the enstatite cell with the following relations corresponding to the Pbca system:

x y z ½ -x ½ + y z x ½ -y ½ + z ½ + x y ½ -z (3) -x -y -z ½ + x ½ -y -z -x ½ +y ½ -z ½ -x -y ½ +z
The lengths of the a, b and c axes are 18.233, 8.8191 and 5.1802 Å respectively and the angles α, β and γ are equal to 90° (orthorhombic system).

The enstatite unit cell was created based on these parameters, and is composed of 16 atoms of magnesium, 16 atoms of silicon and 48 atoms of oxygen. The cell is then repeated following the three space dimensions in order to get a simulated crystal of several tenths to hundreds of nanometers on a side.

Bond-breaking probabilities

Most stochastic simulations performed at the atomic scale on silicate materials consider the cleavage of M-O-M bonds instead of single M-O bonds [START_REF] Kurganskaya | A Comprehensive Stochastic Model of Phyllosilicate Dissolution: Structure and Kinematics of Etch Pits Formed on Muscovite Basal Face[END_REF] . Hence, the oxygen atoms are only considered as bridges between two cations (Obr), and the cations connected to the oxygens of a specific Si or Mg atoms are referred to as their first coordination spheres [START_REF] Kurganskaya | A Comprehensive Stochastic Model of Phyllosilicate Dissolution: Structure and Kinematics of Etch Pits Formed on Muscovite Basal Face[END_REF] . In enstatite, three types of bonds exist: Mg-Obr-Mg, Mg-Obr-Si and Si-Obr-Si bonds.

Two different strategies may be considered to develop stochastic models: the former considers that the release of a cation results from the sequential cleavage of each bond which connects it to the surface, and the latter considers that the cation release results from the simultaneous cleavage of all these bonds at a given iteration step. The difference between the two methods stands in the number of information required to run a simulation (for a single enstatite cell, 160 bonds need to be considered following the first strategy whereas only 32 atoms (16 Mg and 16 Si atoms; O atoms being considered as "bridges") are required for the second strategy). For this reason, the second strategy is generally used and has previously proven successful [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 Approach to Study Carbonate Dissolution[END_REF][START_REF] Zhang | Si Order in Albite and Its Effect on Albite Dissolution Processes: A Monte Carlo Study[END_REF][START_REF] Cailleteau | Insight into Silicate-Glass Corrosion Mechanisms[END_REF][START_REF] Cailleteau | Why Do Certain Glasses with a High Dissolution Rate Undergo a Low Degree of Corrosion?[END_REF] . Accordingly, our simulations were run using the following probability for the release of Mg and Si atoms:

𝑃 = 𝑒 ∆ 𝑒 ∆ (4)
where M = Mg or Si and n and m correspond to the number of Mg and Si atoms linked to the considered M atom, respectively. This equation, together with knowledge of the connectivity of each atom, represents the core description of the physical basis adopted for our stochastic simulations. Importantly, more sophisticated relations may be found in the literature, which attempt to account for various physical processes such as lattice resistance [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] Kurganskaya | Kinetic Monte Carlo Simulations of Silicate Dissolution: Model Complexity and Parametrization[END_REF] , steric constraints [START_REF] Kurganskaya | A Comprehensive Stochastic Model of Phyllosilicate Dissolution: Structure and Kinematics of Etch Pits Formed on Muscovite Basal Face[END_REF] or backward attachment reactions possibly resulting in self-healing of disrupted Si-Obr-Si bonds [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] 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] Cailleteau | Insight into Silicate-Glass Corrosion Mechanisms[END_REF] . Importantly, even if the resulting stochastic models are physically more robust, they also introduce significant degrees of freedom in the parametrization of the resulting bond-breaking probabilities (see bond-breaking probability expressions and parametrization derived in, e.g., [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] Devreux | Leaching of Borosilicate Glasses. Ii. Model and Monte-Carlo Simulations[END_REF] ). In contrast, the present work was intended to assess the extent to which a stochastic approach with the simplest physical basis can satisfactorily account for the main face-specific features of enstatite dissolution, explaining why we used the simple relations provided by Eq. ( 4).

Various algorithms may be used to run stochastic simulations. In the present study, we followed a strategy close to the "random select and test" approach previously described and used in e.g. [START_REF] Devreux | Leaching of Borosilicate Glasses. Ii. Model and Monte-Carlo Simulations[END_REF][START_REF] Meakin | Simple Kinetic Monte Carlo Models for Dissolution Pitting Induced by Crystal Defects[END_REF] (see section 2.3.4). As described in these studies, this approach consists in scaling to unity the most probable event, and scaling accordingly the probabilities of all other events. Equation (4) then simply becomes:

𝑃 = 𝑃 × 𝑃 (5

Dislocations

Etch pits are known to nucleate at dislocations outcropping at crystal surfaces [START_REF] Gilman | Dislocation Etch Pit Formation in Lithium Fluoride[END_REF] . In order to compare the observed and simulated morphology of etch pits on each face, screw dislocations were introduced in the model. The main dislocations for enstatite follow the

[001], [010] and [100] axes [START_REF] Skrotzki | Defect Structure and Deformation Mechanisms in Naturally Deformed Augite and Enstatite[END_REF] . Therefore, screw dislocations are perpendicular to the different studied orientations, except for the (210) face. In order to propagate the dislocation in the volume, the distances between the first atom (located at the center of the reactive surface) and its neighbors in the xOy plan are calculated. The shortest distance defines the atom that will be the next one belonging to the dislocation line. The same process is repeated until reaching the opposite face of the simulated crystal. For the (210) face, a rotation of the volume is made before applying these steps in order to respect the dislocation directions ([100] and [010]).

Finally, when all the atoms that belong to the dislocation line are identified, their bondbreaking probabilities are set to 1 9, 22 .

Dissolution algorithm

The "random select and test" (RST [START_REF] Devreux | Leaching of Borosilicate Glasses. Ii. Model and Monte-Carlo Simulations[END_REF][START_REF] Meakin | Simple Kinetic Monte Carlo Models for Dissolution Pitting Induced by Crystal Defects[END_REF] ) and "divide and conquer" approaches (DaC [START_REF] Meakin | Simple Kinetic Monte Carlo Models for Dissolution Pitting Induced by Crystal Defects[END_REF][START_REF] Bortz | A New Algorithm for Monte Carlo Simulation of Ising Spin Systems[END_REF] ) are two different algorithms commonly used in stochastic simulations. Nowadays, the DaC algorithm is usually preferred because it requires less computer time. Furthermore, it provides a correspondence between iteration and time [START_REF] Kurganskaya | A Comprehensive Stochastic Model of Phyllosilicate Dissolution: Structure and Kinematics of Etch Pits Formed on Muscovite Basal Face[END_REF] . However, DaC algorithms cannot reproduce simultaneous events, and favor a "layer-by-layer" removal of atoms [START_REF] Zhang | Si Order in Albite and Its Effect on Albite Dissolution Processes: A Monte Carlo Study[END_REF] . For these reasons, we preferred using an approach close to the "primitive" RST despite its poor computation efficiency. Moreover, time cannot be associated a priori to a number of iterations. This relation has to be inferred using at least two experimental data sets (e.g., two dissolution rates for two different faces). By measuring mean dissolution rates for 4 different faces, we made sure that the problem was overconstrained, ensuring a well-established relationship between time and simulation steps (iterations).

The crystal is described by the coordinates of the Mg and Si atoms and the connectivity of each atom to its neighbors follows the rules described in Section 2.3.1. Only the cations that have an incomplete first coordination sphere (incomplete connectivity), are considered to be located at the crystal surface, and may experience dissolution. For each Mg and Si atom, an integer number indicates the three possible states of the atom: (i) dissolved,

(ii) at the crystal surface, or (iii) inside the crystal. This indicator allows counting the number of dissolved cations and their nature (Mg or Si) and identifying the atoms located at the surface in an efficient way.

The algorithm for simulating the dissolution process is as follows:

1. At each iteration step k, all cations of the simulated crystal volume are scanned sequentially (Loop 1 in Fig. S2).

2. For each surface cation, a random number Z uniformly distributed between 0 and 1 is generated. The probability of dissolution of that atom PM is computed as the product of the bond-breaking probabilities of the bonds that link it to the crystal following Eq. ( 5). If Z < PM, the atom is released (dissolved). If not, the scan of the atoms is moving on (Loop 2 in Fig. S2).

3. When the scan of all atoms is completed, the first coordination spheres of the remaining atoms are updated if necessary (i.e., when a removed atom was belonging to the coordination spheres of one of the remaining atoms; cf. Loop 3 in Fig. S2). Because the coordination spheres are modified, the probability of dissolution PM of the corresponding atoms are increased.

4. The number of released Si and Mg atoms at iteration k is saved for postprocessing.

5. The iteration number is incremented and the process continues step 1.

Model outputs

The outputs of the simulations were used to document face-specific features of enstatite dissolution, such as (i) the mean dissolution rates at steady-state conditions, (ii) the thickness of the Si-rich surface layer that may result from preferential leaching of Mg (if any) and (iii) the etch pit morphology.

The dissolution rates may be calculated following two different approaches. First, the list of surface cations after a given iteration is used to calculate the corresponding mean becomes constant), this approach can be related to the dissolution rates estimated using VSI (see section 2.2). Of note, the difference in terms of surface areas between the simulated (0.33 µm² for the largest) and measured (dozens of mm²) surfaces preclude any direct use of the rate spectra concept [START_REF] Fischer | How Predictable Are Dissolution Rates of Crystalline Material?[END_REF] , so that in the present study, only the mean surface retreat rates were subsequently used as a criterion to assess the correspondence between simulated and measured dissolution rates.

A second approach considers the amount of Mg and Si atoms released at each iteration, as could be determined experimentally using ICP-AES solution analyses. where it represents the iteration step at which the calculation is performed, NM,i, the number of Mg or Si atoms released after the i th iteration, NM,cell, the number of Mg or Si atoms in the cell (16), Vcell the enstatite cell volume, 𝑉 the molar volume of enstatite and S, the geometric surface area of the reactive surface.

While this approach eventually yielded results very close to the previous method, it further provided information regarding the stoichiometry of the dissolution by comparing the evolution of the departure of the Mg and Si atoms, and eventually, regarding the existence and thickness of a Si-rich surface layer, which was previously documented for inosilicates weathered experimentally [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] Daval | Dynamics of Altered Surface Layer Formation on Dissolving Silicates[END_REF][START_REF] Petit | Mechanism of Diopside Dissolution from Hydrogen Depth Profiling[END_REF][START_REF] Ruiz-Agudo | Mechanism of Leached Layer Formation During Chemical Weathering of Silicate Minerals[END_REF][START_REF] Hellmann | Unifying Natural and Laboratory Chemical Weathering with Interfacial Dissolution-Reprecipitation: A Study Based on the Nanometer-Scale Chemistry of Fluid-Silicate Interfaces[END_REF] . Since the model does not consider the backward reaction of Si attachment to the surface (see discussion in [START_REF] Zhang | Si Order in Albite and Its Effect on Albite Dissolution Processes: A Monte Carlo Study[END_REF][START_REF] Zhang | Theoretical Approach to Evaluating Plagioclase Dissolution Mechanisms[END_REF] ), it only indicates a lower bound for the thickness of the layer resulting from cation depletion, ignoring possible Si re-deposition following an interfacial dissolution-reprecipitation mechanism.

Finally, for simulations run with dislocations, the resulting etch pits and surface retreats were calculated using Surfer©. The coordinates of the atoms at the surface were stored and used as an input file for Surfer. The minimum curvature method was used to calculate the mean surface retreat at a given iteration.

Linking time to iteration

As mentioned above, simulation results are normalized to the number of iterations instead of elapsed time. However, in order to compare modeled and measured dissolution rates, it is necessary to link a single iteration to a physical time increment. The correspondence between time and iteration was determined following Eq. ( 7):

𝑡 = 〈 𝑟 ( ) 𝑟 ( ) 〉 𝑖𝑡 (7)

Results

Face-specific dissolution rates of enstatite revealed experimentally

Fig. 2A shows the temporal evolution of the mean surface retreat measured using VSI.

Whereas the (hk0) faces were reacted for durations of up to three weeks, the rapid increase of the roughness of the (001) face prevented the measurements of surface topography for durations exceeding a couple of days. The main results are as follows:

(i) the dissolution rate is constant through time (at least from the first experimental measurements) for all of the four faces, within experimental uncertainties;

(ii) the dissolution rates of the samples left in the reactor all over the duration of the experiments are similar to those that were periodically retrieved from the reactors;

(iii) enstatite dissolution is an anisotropic process (i.e., the slopes of the linear regressions differ from one face to the other). The face-specific dissolution rates (𝑟

( ) )
observe the following trend:

𝑟 ( ) ≫ 𝑟 ( ) > 𝑟 ( ) ≥ 𝑟 ( ) .
The mean values of dissolution rates are 71.8 ± 13.9, 8.5 ± 1.5, 3.0 ± 1.0 and 2.6 ± 0.6 nmol/m²/s for the (001), (210), (010) and (100) faces, respectively (Fig. 2B).

Face-specific microstructural features revealed experimentally

The etch pits possibly formed on each surface were too small to be properly investigated by VSI, whose resolution is close to 300 nm at the maximum 100xmagnification. Apart from the surface retreat between the masked and unmasked portions of the samples (Fig. 3), no specific topography feature could be evidenced using VSI. SEM imaging of the reacted surfaces revealed that etch pits did nucleate on three of the four investigated faces (Fig. 4), and confirmed that their lateral extension was too modest to be imaged by VSI. The (100) and (010) faces exhibited similar lenticular pits elongated following the c-axis (the major axis is more than 6 times larger than the minor axis). The shape of the pits projected to the horizontal surface shows a central symmetry, consistent with the shape previously described in the literature for pyroxenes in general [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] Ruiz-Agudo | Mechanism of Leached Layer Formation During Chemical Weathering of Silicate Minerals[END_REF][START_REF] Berner | Dissolution of Pyroxenes and Amphiboles During Weathering[END_REF][START_REF] Velbel | Pyroxene Weathering to Smectite: Conventional and Cryo-Field Emission Scanning Electron Microscopy, Koua Bocca Ultramafic Complex, Ivory Coast[END_REF][START_REF] Phillips-Lander | Silicate Dissolution in Las Pailas Thermal Field: Implications for Microbial Weathering in Acidic Volcanic Hydrothermal Spring Systems[END_REF] . Two types of pits are visible on the (210) face: while the former is similar to the ones observed on the (100) and ( 010) faces (with a small minor axis compare to the major axis), the latter, also elongated following the c-axis, presents a wider minor axis (3 times smaller than the major axis). In addition to this difference of minor-axis length, these pits differ from the ones observed on the other faces by their non-symmetrical, half-moon shape. Whereas the lowest points for the pits developed on the (100) and (010) faces are located at the center of the pits, conversely, for pits developed on the (210) face, the lowest point is located at the center of the major-axis but shifted on one edge of the minor-axis (Fig. 4C).

The dislocation density is also a function of crystallographic orientation. The number The dislocation density on the (100) and (210) faces was found to be similar (on the order of 4.9 ± 0.3 × 10 8 dislocations / cm²), and is about two orders of magnitude greater than that measured on the (010) face (9.8 × 10 6 dislocation / cm²). The dislocation density on the (001) face could not be determined because of the rapid roughening of the surface, which impeded the identification of individual etch pits. These dislocation densities are in agreement with other studies [START_REF] Olsen | Analysis of Dislocations in Some Naturally Deformed Plagioclase Feldspars[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] , which reported dislocation densities ranging between 4.0 x 10 7 and 7.0

of
x 10 8 dislocations / cm². The close correspondence between the dislocation density calculated for the (100) and the (210) faces may be explained by the main orientation of screw dislocations for the enstatite structure. As mentioned in section 2.3.3, the main dislocations in enstatite follow the [100], [010] and [001] directions [START_REF] Skrotzki | Defect Structure and Deformation Mechanisms in Naturally Deformed Augite and Enstatite[END_REF] . Therefore, outcropping dislocations at the (210) surface should represent the sum of the dislocation outcropping at the (100) and the (010) surfaces, which is consistent with our observations. TEM investigations performed on FIB thin sections prepared on the reacted samples revealed the occurrence of amorphous surface layers on (010), ( 210) and (001) samples, with a thickness specific to the surface orientations (Fig. 5). The thickest layers were evidenced on (001) samples (15-20 nm; cf. Figs. 5C-D). The layers developed on (210) and (010) samples have similar thickness (≤ 4 nm; cf. Fig. 5B), whereas no surface layer could be clearly evidenced on (100) samples (Fig. 5A). The layers developed on (001) samples were thick enough to analyze their chemical composition with STEM-EDX, which revealed that they are Si-rich and Mg-depleted (Fig. 5D). The similarity of these results with those reported by Daval et al. [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] for diopside dissolved at pH 1 is noteworthy. Finally, it was possible to determine the cross-section morphology of etch pits nucleated on the (100) and (010) faces, as the FIB thin sections were realized in locations where etch pits were particularly abundant.

Regarding the (010) face, TEM images revealed that the pit geometry observes an axial symmetry with respect to the [010] axis, and that the pit walls are particularly steep, forming an angle of 93° ± 2° with the (010) surface (Fig. 5B). Similarly, the walls of the pits developed on the (100) face formed an angle of 97° ± 5° with the (100) surface, although the pits appear to be flat-bottomed (Fig. 5A). Overall, these results confirmed that profilometry techniques such as atomic force microscopy or vertical scanning interferometry would have failed to reveal the 3D-geometry of such etch pits. A common trait to all of the simulation outputs, irrespective to the set of bondbreaking probabilities that have been tested, is the attainment of a steady-state dissolution regime, which is characterized by a congruent release of Mg and Si cations and a constant dissolution rate (Fig. 6). This steady-state regime was found to result from the attainment of a constant number of atoms at the surface of each face, as illustrated in Fig. 6 (color lines).

Model outputs

Dissolution stoichiometry and steady-state dissolution rates

During the initial transient period, Mg is released preferentially, and enstatite dissolution rate based on Si release increases gradually. The duration and the shape of the transient period depend on the set of bond-breaking probabilities that are used to run the simulations and is a function of surface orientation. In a general manner, the steady-state regime is achieved earlier for the (hk0) faces than for the (001) face, for which the initial non-congruence of the dissolution is also more pronounced (Fig. 6). As a consequence, thicker ASSLs are expected to build up on the (001) face, compared to the (hk0) faces. When one considers the set of bond-breaking probabilities that allows for the best agreement between the modeled and measured steady-state face-specific dissolution rates (see below), the modeled thickness of ASSL is on the order of 1-2 atomic monolayers for all (hk0) faces, as opposed to ~ 10 nm for the ASSL developed on the (001) face.

The modeled absolute values of steady-state dissolution rates and the relative differences between the reactivity of the four faces depend on the set of bond-breaking probabilities that were used to conduct the simulations. Considering a range of probabilities that is consistent with the literature (see section 2. Whereas the amplitude of the reactivity measured between the fastest and slowest dissolving face is well reproduced by the model, one can note that it slightly overestimates the reactivity of the (010) face (or slightly underestimates the reactivity of the (210) face).

Morphology of modeled etch pits

The implementation of dislocation lines in the simulated crystals allowed to document the morphology of etch pits for the four investigated faces. Fig. 7 depicts the modeled surface evolution of the four studied faces, centered on an outcropping dislocation line, for simulations ranging from 0 to 2,000 iterations. Screw dislocation orientations were chosen following those determined by Skrotzki [START_REF] Skrotzki | Defect Structure and Deformation Mechanisms in Naturally Deformed Augite and Enstatite[END_REF] for enstatite: [100] for (100) and (210) faces, [010] for the (010) face and [001] for the (001) face. These results were obtained using the set of bond-breaking probabilities corresponding to the best agreement between the modeled and measured steady-state face-specific dissolution rates.

Etch pits growing on the (100) face were found to be lenticular and symmetric with respect to the b-and c-axes, respectively (cross-sections 1 and 2). Their aspect ratio evolves as a function of the reaction progress, as their length following the c-axis increases faster than along the b-axis. In addition, they exhibit a single lower point that may result from a slow step propagation (and thus, difficult pit opening), as suggested by the shape of their bottom along the c-axis. The bottom of the pits developed on the (100) face are surrounded by steep pit walls (θ ≈ 93° ± 1° between the surface and pit walls).

The morphology of etch pits growing on the (010) face displays similar features (elongated following the c-axis, and symmetric with respect to the a-and c-axes, with similar length and width (≈5-10 nm) after 2,000 iterations), although the pits are deeper (≈ 20 nm against ≈ 10-15 nm after 2,000 iterations) and thinner than those developed on the (100) face.

Of note, their growth rate is lower than that of pits developed on the (100) face. The bottom of the pits developed on the (010) face is surrounded by steep pit walls (θ ≈ 95° between the surface and pit walls).

The morphology of etch pits growing on the (210) face is similar to that of etch pits developed on (100) and (010) faces. However, their lenticular shape is not as well defined and approaches that of a half-moon: indeed, one edge of the pit subparallel to the c-axis is straighter than the other, which is more curvilinear. In addition, whereas the deepest point of pits modeled on (100) and (010) faces is located at the center of the pit, the pits developed on the (210) face are asymmetric (see cross-section 7) and the deepest point (3-4 nm) is located at the center of the pit following the c-axis, but shifted on one edge of the minor-axis ([-210]).

Finally, although no etch pit could be evidenced experimentally on (001) faces, the model predicts that small (≈ 10 nm) and deep (> 40 nm) subcircular etch pits may grow on the (001) face if screw dislocations parallel to the c-axis are implemented. However, unlike the other faces, etch pits developed on the (001) face reach a stable radius or tend to vanish, as a result of the increasing roughness of the surrounding surface, possibly explaining why such etch pits remained out of reach of our analytical characterizations.

Discussion

Factors influencing the face-specific dissolution rates of enstatite and comparison with previous studies

As reported above, the face-specific dissolution rates of enstatite range between 2.6 and 71.8 nmol/m²/s at pH 0 and T = 90 °C. These values are between one and two orders of magnitude lower than expected using the dissolution rate law derived from powder dissolution experiments conducted by Oelkers [START_REF] Oelkers | An Experimental Study of Enstatite Dissolution Rates as a Function of Ph, Temperature, and Aqueous Mg and Si Concentration, and the Mechanism of Pyroxene/Pyroxenoid Dissolution[END_REF] (between 730 and 973 nmol/m²/s), which represents the largest enstatite dissolution kinetics dataset that we are aware of. Such discrepancies may either be ascribed to the contribution of edges and corners, whose reactivity is much greater than that of faces [START_REF] Noiriel | Pore-Scale Geochemical Reactivity Associated with Co2 Storage: New Frontiers at the Fluid-Solid Interface[END_REF][START_REF] Briese | The Effect of Crystal Size Variation on the Rate of Dissolution -a Kinetic Monte Carlo Study[END_REF] and which are overrepresented when using powders. It may also be due to the very low pH investigated in the present study, which may result in a saturation of proton surface sites [START_REF] Blum | Role of Surface Speciation in the Low-Temperature Dissolution of Minerals[END_REF] , possibly leading to a weaker pH-dependence of silicate dissolution rate than expected in the extreme acidic pH domain. In any case, the present study illustrates the unsuitability of isotropic dissolution rate laws derived from powder dissolution experiments to predict the amplitude and absolute values of face-specific dissolution rates, as pointed out in previous studies [START_REF] Fischer | How Predictable Are Dissolution Rates of Crystalline Material?[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] Fischer | Beyond the Conventional Understanding of Water-Rock Reactivity[END_REF] .

The results show that the dissolution rate of the ( 001) face is about one order of magnitude greater than that of the other investigated (hk0) faces. Among these latter faces, the dissolution rates of the ( 100) and ( 010) faces are similar within uncertainties, and are about three times slower than that of the (210) face. If one considers the crystallographic similarities existing between the (110) face of diopside and the (210) face of enstatite (both faces correspond to the cleavage plane of these minerals, and intersect the highest density of Mg-O-Si bonds), these results are in excellent agreement with those obtained on diopside (Ca0.5Mg0.5SiO3) [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] , which also belongs to the pyroxene group (Fig. 8B). Note however that the anisotropic reactivity of enstatite seems to be less pronounced than that of diopside. This is likely due to the presence of Ca ions in the diopside structure, which offers more combinations for the bond-breaking probabilities (and therefore, detachment rates) of Mg and Si atoms. Indeed, the hydrolysis rates of Ca-O-Si and Ca-O-Ca bonds are faster than that of Mg-O-Si and Mg-O-Mg bonds, as demonstrated experimentally [START_REF] Casey | Control of Dissolution Rates of Orthosilicate Minerals by Divalent Metal Oxygen Bonds[END_REF] and indirectly confirmed theoretically following ab initio calculations [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] . diopside [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] face-specific dissolution rates. For both graphs, magenta, red, green and blue colors stand for the results obtained for (100), ( 010), ( 001) and ( 210) faces, respectively.

The factors that are most frequently invoked to account for the anisotropy of silicate dissolution include (i) the crystallographic bonding structure and in particular, the connectedness of silicate tetrahedrons [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] Daval | Dynamics of Altered Surface Layer Formation on Dissolving Silicates[END_REF][START_REF] Turpault | The Dissolution of Biotite Single-Crystals in Dilute Hno3 at 24-Degrees-C -Evidence of an Anisotropic Corrosion Process of Micas in Acidic Solutions[END_REF] (ii) the dislocation density and orientation [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] 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] Daval | Dynamics of Altered Surface Layer Formation on Dissolving Silicates[END_REF][START_REF] Burch | Free Energydependence of Albite Dissolution Kinetics at 80°C and Ph 8.8[END_REF][START_REF] Lasaga | Variation of Crystal Dissolution Rate Based on a Dissolution Stepwave Model[END_REF] and (iii) the formation of ASSLs, whose transport properties may differ as a function of the crystallographic orientation [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] Wild | Ph-Dependent Control of Feldspar Dissolution Rate by Altered Surface Layers[END_REF][START_REF] Daval | Dynamics of Altered Surface Layer Formation on Dissolving Silicates[END_REF] .

Considering the conditions that are thought be required for the development of passivating ASSLs (i.e., oxic conditions and mildly acidic pHs, see [START_REF] Wild | Ph-Dependent Control of Feldspar Dissolution Rate by Altered Surface Layers[END_REF][START_REF] Daval | Dynamics of Altered Surface Layer Formation on Dissolving Silicates[END_REF][START_REF] Saldi | Mineralogical Evolution of Fe-Si-Rich Layers at the Olivine-Water Interface During Carbonation Reactions[END_REF][START_REF] Wild | Early Stages of Bacterial Community Adaptation to Silicate Aging[END_REF] ), the latter factor (iii) can be safely ruled out in the present study.

Unravelling the relative contributions of factors (i) and (ii) to the observed anisotropic reactivity may benefit from the comparison between the outputs of the simulations and the experimental results. When the simulations are run without implementing any dislocation in the enstatite structure, the reactivity of the (001) face is strikingly greater than that of the other (hk0) faces, as observed experimentally (Fig. 6). Therefore, this suggests that the connectedness of silicate tetrahedrons dictates the anisotropic reactivity of enstatite at first order: as previously pointed out by Daval et al. [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] , the release of Si atoms from the (001) face to the solution requires the cleavage of one Si-O-Si bond only, while it would require the (less likely) simultaneous cleavage of two Si-O-Si bonds on the (hk0) faces, which are parallel to the silicate tetrahedron chains.

The connectedness of SiO4 tetrahedrons cannot be invoked to distinguish between the reactivity of the three other investigated (hk0) faces, because it is the same for all of them.

However, because the outputs of the simulations show that the dissolution rates of these three faces are different even without implementing any dislocation, this suggests that the crystallographic bonding structure also exerts a strong control on the face-specific reactivity.

As it can be seen in Fig. 1, the Mg and Si atoms are organized following planes that are parallel to the (100) face. If one supposes that the detachment rate of Si atoms connected to two other Si atoms is lower than that of Mg atoms (in agreement with the respective activation energies for the hydrolysis of Si-O-Si, Si-O-Mg and Mg-O-Mg bonds, see section 2.3.1), it can be assumed that the progression of the dissolution in the [100] direction will be hampered by the sluggishness of the dissolution of these Si-rich planes. Conversely, the detachment of accessible Mg atoms located below the first atomic plane of the (010) face may contribute to destroying the first coordination spheres of the topmost Si atoms, which eventually increase their probability of detachment. In turn, this explanation would account for the fact that the slowest dissolving face is (100), and the dissolution rates inferred from the simulations therefore observe the following trend: 𝑟 Finally, these simple crystallographic considerations do not account for the observation that the dissolution rate of the (210) face is greater than that of the (010) face.

Interestingly, running the simulations with the appropriate dislocation density measured experimentally allows to obtain a global trend for the face-specific reactivity of enstatite which matches exactly to that measured experimentally, with 𝑟 , ( )

< 𝑟 , ( ) < 𝑟 , ( ) << 𝑟 , ( )
(see details in section 4.2). These results confirm the critical contribution of dislocations to mineral dissolution rates (e.g., [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 | Does Crystallographic Anisotropy Prevent the Conventional Treatment of Aqueous Mineral Reactivity? A Case Study Based on K-Feldspar Dissolution Kinetics[END_REF][START_REF] Lasaga | Variation of Crystal Dissolution Rate Based on a Dissolution Stepwave Model[END_REF][START_REF] Brantley | Dissolution at Dislocation Etch Pits in Quartz[END_REF][START_REF] Bouissonnié | From Mixed Flow Reactor to Column Experiments and Modeling: Upscaling of Calcite Dissolution Rate[END_REF] ).

To sum up, both the crystallographic bonding structure and the presence of crystal defects (and in particular, dislocations) contribute to the observed anisotropic reactivity of enstatite, in agreement with previous studies dedicated to the face-specific dissolution rates of minerals.

Assessing the agreement between measured and simulated rate data

A total of 45 sets of activation energies for the thermally activated probabilities of breaking Mg-O-Mg, Mg-O-Si and Si-O-Si bonds were tested to cover the possible range of hydrolysis described in section 2.3.2. Three different complementary metrics were used to evaluate quantitatively the agreement between simulation outputs and experimental results: classical metrics such as the relative root mean square error (RRMS) and the relative mean error (RME), and an additional original metrics (SaI) based on a linear combination of the slope and the intercept of the line of best fit, i.e., the line obtained through linear regression between modeled dissolution rates vs. measured dissolution rates. These metrics are defined by: Following this formulation, the smallest the BA, the better is the agreement between experimental dissolution rates and simulated dissolution rates.

𝑅𝑅𝑀𝑆 = 1 𝑁 𝑟 ( ) -𝑟 ( ) 𝑟 ( ) (8)
Based on the BA metric, only a narrow range of probabilities provided a satisfactory agreement between simulation outputs and experimental data (Fig. 9) and one set of probabilities ([PMgOMg = 0.99, PSiOMg = 0.4; PSiOSi = 0.0146]) yielded the lowest value of BA (BA = 0.05). This set has been selected to further discuss the outputs of the simulations below, although it is noteworthy that a couple of sets of probabilities also returns similarly low BA values (e.g., the two sets of probabilities [PMgOMg = 0.99, PSiOMg = 0.25; PSiOSi = 0.0177] and [PMgOMg = 0.99, PSiOMg = 0.3; PSiOSi = 0.0212] correspond to BA values of 0.07 and 0.09, face remains unaffected, probably because of the low dislocation density observed on this face. Finally, the outcropping of dislocations at the surface of the (001) face has not been considered, since it was not possible to evidence etch pits on the corresponding samples. Of note, the implementation of dislocations also has an impact on the agreement between simulation outputs and experimental results. Whereas the implementation of dislocations tend to minimize the RRMS and RME values (RRMS = 22% and RME = 13% with dislocations vs. RRMS = 30% and RME = 22% without dislocations), this also results in a modest increase of the SaI (SaI = 0.11 vs. SaI = 0.05), implying that the dissolution rate anisotropy becomes slightly underestimated. Overall, it can be concluded that adding dislocations to the model slightly improves the general agreement between simulation outputs and experimental results in terms of RRMS and RME but slightly degrades the account for the anisotropy of dissolution rates (increase of the SaI value). However, the BA index remains unaffected (BA = 0.06 in both cases), which indicates that for this set of probabilities, the simulations conducted with dislocations are statistically equivalent to those run without dislocations.

To sum up, the model satisfactorily accounts for enstatite face-specific dissolution rates measured experimentally, using bond-breaking probabilities consistent with what can be inferred from literature. Moreover, adding dislocations to the simulations in proportions that are consistent with the analytical characterizations allows one to reproduce the order of the face-specific dissolution rates observed experimentally.

Comparison between observed and modeled etch pits and related implications for enstatite dissolution mechanisms

The improvement of the agreement between modeled and experimental dissolution rates is not the only interesting output resulting from the implementation of dislocations in the simulated crystals. It also offers the possibility to compare the morphology of simulated etch pits developed at the dislocation outcrop to their counterparts revealed by electron microscopy on the different faces of enstatite (Fig. 4). This represents an aspect of great interest, as it is regularly searched for as a tracer of weathering conditions [START_REF] Bouissonnié | From Mixed Flow Reactor to Column Experiments and Modeling: Upscaling of Calcite Dissolution Rate[END_REF][START_REF] Velbel | Dissolution of Olivine During Natural Weathering[END_REF][START_REF] Phillips-Lander | Can We Use Pyroxene Weathering Textures to Interpret Aqueous Alteration Conditions? Yes and No[END_REF] .

Interestingly, both the lenticular shape, the orientation and the symmetry of the pits generated numerically on (100) and (010) faces using the best set of probabilities (see previous section) are consistent with those observed experimentally. The simulated and experimentally observed pits are very elongated along the c-axis and their minor-axes follow either the [010] (for pits developed on (100) faces) or the [100] (for pits developed on the (010) faces) direction, respectively. These directions also represent symmetry axes for the pits, resulting in the existence of a symmetry point at their center. The TEM characterizations performed on the FIB thin section realized on the altered (010) and (100) faces further indicate that the angle formed by the pit walls with the surface are in good agreement with the outputs of the simulations.

A noticeable difference revealed by the model between the pits developed on the (100) and (010) faces consists in the shape of their bottoms following the c-direction (Fig. 7). These differences may be ascribed to the crystallographic structure of enstatite. Indeed, while the deepening of pits formed on the (100) face requires to sequentially intersect Si-rich and Mgrich planes, the deepening of pits formed on the (010) face entails that the dissolution takes place in a single plane, resulting in a faster opening of the corresponding pits. This may explain why modeled etch pits formed on the (010) face are flat-bottomed in the c-direction, as opposed to the pits formed on the (100) face, since the propagation of the dissolution following the c-axis will be easier in the case of the (010) face than the one of the (100) face. Supporting this assertion, the sequential removal of Mg atoms from Mg-planes, which are readily accessible for pits formed on (010) faces, will result in lowering the coordination of Si atoms in the Si-planes, thereby accelerating the global propagation of the pit along the c-axis. Such a mechanism is not possible for pits developed on the (100) face, where the widening of a pit involves dissolving Mg-and Si-rich planes sequentially. Finally, it is noteworthy that these differences in terms of pit morphology, which follow from the anisotropic structure of enstatite, also result in a difference of reactivity of the pits themselves, with pits developed on the (010) faces being more effective to enhance enstatite dissolution than those developed on the (100) faces. Finally, TEM observations reveal that etch pits nucleated on the (100) face are flat-bottomed following their minor axis ([010]), a feature that is not accurately reproduced by the model. Possible explanations for this observation include the fact that etch pits revealed experimentally may have investigated at a greater reaction progress than the modeled ones, which result in a widening of the pit in the b-direction.

Regarding the modeled pits formed on the (001) face, their main features consist in their narrow subcircular shape (5 nm) and important depth (> 20 nm). As reported above, the simulations reveal that the diameter of the simulated pits never exceeds 5 nm because of the gradual and global roughening of the surface, which may explain why no etch pits were observed by VSI, SEM or TEM on the (001) samples altered experimentally. The absence of a continued lateral growth of the pits in the (001) plane also suggests that dislocation etch pits form on the (001) face (if any) have a negligible contribution to the dissolution rate of this face.

The pits that have been simulated on the (210) face are different from the ones simulated on the (100) and (010) faces. Their shape approaches that of a crescent, with one straight edge along the c-axis on one side and a curvilinear edge on the other side. This shape is consistent with that revealed experimentally for (210) samples (see Fig. 4), and very similar to the one observed for pits developed on the (110) and (11 0) faces of diopside by Daval et al. [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] . The topography profile of the pits following the c-axis is similar to the one observed on the (100) face, with a deepest point (≈ 2.0 -2.5 nm) located at the center of the pit. The main difference between the pits simulated on this face and those simulated on the (100) and (010) faces is revealed when considering the topography profile of the pit following the minor-axis, which is strongly asymmetric. Whereas the angle between the (210) surface and the pit wall at the straight edge is very steep (θ ≈ 95°), the other one is much less sharp (θ ≈ 113°). Although it was not possible to get a FIB foil across a single etch pit nucleated on the (210) face, it is worth mentioning that this asymmetry was also observed for pits developed on the (110) face of diopside (see atomic force microscopy (AFM) measurements in [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] ). This asymmetry is due to the fact that the direction of the dislocation is not orthogonal to the (210) face. An interesting consequence of this specific geometry is that the entire simulated pit is not visible in Fig. 7, since part of the pit is located below the surface of the face. Such features may correspond to the so-called "etch tubes" [START_REF] Brantley | Dissolution at Dislocation Etch Pits in Quartz[END_REF] or "etch tunnels" 81 sometimes reported in the literature, which cannot be imaged using classical surface sensitive techniques such as VSI or AFM.

Finally, the pit geometry revealed by the numerical simulations also has broader consequences for the practical measurement of crystal dissolution rates. As can be seen in Fig. 6D, the dissolution rate of the (210) face based on Mg release exceeds that calculated using the overall surface topography. Such a "hidden" contribution to crystal reactivity may help explain why face-specific mineral dissolution rates determined from solution chemistry measurements sometimes overestimate dissolution rates estimated using surface sensitive techniques, such as VSI [START_REF] Perez | Comparing the Reactivity of Glasses with Their Crystalline Equivalents: The Case Study of Plagioclase Feldspar[END_REF] . Furthermore, considering that the dislocation network in crystals can be connected to the surface as a result of pit opening, intra-crystalline dissolution may occur. A detailed understanding of the solution chemistry within single pits then becomes crucial, since concentration gradients can affect the reactivity of the etch pit [START_REF] Bouissonnié | From Mixed Flow Reactor to Column Experiments and Modeling: Upscaling of Calcite Dissolution Rate[END_REF] .

To sum up, the etch pits modeled on (100), (010) and (210) are very similar to the ones observed on enstatite samples. The orientation, shape, and symmetry of the pits are correctly predicted by the model, and result from the crystallographic structure of enstatite and orientation of dislocation lines.

Relation between the dissolution regime and the formation of Si-rich surface layers

As shown in Fig. 6, the transient regime, which is characterized by a nonstoichiometric release of Mg and Si atoms and a variable dissolution rate, spans over 7,500 to 8,250 iterations (depending on the considered face) corresponding to 2.7 ± 1.0 to 2.9 ± 1.0 days, when the best set of probabilities is selected (section 2.3.6). This result suggests that the dissolution rates estimated experimentally for the (100), ( 010) and (210) faces were actually measured during the steady-state regime, where the dissolution is stoichiometric and the rate is constant, and may explain why the very first VSI measurements of surface retreats conducted on the (100) and (010) faces fall slightly below the regression trend line.

Conversely, part of the measurements carried out on the (001) faces might have been performed during the transient regime, although this assertion is not supported experimentally, since the measured surface retreats observe a linear temporal trend. This result is not totally surprising, considering the expected modest evolution of the dissolution rate with time inferred from the numerical outputs of the simulation for this face (Fig. 6), which would imply that any deviation from the linear trend would fall within experimental uncertainties.

From a physical standpoint, the model outputs suggest that the transient regime corresponds to the period where the number of surface atoms increases. The number of atoms at the surface increases by a 2.5-, 2.8-, 3.5-, and 80 fold-factor for the (100), (210), (010) and (001) faces, respectively. After this transient regime, the number of each surface cation reaches a constant value, resulting in a steady-state dissolution rate. Prior to reaching the steady-state regime, the dissolution rates based on Mg or Si release have opposite trends: while the Mg release rate decreases, enstatite dissolution rate based on Si release increases with time.

The attainment of a steady-state regime thereby appears intimately connected to the dynamics of Si release, which can be explained by the large differences existing between the detachment rates of Mg and Si atoms: the preferential detachment of Mg atoms will result in a depletion of the first coordination sphere of Si atoms, in turn increasing their release rates.

The simultaneous attainment of a constant dissolution rate and a constant number of Si surface atoms for the 4 faces further suggests that the mean coordination of Si surface atoms also reached a steady-state after the transient regime.

This divergence in the dynamics of Si and Mg release has another important consequence: during the transient regime, the simulations predict that the preferential release of Mg atoms results in the formation of an Mg-depleted / Si-rich surface layer. According to the simulation outputs, the thickness of the leached layer is 3.8, 6.0, 6.9 and 107 Å for (100), (210), (010), and (001) faces, respectively. Interestingly, the model therefore predicts that thick surface layers should develop on the (001) face, whereas they should be much thinner on (hk0) faces, which is consistent with the TEM results described in section 3.3, and the results reported by Daval et al. [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] for diopside dissolution. However, the model systematically underestimates the actual thickness of the layer, which may be ascribed to the fact that it neither accounts for (i) the spontaneous reorganization of the Mg-depleted layer, such as the condensation of silanol groups, which may be associated to a global variation of the volume of the surface layers, as suggested in previous studies (e.g. [START_REF] Weissbart | Wollastonite: Incongruent Dissolution and Leached Layer Formation[END_REF][START_REF] Schott | Formation, Growth and Transformation of Leached Layers During Silicate Minerals Dissolution: The Example of Wollastonite[END_REF] ); nor for (ii) the backward reactions of Si attachment to the surface layer, following the alternative interfacial dissolution-reprecipitation mechanism [START_REF] Daval | Mechanism of Wollastonite Carbonation Deduced from Micro-to Nanometer Length Scale Observations[END_REF][START_REF] Ruiz-Agudo | Mechanism of Leached Layer Formation During Chemical Weathering of Silicate Minerals[END_REF][START_REF] Hellmann | Unifying Natural and Laboratory Chemical Weathering with Interfacial Dissolution-Reprecipitation: A Study Based on the Nanometer-Scale Chemistry of Fluid-Silicate Interfaces[END_REF][START_REF] Hellmann | An Eftem/Hrtem High-Resolution Study of the near Surface of Labradorite Feldspar Altered at Acid Ph: Evidence for Interfacial Dissolution-Reprecipitation[END_REF] . Here, and as emphasized in previous kMC studies [START_REF] Zhang | Theoretical Approach to Evaluating Plagioclase Dissolution Mechanisms[END_REF] , the outputs of the model indirectly support that the preferential alkali / alkaline earth depletion during silicate dissolution is limited to the topmost layers of the minerals, and that ASSL formation should, at least partly, result from Si redeposition [START_REF] Daval | Dynamics of Altered Surface Layer Formation on Dissolving Silicates[END_REF] .

To sum up, the comparison between the experimental data and the simulation outputs suggests that the duration of the transient period is reasonably well reproduced by the model.

In addition, the model predicts the formation of Mg-depleted / Si-rich layers, although their actual thickness remains underestimated. This latter result indirectly reveals that certain mechanisms that are not accounted for by the model (e.g. Si attachment reactions; condensation of silanol groups) actually occur over the course of enstatite dissolution.

Conclusions and perspectives

In this study, we focused on the agreement between the outputs of stochastic atomicscale dissolution simulations and the results obtained on face-specific enstatite dissolution experiments carried out at 90 °C and pH 0 (rate data and structural observations of the faces).

Nearly 200 simulations were conducted to demonstrate that the anisotropy of enstatite dissolution rate can be satisfactorily reproduced using bond-breaking probability ratios (and therefore, activation energy differences) that are consistent with the existing literature regarding the hydrolysis of Mg-O-Mg, Mg-O-Si and Si-O-Si bonds. On the same footing, the model also correctly predicts various face-specific features, such as the orientation, symmetry, and shape of dislocation etch pits. From a physical chemistry standpoint, this agreement gives rise to the following two perspectives with respect to silicate dissolution and the silicate-water interface:

(i) Nowadays, ab initio calculations are generally used to feed kMC simulations (e.g. [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 Approach to Study Carbonate Dissolution[END_REF][START_REF] Kurganskaya | Kinetic Monte Carlo Simulations of Silicate Dissolution: Model Complexity and Parametrization[END_REF] ). Conversely, our study suggests that kMC simulations have the ability to be a powerful tool to constrain a range of activation energies for the hydrolysis of bonds that are not

Fig. 1 .

 1 Fig. 1. Projection of enstatite structure following [100] (A), [001] (B), [010] (C) zone axes and a rotation view of the enstatite volume (D) to highlight the chains of silicate tetrahedrons. The Mg, Si and O atoms are represented by red, green and blue spheres respectively. Si tetrahedrons are represented in light green.

) with 0 <

 0 PM-O-Mg, PM-O-Si < 1. According to ab initio and experimental studies, the most probable event is the detachment of an Mg atom bonded to the surface via a single Mg-O-Mg bond. The corresponding probability was then set to 0.99 and can be associated, by identification, to the Mg-O-Mg activation energy term in Eq. (4). The two other elementary probabilities (PMg-O-Si and PSi-O-Si) were calculated in order to observe the theoretical difference of activation energies between the hydrolysis of Mg-O-Mg and that of Mg-O-Si or Si-O-Si, respectively. Considering the ab initio data from 36, 38 for the activation energy of Si-O-Si (∆𝐸 = 69 kJ/mol) and Mg-O-Si (∆𝐸 = 59 kJ/mol) hydrolyses and the experimental data from 55 for Mg-O-Mg hydrolysis (∆𝐸 = 60 ± 12 kJ/mol), the range of difference in activation energy that was explored is [0 ; 10 kJ/mol] for the {Si-O-Si ; Mg-O-Si} pair and [0 ; 4 kJ/mol] for the {Mg-O-Mg; Mg-O-Si} pair.

Fig. 2 .

 2 Fig. 2. A. Face-specific temporal evolution of the mean surface retreat measured with VSI and associated linear regressions. B. Measured face-specific dissolution rates. For both graphs, magenta, red, green and blue colors stand for the results obtained for (100), (010), (001) and (210) faces, respectively.

Fig. 3 .

 3 Fig. 3. A. Surface topography of the (001) surface dissolved for ~ 8h measured using VSI. The red spot represents the mask (non-reacted) area while the blue/green area represents the reacted surface. B. Topography of the surface (straight line in A) after dissolution. The upper part corresponds to mask area while the lower part corresponds to the reacted area.

  dislocations has been calculated on each face by counting the number of etch pits observed by SEM on surface areas of ~ 300 µm². Of note, polishing scratches also represented an important source of etch pits. Because they can contribute to the global dissolution rate by exposing step sites, they have been taken into account in the calculation of the dislocation density.

Fig. 4 .

 4 Fig. 4. SEM images of the surface of (001) (A), (010) (B) and (210) faces reacted at pH 0 for 7 to 21 days. The images reveal the nucleation of lenticular etch pits elongated following the [001] direction. Note the asymmetric shape of the pits nucleated on the (210) face, as opposed to those developed on the other two faces.

Fig. 5 .

 5 Fig. 5. TEM characterizations of FIB thin sections prepared on (100) (A), (010) (B), and (001) (C and D) reacted surfaces. HRTEM images in (B) and (C) reveal the formation of a nm-thick amorphous layer formed on the surface, with a sharp crystallographic boundary between enstatite and the surface layer. For the FIB thin section realized on (001) surface, STEM-EDX profiles measured parallel to the [001] direction reveal that the surface layer is Si-rich and Mg-depleted (D).

Fig. 6 .

 6 Fig. 6. Outputs of the simulations conducted with the bond-breaking probabilities corresponding to the best agreement between measured and modeled dissolution rates (see text for details) for (100) (A), (010) (B), (001) (C) and (210) (D) faces. The evolution of the modeled surface retreat calculated following three different methods (see text) is depicted with black lines. The color lines stand for the evolution of the number of surface

3 . 2 )

 32 and the dislocation density determined experimentally for each face, the modeled face-specific dissolution rates 𝑟 ( ) observe the following trend: relative differences between the modeled dissolution rates for the various (hk0) faces depend on the set of probabilities used to run the simulations. The best agreement between modeled and measured steady-state dissolution rates is obtained for PMg-O-Mg = 0.99; PMg-O-Si = 0.4 and PSi-O-Si = 0.0146 (i.e.,

Fig. 7 .

 7 Fig. 7. Evolution of modeled etch pits and associated cross-sections following the directions depicted by the doted lines. For all pictures, magenta, red, green and blue colors stand for the results obtained for (100), (010), (001) and (210) faces.

Fig. 8 .

 8 Fig. 8. A. Best agreement between modeled and measured face-specific enstatite dissolution rates. The triangles and circles stand for simulations conducted with and without dislocations, respectively. Note that the modeled dissolution rate of the (001) face is also affected because the time for one iteration changes as a function of the mean ratio between experimental and modeled dissolution rates (Eq. 7). B. Comparison between enstatite and

(

  face representing an intermediate configuration between the (100) and (010) faces.

9 )

 9 the dissolution rate for face (hkl) estimated experimentally (resp. by the model), NE is the number of experimental values (4 in our case), a is the slope of the best fit line, b is the intercept of the best fit line and < 𝑟 ( ) > the average of all experimental dissolution rates. The RRMS index indicates how spread out are the differences between the modelled values and the observed values or how close are, in average, these differences around the line of best fit. The RME index allows detecting possible bias, i.e., if RME is significantly greater than 0 (resp. lower than 0), the model overestimates (resp. underestimates) the observed values. The information provided by the best fit line is gathered in the SaI index. Because of differences in units (a is dimensionless, b has the dimension of the dissolution rate), the second term of SaI is divided by < 𝑟 ( ) >. We used this reference for model performance comparisons because it is a common reference for all data. RRMS-RME indices and the SaI index provide different kinds of information. Any set of probabilities (or their corresponding activation energies) that minimizes the RRMS and RME indices does not necessarily minimizes the SaI index as well (e.g., the set of probabilities that minimizes the RRMS and RME indices is [PMgOMg = 0.99, PSiOMg = 0.3 and PSiOSi = 0.0109] (RRMS = 22% and RME = 17%), whereas it overestimates the dissolution rate anisotropy (SaI = 0.52). On the other hand, the set that allows for the best SaI is [PMgOMg = 0.99, PSiOMg = 0.40 and PSiOSi = 0.0146] (SaI = 0.06), corresponding to RRMS and RME of 30% and 22%, respectively). In order to compare the 45 sets of probabilities, we finally defined a composite metric that we called best agreement (BA) for each set i by the following formulation: (𝑅𝑅𝑀𝑆) max(𝑅𝑅𝑀𝑆) -min(𝑅𝑅𝑀𝑆) + 𝑅𝑀𝐸 -min(𝑅𝑀𝐸) max(𝑅𝑀𝐸) -min(𝑅𝑀𝐸) + 𝑆𝑎𝐼 -min (𝑆𝑎𝐼) max(𝑆𝑎𝐼) -min(𝑆𝑎𝐼) (In this formulation, the three metrics are scaled by their max-min values, and we balanced the classical metrics RRMS and RME with the new metric SaI. Other formulations could have been used for BA but we privileged the balanced linear formulation that is more robust than other formulations based on single values like min-max values. Due to the small number of experimental data, more sophisticated composite metrics would have been questionable.

Table 1 .

 1 Atomic coordinates in enstatite cell[START_REF]A Compressional Study of Mgsio3 Orthoenstatite up to 8.5 Gpa[END_REF] used in Eq. (3) to form the entire cell.

	Mg and Si atoms were connected to O atoms according to the following rules: (i) Mg
	atoms are located in octahedral sites and Si atoms in tetrahedral sites (i.e., Mg atoms are
	connected to 6 O atoms and Si atoms are connected to 4 O atoms) and (ii) the closest O atoms
	were considered to create the coordination spheres surrounding each cation. The resulting
	modeled surfaces ((100), (010) and (001)) are shown in Fig. 1.

Acknowledgements

We warmly thank G. Morvan and R. Boutin (LHyGeS, Strasbourg) for their help with EBSD and ICP-AES measurements. M. Cabié (CP2M, Aix-Marseille) and T. Dintzer (ICPEES, Strasbourg) are acknowledged for their help with the preparation of the FIB thin sections and

Abstract graphic

respectively), pointing towards the existence of a narrow range of possible PMgOSi and PSiOSi (and therefore, of ∆𝐸 -∆𝐸 and ∆𝐸 -∆𝐸 ) values. available in the current literature, and/or to validate the values resulting from ab initio calculations. For instance, the value of ∆EMg-O-Mg we used was selected from an experimental work [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] , because to the best of our knowledge, no study provided a theoretical value for this parameter. If the agreement between modeled and simulated rates and etch pit geometry is taken for granted, it could be immediately transposed to the diopside structure, for which farfrom-equilibrium rate data already exist [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] , making possible the determination of a possible range for parameters such as ∆EMg-O-Ca, ∆ECa-O-Ca or ∆ECa-O-Si (note that ab initio data already exist for this latter parameter, see [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] ). This will be the subject of an upcoming paper.

(ii) On the same footing, a similar approach could be usefully transposed to other silicate structures for which rate data already exist (see [START_REF] Casey | Control of Dissolution Rates of Orthosilicate Minerals by Divalent Metal Oxygen Bonds[END_REF] , which report far-from-equilibrium rate data for a series of eight endmember orthosilicate minerals (M2SiO4) which varied only by the nature of the divalent cation M they contain).

Conversely, we showed that the model underestimates the thickness of ASSLs formed on the various faces of enstatite. In itself, this result already provides insights into the dynamics of the dissolution process, as previously pointed out by Zhang and Lüttge [START_REF] Zhang | Theoretical Approach to Evaluating Plagioclase Dissolution Mechanisms[END_REF] : thick (i.e., several 10s of nm) ASSLs cannot be formed by a pure leaching mechanism only. This indicates where future kMC efforts should be directed, as the formation, growth and maturation of such layers are regularly suggested to influence the reactivity of the underlying substrate, at least in SiO2(aq) rich solutions and mildly acidic pHs [START_REF] Wild | Ph-Dependent Control of Feldspar Dissolution Rate by Altered Surface Layers[END_REF][START_REF] Daval | Time-Dependent Feldspar Dissolution Rates Resulting from Surface Passivation: Experimental Evidence and Geochemical Implications[END_REF][START_REF] Saldi | Mineralogical Evolution of Fe-Si-Rich Layers at the Olivine-Water Interface During Carbonation Reactions[END_REF][START_REF] Wild | Early Stages of Bacterial Community Adaptation to Silicate Aging[END_REF] .

for providing high resolution SEM images of etch pits. A.B. thanks the University of Strasbourg for funding his PhD grant.