Luca Stigliano
email: luca.stigliano@mail.polimi.it

Martina Siena

Philippe Ackerer

Alberto Guadagnini

Damien Daval

ACS Paragon Plus Environment ACS Earth and Space Chemistry

Keywords: Calcite, Near equilibrium dissolution, Surface roughness, Spatial statistics, Dissolution models

Linking the evolution of the surface area (as quantified, e.g., through its spatial roughness) of minerals to their dissolution rate is a key aspect of mineral reactivity. Unraveling the nature of its main features requires relying on approaches yielding a quantitative characterization of the temporal evolution of surface topography/roughness. Here, a mechanically-polished {104} calcite surface was dissolved at room temperature and at close-to-equilibrium conditions (Ω=0.6) with an alkaline solution (pH = 8) across a temporal window of 8 days. Surface topography images were acquired daily using vertical scanning interferometry, the ensuing topography data being then embedded within a statistical analysis framework aimed at describing comprehensively the surface roughness evolution. The strongest system variations were observed after one day: the probability density function of surface roughness was observed to transition from being approximately Gaussian to being left-skewed and lepto-kurtic, exhibiting a dramatic increase in the variance and a significant change in the semi-variogram structure. After a relaxation time of approximately two days, the reacting surface appeared to attain a steady-state configuration, being characterized by values of the statistical moments characterizing surface roughness which become virtually independent of time. Attempting to unravel the underlying dissolution mechanism, an original numerical model able to reproduce satisfactorily the statistical behavior observed experimentally was developed and tested. Our results suggest that, under the investigated conditions, dissolution may be characterized as a spatially correlated random process, with the areas most exposed to the flowing fluid being prone to preferential dissolution. The numerical model was also used to obtain insights on the influence of the initial surface roughness and of the fluid composition on the steady-state statistical characterization of the surface roughness. Our results suggest that the influence of the initial surface roughness is limited. The present study suggests that potential empirical relations linking the surface roughness of reacted crystals to the saturation state at which they Page 2 of 42 ACS Paragon Plus Environment ACS Earth and Space Chemistry dissolved may be developed, which would allow to back-estimate the reacting conditions only based on topography data.

Introduction

Chemical weathering is a central part of the Earth system, as it is one of the key processes that contributes to the redistribution of mass on Earth. The dissolution of rock-forming minerals controls the formation of soil and supplies nutrients to sustain ecosystems [START_REF] Beaulieu | Response of a forested catchment over the last 25 years to past acid deposition assessed by biogeochemical cycle modeling[END_REF] , it is responsible for atmospheric CO 2 drawdown, regulating the long-term climate of the Earth 2 , and is involved in various geo-engineered strategies such as enhanced rock weathering for carbon dioxide removal from the atmosphere [START_REF] Beerling | Potential for large-scale CO2 removal via enhanced rock weathering with croplands[END_REF] , radioactive waste storage [START_REF] Frugier | SON68 nuclear glass dissolution kinetics: Current state of knowledge and basis of the new GRAAL model[END_REF] , geothermal energy recovery [START_REF] Alt-Epping | Reactive transport modeling of the geothermal system at Bad Blumau, Austria: Implications of the combined extraction of heat and CO2[END_REF] , or CO 2 sequestration [START_REF] Daval | Carbon dioxide sequestration through silicate degradation and carbon mineralisation: promises and uncertainties[END_REF] . The development of rate laws describing the dissolution kinetics of minerals is critical to an effective strategy for upscaling of mineral reactivity, with the ultimate goal of modeling fluxes associated with chemical weathering across a variety of space and time scales.

The classical hypothesis underlying the rate laws that are employed to interpret mineral reactivity (through modeling frameworks based on either a bottom-up (e.g., [START_REF] Beaulieu | High sensitivity of the continental-weathering carbon dioxide sink to future climate change[END_REF] ) or a top-down (e.g., [START_REF] Berner | A MODEL FOR ATMOSPHERIC CO2 OVER PHANEROZOIC TIME[END_REF] ) approach) rests on a conceptual picture according to which the dissolution flux is proportional to the contact surface area between the fluid and the solid.

While such a simplification rests on the assumption that reactivity is immutable in time and homogeneous at the surface of the minerals, the heterogeneous nature of the surface reactivity appears as a universal feature of solid dissolution, as documented through a variety of materials degradation studies ranging from the corrosion of metals [START_REF] Williams | Stochastic models of pitting corrosion of stainless steels: I. Modeling of the initiation and growth of pits at constant potential[END_REF][START_REF] Ryan | Why stainless steel corrodes[END_REF][START_REF] Mundhenk | Corrosion and scaling as interrelated phenomena in an operating geothermal power plant[END_REF] , glasses [START_REF] Fisk | Experimental abiotic alteration of igneous and manufactured glasses[END_REF] , and ceramics [START_REF] Horlait | Environmental SEM monitoring of Ce1-xLnxO2-x/2 mixed-oxide microstructural evolution during dissolution[END_REF][START_REF] Nicoleau | Analytical model for the alite (C3S) dissolution topography[END_REF][START_REF] Fischer | Direct Measurement of Surface Dissolution Rates in Potential Nuclear Waste Forms: The Example of Pyrochlore[END_REF][START_REF] Robin | Experimental study and numerical simulation of the dissolution anisotropy of tricalcium silicate[END_REF] to the weathering of rocks and minerals (e.g., [START_REF] Gratz | Atomic Force Microscopy of Atomic-Scale Ledges and Etch Pits Formed During Dissolution of Quartz[END_REF][START_REF] Arvidson | Variation in calcite dissolution rates[END_REF][START_REF] Fischer | Converged Surface Roughness Parameters -A New Tool to Quantify Rock Surface Morphology and Reactive Alteration[END_REF][START_REF] Fischer | How predictable are dissolution rates of crystalline material?[END_REF][START_REF] Noiriel | Changes in reactive surface area during limestone dissolution: An experimental and modelling study[END_REF][START_REF] Fischer | Variability of crystal surface reactivity: What do we know?[END_REF][START_REF] Godinho | Importance of surface structure on dissolution of fluorite: Implications for surface dynamics and dissolution rates[END_REF][START_REF] Emmanuel | Short Communication: Evidence for non-Gaussian distribution of rock weathering rates[END_REF][START_REF] Bouissonnié | From mixed flow reactor to column experiments and modeling: Upscaling of calcite dissolution rate[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][START_REF] Noiriel | Examination of crystal dissolution in 3D: A way to reconcile dissolution rates in the laboratory?[END_REF][START_REF] Kurganskaya | Mineral Dissolution Kinetics: Pathways to Equilibrium[END_REF] ). At the heart of the heterogeneous nature of surface reactivity is the concept of surface energy, and its distribution in space and time throughout the evolution of the dissolution process (e.g., [START_REF] Fischer | Direct Measurement of Surface Dissolution Rates in Potential Nuclear Waste Forms: The Example of Pyrochlore[END_REF][START_REF] Perez | Comparing the reactivity of glasses with their crystalline equivalents: The case study of plagioclase feldspar[END_REF], which is then mirrored by the gradual development of surface roughness. The heterogeneous distribution of structural defects and the seemingly random occurrence of hydrolysis events at the mineral surface drive variability of surface reactivity, the latter possibly varying across two orders of magnitude for a single mineral (e.g., [START_REF] Fischer | Direct Measurement of Surface Dissolution Rates in Potential Nuclear Waste Forms: The Example of Pyrochlore[END_REF][START_REF] Perez | Comparing the reactivity of glasses with their crystalline equivalents: The case study of plagioclase feldspar[END_REF] ). These observations thus support the need and opportunity for the development of an alternative approach for the assessment of mineral reactivity focused on the heterogeneous space-time dynamics of surface roughness.

As an example, the rate spectra approach [START_REF] Fischer | How predictable are dissolution rates of crystalline material?[END_REF] , is based on the evaluation of the probability density function of the dissolution rates at each given location of the mineral surface. From a practical standpoint, spatially-resolved dissolution rates at the surface of a given material can be quantified upon monitoring the time-resolved evolution of surface topography resulting from dissolution (e.g., [START_REF] Fischer | Direct Measurement of Surface Dissolution Rates in Potential Nuclear Waste Forms: The Example of Pyrochlore[END_REF][START_REF] Fischer | How predictable are dissolution rates of crystalline material?[END_REF][START_REF] Fischer | Variability of crystal surface reactivity: What do we know?[END_REF][START_REF] Noiriel | Examination of crystal dissolution in 3D: A way to reconcile dissolution rates in the laboratory?[END_REF][START_REF] Perez | Comparing the reactivity of glasses with their crystalline equivalents: The case study of plagioclase feldspar[END_REF][START_REF] Emmanuel | Mechanisms influencing micron and nanometer-scale reaction rate patterns during dolostone dissolution[END_REF][START_REF] Fischer | Inherited control of crystal surface reactivity[END_REF][START_REF] Trindade Pedrosa | A Statistical Approach for Analysis of Dissolution Rates Including Surface Morphology[END_REF][START_REF] Wild | Direct measurement of fungal contribution to silicate weathering rates in soil[END_REF] ). One of the most promising aspects of the rate spectra approach is that the deconvolution of a rate spectrum provides information about the modes (i.e., the main contributors to the overall rates), which can then be compared against outputs of atomic-scale simulations such as kinetic Monte Carlo modeling, and refine our knowledge of the corresponding reaction mechanisms (e.g., [START_REF] Fischer | Beyond the conventional understanding of water-rock reactivity[END_REF][START_REF] Luttge | Kinetic concepts for quantitative prediction of fluid-solid interactions[END_REF] ).

A notable limitation of the rate spectra approach is that determining the distribution of dissolution rates at the mineral surface requires analyzing across time a laterally fixed and motionless surface. Such measurements can only be achieved under well controlled laboratory conditions, using techniques such as in situ atomic force microscopy (AFM) conducted in a fluid reaction cell (e.g., [START_REF] Emmanuel | Mechanisms influencing micron and nanometer-scale reaction rate patterns during dolostone dissolution[END_REF] ) or vertical scanning interferometry (VSI) with a reaction cell fixed to the VSI stage (e.g., [START_REF] Trindade Pedrosa | A Statistical Approach for Analysis of Dissolution Rates Including Surface Morphology[END_REF][START_REF] Arvidson | Mineral dissolution kinetics as a function of distance from equilibrium -New experimental results[END_REF] and references therein), so that a given dissolution rate can be assigned to each pixel (or group of pixels) depicting the topography data collected for a given field of view. When considering crystals subject to weathering in the field (e.g., [START_REF] Wild | In-situ dissolution rates of silicate minerals and associated bacterial communities in the critical zone (Strengbach catchment, France)[END_REF] ) or samples with a low reactivity (which makes them incompatible with measurements of the dissolution progress at room temperature over short time periods), obtaining insights into the temporal evolution of the key statistics of the overall roughness of a given surface remains of significant interest. As documented by Fischer et al. [START_REF] Fischer | Converged Surface Roughness Parameters -A New Tool to Quantify Rock Surface Morphology and Reactive Alteration[END_REF] , the investigation of surface roughness evolution could provide information about the reaction history across the dissolution process. As such, time-resolved statistical analyses of dissolving mineral surfaces are of critical importance, as they underpin the development of empirical and/or theoretically-supported relationships between surface topography and the extent of reaction, or between the surface topography and the reaction conditions.

Pollet-Villard et al. [START_REF] Pollet-Villard | Influence of etch pit development on the surface area and dissolution kinetics of the orthoclase (001) surface[END_REF] showed that a statistical treatment of etch pit nucleation and growth monitored by AFM could be used to successfully develop a model accounting for the long-term evolution of the (001) surface of a dissolving K-feldspar at far-from-equilibrium conditions.

By tuning the model parameters based on the comparison between the statistics of modeled and measured surface topography at the early stage of dissolution, these authors could forecast the subsequent retreat rate of the surface, as well as evaluate, through a sensitivity analysis, the impact of defect density (i.e., screw dislocation) on K-feldspar dissolution rates. To the best of our knowledge, a detailed analysis combining statistical approaches and process models to characterize the evolution of the surface roughness of a crystal dissolved at close-to-equilibrium conditions, where the nucleation of etch pits is no longer spontaneous, is still lacking. The present study aims at bridging this gap by providing a detailed statistically-based analysis of the topography evolution of calcite surface dissolved at close-to-equilibrium conditions. In their preliminary study, Siena et al. [START_REF] Siena | Generalized Sub-Gaussian Processes: Theory and Application to Hydrogeological and Geochemical Data[END_REF] demonstrated that single snapshots of the spatial distribution of surface topography of calcite under such conditions can be successfully characterized by jointly modeling the probability distribution functions of surface heights and their spatial increments as a Generalized Sub-Gaussian (GSG) random field. Since single snapshots of surface roughness are not sufficient to predict or reconstruct reaction rates during dissolution (see, e.g., [START_REF] Fischer | Variability of crystal surface reactivity: What do we know?[END_REF] ), the present study extends their work by providing a time-resolved analysis of the evolution of surface topography. We then assess the ability of three empirical numerical models of dissolution to reproduce the observed behavior of the key statistical descriptors associated with the experimental data. The model which is most successful for reproducing the observed statistical description of the system is then employed to assess the relative impact of initial surface roughness and fluid composition on the statistical evolution of the surface. The outputs of the simulations suggest that a bijective relationship might exist between saturation state and the steady-state statistics of surface topography data.

Materials and methods

Materials and experimental methods

Sample preparation

The fragment used in the present study was obtained from an optically transparent calcite sample originating from the Massif des Ecrins (Alps, France). The chemical composition of the sample was determined, as in Bouissonnié et al. [START_REF] Bouissonnié | From mixed flow reactor to column experiments and modeling: Upscaling of calcite dissolution rate[END_REF] , from 3 g of calcite by loss on ignition and alkaline fusion, and only very few impurities (0.2 wt. % MgO, 0.008% MnO, and trace amounts of Sr) were revealed. The calcite sample was cleaved following the natural {104} calcite plane to obtain a 35 mm² size fragment (Figure 1). The crystal was then embedded in epoxy resin to be subjected to a multi-step abrasive grinding (the finest grit size being 8.4 μm) to establish an initially flat surface and to remove the effects of sample cutting, followed by a polishing sequence (the finest grit size being 0.25 μm), using an aqueous solution saturated with respect to calcite as lubricant to avoid etching of the surface. More details are provided as Supporting Information (Table S1). After polishing, the crystallographic orientation of the sample was verified over the entire surface using electron backscatter diffraction (EBSD) on a TESCAN VEGA II scanning electron microscope (SEM). The sample was analyzed as it was (i.e., without carbon coating) and no grain boundaries were observed. The resin was then removed using a precision saw and the crystal surface was cleaned with ethanol. Finally, a small portion (~ 1-2 mm²) of the surface of the crystal was masked with RTV glue spots (Figure 1) to provide a nonreacting reference surface to enable one to measure the height difference (i.e., surface retreat) between the reference and the unmasked portions using VSI (see Section 3.1). The average roughness (R a ), for a given observation window, is defined as:

𝑅 𝑎 = 1 𝑁 𝑁 ∑ 𝑖 = 1 |𝑍 𝑖 (𝒙) -< 𝑍 > | (1)
where, Z(x) corresponds to the height (evaluated with respect to a given reference) of the cell whose centroid is associated with vector position x, is the average height (evaluated < 𝑍 > with respect to the same reference) and N is the total number of sampled height data. Based on VSI data (50 magnification) collected over 100100 μm² areas and following a second-order degree surface de-trending, which was required to minimize biases stemming from the tilt caused by the VSI acquisition and from the waviness of the sample resulting from polishing (see Section 2.1.5), values of R a were found to range between 1.0 and 1.3 nm.

Figure 1

Intensity map of the polished 35 mm² size fragment employed in the experiment as imaged by VSI at 5 magnification (corresponding to a lateral resolution dl = 2.2 μm) after 188 hours of dissolution. The orange area indicates where the mask was placed during the experiment before its removal. The highlighted rectangle corresponds to the region where VSI data were acquired also at 50 magnification (corresponding to a lateral resolution dl = 0.22 μm; see also Figure 2).

Aqueous solution preparation

Reagent-grade NaCl, NaHCO 3 , and CaCl 2 were added to ultra-pure water equilibrated with atmospheric CO 2 (18MΩ•cm) at concentrations similar to those documented in previous the studies ( 40 and 25 ) to provide a baseline for the analysis of rate data. The concentration of NaCl was set to 5.85 g/L to yield the ionic strength to approximately 0.1 molal, and that of NaHCO 3 was set to 0.095 g/L to fix the alkalinity. The measured pH of the resulting solution was 8.0 ± 0.1. The concentration of CaCl 2 was adjusted to maintain a constant saturation index of the solution with respect to calcite ( ), defined as:

Ω 𝑐𝑎𝑙𝑐𝑖𝑡𝑒 Ω 𝑐𝑎𝑙𝑐𝑖𝑡𝑒 = 𝑎 𝐶𝑎 2 + 𝑎 𝐶𝑂 2 - 3 𝐾 𝑠𝑝,𝑐𝑎𝑙𝑐𝑖𝑡𝑒 (2)
where and are calcium and carbonate ion activities , respectively, and is

𝑎 𝐶𝑎 2 + 𝑎 𝐶𝑂 2 - 3 𝐾 𝑠𝑝,𝑐𝑎𝑙𝑐𝑖𝑡𝑒
the solubility product of calcite. Since the nucleation of etch pits on the {104} calcite faces is no longer spontaneous for  > 0.45 [START_REF] Bouissonnié | From mixed flow reactor to column experiments and modeling: Upscaling of calcite dissolution rate[END_REF] , the experiment was conducted at a saturation index of 0.60 ([Ca 2+ ] = 0.00135 molal).

Experimental approach

The experimental set-up was analogous to the one used in Bouissonnié et al. [START_REF] Bouissonnié | From mixed flow reactor to column experiments and modeling: Upscaling of calcite dissolution rate[END_REF] . In brief, the experiment was conducted in a mixed-flow reactor set-up at room temperature (T = 22 ± 1 °C). A calcite fragment of size ~ 10 mm 3 (which is small in comparison with the volume of the fluid in the reactor, i.e., ~ 100 mL) has been used to maintain a negligible concentration of Ca released throughout calcite dissolution as compared to the background concentration in the inlet solution, thus maintaining a nearly constant saturation index across the reaction progress. The experiment was conducted at a constant flow rate (0.25 mL/min) for a total duration of 188 h (~ 8 days). The peristaltic pump was periodically stopped (every ~ 23 hours) to recover the calcite fragment, clean it with ethanol, and analyze the evolution of the surface roughness as a function of time through VSI. After each analysis, the sample was placed back into the reactor and the pump was restarted for additional 23 hour cycles until the experiment was ultimately stopped. The aqueous solution was also regularly sampled (once a day) to evaluate pH and Ca concentration, the latter being measured with an inductively coupled plasma atomic emission spectroscopy (ICP-AES, Thermo iCAP 6000 Series). Based on these fluid analyses, the saturation index was evaluated to be 0.58 ± 0.01 over the duration of the experiment. ICP-AES analyses are given in the Supporting Information (Table S2).

VSI analyses

As emphasized by Fischer et al. [START_REF] Fischer | Converged Surface Roughness Parameters -A New Tool to Quantify Rock Surface Morphology and Reactive Alteration[END_REF] , a robust statistical analysis of surface roughness requires relying on an experimental method that yields a sufficiently high vertical resolution (similar to that provided by AFM) together with a large field of view (e.g., comparable to that provided by light optical microscopy). The technical specifications of VSI, for which the vertical resolution can range between 0.5 and 2 nm and fields of view can be as large as several tens of mm² in stitching mode, are fully compatible with this goal. Measurements of surface topography, Z(x, y), (x, y) being spatial coordinates in the horizontal plane, were collected using a Zygo NewView 7300 VSI. As roughness is a scale-dependent quantity [START_REF] Fischer | Converged Surface Roughness Parameters -A New Tool to Quantify Rock Surface Morphology and Reactive Alteration[END_REF] , all data were acquired and processed at a fixed magnification (50) over a collection of 9 regions (see Figure 2), each associated with the same planar area and with lateral resolution dl = 0.22 μm. In addition, topography measurements acquired at 5 across the whole crystal surface were used to evaluate the overall dissolution rate of calcite over the duration of the experiment according to (see e.g., [START_REF] Lüttge | An interferometric study of the dissolution kinetics of anorthite; the role of reactive surface area[END_REF] ):

𝑟 = ∆ℎ ∆𝑡 𝑉 𝑚 (3)
where r (mol/m²/s) is the overall dissolution rate, Δh (m) is the spatially-averaged height difference (i.e., surface retreat) between the unreacted reference surface and the reacted mineral surface, Δt (s) is the time interval over which the height difference is evaluated, and 𝑉 𝑚 (m 3 /mol) is the molar volume of calcite.

Topography data processing

The analysis of topography data is focused on the 710325 μm² area delineated in Figure 2. This region has been identified as the most suitable for this study, as it has the lowest density of surface defects (such as macro-steps or cracks), which may lead to optical artifacts that would strongly affect the statistical analyses. Thus, the only features characterizing the initial surface in the selected region were polishing scratches. Nine sub-regions of size 3333 μm² have been selected to provide a representative statistical description of the roughness of this area while attempting to avoid most of the VSI-related artifacts (e.g., artificial spikes arising as a result from a sharp (sub-vertical) height transition between two adjacent locations), which would lead to an overestimation of the average surface roughness parameter R a [START_REF] Spencer | The influence of AFM and VSI techniques on the accurate calculation of tribological surface roughness parameters[END_REF] . As such, for all 9 sampled regions, and at all available times, the total amount of 'dead pixels' (i.e., pixels for which the VSI was not able to acquire a height value) and of VSI-related artifacts was negligible (< 0.1%) as compared to the total number of data. Each of these 'dead pixels' and VSI-related artifacts was replaced by the arithmetic mean of the height values of the 8 neighboring pixels.

If one (or more) neighboring pixel was either a 'dead pixel' or a VSI-related artifact, then these were excluded from the calculation of the mean value. The following two criteria have been considered for the selection of the sample window size: (i) convergence of the average surface roughness parameter R a [START_REF] Fischer | Converged Surface Roughness Parameters -A New Tool to Quantify Rock Surface Morphology and Reactive Alteration[END_REF] , which was observed to remain approximately constant within windows ranging from 2525 to 3535 μm²; and (ii) the presence of a sufficient number of data for a statistical analysis.

We note that the raw topography data in the sampled regions could not be analyzed as acquired. Identification of the actual surface roughness required separating the effects of (i) the waviness generated by the polishing process and (ii) the tilt of the sample support applied during the acquisition phase to adjust the alignment between the sample surface and the VSI microscope objective. For this purpose, the topography data of each sample window at each observation time have been de-trended by removing the best-fitting second order degree surface from all data points (Figure 3). Otherwise, polishing scratches have been considered as part of the initial surface roughness (and, hence, were not removed with data processing) to be able to observe their influence on the surface roughness evolution.

Figure 2

Intensity map of the portion of the crystal surface imaged at 50 magnification (depicted in Figure 1). The 710325 μm² area delineated in the figure (orange color box) represents the region identified as most suitable for the study and contains the 9 square subregions (3333 μm²) investigated across time during dissolution. The initial average surface roughness was evaluated (at 50 magnification) to range between 0.6 and 0.8 nm across the 9 sampling windows.

Figure 3

Exemplary three-dimensional (3D) representation of the mean-removed topography data, , (see Section 2.2) associated with a selected sample window, after replacement of 'dead 𝑍 ′ pixels' and VSI-related artifacts (see Section 2.1.5), and second order degree de-trending at observation times t = 0, 23, and 188 hours.

Statistical approaches

A statistical characterization confined to the univariate analysis of the frequency distribution of heights shadows the information associated with the degree of spatial correlation of the data. We rely on statistical tools to characterize the degree of spatial correlation of surface topography by (i) interpreting the investigated random field, , as the sum of an ensemble 𝑍(𝒙) mean,

, and a zero-mean fluctuation, and (ii) evaluating key statistical features of < 𝑍 > 𝑍 ′ (𝒙); both and spatial increments, , h being separation distance (or

𝑍′(𝒙) 𝛥𝑍(𝒉) = 𝑍(𝒙) -𝑍(𝒙 + 𝒉) lag).
We characterize the temporal evolution of sample probability density functions of 𝑍′ (Figure 4) through the time-resolved analysis of the corresponding sample variance, skewness (Equation 4), and kurtosis (Equation 5). We quantify the temporal variation of the degree of spatial correlation of the roughness field through the analysis of the experimental semivariogram (Equation 6).

We recall that the skewness and the kurtosis of a probability density function are defined as:

skew = ∑ 𝑀 𝑖 ( 𝑧 𝑖 -< 𝑍 > ) 3 (𝑀 -1) 𝑠 3 (4) kurt = ∑ 𝑀 𝑖 ( 𝑧 𝑖 -< 𝑍 > ) 4 (𝑀 -1) 𝑠 4 (5)
where is the i-th observed value of , is the ensemble average, is the number of 𝑧 𝑖 𝑍 < 𝑍 > 𝑀 available observations of , and is the standard deviation of . We recall that skewness is a 𝑍 𝑠 𝑍 measure of the degree of asymmetry of a probability distribution, whereas kurtosis imbues the combined weight of the tails relative to the rest of the distribution. We recall that a Gaussian distribution is characterized by a skewness and kurtosis equal to 0 and 3, respectively. A negative skewness indicates a left skewed distribution; a kurtosis lower or higher than 3 denotes a distribution with tails which are respectively thinner or heavier than Gaussian.

The experimental semi-variogram is defined as half of the sample variance of the spatial increments, i.e.,:

𝛾(ℎ) = 1 2 𝑛(ℎ) 𝑛(ℎ) ∑ 𝑖 = 1 [𝑍(𝑥 𝑖 ) -𝑍(𝑥 𝑖 + ℎ)] 2 (6)
where is the number of data pairs separated by a given omni-directional lag . As 𝑛(ℎ) ℎ = |𝒉| stated above, it is a measure of the spatial dependency of Z.

Numerical modeling

Numerical approach

In addition to providing a statistical description of the evolution of the surface roughness resulting from dissolution, we also aim at comparing such depiction to results associated with three dissolution models, as described in the following. The objective of the development of these models is two-fold: (i) demonstrating that a statistical characterization of surface roughness evolution may help to rule out some of the existing models describing the mechanisms of mineral dissolution; (ii) providing a qualitative assessment of the impact of input parameters such as initial surface topography and fluid composition on the calculated statistical descriptors.

Each of the three dissolution models we analyze is characterized by a given mathematical rendering of dissolution processes, which are viewed in a probabilistic context. The parameters of each of these models (see Section 2.3.2) were adjusted to assess whether spatial statistics based on modeling results were compatible with the corresponding statistics based on the experimentally-observed temporal evolution of calcite roughness. The models considered differ in the way one assigns the probability that dissolution takes place at a given location on the simulated surface, as described in the following (additional details being provided in Section 2.3.2):

(i) Uniform and spatially independent dissolution probability (hereafter termed Model 1): the occurrence of dissolution at each location of the simulated crystal surface and at each time is associated with a given probability, the latter being uniformly distributed across the system, independent of location. Initiation of dissolution at a given location is therefore unpredictable and independent of the neighborhood. Such a conceptual framework is typically common to pitting corrosion models [START_REF] Williams | Stochastic models of pitting corrosion of stainless steels: I. Modeling of the initiation and growth of pits at constant potential[END_REF][START_REF] Ryan | Why stainless steel corrodes[END_REF] .

(ii) Equiprobable dissolution with presence of volume defects dissolving preferentially (hereafter termed Model 2). At close-to-equilibrium conditions, the strain field of a dislocation does not open up hollow cores to form etch pits (e.g., [START_REF] Lasaga | Variation of Crystal Dissolution Rate Based on a Dissolution Stepwave Model[END_REF] ) and dissolution is driven by point defects, pre-existing edges, and corners [START_REF] Arvidson | Mineral dissolution kinetics as a function of distance from equilibrium -New experimental results[END_REF] . As such, whereas the dissolution taking place at defect-free zones is considered to be equally likely, locations termed as 'defects' are associated with a higher probability to dissolve, this being intended to simulate such highly localized point defects.

(iii) Dissolution probability of individual locations in the (x, y) plane controlled by fluid accessibility (hereafter termed Model 3). The model is based on a conceptual picture according to which dissolution is locally proportional to the contact surface between the crystal and the fluid [START_REF] Jeschke | Dissolution rates of minerals and their relation to surface morphology[END_REF] , thus rendering the sites that are more accessible to the fluid more likely to be dissolved.

Details of the Numerical models

For all models, the simulated crystal can be interpreted as a three-dimensional collection of voxels, each with a lateral length equal to the VSI lateral resolution (dl = 0.22 μm), their vertical size being assigned a value of dv = 2 nm, consistent with the VSI vertical resolution range (0.5 -2 nm; note that we selected the largest value to reduce computational time).

To ensure that the outputs of the numerical simulations were representative of the experimental scenarios, several settings were considered for all models: (i) the size of the simulation window was taken equal to the size of a single experimental sample window (150150 in terms of pixels, corresponding to a crystal surface area of 3333 μm²); (ii) model evaluations were stopped when the simulated surface retreat, averaged over the whole surface, was equal to its experimental counterpart at the end of the experiment (at 5 magnification); and (iii) the initial topography employed in the simulations coincided with one of the experimental sample windows at t = 0 (as the sample windows were characterized by similar initial average surface roughness R a , as described in Section 2.1.5). Finally, the outputs of the simulations were compared with the experimental data in terms of the statistics of and of the 𝑍′ associated semi-variogram function to identify the numerical model which was most compatible with the statistical results stemming from the experimental observations.

The simulation of the dissolution process at a given time is based on the generation of a random number sampled from a uniform distribution at each (x, y) location of the simulated solid matrix. Such number was then compared to the probability value assigned to the corresponding voxel: if the generated number was lower than the corresponding probability, the voxel was released from the crystal volume; otherwise, it remained in place.

Model 1 involved the presence of 5 types of surface voxels assigned at the beginning of the simulation, each one associated with a given probability of dissolution, being randomly distributed over the entire simulated crystal volume (i.e., the frequency of each voxel type was 20% and they were distributed without any spatial dependency; see Figure 5a).

Model 2 was initiated upon assigning to the whole crystal volume the same probability of dissolution (P min ). Then, a given number of voxel-defects (characterized by a given probability P max > P min ) were randomly and uniformly distributed over the whole crystal volume. Each of these defects was associated with a degree of spatial influence, i.e., all neighboring voxels up to a given distance (L, which is intended as an adjustable model parameter) were characterized by a probability of dissolution (P min < P < P max , P decreasing linearly with distance from the voxel-defect, until the non-disturbed region was reached; see Figure 5b for an example of a voxel-defect).

Model 3 was characterized by 2 parameters (here denoted as P 1 and ), and is based on 𝑃 * 1 the concept that the voxels having the largest number of faces exposed to the fluid are the most reactive ones. All voxels having only their top faces exposed to the fluid are characterized by a probability of dissolution . All voxels having also one or more lateral faces exposed to the 𝑃 * 1 fluid are characterized by a higher probability of being associated with dissolution: the higher the number of faces exposed to the fluid, the higher the probability that dissolution takes place at that voxel. Model parameter P 1 corresponds to the contribution to the probability of dissolution of a single voxel given by a single lateral face exposed to the fluid. Hence, the probability of dissolution linked to a voxel having j lateral faces in contact with the fluid is P 1  j. This model also includes the possibility that various voxels which are superimposed (henceforth termed voxel spines, see Figure 5c) dissolve simultaneously. Consistently, the larger the voxel spines lateral area in contact with the fluid, the higher the probability that one or more voxels are dissolved: for a given voxel within a voxel spine, the probability that the voxel considered is dissolved along with all the upper voxels is P 1 

, where ∑ to dissolution (and hence are released from the solid matrix) at the same time is zero, as it is given by the probability of voxel F (with no lateral faces exposed to the fluid) to be released, 𝑃 𝐹 which is null (note that the configuration of the voxel layers along the y-axis is the same as the one here depicted). Following a similar reasoning:

𝑃 𝐴 + 𝐵 + 𝐶 + 𝐷 + 𝐸 = 𝑃 𝐴 + 𝐵 + 𝐶 + 𝐷 + 𝐸 + 𝐹 + 𝑃 𝐸 ; ; = 0 + 𝑃 1 × 1 𝑃 𝐴 + 𝐵 + 𝐶 + 𝐷 = 𝑃 𝐴 + 𝐵 + 𝐶 + 𝐷 + 𝐸 + 𝑃 𝐷 = 𝑃 1 × 1 + 𝑃 1 × 1 𝑃 𝐴 + 𝐵 + 𝐶 = ; ; 𝑃 𝐴 + 𝐵 + 𝐶 + 𝐷 + 𝑃 𝐶 = 𝑃 1 × 2 + 𝑃 1 × 2 𝑃 𝐴 + 𝐵 = 𝑃 𝐴 + 𝐵 + 𝐶 + 𝑃 𝐵 = 𝑃 1 × 4 + 𝑃 1 × 2 𝑃 𝐴 = 𝑃 1 . × 6 + 𝑃 1 × 2 + 𝑃 * 1

Results

Calcite dissolution rate

The average surface retreat (Δh) measured after 188 hours of dissolution at Ω = 0.60 was determined to be (0.20 ± 0.01) μm (mean value ± standard deviation), which resulted in an estimate of the {104} overall dissolution rate equal to (8.0 ± 0.4) × 10 -9 mol m -2 s -1 . We note that this value is in good agreement with close-to-equilibrium values documented in previous studies ( 40 and 25 ) using similar experimental and analytical protocols (these values ranging from 1.8×10 -8 mol m -2 s -1 at Ω = 0.45 to 5.6×10 -9 mol m -2 s -1 at Ω = 0.80). This imbues us with confidence that the desired close-to-equilibrium thermodynamic conditions were attained in our experimental setting.

Statistical analyses

Time-resolved statistical analysis of surface roughness

Key results of the time-resolved statistical analysis of the mean-removed surface topography (Z') are collected in Figure 6, which depicts the temporal dependence of variance (Fig. 6a), skewness (Fig. 6b), and kurtosis (Fig. 6c) of Z', as obtained upon averaging the corresponding quantities evaluated for each of the 9 samples analyzed. Intervals of width equal to ± one standard deviation are also depicted to provide an appraisal of the variability of such statistical moments among the samples. These plots reveal that temporal variations of the main statistics of the surface topography take place mostly at early times (reaction times shorter than 46 hours). All of the investigated statistical moments are seen to exhibit a plateau after such a relaxation time. Fig. 6a clearly shows a steep increase of the variance for time t < 46 h, until an asymptotic value of (2.85 ± 0.6) × 10 -5 is attained. The sharp variation observed at the µm 2 very beginning of the experiment is possibly related to the initial condition rendered by the polishing process which might have locally weakened the surface. Polishing-resulting scratches are likely to induce the formation of surface defects which are in turn prone to preferential dissolution (see e.g., [START_REF] Perez | Comparing the reactivity of glasses with their crystalline equivalents: The case study of plagioclase feldspar[END_REF][START_REF] Trindade Pedrosa | A Statistical Approach for Analysis of Dissolution Rates Including Surface Morphology[END_REF] ). Attainment of constant values of these statistical moments suggests that the rate at which the surface topography dissolves tends to become spatially constant as time progresses. Analysis of the third and fourth moments (i.e., skewness and kurtosis) of Z' reveal that the initial surface configuration is characterized by an approximately Gaussian probability density function, whereas the surface configuration at late times emerges as leftskewed (skew < 0) and lepto-kurtic (kurt > 3). Fig. 6b shows a transition from a Gaussian pattern (skew ≈ 0) to an asymptotic negative value (skew ≈ 0.2), displaying a minimum value at early times (t ≈ 23h). Similarly, Fig. 6c illustrates a deviation from values typical of a Gaussian distribution at t = 0h to a leptokurtic behavior at later times . This implies that dissolution leads to the formation of more pronounced peaks and valleys as the reactions starts, the valleys being predominant. This behavior may also provide a quantitative criterion to distinguish, at steadystate, between reacted and unreacted (mechanically-polished) calcite surfaces based only on surface topography measurements. Consistent with the observation of Siena et al. [START_REF] Siena | Generalized Sub-Gaussian Processes: Theory and Application to Hydrogeological and Geochemical Data[END_REF] , (mechanically-polished) dissolved calcite surfaces are characterized, at late times, by leftskewed topography distributions. As such, while the probability distribution of the surface roughness across the mechanically-polished unreacted surface is symmetric, reaction is seen to be associated with a larger frequency of low rather than high heights values, implying that dissolution was very active in the majority of the area. As long as initial symmetric unreacted surfaces are concerned, the asymmetry of Z' distributions can then be considered as a distinctive feature of dissolution, being also consistent with the right-skewed rate-spectra observed in several studies related to a variety of minerals altered through various dissolution regimes (e.g., starts, may provide some insights regarding surface reactivity. Respectively: (i) a skewness value that decreases from 0 at t = 0h to a negative value at t = 23h might be the result of a nonuniform distribution of dissolution rates; (ii) a kurtosis value that increase from 3 at t = 0h to a larger value at t = 23h might be related to dissolution rates that are, respectively, extremely lower and extremely higher with respect to the average dissolution rate. Consistent with the results of Section 3.2.1, the most pronounced variations take place at early times. Three regimes can be identified. The first regime (corresponding to 0 h  t  23 h) is characterized by a strong increase of the experimental semi-variogram sill with time, consistent with the sample variance rise (Figure 6a and Figure 7a). The semi-variogram at t = 0h exhibits higher values at shorts lags and a more persistent spatial correlation than its counterpart at t = 23 h (Figure 7b). This feature, combined with the (previously discussed) observation that the distributions of Z' are left-skewed, may be employed to quantitatively distinguish between reacted and unreacted surfaces. Semi-variograms evaluated for 23 h < t  69 h are characterized by a mild temporal increase of their sill (see also Fig. 6a) and a virtually constant level of spatial correlation across time (Figure 7c). Our results show that differences between semi-variograms evaluated at subsequent times tend to diminish as time increases (see also Fig. 7d), consistent with the observation that the system tends to attain a certain degree of steady-state, in terms of the statistical characterization of surface roughness.

Time-resolved semi-variogram analyses

Numerical model selection

As discussed in Section 2.3.1, model parameters are tuned until the best possible agreement with the statistics of the experimental data is attained. The least stringent condition to be met to this end was the constraint on the observed average surface retreat. Model-based results correspond to averages across a collection of 100 simulation runs for each set of parameters tested. Figure 8a depicts the behavior of the variance of surface roughness associated with (i) the experimental observations (and obtained upon averaging results associated with all of the spatial sub-regions investigated) and (ii) their model-based counterparts corresponding to the optimal parameter value sets versus the extent of reaction (i.e., the fraction of total iterations, equivalent to time). Figure 8b shows a corresponding depiction of sample semi-variograms associated with the longest observation time (t = 188 h, corresponding to 100% extent of the reaction).

Results from Model 1 were deemed as not consistent with the documented experimental evidences because the model-based semi-variogram in Figure 8b displayed a lack of spatial correlation, in line with the nature of the underlying physical model. A large number of parameter value combinations have been investigated for Model 2. It was found that, even though the model does imprint a certain degree of correlation to the spatial field of roughness, none of the analyzed parameter sets allowed attaining a variance plateau over time. The latter result is related to the observation that the considered dissolution mechanism implied the variance to continuously increase with time. Otherwise, results from Model 3 are consistent with the experimental observations with reference to both the attainment of the variance plateau at long times (see Figure 8a) as well as the ability to reproduce the semi-variogram structure after 188 hours of dissolution (see Figure 8b). As an additional advantage, we note that Model 3 requires only two model parameters. All of these features contribute to identify Model 3 as the most successful, in terms of consistency with the experimental results, across the collection of models tested.

Parametrization and validation of the selected model

The set of parameters for Model 3 that results in the best match (based on the minimization of MSE, as described in Section 3. This model is able to capture the observed behavior of the variance, as it causes the variance to steeply increase at early times through the generation of areas dissolving preferentially, until a plateau is reached as these areas merge at long times, while the average retreat increases linearly. The implementation of the dissolution mechanism associated with Model 3 imprints a certain degree of spatial correlation due to the assumption that the voxels associated with a larger number of faces exposed to the fluid are more reactive (i.e., more likely to be dissolved).

As a consequence, the simulated dissolution process advances through time by propagating a number of randomly generated valleys (i.e., areas resulting from dissolution being locally more active than in the remaining portions of the surface). The temporal evolution of the simulated surface is probabilistically related to the locations where these valleys are randomly generated at early iterations. As a result, diverse simulation runs may lead to moderately different statistical evolutions of the simulated surface. We investigate this aspect by performing 100 runs of Model 3 with the optimal parameter set. Figure 9 depicts a comparison of Z' statistics between ± standard deviation resulting from 100 simulation runs of Model 3 and ± standard deviation evaluated across the 9 sub-regions investigated.

We observe that the steady-state model-based values of Z' variance are in good agreement with their experimental counterparts (see Figure 9a), whereas mild discrepancies can be observed as the order of the statistical moments considered increases (see Figures 9b and9c).

The main differences in the behavior of the variance of Z' are detected at early times, the model plateau being attained later than the experimental one. Whereas the early time system variation observed experimentally is very sharp (and possibly enhanced by the weakening effect of polishing scratches), numerical results show that merging of areas dissolving preferentially is more gradual. Nevertheless, the main temporal trend of the modeling results is satisfactorily consistent with the experimental observations, thus supporting the selection of Model 3 as capable of embedding the main statistical traits associated with the observed system dynamics.

Discussion

Interest of a statistical approach to surface topography of dissolved mineral surfaces

Attempting to take into account the spatial heterogeneity of mineral dissolution, Fischer et al. [START_REF] Fischer | How predictable are dissolution rates of crystalline material?[END_REF] suggested the rate spectra approach, which was then largely applied in several subsequent works (e.g., [START_REF] Fischer | Direct Measurement of Surface Dissolution Rates in Potential Nuclear Waste Forms: The Example of Pyrochlore[END_REF][START_REF] Fischer | How predictable are dissolution rates of crystalline material?[END_REF][START_REF] Fischer | Variability of crystal surface reactivity: What do we know?[END_REF][START_REF] Noiriel | Examination of crystal dissolution in 3D: A way to reconcile dissolution rates in the laboratory?[END_REF][START_REF] Perez | Comparing the reactivity of glasses with their crystalline equivalents: The case study of plagioclase feldspar[END_REF][START_REF] Emmanuel | Mechanisms influencing micron and nanometer-scale reaction rate patterns during dolostone dissolution[END_REF][START_REF] Fischer | Inherited control of crystal surface reactivity[END_REF][START_REF] Trindade Pedrosa | A Statistical Approach for Analysis of Dissolution Rates Including Surface Morphology[END_REF][START_REF] Wild | Direct measurement of fungal contribution to silicate weathering rates in soil[END_REF] ). This strategy has opened the possibility for new developments in the field of mineral reactivity, as the deconvolution of rate spectra allows for the determination of the main contributors to a given rate distribution and can shed some light on our ability to identify dissolution mechanisms [START_REF] Fischer | Beyond the conventional understanding of water-rock reactivity[END_REF][START_REF] Luttge | Kinetic concepts for quantitative prediction of fluid-solid interactions[END_REF][START_REF] Emmanuel | Limestone weathering rates accelerated by micron-scale grain detachment[END_REF] or the specific contribution of microbes to chemical weathering [START_REF] Wild | Direct measurement of fungal contribution to silicate weathering rates in soil[END_REF] . However, this approach suffers from two main limitations: (i) it requires the presence of a laterally fixed (motionless through time) surface and/or reference points at the crystal surface, which may not be compatible with a number of field applications; and (ii) it does not contain information about the spatial correlation of the surface roughness distribution, as the evolution of dissolved surfaces characterized by diverse extents of spatial correlation might be associated with the same rate spectra.

Our statistically-based approach, targeting close-to-equilibrium conditions (which are characteristic of most natural settings) can be seen as complementary to the rate spectra approach.

Our observations that the initial surface configuration of a mechanically-polished calcite sample is characterized by a Gaussian probability density function, whereas the surface equilibrium configuration is left-skewed and lepto-kurtic, may provide a quantitative criterion to distinguish between reacted and unreacted (mechanically-polished) calcite surfaces based only on surface topography measurements. More specifically, the left-skewness of dissolved calcite surfaces (already observed by Siena et al. [START_REF] Siena | Generalized Sub-Gaussian Processes: Theory and Application to Hydrogeological and Geochemical Data[END_REF] ), which intrinsically implies that the surface is characterized by a larger number of low heights values rather than high heights values, could be viewed as a specific feature of dissolution, as it results in a right-skewed rate spectrum consistent with the observations of numerous studies related to various minerals altered through various dissolution regimes (e.g., [START_REF] Fischer | How predictable are dissolution rates of crystalline material?[END_REF] or 30 ).

Insights on surface dynamics of calcite dissolved at close-to-equilibrium conditions

Our results show that the surface roughness variance increases with time until a plateau is reached, indicating the achievement of a steady-state surface configuration, as supported by the attainment of a virtually invariant semi-variogram. We note that the numerical simulations and scanning tunneling microscopy observations of silver electro dissolution by Hernández et al. [START_REF] Hernández Creus | Three-Dimensional Monte Carlo Simulations of Roughness Development from Different Mechanisms Applicable to the Dissolution of a Pure Solid[END_REF] yielded similar trends. A qualitatively similar behavior was also observed in previous works ( 23 and 47 ) at far-from-equilibrium conditions for fluorite, these authors noting that the value of the variance plateau may be specific to the considered crystallographic orientation.

With the aim of obtaining further insights on the underlying dissolution mechanism, three diverse numerical models have been examined (see Section 2.3). Model 3 (i.e., Dynamic dissolution probability based on fluid accessibility) has been identified as the one that best captures the results of our statistical analyses on the evolution of the surface topography. Model 1 (i.e., Uniform and spatially independent dissolution probability similar to pitting corrosion models) failed to account for spatial correlation, whereas Model 2 led to a monotonous temporal increase of the variance of the surface topography and the formation of geometric etch pits, these features being incompatible with dissolution at close-to-equilibrium conditions.

Otherwise, the results of Model 3 (i.e., Equiprobable dissolution with presence of volume defects dissolving preferentially), by involving a dynamic set of dissolution probabilities, were fully consistent with the main statistics and semi-variogram structure of surface roughness.

Pollet-Villard et al. [START_REF] Pollet-Villard | Influence of etch pit development on the surface area and dissolution kinetics of the orthoclase (001) surface[END_REF] developed a model accounting for the long-term evolution of the (001) surface of a dissolving K-feldspar at far-from-equilibrium conditions, based on the correspondence between modeling and experimentally based statistical behaviors. Similarly, the results of our study enable us to conclude that the dissolution of a {104} calcite surface at close-to-equilibrium conditions progresses in a way which is consistent with a mechanism based on a random generation of a number of valleys, which tend to increase with time, and grow until they finally merge to attain a steady-state surface area. We note that, while the temporal evolution of the roughness variance reaches a plateau, the average retreat increases according to a linear trend, which is also consistent with other works (e.g., [START_REF] Pollet-Villard | Influence of etch pit development on the surface area and dissolution kinetics of the orthoclase (001) surface[END_REF] ).

One should emphasize that our results are not intended to imply that Model 3 can be considered as a universal descriptor of mineral dissolution that would apply to any mineral reacted at close-to-equilibrium conditions. As opposed to mechanistic models such as kinetic Monte Carlo simulations (e.g., [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][START_REF] Kurganskaya | Mineral Dissolution Kinetics: Pathways to Equilibrium[END_REF][START_REF] Kurganskaya | A comprehensive stochastic model of phyllosilicate dissolution: Structure and kinematics of etch pits formed on muscovite basal face[END_REF] ), Model 3 is essentially phenomenological and its parameters need to be estimated against experimental observations. Nonetheless, and more importantly, the development and analysis of these various models underlined the need to consider some key elements that have to be implemented into any future model to ensure its ability to capture major features of the statistics of surface roughness, as described in the following:

(i) The dissolution probability of a given voxel is dynamic, as it depends on its immediate neighboring environment. Interestingly, this feature is also observed at a smaller (i.e., atomic) scale (e.g., [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][START_REF] Kurganskaya | A comprehensive stochastic model of phyllosilicate dissolution: Structure and kinematics of etch pits formed on muscovite basal face[END_REF][START_REF] Cailleteau | Insight into silicate-glass corrosion mechanisms[END_REF][START_REF] Kurganskaya | Kinetic Monte Carlo Approach To Study Carbonate Dissolution[END_REF][START_REF] Bouissonnié | Dissolution Anisotropy of Pyroxenes: Role of Edges and Corners Inferred from Stochastic Simulations of Enstatite Dissolution[END_REF] ), although in this latter case, the relationship that links the dissolution probability of an atom to its environment corresponds to the product of the probabilities to break the bonds linking the atom to the surface (thus corresponding to the hydrolysis of all bonds simultaneously). Applying the exact same rule also at the nm-scale would most likely be unreasonable without the development of upscaling laws, as the detachment of a nm-scale voxel already involves the dissolution of hundreds of thousands of atoms. Here, we show that the detachment of a voxel is related, in a statistical sense, to the number of faces of the voxel exposed to the solution, consistent with a macroscopic description of the dissolution process where the dissolution rate is scaled to the surface area in contact with the solution.

(ii) Specific features that are shown to be consistent with the observation that the variance of the surface roughness reaches a plateau include (a) considering that columns of voxels (i.e., voxel spines) can dissolve at the same time, as well as (b) assigning a higher dissolution probability to voxels located at the top of a given column (see Section 2.3.2). Interestingly, such a condition is similar to that implemented in one of the models developed to describe the dissolution of silver [START_REF] Hernández Creus | Three-Dimensional Monte Carlo Simulations of Roughness Development from Different Mechanisms Applicable to the Dissolution of a Pure Solid[END_REF] , where a given probability is assigned to each voxel depending on its distance to the basal plane.

(iii) Finally, we showed that the strength of the spatial correlation associated with the semivariogram for late times is not necessarily related to a specific distance between preexisting defects in the simulated crystal (as in Model 2 or in Pollet-Villard et al. [START_REF] Pollet-Villard | Influence of etch pit development on the surface area and dissolution kinetics of the orthoclase (001) surface[END_REF] ), and can be otherwise associated with the specific choice of a dynamic set of dissolution probabilities.

Statistical characterizations of surface topography as a metric for identification of reacting conditions?

Dissolution at far-from-equilibrium conditions results in a final topography which is markedly different from the steady-state surface observed in our study at close-to-equilibrium conditions [START_REF] Pollet-Villard | Influence of etch pit development on the surface area and dissolution kinetics of the orthoclase (001) surface[END_REF][START_REF] Beig | Albite dissolution kinetics as a function of distance from equilibrium: Implications for natural feldspar weathering[END_REF][START_REF] Bandstra | Surface evolution of dissolving minerals investigated with a kinetic Ising model[END_REF] . As suggested by Godinho et al. [START_REF] Godinho | Importance of surface structure on dissolution of fluorite: Implications for surface dynamics and dissolution rates[END_REF] , the temporal evolution of surface topography may be associated with decreasing dissolution rates. This opens the possibility that the observed steady-state surface topography is closely related to the fluid composition. The experimental investigation of such conceptual element remains the focus of ongoing research.

The calibrated numerical model developed in this work which has been shown to satisfactorily reproduce the experimentally observed statistical behavior at given conditions of initial roughness and fluid composition (i.e., Model 3), is then employed to explore additional scenarios with the aim of discriminating the contribution to the surface equilibrium configuration given by the initial roughness from that resulting from the fluid composition.

Results obtained considering various initial topographies in the absence of macrofeatures suggest that the initial surface roughness has only a negligible impact on the simulated steady-state surface. Our simulation results obtained upon considering the effect of a random noise associated with the initial surface document that there are no significant variations in the fluid-mineral contact surface (and, hence, in preferential dissolution), whereas considering the presence of surface macro-features (e.g., cracks or fractures) would possibly lead to different results. As an example, Figure 10 a, b depict the temporal evolutions of variance and the steadystate semi-variograms resulting from 4 different initial surface topographies: (i) a perfectly flat surface, (ii) the experimental topography measured at t = 0 hours for one of the sample windows, (iii + iv) the topography associated with (ii) which is then corrupted through a white noise with standard deviation equal to: (iii) 1 standard deviation and (iv) 2 times the standard deviation of the experimental topography considered. If supported by experimental data, these results would enable one to rule out the starting topography as a major parameter controlling the steady-state variance of the surface topography resulting from mineral dissolution. Conversely, differences in "fluid compositions" were shown to have a major impact on the steady-state topography. The impact of fluid composition was simulated based on the observation that and are proportional to the dissolution rates, hence inversely 𝑃 1 𝑃 * 1 proportional to the saturation state. Here, we limit our analysis to conditions where no etch pits would actually be observed, so that only dissolution rates corresponding to saturation states ranging between Ω=0.45 and Ω=0.8 are simulated [START_REF] Bouissonnié | From mixed flow reactor to column experiments and modeling: Upscaling of calcite dissolution rate[END_REF] . As an example, Figure 10c,d model parameters, which resulted respectively in higher (Δh = 0.26 μm) and lower (Δh = 0.14 μm) surface retreats. Model results suggest that higher saturation states (i.e., closer to equilibrium conditions) result in higher values of the variance plateau due to the delayed merging of the areas subject to preferential dissolution. This observation may seem counterintuitive: Fischer et al. [START_REF] Fischer | Variability of crystal surface reactivity: What do we know?[END_REF] showed that kink site density at the atomic scale is positively correlated with an increase in surface roughness, the latter being in turn strongly correlated to the dissolution rate (see their Fig. 6). On the other hand, nm-scale experimental and modeling results of previous relevant studies ( 23 and 47 ) showed that an increase in roughness can be associated with decreasing dissolution rates, which would be in close agreement with our overall analysis. We note that the results from our simulations, albeit constrained to the available experimental settings, suggest that the steady-state variance of the topography and the fluid composition may be related. From a practical aspect, generalizing this result would rest on the assumption that steady-state surface configuration would be achieved also at saturation states differing from the one here investigated: if supported through additional experimental observations, direct measurements of surface roughness would allow to unambiguously distinguish surfaces that have reacted from surfaces that have not, while possibly getting insights into the saturation state of the fluid at which dissolution took place. This justifies the need for further experimental investigation on the relationship between mineral dissolution rate and surface roughness, to enhance the ability of our model to characterize these processes.

While further experimental studies are required to explore the behavior seen in our numerical simulations, our results suggest that the reacting conditions could be back-estimated through a simple statistical characterization of topography data of a reacted surface, paving the way to the development of empirical transfer functions linking saturation states and reacted surface topographies.

Conclusions

Time-resolved statistical analyses of a mechanically polished calcite surface reacting at close-to-equilibrium conditions were demonstrated to be a powerful tool to assess attainment of steady-state configuration for the reacting surface. This is reflected through the temporal dynamics of the key statistics and spatial semi-variograms of surface roughness which tend to stabilize after a given relaxation time (approximately two days in our experiments). Our approach can be used as a quantitative reference to compare (i) differently polished (or unpolished, e.g., cleavage) calcite surfaces reacting at the same physico-chemical conditions, or (ii) mechanically polished calcite surfaces reacting at different physico-chemical conditions.

In both cases, spatial analysis of semi-variogram functions can provide a quantitative appraisal of the strength of the space-time correlation of surface roughness as it evolves subject to the progress of the reaction. In this context, additional investigations may contribute to possibly distinguish between the individual contribution given by the initial surface roughness and the reacting conditions to the plateau values of the statistical moments characterizing the dissolving surface area at steady-state. As a consequence, a thorough characterization and understanding of the statistical evolution of surface roughness could possibly be key to back-estimate the reacting conditions (e.g., saturation state of the solution) by simply measuring the topography of a reacted surface after attainment of steady-state.

We developed and implemented an original conceptual model which enabled us to reproduce the experimentally observed space-time evolution of the key statistical descriptors of surface topography. The documented consistency between modeling-and experimentally-based statistical analyses may imply that, under the investigated conditions, the dissolution mechanism is characterized by an overall space-time organization, as reflected by the semivariogram structure, with the areas most exposed to the fluid being prone to preferential dissolution. The results of our statistical analyses document that such a dissolution mode is fully compatible with the achievement of a steady-state surface area observed through the experiments in this work.

Ongoing work is geared towards analyzing whether this model, while being conceptually and mathematically simple, can be a viable tool to reliably predict the statistical evolution of calcite surface roughness as a function of the fluid composition.

Figure 4

 4 Figure 4 Comparison of the sample distributions probability densities of the variable Z' at observation times t = 0, 23, and 188 hours (open circles) with Gaussian distributions densities (dashed lines) characterized by the same mean and standard deviation of Z'. The distribution of Z' appears to be symmetric at t = 0 hours and slightly asymmetric (left-skewed) at longer times.

4 𝑗 = 1 𝑗

 1  𝑛 𝑗 + 𝑃 * 1 𝑛 𝑗 is the number of lower voxels (the one considered being included) having j faces exposed to the fluid and is nonzero only when considering the top voxel (see Figure 5c for a schematic 𝑃 * 1 example).

Figure 5

 5 Figure 5 (a-b) Conceptual representation of the main features of Model 1 (a) and Model 2 (b) in a horizontal section. In these examples, values of dissolution probabilities are represented as Pmin<P2<P3<P4<Pmax. (c) Schematic representation of a voxel spine in a vertical section. As an example, the probability that voxels A, B, C, D, E, and F are all subject 𝑃 𝐴 + 𝐵 + 𝐶 + 𝐷 + 𝐸 + 𝐹to dissolution (and hence are released from the solid matrix) at the same time is zero, as it is given by the probability of voxel F (with no lateral faces exposed to the fluid) to be released, 𝑃 𝐹 which is null (note that the configuration of the voxel layers along the y-axis is the same as the

Figure 6

 6 Figure 6 Temporal evolution of the key higher-order statistics of Z', i.e., variance (a), skewness (b) and kurtosis (c). The figure depicts mean values (circles) evaluated across the 9 sample windows identified in Figure 2 and intervals of width corresponding to ± one standard deviation (SD).

Figure 7

 7 Figure 7 depicts sample omni-directional semi-variograms of surface roughness evaluated

Figure 7

 7 Figure 7 Evolution of experimental semi-variogram in time for a selected sample window (data associated with the remaining windows exhibit a qualitatively similar behavior). (a) Experimental semi-variogram (log-log scale) at observation time t = 0 h and 23 h; (b) Normal scored transformed variogram (i.e., semi-variogram normalized to the sample variance) at 0 and 23 hours. The main advantage of this visualization is that it makes it possible to compare the structure of semi-variograms characterized by strongly different sample variances; (c) Experimental semi-variograms for 23 h  t  69 h; (d) Experimental semi-variograms for 69 h  t  188 h.

  Optimal parameter values correspond to minimization of the mean square error (MSE) between (i) the average variance values depicted in Fig.6aand their model-based counterparts and (ii) the late-time model-based semi-variogram and the corresponding averaged experimental semi-variogram. We recall that for a stationary random field, as is, the semi-𝑍′(𝒙) variogram sill corresponds to the variance.

Figure 8

 8 Figure 8 (a) Variance of surface roughness (Z') associated with (i) the experimental observations (symbols; ± standard deviation evaluated across the 9 spatial sub-regions investigated) and (ii) their model-based counterparts corresponding to the optimal parameter value sets (solid lines) versus the extent of reaction (i.e., the fraction of total iterations, equivalent to time). (b) Experimental and numerically-based sample semi-variograms associated with the longest observation time (t = 188 h, corresponding to 100% extent of the reaction).

1 and

 1 3.1) between numerical and experimental data is: =0.05 𝑃

Figure 9

 9 Figure 9 Temporal evolution of the key statistics of Z' evaluated over 100 realizations of Model 3. Envelopes of width equal to ± one standard deviation around the corresponding mean are depicted (shaded areas). Symbols denote the corresponding intervals associated with the experimental data monitored across the 9 sample windows investigated.

Figure 10

 10 Figure 10 Temporal evolution of the variance of surface roughness as a function of the initial surface roughness (a) and of the fluid composition (i.e., saturation state) (c) resulting from numerical simulations associated with Model 3. Corresponding results for the normal scored transformed experimental semi-variogram evaluated at the last observation time (b, d) are also included. No significant variation in the semi-variogram was observed in both cases, suggesting that the semi-variogram structure may significantly change only at far from equilibrium conditions regardless of the initial surface roughness.

or[START_REF] Emmanuel | Mechanisms influencing micron and nanometer-scale reaction rate patterns during dolostone dissolution[END_REF] ). Furthermore, the decrease in skewness and the increase in kurtosis as the reaction

Acknowledgements

This work was funded by the project ANR-20-ERC8-0006-01_t erc Mobidic awarded to D.D.

We thank Gilles Morvan and René Boutin for their help with EBSD and ICP-AES measurements. Finally, L.S. thanks Arnaud Bouissonnié for the assistance in the implementation of experimental protocols.

Data Availability

Data are available from the corresponding author upon request.

Supporting Information

The Supporting Information is available at: … (1) Grit sizes of the multi-step grinding and polishing sequence.

(2) ICP-AES analyses of the aqueous solution at different times.

ACS Paragon Plus Environment ACS Earth and Space Chemistry