Process Modeling of Mineral Dissolution From Nano‐Scale Surface Topography Observations

,


Introduction
Chemical weathering is a critical and ubiquitous element of the evolution of the Earth system.It drives a variety of processes with remarkable feedbacks on human life such as, for example, formation of soils (Yu & Hunt, 2018), nutrient supply for ecosystems (Bañuelos et al., 2023, and references therein), and trace metal incorporation in rocks (Heberling et al., 2014;Hövelmann et al., 2018;Renard et al., 2018Renard et al., , 2019)).It is also at the core of various geo-engineered strategies for sustainable development such as carbon dioxide removal and natural underground storage (Daval, 2018;Fitts & Peters, 2013;Jun et al., 2017;Noiriel & Daval, 2017) or geothermal energy exploitation (Erol et al., 2019).Providing robust and physically based fundamental understanding on the kinetics of mineral dissolution is key to effectively quantify and transfer chemical weathering patterns across diverse spatial and temporal scales.These elements provide the core motivation for our study, together with the observation that there is a dire need to pair process model development with highly accurate experimental observations at such scales in a well-controlled environment (as also evidenced by Soulaine (2024)).
Classical interpretive approaches based on the assumption of spatially homogeneous reactivity do not align with the well-documented presence of significant heterogeneity in mineral surface reactivity (Fischer et al., 2014;Lüttge et al., 2013).Evidences comprise detailed micro-and nano-scale direct observations of time-resolved evolution of mineral surface topography resulting from dissolution.These types of experimental data are obtained under controlled laboratory conditions upon relying on techniques such as, for example, in situ Atomic Force Microscopy (AFM) and Vertical Scanning Interferometry (VSI).These types of observations reveal that the exceptional spatial variability of reaction rates is linked to the inherent presence of inhomogeneities and/or discontinuities at the level of the mineral lattice.Such irregularities in the crystal structure give rise to uneven distributions of surface energy driving dissolution kinetics.Given the random nature of defect distributions, the associated spatial variability of rates can be interpreted in a stochastic framework.A first approach is based on the analysis of sample probability densities, also termed as rate spectra (Fischer et al., 2012), of spatially and temporally resolved reaction rates (Bibi et al., 2018;Bollermann & Fischer, 2020;Brand et al., 2017).Some recent studies are also focused on providing a geostatistically based analysis of the evolution of the topography of a mineral surface.These consider the characterization of the spatial heterogeneity of the topography through the assessment of its spatial variogram as an indicator of spatial correlation (Pollet-Villard et al., 2016;Stigliano et al., 2021).With reference to this modeling approach, Siena et al. (2021Siena et al. ( , 2023) ) provide a significant advancement upon proposing a unified theoretical framework for the stochastic analysis and interpretation of mineral reaction rates and their spatial increments.In this sense, their theoretical framework enables one to infer distributions of quantities of interest through a joint analysis of data about values of a target quantity (i.e., mineral surface topography roughness or reaction rates) and its (spatial/temporal) increments.Relying on such an approach has been shown to ensure consistency between these two sets of data, thus going beyond the use of rate spectra of the target variable and/or variograms.
A parallel line of research entails development of process models that can provide enhanced understanding of the role of the mechanistic processes underpinning the evolution of a mineral surface subject to chemical weathering.These include, for example, (a) kinetic Monte Carlo (kMC) models that are tailored to simulate dissolution at the atomic scale upon considering the reactivity of individual surface sites on the crystal and explicitly embedding it in the evaluation of the probability of dissolution events occurring at these sites (see, e.g., Kurganskaya & Lüttge, 2016;Lüttge et al., 2019, and references therein); (b) simplified conceptual models grounded on a set of probability-based rules according to which dissolution does (or does not) take place at a given location on the crystal surface (Stigliano et al., 2021); and (c) reactive transport models formulated at the continuum level that involve tracking the displacement of a fluid-solid interface due to geochemical processes that consider the intrinsic variability of surface reactions (Schabernack & Fischer, 2024;Soulaine et al., 2017).
Reactive transport models have been successfully used to simulate mineral dissolution with evolving fluid-solid interfaces at the scale of single pores across a porous system.For example, Soulaine et al. (2017) rely on the micro-continuum approach to simulate dissolution of a calcite crystal in a flowing acidic solution, showing close agreement with high-fidelity microfluidic experiments as well as with other numerical approaches (Molins et al., 2021).An extension of this model is then obtained by Soulaine et al. (2021) through a coupling with the geochemical package PHREEQC.Schabernack andFischer (2022, 2024) provide a first example of investigation of the mechanisms dominating calcite dissolution at the micro-scale through a reactive transport model.Their model is implemented in COMSOL Multiphysics software using VSI data of polished calcite from Fischer and Lüttge (2018) as targets, and is grounded on an empirical parametrization of the surface reactivity through the local topographic gradient.Notably, none of the above-mentioned studies provide quantitative comparison against experimental spatial maps of reaction rates or surface topography, nor provide analysis of the convergence properties of the dynamic models down to a sub-nm vertical resolution, these elements being critical to gain fundamental understanding of the dynamics underpinning such reactive processes.
In this Letter, we show an innovative approach where we combine for the first time a unique real-time (and openly available) data set documenting absolute dissolution rates of a calcite crystal with sub-nm vertical resolution with an original reactive transport model that has been specifically developed to be compatible with the experimental workflow.Experimental observations correspond to high-resolution (i.e., horizontal and vertical resolution of 19.5 and ∼0.1 nm, respectively) in-situ AFM data obtained across a cleaved (non-polished) calcite sample subject to dissolution.In this context, we provide a mathematical framework for the representation of the evolution (in space and time) of the mineral surface topography, and analyze convergence of the numerical approach for vertical grid spacing down to sub-nm resolution.

Experimental Observations
Spatial distributions of the topography of a calcite crystal are obtained from nano-scale (in situ and real-time) measurements collected through imaging of the crystal surface upon relying on an AFM technology.The sample consists of a calcite fragment cleaved along the {104} plane (these correspond to the so-called Miller indices, {hkl}, see, e.g., De Graef and McHenry (2012), that are widely used to identify lattice planes in crystal lattices).The surface is partly covered with a metallic (Ti) mask that serves as a reference inert layer.The latter is a necessary requirement to obtain absolute reaction rates (Emmanuel, 2014;Fischer & Lüttge, 2017;Recalcati et al., 2024b).Key elements of the experimental setting are provided in the Supplementary Information S1, its exhaustive description being available in Recalcati et al. (2024b), together with details on the mask fabrication workflow.Dissolution is induced at ambient conditions (temperature T = 22°C, pressure p = p atm ).The sample is in contact with a continuously flowing fluid (MilliQ water is employed in the experiments).These settings ensure that chemical conditions are consistent throughout the entire duration of the experiment (Recalcati et al., 2023).A 20 × 20 μm 2 portion of the crystal surface including an area ≥40% covered by the Ti mask is imaged across a 1,024 × 1,024 horizontal grid (pixel size equal to 19.5 nm).The scanning frequency is set to ∼1 Hz, yielding an acquisition time Δt = 16.43 min.
Under the conditions set in the experiments, calcite dissolution involves the formation of rhombohedral monoand multilayer etch pits (see, e.g., Bouissonnié et al., 2018;Teng, 2004).The former are shallow ephemeral features whose nucleation can take place randomly on crystal terraces or at local crystal defects.Their depth is typically corresponding to a single crystal layer (Harstad & Stipp, 2007).Multilayer pits are formed when screw dislocations are present.They propagate horizontally and vertically and remove multiple layers of the lattice.The occurrence of these pits is related to an excess of strain in the crystal structure at the molecular level, enhancing local material flux and resulting in the emission of trains of steps (also termed stepwaves (see also Lasaga & Lüttge, 2001;Lasaga & Lüttge, 2003;Lüttge et al., 2019)).The space-time dynamics of these mechanistic processes yields highly heterogeneous crystal topography maps.A preliminary pre-processing phase of AFM raw data is required to retrieve actual topography maps, η(x ‖ , t), for points x ‖ on the (x, y) plane.These are obtained through subtraction of the underlying background.The latter is modeled as a first order polynomial that is assessed upon relying on the topographic data corresponding to the non-reacted portion of the surface.Additional details are offered in Recalcati et al. (2024b).Spatial distributions of rates, R (x ‖ , t) [mmol m 2 s 1 ], are then evaluated as where V m = 36.9× 10 9 m 3 mmol 1 is calcite molar volume, with the sign convention that reaction rates are positive for dissolution processes, that is, dissolution proceeds from higher to lower local values of η.

Mathematical Model
The mathematical model considers a system comprising two phases, a solid and a fluid phase, separated by an evolving surface of discontinuity S(t).An indicator function γ s (x, t) is defined for the solid phase at any time t and any position x = (x ‖ , z) within the system as: Here, H is the Heaviside step function and F (x, t) = z η (x ‖ , t) is the surface topography function (F = 0 on S).A smooth variation of γ s over a length corresponding to the vertical grid spacing Δz, ϕ s , is then introduced with the property that ϕ s approaches γ s in the limit of Δz → 0. On these bases, following the continuous surface approach of Brackbill et al. (1992), the crystal surface is replaced by nested surfaces whose local normal vectors correspond to the gradient of ϕ s : for points x σ at the crystal surface.Since ∇ϕ s is nonzero only in the transition region, a delta function δ σ is defined through the limit for Δz → 0 of ∇ϕ s as: Geophysical Research Letters Relying on such a formulation yields: where r (x σ , t) is a reaction rate (expressed per unit interfacial area) at the surface S, V corresponding to a volume containing S. The delta function enables us to convert the surface integral over S of r (x σ , t) evaluated at that surface to a volume integral of r (x, t) over V. Hence, the equation governing the evolution of the mineral surface topography can be obtained from the mass balance of the solid phase according to: The chemical problem is solved according to the following stoichiometric reaction: For simplicity, in this work we rely on a broadly employed (e.g., Lasaga, 1998) linear relationship between the reaction rate and the H + activity, a, that is, where k is a reaction term embodying the heterogeneous surface reactivity.We defer the analysis of diverse formulations for r to future studies.Upon considering the molecular diffusion of H + ions in water (Parkhurst and Appelo, 2013) as D = 9.31 × 10 9 m 2 s 1 , a characteristic length scale of 20 μm (corresponding to the observation window size), and an average fluid velocity V = 0.09 mm s 1 (see also Supporting Information S1), one obtains Pe = 0.2 (here, Pe corresponds to the Péclet number characterizing the system), which is indicative of a purely diffusive setting.On these bases, and taking advantage of the time scale separation between boundary movement (dissolution) and transport processes (diffusion), the steady-state diffusion-reaction equation for H + concentration, c, reads: where ϕ = 1 ϕ s is the fluid volume fraction.The complete set of details of the mathematical model is offered in the SI file.The code is open-source, its exhaustive description and validation being available in Starnoni and Sanchez-Vila (2024) and Starnoni, Dawi, and Sanchez-Vila (2024).

Results
We perform our analyses upon considering two scenarios: (scenario 1) we only rely on the solution of the surface topography (Equation 5where we set dissolution rates r at the value of their experimental counterparts R) and (scenario 2) we solve the full reactive transport problem (Equations 5 and 7 combined upon relying on Equation 6for the dissolution rates).This strategy enables us to: (a) verify the convergence properties of our continuous delta function-based formulation (in both scenarios), and (b) quantify the effects of transport processes on the spatially heterogeneous surface reactivity (in scenario 2).
In scenario 1, we consider a ∼2.5 × 2.5 × 0.4 μm 3 computational domain extracted in the middle of the AFM sample, which corresponds to the location of the most pronounced etch pit observed (see Box 1 in Figure 1a).The in-plane grid spacing is fixed and set to the AFM pixel resolution of 19.5 nm.The vertical grid spacing Δz is varied down to a sub-nm scale, that is, Δz = {16, 8, …, 0.125} nm, the finest grid spacing being consistent with the AFM vertical resolution (i.e., ∼0.1 nm).The solid volume fraction ϕ s is initialized from the fully resolved three-dimensional surface topography extracted from the AFM measurements at time t 0 .The convergence study employs the same spatially heterogeneous field of dissolution rates R x ‖ ) obtained from AFM (see Figure 1b).
Simulations are terminated at t end = 16.43 min, corresponding to the time interval between two successive AFM topography measurements.
We assess the error ɛ between model-based and experimental surface topography variation (Δη) through: Here, the sum is over all the computational cells, and Δη and Δη AFM are the computed and measured quantities, respectively defined as: where η 0 (x ‖ ) is the initial surface topography and N k is the number of grid elements along the vertical direction.
Results for scenario 1 are depicted in Figures 1c (surface profiles) and 1e (convergence).The numerical solution approaches nearly first-order convergence.The order of convergence is documented to increase with vertical resolution (Figure 1e; the solution exhibits order 0.66 on average, converging to an error below 0.01 with order 0.97 for the last grid refinement).Figure 1c complements these results with a comparison between the experimental and numerical crystal surface profiles along a given cross-section of the domain (AA' in Figure 1a).The profile obtained with the coarsest resolution (i.e., Δz = 16 nm) exhibits significant deviations from its experimental counterpart (ɛ = 0.2), these deviations vanishing with increasing vertical resolution (ɛ = 0.045 and 0.008 for Δz = 1 and 0.125 nm, respectively).We conclude that, for the type of problem here analyzed, where there is a sharp transition between the solid and fluid region, our delta function based discretization provides an optimal trade-off between highest quality numerical results and computational efficiency.
For the second scenario, the reactive transport model is applied to five distinct regions (see SI) across the calcite sample to limit possible dependencies of the results on the specific surface topography.To assess the influence of the transport processes on surface reactivity, we consider two cases comprised in this scenario.These correspond to considering the surface reactivity term k in Equation 6as either (a) homogeneous over the entire computational domain (hereafter termed scenario 2a) or (b) spatially heterogeneous and modulated on the basis of the experimental rates (hereafter termed scenario 2b), that is, where R is a mean experimental rate and a is H + activity at conditions of deionized water used in the experiments (see Supplementary Information S1).
Convergence is depicted in Figure 1e, where numerical values correspond to averages of the numerical results obtained on the five regions of the surface described above (individual curves are included in the SI).Convergence is not attained for the case with homogeneous reactivity.This is somehow consistent with the results obtained by Schabernack and Fischer (2022).Scenario 2b exhibits better convergence features, even as at a reduced order (about 0.3 on average) as compared to scenario 1.This is possibly related to the explicit dependence therein of the rates on the H + activity (the latter varying both in space and time) that contributes to increase complexity of model formulation and convergence rate.A comparison between experimental and numerical results calculated for the finest resolution (Δz = 0.125 nm) is depicted in Figure 1d.Numerical results for scenario 2b remain close to their experimental counterparts (ɛ = 0.08 for this region).Otherwise, the agreement between numerical and experimental quantities deteriorates for scenario 2a (ɛ = 0.43), where spatial variability of reaction rates is not included in the model.
Notably, the complete set of our results (not shown here for brevity) suggests that the magnitude of the error in scenario 2b is sensitive to molecular diffusion D, an analysis revealing that errors decrease with increasing D. The effect of an enhanced diffusion could be consistent with the spatially varying surface roughness promoting various degree of the strength of local diffusive processes.We further note that a reduction in the error could also be obtained upon adopting a strategy akin to the work of Schabernack and Fischer (2024), that is by introducing an empirical (tunable) corrective factor β as a multiplier of the rates in Equation 6.This factor would offset the local reduction in the calculated H + activity in the proximity of the crystal surface (with respect to its reference value a) as dissolution progresses, the exact experimental rates being retrieved for spatially heterogeneous (and evolving in time) values of β equal to the local activity ratio a/a.Here, we do not pursue further this strategy as no unambiguous physical meaning can be associated with such an empirical tuning factor.

Perspectives
We combine a unique real-time data set documenting absolute dissolution rates of a calcite crystal with sub-nm vertical resolution with an original reactive transport model tailored to the analysis of the dynamics of nano-scale mineral dissolution processes.Data are obtained through Atomic Force Microscopy (Recalcati et al., 2024b).Our numerical results enable us to demonstrate for the first time convergence of a continuum-scale model against high resolution experimental data, bridging an important gap in geosciences and physics of fluids.This study is considered fundamental to enhance our understanding of the dynamics of mineral dissolution across diverse spatial and temporal scales.
We envision two possible avenues along which future research can be focused.A stream of research is directed toward incorporating the theoretical framework of Siena et al. (2023) and Recalcati et al. (2024b) according to which reaction rates are viewed as multimodal spatially correlated random fields.The theoretical framework of these authors explicitly embeds the scaling features experimentally documented by spatial increments of reaction rates evaluated on the mineral surface.Upon leveraging on such a framework, one can then promptly view Equation 5 in a stochastic context, due to randomness of the rates, which is in turn consistent with the random nature of the distribution of structural defects in the mineral lattice (e.g., screw dislocations).As such, one could then consider the solution of Equation 5 in a typical Monte Carlo framework, upon relying on a collection of aptly generated random rate fields.
A parallel set of efforts should also be devoted to derive a theoretically robust relationship between local rates and spatial gradients of crystal surface topography, thus going beyond the formulation of Schabernack and Fischer (2024) based on a so-called empirical and tunable surface slope factor.

Figure 1 .
Figure 1.(a) Initial surface topography, η 0 (x ‖ ), from AFM measurements; (b) spatial field of reaction rates, R (x ‖ ), associated with Δt = 16.43 min; (c) comparison between experimental and numerical profiles taken along cross-section AA' for scenario 1 (surface topography evolution); (d) comparison between experimental and numerical profiles taken along cross-section BB' for scenario 2 (full reactive transport model); (e) convergence of the error, ɛ, with vertical grid spacing, Δz, for all scenarios.