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

Damien Daval

Philippe Ackerer

Dissolution Anisotropy of Pyroxenes: Role of Edges and Corners Inferred from Stochastic Simulations of Enstatite Dissolution

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Dissolution kinetics is a key parameter to understand the Earth's surface evolution and is crucial for several industrial and environmental problems (cement dissolution, nutrient depletion in soils, etc.). Over the last decades, numerous studies have helped understand the processes that occur at the mineralwater interface, as well as developing different methods to measure their kinetics for a wide variety of minerals. [1][2][3][4][5][6][7][8][9] While most of these studies are based on powder dissolution experiments, 2,10-12 fewer studies succeeded to measure dissolution rates on different faces of a given mineral. 3,4,6,7,9,[START_REF] Salehikhoo | Magnesite Dissolution Rates 1082 at Different Spatial Scales: The Role of Mineral Spatial Distribution and 1083 Flow Velocity[END_REF][START_REF] Perez | Comparing the Reactivity of Glasses with Their Crystalline Equivalents: The Case Study of Plagioclase Feldspar[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] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] This latter kind of studies has shown the importance of the crystallographic structure (i.e., anisotropy) on mineral dissolution kinetics. As a relevant example, pyroxenes, which are strongly anisotropic due to their silica tetrahedron chains, exhibit a strong dissolution anisotropy consistent with the periodic bond chain theory. [START_REF] Hartman | On the relations between structure and morphology of crystals[END_REF] However, the comparison between the dissolution rates measured on single faces and those measured on powder highlights a significant difference (powder dissolution rates are generally higher than face-specific dissolution rates [START_REF] Arvidson | Variation in Calcite Dissolution Rates[END_REF] ). This adds complexity to our experimental setups have emerged, using analytical techniques such as X-ray tomography [START_REF] Noiriel | Direct Determination of Dissolution Rates at Crystal Surfaces Using 3D X-Ray Microtomography[END_REF][START_REF] Noiriel | Examination of Crystal Dissolution in 3D: A Way to Reconcile Dissolution Rates in the Laboratory? Geochim[END_REF][START_REF] Kahl | Crystal Surface Reactivity Analysis Using a Combined Approach of X-Ray Micro-Computed Tomography and Vertical Scanning Interferometry[END_REF] or vertical scanning interferometry. [START_REF] Luettge | An Interferometric Study of the Dissolution Kinetics of Anorthite: The Role of Reactive Surface Area[END_REF] However, despite their increasing efficiency to study the processes that occur at finer scales (micrometer to nanometer scale), they often need to be combined with numerical studies to unravel what happens at the solid-fluid interface.

Numerical experiments have become a powerful tool to improve our understanding of processes occurring at the watermineral interface, ranging from molecular dynamics to reactivetransport model used at pore scale. [START_REF] Bourg | Molecular Dynamics Simulations of Water Structure and Diffusion in Silica Nanopores[END_REF][START_REF] Molins | Pore-Scale Controls on Calcite Dissolution Rates from Flow-through Laboratory and Numerical Experiments[END_REF][START_REF] Steefel | Reactive Transport Codes for Subsurface Environmental Simulation[END_REF] The last decades have also seen the development of stochastic dissolution models, such as kinetics Monte Carlo models or Voronoi ̈method. [START_REF] Salehikhoo | Magnesite Dissolution Rates 1082 at Different Spatial Scales: The Role of Mineral Spatial Distribution and 1083 Flow Velocity[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] Luẗtge | A Stochastic Treatment of Crystal Dissolution Kinetics[END_REF][START_REF] Kurganskaya | Kinetic Monte Carlo Simulations of Silicate Dissolution: Model Complexity and Parametrization[END_REF][START_REF] Kurganskaya | Kinetic Monte Carlo Approach To Study Carbonate Dissolution[END_REF][START_REF] Rohlfs | Crystal Dissolution Kinetics Studied by a Combination of Monte Carlo and Voronoi Methods[END_REF][START_REF] Kurganskaya | Atomistic to Meso-Scale Modeling of Mineral Dissolution: Methods, Challenges and Prospects[END_REF][START_REF] Zhang | Si Order in Albite and Its Effect on Albite Dissolution Processes: A Monte Carlo Study[END_REF][36][START_REF] De Assis | Dissolution of Minerals with Rough Surfaces[END_REF][START_REF] Lasaga | Surface Chemistry, Etch Pits and Mineral-Water Reactions[END_REF] These models, based on a relation between the activation energies of individual processes (bond hydrolysis, for example) and their occurrence probability, allow for fast computation and the simulation of larger objects, which classical molecular dynamics is not able to deal with. These models, first applied on Kossel crystals, have brought or validated many important results. As an example, numerous studies applied on Kossel crystals validate the difference of reactivity between the crystal's edges, corners, and surfaces. [START_REF] Briese | The Effect of Crystal Size Variation on the Rate of DissolutionA Kinetic Monte Carlo Study[END_REF] The work of Zhang and Luttge [START_REF] Zhang | Si Order in Albite and Its Effect on Albite Dissolution Processes: A Monte Carlo Study[END_REF] also mentioned these differences for feldspar dissolution. However, providing a generalization to other minerals, which is an intrinsic issue belonging to the field of physical chemistry, is important and may present a great interest to link laboratory and field observations, as mentioned by Noiriel et al. [START_REF] Noiriel | Examination of Crystal Dissolution in 3D: A Way to Reconcile Dissolution Rates in the Laboratory? Geochim[END_REF] For several reasons, pyroxenes represent a relevant target to address these aspects. Indeed, not only do they have a strongly anisotropic structure (silica chains running parallel to a given crystallographic axis), but they are also present worldwide and represent a major group of minerals for reactive-transport models. [START_REF] Berner | GEOCARB II: A Revised Model of Atmospheric CO[Sub 2] over Phanerozoic Time[END_REF]40 In two previous studies, [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF][START_REF] Bouissonnié | Dissolution Anisotropy of Pyroxenes: A Surrogate Model for Steady-State Enstatite Dissolution Resulting from Stochastic Simulations of the Hydrolysis Process[END_REF] we described a stochastic dissolution model operating at atomic scale and applied to enstatite, the magnesian endmember of pyroxenes. This model, validated by a comparison between simulations and experimental results, has shown its ability to reproduce the dissolution anisotropy of enstatite as well as the surface features observed at the end of the experiments (etch pit morphology, and relative thickness of the silica layers developed on each face). [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] However, in these two studies, we focused on facespecific dissolution, neglecting the contribution of "macrocorners" and "macro-edges", which can represent, for some minerals such as calcite, the most reactive sites. [START_REF] Noiriel | Examination of Crystal Dissolution in 3D: A Way to Reconcile Dissolution Rates in the Laboratory? Geochim[END_REF] In the present study, we propose to simulate the dissolution of enstatite rectangular parallelepipeds. Different sizes and aspect ratios are summarized in Figure 1. [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] Briefly, atoms contained in 120 the enstatite cell are placed using Pbca space group symmetry 121 elements, using the cell parameters (a = 18.233 Å, b = 8.8191 Å, 122 and c = 5.1802 Å) and the coordinates given by Hugh-Jones and 123 Angel. [START_REF] Angel | A Compressional Study of MgSiO3 Orthoenstatite up to 8.5 GPa[END_REF] Each magnesium and silicon atom is linked to oxygen 124 according to their coordination. However, the dissolution of 125 oxygen atoms is never considered in the dissolution model. 126 Indeed, they are considered as bridging atoms [START_REF] Salehikhoo | Magnesite Dissolution Rates 1082 at Different Spatial Scales: The Role of Mineral Spatial Distribution and 1083 Flow Velocity[END_REF] and the 127 dissolution model is applied on M-O br -M bonds (M being a 128 Mg or Si atom). All atoms linked to another one via a bridging 129 oxygen are part of its first coordination sphere. [START_REF] Salehikhoo | Magnesite Dissolution Rates 1082 at Different Spatial Scales: The Role of Mineral Spatial Distribution and 1083 Flow Velocity[END_REF] Only the first 130 coordination sphere is taken into account to compute the bond-131 breaking probabilities, since previous results showed that this 132 strategy results in a satisfactory agreement between modeled and 133 experimental results. [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] 134 Usually, the probability attributed with one event (i.e., the 135 bond-breaking probability) is written as follows [START_REF] Lasaga | Surface Chemistry, Etch Pits and Mineral-Water Reactions[END_REF] = - P e E k T / a B

METHODS

136

(1) 137 where P stands for the bond-breaking probability, E a is the 138 activation energy of the bond hydrolysis (J), k B is the Boltzmann 139 constant (J/K), and T is the temperature (K). In the present 140 study, the greatest bond-breaking probability (i.e., the bond 141 hydrolysis associated with the lowest activation energy, Mg-O-142 Mg) was arbitrarily set to 0.99 and the two other probabilities 143 were scaled according to the differences between their 144 corresponding activation energies following It has been demonstrated on Kossel crystal and calcite that, during the dissolution process, the shape of the crystal evolves, and edges and corners are rounded, leading to an almost spherical shape (for Kossel crystal, depending on the activation energy of the hydrolysis) when the primal material is a cube. In this study, we aim to determine whether these considerations are applicable on inosilicates, with a specific focus on enstatite. The study is based on simulations using a stochastic dissolution model applied at the atomic bond level and aims to explore different aspect ratios of the initial mineral parallelepipedal shape.

The Journal of Physical Chemistry C

= -- P P e E E k T A B ( ) / a A a B B

145

(2)

146 where A and B represent two different bonds. According to the 147 best agreement between modeled and experimental results 148 developed 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] the other two probabilities 149 were calculated as 0.4 and 0.0146 for Mg-O-Si and Si-O-Si, 150 respectively. The probability of backward reactions (i.e., the 151 creation of a new bond, previously destroyed or not) is not 152 considered here since the model is applied to reproduce 153 experiments run at pH 0, [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] where these reactions can be 154 neglected.

155

As it has been shown to successfully reproduce experimental 156 data and observations, an "all-or-none" approach was used in 157 this study. [START_REF] Salehikhoo | Magnesite Dissolution Rates 1082 at Different Spatial Scales: The Role of Mineral Spatial Distribution and 1083 Flow Velocity[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] Zhang | Si Order in Albite and Its Effect on Albite Dissolution Processes: A Monte Carlo Study[END_REF][36][START_REF] De Assis | Dissolution of Minerals with Rough Surfaces[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] This approach considers that an atom 158 is released only if all bonds that connect it to the surface are 159 simultaneously broken during a single iteration step. Hence, the 160 numerical model works directly on the atoms, and not on the 161 bonds. Then, the overall probability of an atom to be released 162 into solution during one iteration is given by the product of all of 163 the probabilities corresponding to the different bonds it shares 164 with other atoms

= = - - -- -- -- -- P P P e e nE k T mE k T n m M / / M O Mg M O Si M O Mg B M O Si B

165

(3) 166 where n and m stand for the number of bonds that the 167 considered atom shares with neighboring Mg and Si atoms, 168 respectively. 169 2.2. Dissolution Algorithm. At each iteration step, all 170 atoms located at the surface are scanned. An atom is considered 171 at the surface if its first coordination sphere is incomplete. A 172 random number Z is generated for each atom; if Z < P M , the 173 atom is released into the solution, if not, the scan of the atoms 174 moves on. At the end of the iteration, i.e., after scanning all atoms 175 located at the surface, the attributes of the atoms contained in 176 the first coordination sphere of those which have been dissolved during the iteration, are updated to consider this change, resulting in an increase of their departure probability at the next iteration. It is essential to update the attributes of the considered atoms only at the end of an iteration to avoid an "in cascade" dissolution, where the atoms would be released one after the other due to a direct change of the probabilities. A scheme of the f2 algorithm is shown in Figure 2.

Output of the Model. The dissolution was monitored

based on the release of Si atoms, as in classical studies dedicated to silicate dissolution kinetics. Note that because the dissolution at steady state (i.e., the dissolution rate can be considered to be constant) is stoichiometric, [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] this choice has no consequence on what will be described in the present study. However, the behavior of Mg release is important at the beginning of the dissolution process [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] and will be discussed in the section corresponding to this part of the dissolution timeline.

In the literature, the term "mineral dissolution rate" is usually used to refer to the release rate of atoms in the solution from a mineral surface, whether or not it is normalized to the contact surface between the fluid and the mineral. To avoid any confusion, in the following, we will use the term "dissolution rate" to refer to a velocity that is not normalized to the surface area (unit: mol/iteration) and "dissolution flux" to refer to dissolution rates normalized with respect to the surface area (unit: mol/m 2 /iteration). Both quantities are calculated at each iteration for any numerical simulation.

The cumulated number of released atoms has been used to calculate the total surface retreat during the simulations corresponding to the dissolution of single faces (i.e., only one face is considered to be in contact with the solution). The surface retreat is defined as the average height of material lost during the dissolution. It is calculated as follows

∑ Δ = * * * = h a b c S nb 16 (Si ) i hkl i n d n ( ) f 1 , (4)
Figure 2. Scheme of the dissolution algorithm used in the dissolution model. This scheme is adapted from Bouissonniéet al. [START_REF] Bouissonnié | Dissolution Anisotropy of Pyroxenes: A Surrogate Model for Steady-State Enstatite Dissolution Resulting from Stochastic Simulations of the Hydrolysis Process[END_REF] The Journal of Physical Chemistry C

where Δh i (hkl) is the surface retreat at the ith iteration (m); a, b, 211 and c are the cell parameters (m); S f is the surface area of the face 212 (m 2 ); and nb(Si d,i ) is the number of Si atoms dissolved at the ith 213 iteration. The number 16 corresponds to the number of Si atoms 214 in the enstatite cell. This method has been validated in 215 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 gives identical results to those obtained 216 if the mean coordinates in the direction normal to the surface 217 were used.

218

Finally, the last parameter used in this study is the 219 advancement of the dissolution reaction. It is defined as the 220 ratio between the number of Si atoms dissolved from the 221 beginning of the simulation and the initial number of Si atoms at 222 the beginning of the simulation. • The dissolution algorithm is applied on atoms rather than atomic bonds. This method is commonly used in stochastic dissolution studies, however, combined with the first assumption, this precludes the formation of amorphous silica layers resulting from interfacial dissolution-reprecipitation [START_REF] Hellmann | Nanometre-Scale Evidence for Interfacial Dissolution-Reprecipitation Control of Silicate Glass Corrosion[END_REF] or reorganization of dangling bonds. [START_REF] Gin | Dynamics of Self-Reorganization Explains Passivation of Silicate Glasses[END_REF] In the present study, such layers are suspected to have no impact on the dissolution rate because their passivating ability is thought to be pH-dependent [START_REF] Perez | Comparing the Reactivity of Glasses with Their Crystalline Equivalents: The Case Study of Plagioclase Feldspar[END_REF][START_REF] Wild | PH-Dependent Control of Feldspar Dissolution Rate by Altered Surface Layers[END_REF] and passivating layers are unlikely to form in strongly acidic environments. Such conditions may limit the extension of the conclusions of the study to less acidic environments. Contrary to the Si curve, which increases continuously at the beginning of the simulation, the Mg curve is marked by a sharp decrease of the number of atoms released per iteration. This is explained by the fast depletion of superficial Mg, which have a high probability of departure. After reaching the maximum value of the amount of Si and Mg released per iteration (corresponding to the iteration where the silica layer reaches a constant thickness), the dissolution rate decreases continuously and the dissolution is stoichiometric. Legends indicate how many times the crystal cell is multiplied in x, y, and z directions or exact dimension (black line). Each panel represents Si dissolution rate as a function of the reaction progress for the set of simulations that observe realistic aspect ratios for the grains, plotted in standard and logarithmic scales, respectively. Black triangles indicate the reaction progress where the dissolution rate reaches its maximum. Black circles indicate the reaction progress where the volume of the enstatite grain is equivalent to that of the initial smallest simulated rectangular parallelepiped (14 × 29 × 250). For each simulation, the dissolution rate increases at the beginning of the process, until reaching Aconitum maximum value. After the maximum value is reached, the dissolution rate decreases linearly until the reaction progress reaches ∼95%. The duration of the linear decreasing phase depends on the crystal size: the larger the crystal, the longer the linear phase.

The Journal of Physical Chemistry C

The second subgroup of simulations is composed of six 295 different volumes. While these volumes differ in height (length 296 in the z direction parallel to the c axis), the lengths in the x and y 297 directions were kept constant. Hence, the surface area of (100) 298 and (010) faces varied, unlike that of the (001) faces. The onset 299 of the simulations exhibits a similar behavior to that described 300 for the first subgroup: the first 3% of the reaction is characterized 301 by a sharp increase (from 0 to 350-600 Si atoms/iteration) of 302 the release rate of Si. This step is followed by a linear increase of 303 the dissolution rate until it reaches its maximum. Contrary to the 304 first subgroup of simulations, the advancement corresponding to 305 this maximum depends on the size of the initial volume. Indeed, 306 if the maximal dissolution rate is greater for the largest volume 307 (1390 Si atoms/iteration for the volume 40 × 75 × 500) 308 compared to the smallest one (1190 Si atoms/iteration for the 309 volume 40 × 75 × 250), the advancement at which the 310 maximum is reached is lower (0.30 against 0.45, respectively). 311 Another difference between the first and second subgroups is the behavior of the dissolution after the maximum is reached. While, in the first case, the dissolution rate decreases sharply until the end of the reaction, in the second case, a slow and linear decrease is observed from an advancement of 30-45 to 85-90% depending on the considered volume. During this linear decrease, regardless of the volume, all simulation results are almost identical, showing the same dissolution rate of Si atoms at the same advancement of reaction. The last 10-15% of the reaction is characterized by a sharp decrease of the dissolution rate of Si atoms to finally reach 0 when the reaction is complete f4 (Figure 4). The second group of eleven simulations is shown in Figure 5.

Contrary to the previous group of simulations, the aspect ratios which were used are representative of natural enstatite crystal habit. In spite of these differences, the evolution of the dissolution rate of Si atoms is almost identical to that described for the second subgroup of the first group described above, with the following steps: (1) a sharp increase of the Si dissolution rate during the first 3% of the dissolution, (2) a linear increase until the maximum is reached (from 145 to 3660 Si atoms released per iteration, reached after a reaction progress of 40-20% for the smallest and the largest volumes, respectively), (3) a linear decrease in the release rate of Si until 90-95% of the reaction, and finally (4) a sharp decrease in the release rate of Si over the last 5-10% of the reaction.

Evolution of Mg Dissolution. The behavior of Mg

dissolution is similar to the one described above for Si dissolution (Figure 4). The main differences are observed at the very beginning and at the end of the simulations. Mg dissolution evolution has been studied only for simulated volumes representative of natural aspect ratios of enstatite. Because all simulations, as for Si atoms, exhibit a similar behavior, only one case is shown in Figure 4 (40 × 82 × 700).

The global trend of Mg dissolution evolution is as follows:

• A dramatic decrease of the dissolution rate at the beginning of the reaction (325 iterations for the volume shown in Figure 4).

• An increase of the dissolution rate from the 315th to the 8000th iterations. This increase is marked by two pulses for the volume presented in Figure 4.

• A linear, slow decrease of the dissolution rate, similar to the one observed in the case of Si dissolution. In the 40 × 82 × 700 volume, this step extends from the 8000th to the 26 990th iterations.

• A more pronounced decrease in the dissolution rate for the last iterations until all of the Mg atoms have been dissolved.

Dissolution Normalized Per Surface Unit.

Because the most realistic cases treated in this study are those that observe the aspect ratio of pyroxenes, only the corresponding simulations were investigated with respect to a normalization per surface unit. The evolution of the surface area is not derived from direct outputs of the simulations. Instead, it is assumed that the surface decreases according to the retreat of each single face of the simulated volumes. 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] we investigated the surface retreat of (100), (010), and (001) faces, demonstrating that our model was able to reproduce observed experimental data. Similar simulations were then used to calculate the evolution of the surface area of the simulated crystals according to the following equation (5) 373 where A(i) is the surface area at the ith iteration; x 0 , y 0 , and z 0 are 374 the initial lengths in the x, y, and z directions, respectively; and 375 Δh i (hkl) is the surface retreat of the corresponding (hkl) face at 376 the ith iteration, derived from the simulations performed on 377 single faces where the edges were constrained to be nonreactive. 378 This method, however, does not take into account a possible 379 rounding of the edges, often observed in mineral dissolution 380 experiments and modeling. [START_REF] Luẗtge | A Stochastic Treatment of Crystal Dissolution Kinetics[END_REF][START_REF] Noiriel | Examination of Crystal Dissolution in 3D: A Way to Reconcile Dissolution Rates in the Laboratory? Geochim[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] The number of Si atoms 381 released at each iteration is then divided by the surface area at 382 the corresponding iteration. Results are shown in Figure 6. The trends followed by the 384 dissolution fluxes differ from those obtained for dissolution 385 rates. The first 3.5-1% of reaction is characterized by a sharp 386 increase in the dissolution flux from 0 to 5.4-6.2 × 10 15 Si 387 atoms/(m 2 iteration) for the smallest and largest volumes, 388 respectively. The increase in the dissolution flux is then less 389 pronounced over the next 42.5-19% of the reaction. The Si 390 dissolution flux reaches 1.6 × 10 16 Si atoms/(m 2 iteration) at 391 47% of the reaction progress for the smallest volume against 1.8 392 × 10 16 Si atoms/(m 2 iteration) at 20% for the largest one. After this period, the dissolution flux tends to stabilize for the three largest volumes but only over a short reaction progress (approximately 10% of the total reaction). After this phase, the dissolution flux of the largest volumes follows that of the smallest ones. The dissolution flux still increases, but the difference between each iteration is less pronounced until the reaction progress reaches 90-96% (where the dissolution flux reaches its maximum: (2.2-3.5) × 10 16 Si atoms/(m 2 iteration) for the smallest and largest volumes, respectively). Finally, from 90 to 96% of the dissolution, the Si dissolution flux decreases sharply to reach 0 when the volume has been completely consumed. However, a clear pattern is highlighted when this evolution is considered as a function of the number of iteration steps (i.e., f7 time) (Figure 7a). The dissolution rate is always at maximum around the same iteration (≈8000th) for the different volumes. In each case, the dissolution flux increases at the beginning of the simulations. After reaching a given value of the reaction progress (corresponding to the 8000th iteration), the increase slows down and the dissolution flux can be considered as constant for few iterations (the larger the crystal, the longer the "steady-state" regime). Then, a moderate increase is observed until the dissolution reaches its maximum. After reaching its maximum, the dissolution flux decreases dramatically and reaches zero when all of the atoms have been dissolved.

The Journal of Physical Chemistry C F = [ -Δ -Δ + -Δ -Δ + -Δ -Δ ] A i x h y h y h z h x h z h ( ) 2 ( 2 )( 2 ) ( 2 )( 2 ) ( 2 )( 2 ) i i i i i i 0 ( 100

The Journal of Physical Chemistry C

417 Interestingly, this value corresponds to the one where the surface 418 retreat of the (001) face reaches the steady state when the set of 419 probabilities used is [0.99, 0.4, 0.0146] for [Mg-O-Mg, Si-420 O-Mg, Si-O-Si] bonds respectively. [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] It was interpreted as 421 the moment when the dissolution becomes congruent. 422 Considering that the other faces reach their steady-state 423 dissolution regime before the (001) face, [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] the fact that the 424 dissolution rate of the simulated rectangular parallelepipeds is at 425 maximum at the iteration corresponding to the attainment of a 426 steady state for face (001) may mean that the evolution of the 427 dissolution regime of a simulated rectangular parallelepiped 428 enstatite is controlled by the face for which this steady state is 429 reached the latest (after the largest number of iterations). 430 Alternatively, another explanation could be that the face (001) 431 being the fastest to dissolve, the amount of Si released from this 432 face is greater than those of the other faces. Hence, this relation 433 between a simulated rectangular parallelepiped enstatite and the 434 (001) face reaching their steady state at the same iteration may 435 come from the fact that a significant proportion of their 436 dissolution is supported by the (001) face.

437

Regarding the dissolution rate of Mg atoms, Figure 4 shows 438 that it decreases dramatically during the first few hundred 439 iterations (325 for the volume 40 × 82 × 700). This can be 440 explained by the fast departure of Mg atoms from the reactive 441 surface at the beginning of the simulations because of their greater probability to be dissolved. This results in a fast leaching step at the beginning of the reaction, when these atoms are numerous. After a few iterations, the main coordination of Mg at the surface increases prior reaching a constant value, [START_REF] Bouissonnié | Dissolution Anisotropy of Pyroxenes: A Surrogate Model for Steady-State Enstatite Dissolution Resulting from Stochastic Simulations of the Hydrolysis Process[END_REF] making them more difficult to dissolve.

The fast decrease of Mg dissolution rate can also be linked to the increase of Si dissolution rate. Indeed, the fewer the Mg-O-Si bonds, the higher the Si departure probability. In other words, when numerous Mg atoms are released into solution, the dissolution of Si atoms is enhanced. This departure of Si atoms resulting from the fast leaching of Mg atoms at the surface also makes it easier for the "solution" to reach the Mg atoms located deeper in the simulated crystals. These combined effects possibly explain why, after a dramatic decrease, the Mg dissolution rate increases with Si dissolution, until the maximum is reached. After the maximum is reached, both Si and Mg dissolution rates decrease. More importantly, these decreases are almost identical (same slope), indicating that the dissolution is stoichiometric.

4.1.2. Parameters Controlling the Dissolution. Although the first set of simulated rectangular parallelepipeds was not representative of the aspect ratio of natural enstatite grains, [START_REF] Dana | The System of Mineralogy[END_REF] it allows highlighting some behaviors and parameters that may control the dissolution. Indeed, the two subgroups presented in Figure 3 have their own patterns that can help understand the 473 While the beginning of the dissolution differs depending on the 474 volume, the stage depicting a decreasing reactivity is 475 characterized by a similar trend and value of Si dissolution 476 rates for all simulations (Figure 3b). Therefore, it can be 477 assumed that this stage is controlled mostly by the surface area of 478 the (001) faces. Then, the surface area of the (100) and (010) 479 faces would only affect the reaction progress required to reach 480 this stage (the greater these surfaces, the longer the preliminary 481 stage). 3a). While this similarity supports that the dissolution is primarily controlled by the face (001), it also suggests that the other two faces have only a minor contribution to the overall dissolution rate of enstatite grains. This is further supported by the differences in reaction rates between the three faces ((100), (010), and (001)). The dissolution rate of the (001) face is more than 1 order of magnitude higher than those of faces (100) and (010). To have a significant impact on the total dissolution rates, the surface area of the (100) and (010) faces should be at least The Journal of Physical Chemistry C I 498 10 times greater than that of the (001) face, which was never 499 investigated in this study due to computational capacities and 500 our commitment to perform simulations that observe as much as 501 possible the aspect ratios of enstatite grains encountered in 502 natural environments.

503

The difference in shape between the released Si curves of the 504 two subgroups (Figure 3a vs Figure 3b) is more difficult to 505 explain. The main difference in the simulation parameters 506 resides in the dimensions of the simulated volumes and mainly, 507 their heights. Indeed, the heights (length in the z direction) of 508 the subgroup composed of the 50 × 103 × 176, 75 × 150 × 75, 509 150 × 75 × 75, 90 × 90 × 180, and 120 × 120 × 120 rectangular 510 parallelepipeds are smaller than those of the volumes of the 511 second subgroup. It is therefore possible that the former 512 volumes are consumed before reaching a stage where the Si 513 dissolution rate decreases linearly, which would explain this 514 difference in behavior. If the overall dissolution is mainly 515 controlled by the dissolution of the (001) faces, it is indeed 516 possible that these volumes cannot reach this dissolution stage 517 because their height is too short. This is consistent with the 518 evolution of the dissolution rate of the other subgroup for which 519 the smaller the volume, the shorter the linear stage. It can be 520 reasonably assumed that rectangular parallelepipeds with a very 521 short height never reach a stage where the dissolution rate 522 decreases linearly.

523

To conclude, the simulations shown in Figure 3 provide 524 important information regarding the global dissolution of a 525 rectangular parallelepiped enstatite: the dissolution rate seems 526 mainly controlled by that of the (001) faces.

527

The simulations conducted with the second set of rectangular 528 parallelepipeds, more representative of the natural crystal habit 529 of pyroxenes, also exhibit specific patterns. Indeed, while the 530 amount of Si atoms released at each iteration varies from one 531 volume to another, trends and shapes are identical (Figure 5). 532 This could indicate that the overall evolution of the dissolution is 533 directly controlled by the shape of an enstatite grain. For the 534 same aspect ratio, the initial volume only influences the absolute 535 amount of released atoms.

536

This simulation set allows focusing on a quantity that 537 delineates a transition in the evolution of the dissolution rate 538 of a rectangular parallelepiped enstatite, namely, the maximum 539 value this latter can reach (black triangles in Figure 5a). While it 540 is obvious that the bigger the volume of a rectangular 541 parallelepiped, the greater the instantaneous Si dissolution 542 rate, interestingly, it is possible to link these two quantities f8 543 mathematically (Figure 8b). In the same way, it is possible to 544 link this maximum rate to the reaction progress (Figure 8a). The 545 corresponding relations are given by Based on this equation, it is theoretically possible to estimate 553 the time at which the release of Si is maximal when knowing the 554 initial volume of an enstatite grain. Such empirical relations are 555 very relevant for experimental works and practical applications because it corresponds to the attainment of the steady-state regime of the dissolution (see Section 4.5).

Variability of Dissolution Rates for a Given Grain

Size. A long-standing concern in mineral/water interaction deals with the uniqueness of the dissolution rate. [START_REF] Arvidson | Variation in Calcite Dissolution Rates[END_REF][START_REF] Fischer | How Predictable Are Dissolution Rates of Crystalline Material? Geochim[END_REF] While the variability of the dissolution rate is known to be significant when the geochemical conditions are changed, it is often assumed, particularly in reactive-transport studies, that the intrinsic dissolution rate constant remains unchanged when the fluid chemistry and the temperature are constant. However, several studies tend to show that this assumption may be too simplistic, as the dissolution rate may vary with time. While some studies highlight a monotonic decrease of the dissolution rate with time which can be associated with several factors such as the development of amorphous silica-rich layers, [START_REF] Daval | Time-Dependent Feldspar Dissolution Rates Resulting from Surface Passivation: Experimental Evidence and Geochemical Implications[END_REF] others emphasize a complex evolution of the dissolution rate resulting from the gradual modification of the crystal habit, [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] Luẗtge | A Stochastic Treatment of Crystal Dissolution Kinetics[END_REF][START_REF] Carrasco | Modeling the Kinetics of Calcite Dissolution in Neutral and Alkaline Solutions[END_REF] the inherited history of the crystal reactivity [START_REF] Fischer | Inherited Control of Crystal Surface Reactivity[END_REF][START_REF] Beig | Albite Dissolution Kinetics as a Function of Distance from Equilibrium: Implications for Natural Feldspar Weathering[END_REF] or the microstructural differences (e.g., defect density) between two individual crystals of a same mineral. 3,5,[START_REF] Arvidson | Variation in Calcite Dissolution Rates[END_REF][START_REF] Luẗtge | A Stochastic Treatment of Crystal Dissolution Kinetics[END_REF][56][57][58] Therefore, it is central to understand if two minerals placed under the same biogeochemical conditions but differing in their history (one being fresh and the other one being aged for example) will have the same dissolution rate; hence, if two minerals of the same volume have the same unique dissolution rate.

To investigate this question using the dissolution model, it has been decided to focus on the dissolution rate of the rectangular parallelepipeds when their volume reaches a volume equivalent to that of the smallest simulated crystal ( 14 Interestingly, a relation between the initial volume and the dissolution rates for V = V 14×29×250 exists

= × = × × V Si 1 0 V V released,
12.0 0 0.44 14 29 250 (9)

This kind of relationships means that for a given volume, different values of dissolution rates exist depending on the initial size of the minerals.

Simple Mathematical Model of Si Dissolution

Rate. The evolution of Si dissolution rates and fluxes has a specific behavior as described above. If this behavior seems difficult to reproduce with a theoretical model, it is possible to test hypotheses on how the dissolution evolves: here, we suggest testing the assumption that the overall enstatite grain dissolution simply corresponds to the sum of the contribution of the various single faces delineating the enstatite volume. This assumption may be true if, contrary to what is often proposed in the literature, the edges and corners do not have a strong impact on

The Journal of Physical Chemistry C 614 the dissolution rate/flux and thus, if these sites do not lead to an 615 excessive rounding of the crystal. In this case, a simple cross-616 multiplication would be able to reproduce the evolution of the 617 dissolution rate/flux.

618

To model the dissolution of Si atoms, the first step is to 619 consider the volume (ΔV i ) that is dissolved during an iteration . Scheme of the mechanism driving the dissolution of Kossel and enstatite crystals. Kossel: Black cubes represent atoms that are not at the mineral reactive surface. Purple, orange, and yellow cubes represent surface-, edge-, and corner-site atoms, respectively. During the dissolution, the atoms with fewer neighbors have a higher probability to be dissolved. This leads to the formation of rounded corners and edges. The atoms at the surface with only one neighbor missing (purple) are particularly difficult to dissolve compared to the others. Enstatite: at the beginning of the dissolution, Mg atoms are dissolved first. This is particularly visible on the (001) face where they exhibit a less populated first coordination sphere, leading to Mg depletion of the surface and the formation of a Si-enriched layer with a thickness not exceeding a few atomic layers. Because silica chains are not connected to each other, the probability of departure of Si atoms is the same on the entire (001) face, making it unlikely to exhibit rounded edges and corners.

Δ = - - V V V i i i 1 620 (10a) = -Δ -Δ -Δ V x h y h z h ( 2 )( 2 )( 2
) i i i i 0 (100) 0 (
The Journal of Physical Chemistry C 663 observation stems from the fact that stochastic variations in the 664 instantaneous dissolution rate of each single face are cumulated, 665 increasing the overall standard variation artificially. Interestingly, 666 when one looks at the variation of the dissolution rate as a 667 function of the number of elapsed iterations (Figure 7a), the 668 mean value of the results derived from the face-specific approach 669 is slightly lower than the one derived from the simulations 670 performed with rectangular parallelepipeds. This slight differ-671 ence may be explained by the contribution of corners and edges 672 to the dissolution rate/flux. However, contrary to other 673 minerals, the dissolution of these surface sites seems to have a 674 negligible impact. A scheme exhibiting the difference between 675 enstatite and a Kossel crystal 9,19,36,37,59-61 dissolution, for which 676 the contribution of edges and corners was shown in previous f9f10 677 studies to be significant, is depicted in Figures 9 and10 to 678 explain this counterintuitive result. These schemes are based on 679 the periodic bond chain theory, [START_REF] Hartman | On the relations between structure and morphology of crystals[END_REF] which governs silicate 680 dissolution. In a Kossel crystal, the volume is composed of a 681 unique component ("atom"), leading to the existence of a single 682 type of bond. By nature, a Kossel crystal will have specific surface 683 sites (adatoms, steps, and kinks) with a high probability of 684 departure, leading to the rounding of the simulated crystals over 685 time. In the case of enstatite, the anisotropic nature of pyroxenes 686 makes this morphological evolution simpler. As opposed to, e.g., 687 tectosilicates, which have a three-dimensional framework of 688 silica tetrahedrons, inosilicates are composed of silica chains 689 oriented following a specific direction (parallel to the c-axis in 690 enstatite). In enstatite, Mg and Si atoms are also organized in 691 planes parallel to the (100) face. This particular organization 692 explains the anisotropy of dissolution of enstatite, as discussed in 693 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 likely explains why macro steps and 694 corners do not play an essential role in the evolution of the dissolution rate with time. Indeed, as shown in Figure 9, most of enstatite dissolution flux is supported by (001) faces. This is due to the fact that the release of Si from (hk0) faces requires the simultaneous cleavage of two Si-O-Si bonds, while the release of Si from (001) faces requires the cleavage of one Si-O-Si bond only. This specific anisotropic structure also explains why the model depicted in Section 4.2 successfully matches the outputs of the simulations run with enstatite parallelepipeds, indirectly suggesting that no rounding of the crystals occurs with time. The greater probability to release Mg atoms to the solution results in the formation of a Mg-depleted zone at the mineral surfaces during the first iterations (Figure 4), where the mineralsolution interface is only composed of chains of Si atoms (Figure 9). This process is more pronounced at the surface of the (001) faces, where Mg atoms have the greater departure probabilities. This is fully consistent with the experimentally observed thicker amorphous silica layers on (001) faces. If the (001) faces exhibit a Mg-depleted mineral-solution interface, then it means that it is composed essentially on silica chains, not connected between them (Figure 9). These silica chains will be present throughout the entire surface of the (001) face but, as their topmost Si atoms all exhibit the same departure probability, a rounding of the edges and kinks is unlikely to happen.

In summary, the dissolution of a Kossel crystal can be with only disconnected silica chains. This may present a limitation to the mechanisms developed above as, in natural environment, a thicker amorphous silica layer may impact the crystal shape evolution. Extending this discussion to other minerals requires a specific treatment of other classes of minerals, which is out of the scope of this study. However, it is likely that many minerals, which highlight an anisotropy similar to pyroxenes could observe a similar dissolution mechanism. One could argue that many minerals are more anisotropic than enstatite because of its orthorhombic Bravais lattice, which has a high symmetry level. However, the notion of anisotropy here refers to the organization of the Si tetrahedrons.

In this regard, orthorhombic inosilicates are less isotropic than, e.g., triclinic tectosilicates, which have a three-dimensional framework of silica tetrahedrons, even though the symmetry level of the orthorhombic system is higher than that of the triclinic system.

Finally, the main discrepancy between model's outputs and theoretical results is observed at the end of the simulations. While the dissolution rates calculated using the simulations run with rectangular parallelepipeds reach 0 at the end of the reaction, indicating that all of the material is consumed, it is not the case for the rates calculated based on face-specific dissolution rates and using eqs 10a and 10b. This observation is an artifact resulting from the mathematical construction of the model. While faces (100), (010), and (001) are considered, the Figure 10. Scheme of the effect of dissolution through time on Kossel and enstatite crystals. Kossel: the high probability of departure of edge and corner sites leads to a rounding of the corresponding locations. This modifies significantly the crystal shape during the dissolution, which evolves from a cubic to a spherical shape. Black squares represent the crystal shape at t = 0. The gray area represents a section of the crystal. Enstatite: At the beginning of the dissolution, Mg atoms are preferentially dissolved due to their particularly incomplete first coordination sphere, leading to the formation of a Si-enriched layer (orange area) on the (001) face (a-c). The Si-enriched layer reaches its constant depth (b, c) and the dissolution becomes stoichiometric. Both Si and Mg atoms are dissolved preferentially on the (001) face, but a few are dissolved on the (100) and ( 010) faces (c-e). This dissolution behavior makes the rounding of the crystals unlikely.

The Journal of Physical Chemistry C dissolution is mainly controlled by the sole (001) face. Equation 10b implies that all faces dissolve independently so that when the length of the parallelepiped following the z-axis reaches 0, it is not the case for the two other axes. Therefore, following eq 10b, a plan of enstatite contained in the two-dimensional (2D) (xOy) plan still exists and continues to fuel the dissolution. Conversely, in the real world (and in the simulations run with rectangular parallelepipeds), the dissolution stops when the crystal is completely dissolved. In the previous section, it has been shown that it is possible to reproduce the grain dissolution behavior with a simple mathematical model considering the dissolution fluxes of the faces in contact with the fluid. However, this good agreement was obtained using simulation results only. While, theoretically, the model used in the present study is also able to reproduce experimental data, [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] the corresponding validation has been previously based on a comparison with measured face-specific dissolution fluxes only, and not with results from powder experiments. This missing comparison is therefore attempted hereafter. As a first step, we use below the additivity property of single-face dissolution fluxes illustrated in Section 4.3 to compare the results of dissolution experiments previously conducted with single faces with those conducted with powdered enstatite.

In their study, Oelkers and Schott 11 have determined the dissolution flux of enstatite powder over a wide range of pH and temperature. Using the reaction order with respect to [H + ] and the activation energy derived from their study, the dissolution flux of enstatite powder at 90 °C and pH 0 (the chemical conditions simulated in this study) should be around 2.35 × 10 -8 mol/(m 2 s). In our previous study, the different dissolution fluxes of the four main faces naturally exposed in enstatite have been determined to yield 5.22 × 10 -9 , 5. 92 × 10 -9 , 1.7 × 10 -8 , and 1.44 × 10 -7 mol/(m 2 s) for the faces (100), (010), (210), and (001) respectively [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] (note that the results of Bouissonniéet al. [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] were corrected by a 2-fold factor, as they considered Mg 2 Si 2 O 6 as the chemical formula for enstatite, whereas Oelkers and Schott used MgSiO 3 ). The main cleavage plan of enstatite is {210} 62 and is likely to be the most abundant when enstatite is crushed prior to run powder dissolution experiments. If one assumes that, in enstatite powders, mostly (210) and (001) faces are present, then based on the results of the previous section, the dissolution flux of the powder can be calculated following:

= * * + * * r S r S r S 4 2 (210) (210) (001) (001) 0 ( 12
)
where S (hkl) is the surface area of the (hkl) face (m 2 ), r (hkl) is the dissolution flux (mol/(m 2 s)) of the (hkl) face, and S 0 is the initial surface area (m 2 ). Depending on the size of the grains used in this calculation (10 × 10 × 100 or 20 × 20 × 100 μm 3 , consistent with the grain size range reported in the Oelkers and Schott study), the dissolution flux is estimated to be 2.85 × 10 -8 and 2.31 × 10 -8 mol/(m 2 s) for 20 × 20 × 100 and 10 × 10 × 100 μm 3 crystals, respectively. These values are in excellent agreement with the dissolutions flux calculated based on the study of Oelkers and Schott (hereafter referred to as r OS ) (2.35 × 10 -8 mol(m 2 s)). This agreement further supports that edges and corners must negligibly contribute to the dissolution flux of enstatite grains, consistent with the previous section.

Δ = - -Δ -Δ -Δ Δ | Δ | Δ r V t a b c t S N V t V x h y h z h 16 ( ) ( ) ( 2 )( 2 )
( 2 )

t
where r is the dissolution flux (mol/(m 2 s)), Δt is the time interval considered for the calculation (s), S 0 is the initial surface (m 2 ), N A is the Avogadro number (mol -1 ), ΔV(Δt) is the difference of volumes over the time interval Δt (m 3 ), and Δh Δt

(hkl)
is the surface retreat of the (hkl) face over the time interval Δt.

The second method uses the surrogate relation between the dissolution flux and the probabilities (i.e., activation energies) determined in Bouissonniéet al. [START_REF] Bouissonnié | Dissolution Anisotropy of Pyroxenes: A Surrogate Model for Steady-State Enstatite Dissolution Resulting from Stochastic Simulations of the Hydrolysis Process[END_REF]

Δ

Δ + - α β γ α β γ ρ ω -- -- -- -- -- -- - h
Ä Ç Å Å Å Å Å Å Å Å Å Å Å Å Å Å É Ö Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ (14)
with k, k dis , α, β, γ, α dis , β dis , γ dis , and ω being empirical parameters and ρ d the dislocation density. The left-hand side of the sum referred to the "bulk" dissolution flux, while the righthand side referred to the contribution etch pits at dislocation Bouissonniéet al. [START_REF] Bouissonnié | Dissolution Anisotropy of Pyroxenes: A Surrogate Model for Steady-State Enstatite Dissolution Resulting from Stochastic Simulations of the Hydrolysis Process[END_REF] and are listed in Table 1.

891

The last method considers the modeled face-specific [START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF][START_REF] Bouissonnié | Dissolution Anisotropy of Pyroxenes: A Surrogate Model for Steady-State Enstatite Dissolution Resulting from Stochastic Simulations of the Hydrolysis Process[END_REF] are reported in the last column. Crystal size has been chosen according to that reported in Oelkers and Schott, and an aspect ratio consistent with Dana. [START_REF] Dana | The System of Mineralogy[END_REF] The Journal of Physical Chemistry C by Oelkers and Schott, which is consistent with the fact that {210} represents the main cleavage plan of enstatite.

4.5. Relevance for Mineral Dissolution Studies and Limitations. Both experimental results obtained with facespecific dissolution measurements and those obtained with the stochastic dissolution model run with rectangular parallelepipeds succeed to reproduce the dissolution rates/fluxes measured during powder dissolution experiment. Here, we showed that face-specific and powder dissolution rates can be linked using a simple relation. However, this may be only true in the case of very anisotropic minerals such as enstatite (and pyroxenes) because edges and corners contribute negligibly to the dissolution flux. This conclusion does not hold for calcite, for example, as Noiriel et al. [START_REF] Noiriel | Examination of Crystal Dissolution in 3D: A Way to Reconcile Dissolution Rates in the Laboratory? Geochim[END_REF] showed that corners dissolve faster than (104) cleavage plans. In our opinion, our results would justify investigating the dissolution of other minerals according to the methods reported here to determine the crystallographic structures for which dissolution obeys relations similar to those evidenced for enstatite.

The results presented here may also be useful to model silicate weathering in natural environments. Subsurface environments present a wide diversity of textures, from porous to fractured media. Although fractured media should involve a preferential water/solid interaction of a limited number of mineral cleavage planes, porous media can exhibit a greater variety of mineral faces in contact with water. Actually, as a first approximation, micrometer-wide enstatite grains may be simplified as rectangular parallelepipeds made of (210) and (001) faces only, for which the length of the crystal in the z-direction is 5-10 times greater than in the other directions leading to consider the dissolution rate given by eq 12. This result may still hold true in less acidic environment as long as the fluid is far from equilibrium with respect to SiO 2(am) . [START_REF] Daval | Influence of Amorphous Silica Layer Formation on the Dissolution Rate of Olivine at 90 °C and Elevated PCO(2)[END_REF] The fact that facespecific and powder dissolution studies can be linked through a simple relation for anisotropic minerals may imply that facespecific dissolution studies offer a greater interest than those conducted with powders for minerals presenting an anisotropy similar to the one of enstatite. In addition, our study shows the prevalence of fast dissolving faces over the others, suggesting that face-specific dissolution studies should target in priority such faces. Overall, the agreement between the outputs of the model and the theoretical results detailed in Section 4.2 shows how important it is to consider the evolution of the reactive surface area in geochemical/reactive-transport studies/models. The evolution of the aspect ratio was shown to be important in both cases when the reaction progress reaches approximately 80%. Indeed, at this step, the aspect ratio has evolved enough so that the reaction flux is characterized by a sharp increase. This may be of great interest to model chemical weathering in the regolith, where primary materials are altered into secondary minerals and may be the main contributors of dissolved elements in springs and rivers.

Interestingly, the simulated instantaneous release rates of Mg and Si atoms as a function of time (Figure 4) exhibit a specific shape that corresponds in every way to the calorimetry curves depicting the progress of cement hydration (cf. Figures 2-4 in Luc et al. [START_REF] Skrotzki | Defect Structure and Deformation Mechanisms in 1239 Naturally Deformed Hornblende[END_REF] ), for which the rate-limiting step is supposed to be the dissolution of silicate minerals such as alite. It is noteworthy that such a shape has puzzled the materials sciences community for a couple of decades, since it reveals a strong nonlinearity of the dissolution rate with time, which is at odds with existing models of mineral reactivity. 64 While previous attempts did not succeed to fully account for such observations, 65,66 the present

  Model Description. A complete description of the 119 model is given in Bouissonniéet al.

Figure 1 .

 1 Figure1. Schematic representation of the objectives of this paper. Dissolution is often assumed to take place on reactive sites such as edges and corners. It has been demonstrated on Kossel crystal and calcite that, during the dissolution process, the shape of the crystal evolves, and edges and corners are rounded, leading to an almost spherical shape (for Kossel crystal, depending on the activation energy of the hydrolysis) when the primal material is a cube. In this study, we aim to determine whether these considerations are applicable on inosilicates, with a specific focus on enstatite. The study is based on simulations using a stochastic dissolution model applied at the atomic bond level and aims to explore different aspect ratios of the initial mineral parallelepipedal shape.

  223 2.4. Assumptions and Limitations of the Model. The 224 model has proven its ability to reproduce both experimental 225 dissolution rates and observed surface features evolutions. 226 However, several assumptions (all based on the studied 227 chemical condition) are considered and may impact the 228 model outputs and their upscaling to study natural environment: • No backward reactions (i.e., formation of atomic bonds) are implemented. This simplification was made because of the very low pH of the studied solution (pH 0). This pH makes the solution to be very far from the thermodynamic equilibrium with respect to enstatite. Therefore, backward reactions (i.e., formation of atomic bonds) are considered unlikely to happen.

Figure 3 .

 3 Figure 3. Results of the dissolution model showing the dissolution rate (defined as the number of Si atoms released per iteration) as a function of the reaction progress. Legends indicate how many times the crystal cell is multiplied in x, y, and z directions. The results reported here represent sets of simulations applied on arbitrary parallelepiped sizes (i.e., not representative of the aspect ratio of natural enstatite). (a) First subgroup of simulations, with arbitrary rectangular parallelepiped dimensions. (b) Second subgroup: Length and width (parallel to the a and b axes of the enstatite cell) of the rectangular parallelepiped were kept constant while the height (parallel to the c axis) was modified in the different simulations. In both cases, a similar behavior is observed. (a) Simulations highlight a parabolic behavior of the dissolution. The maximum seems to be observed at similar values of the reaction progress. (b) Contrary to the first set, the simulations of the second set are different at the beginning of the dissolution.Here, the larger the crystal, the higher the maximum. After reaching the maximum, the number of Si atoms dissolved decreases linearly. The slope is the same for each simulation.

Figure 4 .

 4 Figure 4. Example of the results of a simulation conducted on a 40a × 82b × 700c volume. This figure shows the dissolution rate (number of dissolved atoms per iteration) as a function of the number of iteration steps. The cyan curve represents Si atoms, and the dark yellow curve represents Mg atoms.Contrary to the Si curve, which increases continuously at the beginning of the simulation, the Mg curve is marked by a sharp decrease of the number of atoms released per iteration. This is explained by the fast depletion of superficial Mg, which have a high probability of departure. After reaching the maximum value of the amount of Si and Mg released per iteration (corresponding to the iteration where the silica layer reaches a constant thickness), the dissolution rate decreases continuously and the dissolution is stoichiometric.

Figure 5 .

 5 Figure 5. (a, b) Results of the dissolution model showing the dissolution rate (defined as the number of Si atoms released per iteration) as a function of the reaction progress.Legends indicate how many times the crystal cell is multiplied in x, y, and z directions or exact dimension (black line). Each panel represents Si dissolution rate as a function of the reaction progress for the set of simulations that observe realistic aspect ratios for the grains, plotted in standard and logarithmic scales, respectively. Black triangles indicate the reaction progress where the dissolution rate reaches its maximum. Black circles indicate the reaction progress where the volume of the enstatite grain is equivalent to that of the initial smallest simulated rectangular parallelepiped (14 × 29 × 250). For each simulation, the dissolution rate increases at the beginning of the process, until reaching Aconitum maximum value. After the maximum value is reached, the dissolution rate decreases linearly until the reaction progress reaches ∼95%. The duration of the linear decreasing phase depends on the crystal size: the larger the crystal, the longer the linear phase.

  f5

  Evolution of Mg and Si Dissolution Rates. 4.1.1. Evolution during the First 8000 Iterations. Despite the differences in volumes used in the simulations, all simulation outputs are characterized by a sharp increase of the number of Si atoms released at the beginning of the dissolution. Figures 3a and 5b show this increase as a function of the advancement of the dissolution reaction. Depending on the volumes, the advancements corresponding to this stage are different.

Figure 6 .

 6 Figure 6. Dissolution flux as a function of the reaction progress. Legends indicate how many times the crystal cell is multiplied in the x, y, and z directions. (a) Results of the simulations. (b) Results of the simulations and modeled theoretical result for enstatite grain with realistic grain size (black curve). The black patch at the beginning of the realistic grain size dissolution corresponds to the stages 1 and 2 of Figure 9a,b.In each case, the dissolution flux increases at the beginning of the simulations. After reaching a given value of the reaction progress (corresponding to the 8000th iteration), the increase slows down and the dissolution flux can be considered as constant for few iterations (the larger the crystal, the longer the "steady-state" regime). Then, a moderate increase is observed until the dissolution reaches its maximum. After reaching its maximum, the dissolution flux decreases dramatically and reaches zero when all of the atoms have been dissolved.

Figure 7 .

 7 Figure 7. (a) Dissolution rate as a function of the total number of iteration steps (equivalent to time). Legends indicate how many times the crystal cell is multiplied in the x, y, and z directions. Light-colored curves represent the evolution of the number of Si atoms released per iteration based on theoretical considerations developed in Section 4.2. Dark-colored curves represent the simulation results of the stochastic dissolution of rectangular parallelepipeds. (b) Dissolution flux as a function of the total number of iteration steps. Light-colored curves represent the evolution of the number of Si released per iteration and per surface unit based on theoretical considerations developed in Section 4.2. Dark-colored curves represent the simulation results of the stochastic dissolution of rectangular parallelepipeds. In each panel, the agreement between theory and simulation outputs highlights the slight effects of edges and corners on the global rate/flux. The difference between theory and simulation outputs at the end of the simulations is due to the assumptions underpinning the theoretical model (see text for details).

  This assumption can be verified with the second subgroup of simulated rectangular parallelepipeds (50 × 103 × 176, 75 × 150 × 75, 150 × 75 × 75, 90 × 90 × 180, and 120 × 120 × 120). Indeed, the only simulations presenting the same (001) surface area are those conducted with the 75 × 150 × 75 and 150 × 75 × 75 volumes. These two rectangular parallelepipeds have the same dissolution behavior throughout the entire reaction (Figure

Figure 8 .

 8 Figure 8. (a) Maximum dissolution rate as a function of the progress of the dissolution process for various initial volumes. Each triangle represents a different initial volume. The maximal value corresponds to the largest volume and the minimal to the smallest. (b) Maximum dissolution rate as a function of the initial volume. (c, d) Dissolution rate when the volume of the enstatite grain reaches the initial volume of the smallest simulated rectangular parallelepiped (14 × 29 × 250) as a function of the initial simulated enstatite rectangular parallelepiped volume represented in standard and logarithmic scales, respectively.

  V 0 is the initial volume (m 3 ) and ξ is the reaction progress.548Combining these two relations allows predicting the reaction 549 progress corresponding to the maximal dissolution rate as a 550 function of the initial volume.

× 29 ×

 29 250). The remaining volume is calculated with the number of atoms that are still in the rectangular parallelepiped at each iteration following the volume at the ith iteration; a, b, and c the cell parameters (m); N 0 the number of Si atoms in the enstatite parallelepiped at the beginning of the simulation; and N n , the total number of Si released until the ith iteration. This mass balance is divided by 16 since the enstatite cell is formed by 16 atoms of Si. The attainment of the volume corresponding to the smallest crystal is indicated by the black circles in Figure 5, and the dissolution rates at V = V 14×29×250 (named Si released, V=V 14×29×250 ) as a function of the initial volume are shown in Figure 8c,d.

  volume is known, the number of Si atoms 623 released per iteration can easily be calculated knowing the 624 volume of the enstatite cell and the number of Si atoms it 625 contains (16) rate of Si requires to know the values of 628 Δh i (hkl) at each iteration step i, which requires to know the face-629 specific dissolution rate of enstatite both at steady state and 630 during the transient states, i.e., at the beginning and at the end of 631 the simulations. To the best of our knowledge, while an 632 analytical expression has been developed in Bouissonniéet al.[START_REF] Bouissonnié | Dissolution Anisotropy of Pyroxenes: A Surrogate Model for Steady-State Enstatite Dissolution Resulting from Stochastic Simulations of the Hydrolysis Process[END_REF] 633 for the steady state, a similar expression does not exist for 634 transient phases. To model the dissolution of the rectangular 635 parallelepipeds used in the different simulations, we first performed three additional face-specific simulations (i.e., simulations where only the atoms belonging to the face (hkl) are allowed to react, thus excluding atoms belonging to edges, corners, and all other (h′k′l') faces), with a depth equivalent to half the depth of the volume in a given direction. The evolution of the surface retreats as a function of time for the (100), (010), and (001) faces were then retrieved and implemented in eq 10b to model the dissolution of a rectangular parallelepiped and determine the evolution of the release rate of Si atoms as a function of time. The results of this dissolution model are presented in Figure7(light-colored curves) and compared to the outputs of the simulations conducted with rectangular parallelepipeds (dark-colored curves). The agreement between both models is striking, with a sharp increase of the dissolution rate observed during the first 8000 iterations, followed by a slow decrease of the dissolution rate and a sharper decrease at the end of simulations.4.3. Comparison between Results of the Simulationand Modeled Dissolution.

  Figure 7 shows that the results derived from theoretical dissolution and modeled dissolution are almost identical. With the exception of the very end of the dissolution, the results of the simulations performed with rectangular parallelepipeds and those derived from face-specific dissolution are superimposed. However, contrary to the results of the simulations performed with rectangular parallelepipeds, those derived from face-specific dissolution exhibit a greater standard variation around the mean value of the curves. This

Figure 9

 9 Figure9. Scheme of the mechanism driving the dissolution of Kossel and enstatite crystals. Kossel: Black cubes represent atoms that are not at the mineral reactive surface. Purple, orange, and yellow cubes represent surface-, edge-, and corner-site atoms, respectively. During the dissolution, the atoms with fewer neighbors have a higher probability to be dissolved. This leads to the formation of rounded corners and edges. The atoms at the surface with only one neighbor missing (purple) are particularly difficult to dissolve compared to the others. Enstatite: at the beginning of the dissolution, Mg atoms are dissolved first. This is particularly visible on the (001) face where they exhibit a less populated first coordination sphere, leading to Mg depletion of the surface and the formation of a Si-enriched layer with a thickness not exceeding a few atomic layers. Because silica chains are not connected to each other, the probability of departure of Si atoms is the same on the entire (001) face, making it unlikely to exhibit rounded edges and corners.

  Figures 9 and 10. Also, as mentioned in Section 2.4, the model is not able to reproduce the amorphous silica layer. Then, the interface mainly composed of Si atoms remains "crystalline"

4. 4 .

 4 Comparing Modeled Dissolution Fluxes and Experimental Results. 4.4.1. Comparison Based on Dissolution Experiments Conducted with Powdered vs Single-Face Enstatite.

r

  OS ). 869 Three additional upscaling methods are tested and compared 870 with r OS . 871 The first method consists of considering the dissolution fluxes 872 determined with single-face dissolution experiments: 16 (a) with 873 faces (100), (010), and (001) and (b) with faces (210)

  889outcrops. The values of these parameters were determined in 890 t1

892 dissolution fluxes given in Bouissonniéet al. 16 893

 16 The three methods were applied on 20 × 20 × 100 and 10 × 894 10 × 100 μm 3 rectangular parallelepipeds (considered as a 895 model shape and size for the grains used in Oelkers and 896 t2 Schott 11,50 ), and the results are summarized in

Figure 11 .

 11 Figure 11. Scheme of the method followed to get the dissolution flux of enstatite rectangular parallelepipeds with a size comparable to those used in Oelkers and Schott. 11 (a) Surface retreat simulated as a function of the iterations is used for each considered face (e.g., (100), (010), and (001)). It is divided into three different phases: (1) the transient regime, (2) the steady-state regime, and (3) the terminal regime, corresponding to the consumption of the remaining atoms in the materials. A cut is operated between the second and third phases. (b) A linear regression (dashed red line) is applied, based on the results of the steady-state phase starting from the transition between the first and second phases and ended between the second and third phases. (c) Linear regression is used to increase the size of the parallelepiped to reach the size used by Oelkers and Schott. The third phase is then added at the end of the extended second phase.

  965study suggests that the anisotropic dissolution of silicates may 966 indeed result in a general dissolution pattern consistent with 967 previous measurements of cement hydration, representing 968 relevant new avenues to be investigated in greater details.969 Finally, among the main limitations of our study, the 970 development of amorphous silica-rich surface layers through 971 silica redeposition or surface network reorganization was not 972 considered, as the impact of such layers is suspected to be 973 nonsignificant in the experimental conditions investi-974 gated 3,47,67,68 by Bouissonniéet al.[START_REF] Bouissonnié | The Dissolution Anisotropy of Pyroxenes: Experimental Validation of a Stochastic Dissolution Model Based on Enstatite Dissolution[END_REF] (i.e., very acidic pH and 975 low concentrations of aqueous silica). The implementation in 976 the stochastic dissolution code of backward reactions of Si 977 attachment, as described in, e.g., Zhang and Luẗtge 69 or the 978 consideration of M-O hydrolysis instead of M-O-M 979 hydrolysis would be required to ultimately get closer to 980 conditions operating in natural settings. 5. CONCLUSIONS 981 In this work, we simulated the dissolution of enstatite grains with 982 a stochastic dissolution model developed at the atomic scale. By 983 varying the size of the simulated grains, we showed that the main 984 factor impacting the temporal evolution of the dissolution rate is 985 the aspect ratio. A similar evolution of the dissolution rate 986 through time is observed for grains with the same aspect ratio. 987 Furthermore, different parameters such as the maximum 988 dissolution rate, initial volume, dissolution at a given volume, 989 and reaction progress were linked together through empirical 990 functions. 991 One of the main outcomes of this study was provided by 992 comparing the outputs of the stochastic grain dissolution with a 993 theoretical dissolution model based on the sole contribution of 994 face-specific dissolution, neglecting the contribution of edges 995 and corners. The very good agreement between the two 996 methods highlights the modest effects of edges and corners 997 during enstatite dissolution, and probably, during the dis-998 solution of other very anisotropic minerals as well (i.e., minerals 999 where the stronger atomic bonds are aligned in a specific 1000 direction). 1001 The comparison between modeled and experimental results 1002 obtained by Oelkers and Schott 11 on powder dissolution also 1003 revealed an excellent agreement. This led to the conclusion that 1004 powder dissolution results from the specific contribution of each 1005 face in a very simple way, making it possible to link the different 1006 experimental works performed in a laboratory. 1007 Finally, this study shows the importance of fast dissolving 1008 faces. Indeed, almost the entire dissolution rate relies on the 1009 dissolution of the face (001). However, if the dissolution rates of 1010 the fastest dissolving faces are prominent, the overall dissolution 1011 rate is strongly correlated to the surface area of the other slower 1012 faces. This result supports the idea that reactive surface is one of 1013 the most important parameters in water/mineral interactions 1014 studies and that one must be aware, particularly in reactive-1015 transport studies, of what is the reactive surface and how it 1016 should be taken into account. 1017 ■ AUTHOR INFORMATION

Table 2 .

 2 Overall,

	897

Table 1 .

 1 Numerical Values of the Parameters in Equation12and Taken from Bouissonniéet al.41 a

	face	(100)	(010)	(001)
	k (Å/it)	8.87	7.8	4.72
	α	2.46	2.22	6.36
	β	4.62	4.65	4.19
	γ	0.94	0.75	0.14
	k dis (Å/it)	44.7	44.7*
	α dis	2.09	2.09*
	β dis	4.09	4.09*
	γ dis	0.67	0.67*
	ω	5.51 × 10 10	5.51 × 10 10 *
	ρ d (/cm 2 )	4.9 × 10 8	9.8 × 10 6	0

a α, β, and γ represent fitting parameters for bulk dissolution rate (without the presence of dislocation). α dis , β dis , and γ dis represent fitting parameters for dissolution rate impacted by the presence of dislocations. The parameters k and k dis represent the dissolution constant for bulk dissolution rate and dissolution rate impacted by the dislocation, respectively. ω is a fitting parameter, and ρ d is the dislocation density measured 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]

(*) Parameters corresponding to the presence of dislocation on the (010) face were not studied. For this study, it has been assumed that they were equal to those of the (100) face.

Table 2 .

 2 Numerical Values of the Steady-State Dissolution Flux of an Enstatite Grain Considering the Simple Relation between Face-Specific Dissolution Flux and Powder Dissolution Flux Following Several Methods a

	method	crystal size (m 3 )	dissolution flux (r = mol/(m 2 s)) (×10 -8 ) r/r OS
	surface retreat	20 × 20 × 100 μm 3	1.78	0.76
	of (100), (001) 16 (010), and	10 × 10 × 100 μm 3	1.18	0.50
	surface retreat	20 × 20 × 100 μm 3	2.75	1.17
	of (210) and (001) 16	10 × 10 × 100 μm 3	2.17	0.92
	surrogate model 41	20 × 20 × 100 μm 3 10 × 10 × 100 μm 3	2.22 1.62	0.94 0.69
	this study	20 × 20 × 100 μm 3	2.62	1.11
	single-face-	20 × 20 × 100 μm 3	1.88	0.8
	modeled results 16	10 × 10 × 100 μm 3	1.35	0.57

a The agreement (evaluated by the r/r OS ratio) between powder dissolution experiments of Oelkers and Schott 11 and results obtained from face-specific dissolution experiments and modeling taken from Bouissonniéet al.

https://doi.org/10.1021/acs.jpcc.0c11416 J. Phys. Chem. C XXXX, XXX, XXX-XXX

https://doi.org/10.1021/acs.jpcc.0c11416 J. Phys. Chem. C XXXX, XXX, XXX-XXX L

The first important highlight is that no real steady state is 852 reached for dissolution rate/flux (Figures 5b and6b). Indeed, 853 the dissolution rate and flux steadily increase over the course of 854 the simulation. This may be due to the evolution of the surface 855 area that is not correctly captured, as the calculation does not 856 consider its intrinsic variability over the entire surface (intrinsic 857 reactivity and etch pits formation for example). However, the 858 increase is slow and not significant (below 14 nmol/s, which