Implementing the Variability of Crystal Surface Reactivity in Reactive Transport Modeling

Current reactive transport model (RTM) uses transport control as the sole arbiter of differences in reactivity. For the simulation of crystal dissolution, a constant reaction rate is assumed for the entire crystal surface as a function of chemical parameters. However, multiple dissolution experiments confirmed the existence of an intrinsic variability of reaction rates, spanning two to three orders of magnitude. Modeling this variance in the dissolution process is vital for predicting the dissolution of minerals in multiple systems. Novel approaches to solve this problem are currently under discussion. Critical applications include reactions in reservoir rocks, corrosion of materials, or contaminated soils. The goal of this study is to provide an algorithm for multi-rate dissolution of single crystals, to discuss its software implementation, and to present case studies illustrating the difference between the single rate and multi-rate dissolution models. This improved model approach is applied to a set of test cases in order to illustrate the difference between the new model and the standard approach. First, a Kossel crystal is utilized to illustrate the existence of critical rate modes of crystal faces, edges, and corners. A second system exemplifies the effect of multiple rate modes in a reservoir rock system during calcite cement dissolution in a sandstone. The results suggest that reported variations in average dissolution rates can be explained by the multi-rate model, depending on the geometric configurations of the crystal surfaces.


Introduction
Current reactive transport codes rely on a constant surface reactivity for a given set of chemical parameters. Variability of surface reaction rates is tied to changes of any of those extrinsic parameters (Steefel et al. 2015). Constraints, such as pH, temperature, ionic strength, saturation index, etc., may vary as a function of local differences of solute transport conditions in complex porous media (Molins et al. 2014). The variability of effective surface rates of reacting solid material is thus controlled solely via transport heterogeneities in RTMs (Menke et al. 2018). This exclusively extrinsic view of surface reactivity and its variability is in contrast to multiple studies confirming the importance of intrinsic variability on reaction rates (Emmanuel and Levenson 2014;Fischer et al. 2012). Also, in Gray et al. (2021) it has been found that under far-fromequilibrium conditions, this approach underpredicts the dissolution rate of calcite crystals. The authors speculate that this discrepancy mighty come from a dependence of the dissolution rate on the local orientation of the cutting plane to the crystal lattice directions.
More specifically, crystal surface reactivity shows a variability of two to three orders of magnitude at constant chemical conditions. Such variability is dictated by specific material properties, such as the density of structural defects, crystal orientation, and the surface topography at the nanometer scale ( de Assis and Reis 2018; Fischer and Luttge 2017; Godinho et al. 2014;Saldi et al. 2017;Tang et al. 2004). The quantification of the surface reactivity is possible via surface direct measurements of the material flux (Fischer et al. 2012). Resulting rate maps do not rely on assumptions regarding the size of a so-called reactive surface area. Instead, the analysis of the local material flux provides insight into the spatial pattern of contrasting surface reactivity. The existence of multiple rate modes (instead of a simple mean rate) is identified by the statistical analysis of rate datasets in the frequency domain. Such rate modes are utilized to quantify and localize rate contributions that combine to an overall rate of a reacting complex material, e.g., in Noiriel et al. (2018). Thus, the rate spectrum and rate mode approaches offer an opportunity to parameterize the reactivity of solid materials in numerical simulations . In applying this strategy to reactive transport modelling, RTM concepts are able to address the intrinsic variability of surface reactivity. This is crucial for the analysis of reactivity pattern at the pore scale with feedback to local transport heterogeneities (Deng et al. 2018). In this study, we suggest a parametrization of surface reactivity via multiple rate components. The differences in surface reactivity and, consequently, the different rate modes are based on the coordination number (number of next neighbors) of reacting surface building blocks considered (Lasaga 1998). Hence, the reactivity is implemented via geometric constraints in a sense that the kink site reactivity is greater compared to the reactivity of edge sites or surface sites.
Recently, heterogeneous surface reactivity has been parameterized in reactive transport models by using the surface slope factor as a proxy value for the intrinsic reactivity of reacting materials (Karimzadeh and Fischer 2021). Here, as a special case of single crystals, we categorize specific crystals portions and parametrize them with measured rate data. In the following study, the impact of the variability of the dissolution rate is simulated on a single crystal exposed in a channel flow and a single crystal grown into a cylindrical pore. These test cases exemplify how the variability of the dissolution rate influences the geometrical evolution of the crystal. The impact of the multi-rate model is quantified, and the results are compared with a single rate standard model based on measured calcite dissolution rates. As a second model system, the dissolution of calcite crystals in a porous matrix based on a micro-CT image is simulated using the new model and standard reactive transport modeling. Finally, the multi-rate and the standard RTM model are compared and the impact of the dissolution rate heterogeneity on the permeability evolution during dissolution is investigated.

Dissolution Modeling
Three models are used in this paper to simulate crystal dissolution. In the single-rate model and multi-rate model, it is assumed that the dissolution occurs under far-from-equilibrium conditions and that the dissolution reaction is purely kinetically limited. This assumption is made to study and model the impact of surface-controlled conditions, as experimentally and analytically shown in previous studies. Thus, it is assumed that the residence time of the fluid is short enough to allow for far-from-equilibrium conditions. Under such conditions, fluctuations of fluid saturation play no role on the effective reaction rate due to the plateau-shaped rate curve as a function of G (Lasaga 1998) and the problem reduces to a time-dependent surface evolution problem. We are aware of the competing impact of both factors, i.e., surface and transport control toward chemical equilibrium conditions, where this superimposition of controlling factors is hard to deconvolve. For further simplification, we utilize a Kossel-type cubic crystal and parameterize specific surface sites, e.g., corners, edges and faces with varying kink-site densities according to literature data based on calcite dissolution experiments. We are aware of the different crystal symmetries. Nevertheless, we wanted to reduce additional complexity of the geometry under consideration.

Numerical Simulations
Simulation of crystalline dissolution is a long-standing problem in computational geoscience, and various different approaches have been developed, for an overview see Molins et al. (2020). Further, numerous approaches to simulate the crystal evolution on Cartesian grids have been proposed, such as immersed interface or diffuse interface methods. For this study, we adopt a volume of solid method similar to the ones described in Molins et al. (2012), Gray et al. (2016), Prasianakis et al. (2017) and, by extension Soulaine et al. (2018). This formulation or variations thereof are used when the grid is not aligned with the phase boundaries in the domain, as in this study. In all cases considered in this paper, the computational domain is discretized with a structured hexahedral grid or voxel grid.
To simulate the dissolution of the crystal, each voxel is assigned a solid mass fraction m given by the volume fraction of the crystal in this respective voxel.
Each voxel is assigned either to the fluid, interphase, or solid phase, where the solid phase is defined as the part of the domain having a solid mass fraction of 1, the interface phase as having a nonzero solid mass concentration 1 > m > 0 and the fluid domain as having m = 0 (see Fig. 1). Generally, flow and transport problems are solved only in the fluid domain.
(1) m(x, t) ∈ [0, 1]. The time-dependent surface evolution (dissolution) problem is then discretized using the explicit Euler method, meaning that the solid mass fraction is updated according to: with the volumetric dissolution rate K depending on the local topography of the crystal surface and the time step of the model t.
The volumetric dissolution rates in the interface voxels are related to the surface dissolution rate R given in mol m 2 s , which can be determined experimentally. To make the model comparable with measurements of surface dissolution rates, it must hold that the flux from the surface into the fluid domain is given by: with a suitable approximation of the area of the reactive surface exposed to the fluid domain A V , depending on the model assumptions.

Single Rate Model
In standard reactive transport modeling, it is assumed that the dissolution reaction proceeds with a fixed rate on the surface of the crystal, i.e., the normal flux from the surface is given by a function of the concentration on the face. If it is assumed that the solution near the face happens under far-from-equilibrium conditions, this reduces the dissolution rate to a single constant, i.e., the normal fluxes given by: where R is the dissolution rate given in mol m 2 s , D is the diffusion coefficient of the dissolving species in the fluid and is the surface normal of the interface. Hence, the volume of solid function evolves according to: where A V is the area of the interface between the crystal and the fluid for the specific voxel.

Numerical Estimation of the Surface Area
The surface area A V of an interface voxel is approximated using a method similar to the scheme described in Gray et al. (2016). The local surface area for each voxel is estimated by summing the area of the surfaces exposed to the fluid phase. To improve the estimation, each face 's surface area is multiplied by a correction factor , which is determined from the gradient of the volume of solid function V (m) estimated for each interface voxel. Hence, the solid-fluid interface area A V for each voxel is given by: where | | A n | | is the area of the face shared with the neighbor n in the six-neighborhood N 6 of the voxel, and is the Heaviside-step-function selecting only fluid voxels.
The gradient of the volume of solid function is estimated using Young's method as described in Youngs (1982). Given the gradient, the surface area correction factors are i.e., | | A n | | n ( V ) is the projection of the voxel surface area in the direction of the estimated surface normal. This method yields a better estimate of the surface area then adding the voxel surface areas.
In order to assess the accuracy of the method, the surface area of different sample geometries was estimated, where the analytical solution could be computed. The test cases include a sphere, a cylinder and two cases of an inclined plain. All geometries are discretized in a unit cube. The sample geometries used are shown in Fig. 2 , while Table 1 shows the relative error of the estimated surface area for the different geometries. For comparison, the estimations using three different numerical schemes are shown: first, an area estimation using the sum of the voxel surface areas, second, an estimation of the surface area by the norm of the gradient of the volume of solid functions, estimated at the interface voxels with the method in Youngs (1982), and finally, the estimation with our scheme, using correction factors. Evidently, the corrected voxel surface areas estimate the surface area with a relative error smaller than 4%, while other methods lead to a larger discretization error of at maximum 40% and 70%, respectively.

Multiple Rate Model
As multiple experiments have shown, the dissolution rate of crystals varies locally throughout the crystal surface. A common model for this is to propose a number of dissolution rate modes, i.e., a discrete set of dissolution rates utilized for parametrization of surface reactivity. The implementation is provided via a set of surface dissolution rates R which are assigned to the voxels on the interface of the crystal surface with the fluid. Most experimental rate mode data are calculated via surface normal retreat. These are thus a projection of dissolution rates with respect to the normal direction of the macroscopic crystal surface. In order to make our 3D model comparable to the twodimensional measurements, this projection is reproduced by assuming that the dissolution only happens on the surface area of the control volumes facing in normal direction (as illustrated in Fig. 1). Using Eq. 3, and assuming that only the surface normal to some direction is active, this means that the volumetric dissolution rates relate to the projection surface dissolution rates by: where l is the edge length of the cubic control volume.
To implement the different rate modes, we parameterized the surface rates, according to literature. Thus, surface rates at edges or corners of dissolving crystals are higher than reaction rates of crystal faces (see Sect. 3). To model this behavior, an heuristic model was implemented using the coordination number of the control volumes in the structured grid, which has the same symmetries as the Kossel crystal structure. More specifically, a set of dissolution rates is introduced, each of which depends on the coordination number of the cell in the grid. The coordination number is defined for every solid voxel at each time step as the number of neighboring voxels belonging to the fluid phase.
A parameterization of the dissolution rates, depending on the coordination number has two essential properties: • Locality: It is computed on a local neighborhood, i.e., where N i is some neighborhood of the voxel at x. Hence, it is transferred to any crystal dissolution state. • Invariance: The coordination number is invariant under the crystal symmetry transformations, i.e., where is a coordinate-transformation in the crystals symmetry group. In the model, the dependency of the dissolution rate on the coordination number distinguishes three cases: solid voxels with only one face having contact to the fluid, called here face voxels, i.e., N C = 5 , solid voxels with two faces having contact with the fluid, called here edge voxels ( N C = 4 ) and solid voxels having more than two faces in contact with the fluid, called here corner voxels with N C ≤ 4 . Sample configurations with the corresponding coordination number are shown in Fig. 3b. In the model, face voxels are dissolved with the lowest rate, while edge and corner voxels are dissolved with an elevated dissolution rate, where each of these rates is determined from experimental measurements. The dissolution rates for the face, edge and corner sites are given taken from Fischer and Luttge (2017) in Table 2.
We conclude that this model is a representative of a wider class of models which use local geometric descriptors of the crystal surface which are invariant under the crystal   (2017) far-from-equilibrium conditions, flow velocity high enough to prevent any transport-controlled conditions, (2) Bollermann and Fischer (2020) symmetries. Enhanced functionality might include a larger number of rates, additional configurations of the crystal's surface building block number and type of neighbors, e.g., distinguishing between acute and obtuse steps in the edge dissolution rate. Further, to account for the variability in the dissolution rate due to structural defects in the crystal, the rates can be varied across the volume of the crystal.

Flow Model
To estimate the impact of the multi-rate dissolution model on the evolution of the hydrodynamic properties during the dissolution process of reservoir rocks the effective permeability of the sample geometries was computed following the method in Wiegmann (2007). Further on, as it is assumed that the typical time scale of the crystal dissolution is much longer than the residence time of the fluid and the quasistationary flow equations can be solved to compute the pressure drop across the domain. This shall give an indication of how the flow permeability of a rock sample containing the dissolving crystal structure might change depending on the dissolution model. To compute the permeability, the pressure and velocity fields in the pore space are computed at selected time steps by solving the steady state Stokes equations in the evolving fluid domain F(t). That means solving the problem: for the velocity (x, t) and pressure p(x, t) fields in the pore-space with the no-slip interface condition on the crystal surface, i.e., (x) = 0 . For all flow computations, symmetric boundary conditions were implemented at the boundaries normal to the mean flow direction, while symmetric boundary conditions were enforced at the boundaries transversal to the flow direction. The value for the dynamic viscosity of the fluid was chosen for water at 25• C; however, the choice of does not influence the result. The flow problems were solved using the LIR-Stokes solver in the software GeoDict (Linden et al. 2015). The effective permeability can be computed from the solution of Eq. (12) by relating the average pressure gradient through the domain to the Darcy velocity, i.e.,

Transport Model
To compute the species concentration needed for the mean rate model, a steady-state transport equation for a dimensionless ion concentration c(x), dissolving the calcite, was solved in the fluid-filled pore space: At the inlet and outlet at the top and the bottom of the domain, respectively, the corresponding boundary conditions are applied, i.e., Dirichlet boundary conditions at the inlet with an inlet concentration of c I = 1 and outlet boundary conditions at the outlet, i.e., ∇D c ⋅ = 0 , where is normal to the boundary. The same boundary conditions are imposed as no flux boundary conditions on the fluid solid interface as well as on the exterior boundaries perpendicular to the flow direction as symmetry boundary conditions.
The boundary conditions on the dissolving, i.e., reactive crystal boundary walls, are given by: where k is the reaction rate constant and c F is the ion concentration in the fluid voxel adjacent to the dissolving boundary voxel. The dissolution rate R is coupled to the dissolution rate in Eq. 5.
The transport and the crystal dissolution simulation were implemented in the software package PoreChem (Iliev et al. 2017;Greiner et al. 2019). The discretization of the fluxes was done using the finite volume method, with the concentrations of the mobile and immobile species taken at the cell centers and the velocity components at the respective cell walls (staggered grid). After each two domain evolution time steps, the flow was recomputed by solving the above steady-state flow problem.

Characterization of the Dissolution Process
Three different metrics are employed, to characterize and compare the different dissolution processes for different geometries.
Undissolved Mass Fraction: To characterize how far the dissolution has progressed, we compute the fraction of the undissolved mass remaining with respect to the initial mass of the crystal, i.e., we compute: This measure quantifies how far the dissolution process has progressed up to a time t.
Integrated Dissolution Rate: To estimate how fast the dissolution progresses at the current time t, we compute the overall flux from the dissolving crystal surface: where local rate R(x, t) is given by the local rate modes in the control volume at time t.
Contribution of Each Rate Mode: Further, the contribution of each rate mode to the overall dissolution process was computed. This means that the number of voxels participating in the dissolution process and the number of face, edge and corner voxels are counted and then the fraction of voxels taking part in the dissolution in each rate mode is computed. That is, the relative contribution of different rate modes p(R = F, E, C) is given by: This will yield information on the activity in each rate mode over time during the dissolution process and will allow to compare and characterize the temporal evolution of the dissolution process for different crystal geometries.

Rate Components
Parametrization of the new model is according to literature data and utilizes several sitespecific rate contributors. These rate contributors are determined by statistical analysis of dissolution rate maps in VSI and -CT experiments . More specifically, local dissolution rates are measured at the micrometer-scale and the site-specific rate contributors are identified by clustering, where the rate contributors form centers of rate-clusters. As these rate-modes are derived from empirical data analysis at the micrometer-scale, they are "effective" rates including nanometer-scale effects. Hence, nanometer-scale information, e.g., step velocities, cannot be directly transferred to the micrometer-scale; nevertheless, this information is included in the effective micrometerscale rates. The rate contributors in the literature vary depending on the type of crystal surface site and nanotopography as well as defect density (Fischer and Luttge 2017). Noiriel et al. (2018) showed the existence of specific rate modes controlling the dissolution rate at crystal corners, edges, and faces under pH = 4. Additional data (pH = 2.8) similarly calculated from synchrotron-based CT imaging are available (Yuan et al. 2019). Multiple dissolution studies were performed at elevated pH conditions (pH ∼ 9 ). Site-specific rate data (crystal corner, edge) are available based on interferometry measurements (Bollermann and Fischer 2020). The rate distribution of dissolving (10-14) faces is more complex ( Brand et al. 2017;Fischer et al. 2012;Levenson and Emmanuel 2013). Here, comparatively low rates are locally superimposed with rate components of single etch pits as well as coalescing pits ( Fischer et al. (2014); Kahl et al. (2020)) that may form very reactive "super" steps (Bibi et al. 2018). Consequently, we parameterize our kinetic dissolution model based on a compilation of simplified literature data, summarized in Table 2. Kinetic data summarize results for both elevated ( ∼ 9 ) and low (2.8 and 4) pH values. Thus, the available rate components cover in part the pH-depending calcite dissolution rate function. Highest rates at crystal corners [ Table 2 (1)] are based on surface data by Bollermann and Fischer (2020) and powder data measured by Rickard and Sjoeberg (1983). The specific material preparation of powder for dissolution experiments results in a high kink site density of such material; thus, the rate range is similar to that of the crystal corners. Rate data of crystal edges are rare. Combined single crystal dissolution rate data are available (e.g., Busenberg et al. (1986)). However, such combined data, i.e., averaged data of the complex multimodal rate spectrum, provide no site-specific information. Thus, Table 2 (2) provides crystal edge data based on a single dataset only . Crystal facets were analyzed using multiple microscopic methods for direct rate measurements, e.g., in Brand et al. (2017). Surface step velocities are available from several datasets ( Duckworth and Martin 2004;Jordan and Rammensee 1998;Liang and Baer 1997;Shiraki et al. 2000). Surface data based on step velocities integrated over larger field-of-views are provided in Table 2 (3), e.g., Bibi et al. (2018), and Arvidson et al. (2003).

Results
To investigate and quantify the impact of the extended dissolution model, several simple test cases were simulated which should reveal differences between the single and multirate model. First, the dissolution of a single cubic Kossel crystal was simulated as a simple example which illustrates the differences in the modeling approaches, but also serves as a test case for a grid convergence study. We are aware of the differences in the structure of calcite vs. cubic crystals; thus, we are assuming the simplified case of no anisotropy of crystal facet evolution due to the existence of acute vs. obtuse surface steps. A further test case is investigated, of a crystal within inert material, e.g., grown into a single cylindrical pore which should demonstrate a dissolution process a convex pore. Finally, a test case where reactive crystals are found in an inert complex porous sandstone material is shown. The geometry of this dataset is based on segmented micro-CT data. This final example shall demonstrate the applicability of the new numerical approach to study the competing impact of surface vs. transport control during dissolution processes in sandstone reservoir rocks.

Dissolution of a Single Kossel Crystal
As a first test case for the multi-rate model, the dissolution of a single crystal in a channel flow was simulated, as shown in Fig. 3a. The computational domain was a cube with edge length of 0.9 mm and containing a cubic single Kossel crystal with an edge length L = 0.6 mm of surrounded by pore space. As a reference case, the crystal dissolution was simulated using the multi-rate model with different rates for face, edge and corner voxels, with the dissolution rates from Fischer and Luttge (2017) in Table 1, i.e., the case of pH = 9.2 . A small time step of 1 hour was chosen for the dissolution modeling as it is not computationally intensive.

Convergence Study
Before further analysis, a brief convergence study was done to test whether the model converges to a continuous process under refinement of the computational grid and to estimate the error due to the finite resolution of the dissolution simulation. For this purpose, a cubic Kossel crystal with an edge length of L = 60 m was discretized into 50 3 , 100 3 and 200 3 voxels and the dissolution was simulated. Then, the integrated remaining mass, and the dissolution rate over time were computed for each resolution. The remaining mass over time is shown in Fig. 4a and the dissolution rate over time is shown in Fig. 4b. Accordingly, general convergence can be observed for both the dissolved mass and the dissolution rate.
As there is only a small difference between a resolution of 100 3 and 200 3 , in both measures, we restricted the system to a resolution of dx = L∕100.

Comparison Between the Single-and Multi-Rate-Model
The simulated state of dissolution of the crystal is shown in Fig. 5 for four different reaction periods. The dissolved mass and the total dissolution rate are computed over time and plotted in Fig. 6. For comparison, the same dissolution process was simulated using a single-rate model with a dissolution rate of R = 2.25 * 10 −6 mol m 2 s . The results show that the multi-rate model results in a slower dissolution during the first 4000 hours. Subsequently, the overall dissolution rate increases and surpasses the rate of the single rate model and keeps dissolving faster until the crystal is completely dissolved. This behavior can be explained by varying contributions of the individual rate modes as shown in Fig. 7, which shows the relative contribution of the rates of face, edge and corner voxels. It shows that the faces contribute mostly in the beginning, while subsequently the edge rate contribution increases, which again drops in favor of the corner rates, which dominate the final stages of the dissolution process. This is also reflected in Fig. 5, where the crystal assumes a shape where the corner rates dominate. Hence, the slowest rate dominates the beginning of the dissolution process, while the fastest rate dominates the end. As for the single-rate model, the rate decreases monotonically from the beginning, due to the monotonically decreasing surface area. Both models lead to qualitatively different results.

Dissolution of a Kossel Crystal in a Cylindrical Pore
As a second test case, moving toward dissolution in porous media, a dissolving crystal grown into cylindrical pore was investigated. The geometric setup as well as different stages of the dissolution process are shown in Fig. 8. The cylindrical pore has a radius r = 0.15 mm and length d = 0.3 mm. For analysis, the dissolved crystal mass and integrated dissolution rate are shown in Fig. 9. Again, the multirate model was compared with the single rate model where the dissolution rate was manually adjusted, such that the remaining mass curve most closely matches the curve of the multirate model, again with the dissolution rate of R = 2.25 * 10 −6 mol m 2 s . However, both models cannot be brought in agreement, since the multirate model shows a large contrast between the stages of dissolution taking place outside the pore and inside, where the dissolution inside the pore progresses only with the very small face dissolution rate. This is also visible in Fig. 10  showing the contribution of each rate mode over time. The dissolution starts predominantly over the faces, and then between 2500 h and 7500 h, the process is dominated by edge and corner rates and in the final stage, as the dissolution front reaches the pore throat, only dissolution over the face mode remains. As in most porous media, crystals are formed inside convex pores, we expect this behavior to be dominant in real porous media.

Application of the New Approach to RTM
Finally, both models were compared in an artificially generated computational domain designed to mimic calcite crystals in a sandstone reservoir. To exemplify this case, simulations of the crystal dissolution models were done on a micro-CT image from Chagneau et al. (2015), to which artificially generated calcite crystals were added in the pore space. The investigated example serves as an example sandstone reservoir, where early diagenetic calcite precipitation occupies a large portion of the initial porosity (remaining porosity + cement volume = intergranular volume, IGV). Such high concentrations of early cements prevent strong compaction of the sandstone, thus preserving a high intergranular volume during burial. Subsequent flushing with fluid undersaturated with respect to CaCO 3 results in dissolution of calcite. Here, a situation with constant IGV and a cement volume of about 35% is investigated.
For the model comparison, a sub-domain was selected from data from Chagneau et al. (2015), with a dimension of 300 x 300 x 320 voxels with a resolution of 6 m 3 , leading to physical dimensions of L X x L Y x L Z = 1.8 mm x 1.8 mm x 1.92mm. The investigated region of interest had a porosity of roughly 45% and a permeability of K = 1.61 * 10 −11 m 2 . The solid phase consists of solid grains with a typical grain size of roughly 100 m . Further geometrical characteristics are shown in Table 3.

Addition of Reactive Crystals
Calcite crystals were added to the micro-CT data by artificially growing crystals into the pore space. This was done by using a simple morphological method described in Algorithm 1, the geometric set up of which is illustrated in Fig. 11. In the first step, nucleation points (G) for the growing crystals are placed into the pore space (F). Then, by morphological dilation [denoted with ⊕ , for a definition see Soille (2013) or Serra (1984)] with a structuring element given by a sphere of radius r , a layer of solid voxels is grown around the seed points. To ensure that the crystals grow only in the pore space, the added crystals are intersected with the pore space, and portions which grow into the grain structure are removed. This process is repeated until the crystals have grown to a Table 3 Physical and hydrodynamic properties of the sandstone sample. Remarks: (1) Pore size D 50 as measured by pore size distribution by granulometry with GeoDict, i.e., 50% of the grains have a diameter smaller than D 50 (2)  given size d max . Further, a flow channel (C) with radius (r) was introduced into the pore space along the center axis of the domain with dimension (L) in flow direction. This is done by prohibiting seed points being placed in a cylinder of radius r around the center axis.
The resulting computational domain includes a volume fraction of 35% of calcite crystals in the segmented geometry, as is shown in Fig. 12. The geometric characteristics of the new domain are summarized in Table 3. The domain including the artificially grown crystals has a porosity of less than 10% and also a significantly reduced permeability of 1.27 * 10 −13 m 2 , a reduction by more than two orders of magnitude. A flow field through the artificially generated computational domain is shown in Fig. 12c , where the fluid predominantly flows through a channel across the center of the domain induced by the cylinder in the seed algorithm.

Simulation Results in Artificial Sandstone
In this final test case, the multi-rate dissolution model without transport control is compared to the RTM model given by Eqs. (14)-(15) for the transport and Eq. (12) for the flow computation. Using both models, the dissolution of roughly 70% of the calcite in the artificial sandstone domain was simulated, and then, a quantitative comparison of the dissolution processes was done by comparing average dissolution rates and the computed permeabilities of the sandstone, sample over time.

Multi-Rate Model
As before, the dissolution of the calcite cement was simulated using the multi-rate model with the dissolution rates from Reference (1) in Table 2. The resulting dissolution process is shown in Fig. 13. Figure 13a shows the dissolving calcite surface with the associated local rate component as color coding. At this (later) stage of the process, the dissolution progresses predominantly via the surface mode and the dissolving crystal surface is mainly aligned with the crystal planes. This is due to the fact that the Fig. 13 (a) Section from a centerpiece of the computational domain: Dissolving calcite crystals (green) as simulated with the multi-rate model (sand grains were set to be transparent) and the pore space. The formation of etch pits is clearly visible and also that the dissolution takes place mainly at the crystal planes. The color coding on the crystal surface shows the active rate component: corner (red), edge (white), surface (blue). (b) The same section as in Fig. 13a, showing the calcite dissolution using the RTM-model with the dissolution rate color-coded from blue to red nonaligned surfaces, which contain a higher density of edge and corner modes, are dissolved in the beginning, leaving the slower dissolving face mode. This is also clearly visible in Fig. 14, which shows the fraction of active rate modes. While between zero and 3000 h, corner and edge modes are active and make up to up more than 20% of the active voxels. Subsequently, this drops to below 10% in the later stages of the process, as only surface modes survive. This effect is also caused by the generally convex shape of the pores, similar to the behavior observed in the cylindrical pore in Sect. 4.1.3. The values of the dimensionless numbers indicate a strongly convection dominated and reaction limited transport regime, which is consistent with the fact that only a very small variation in the concentration of the dissolving species is observed and which is consistent with the far-from-equilibrium assumption. This leads to a nearly uniform concentration on the pore walls and thus a nearly uniform dissolution rate, showing that transport control can be neglected in the present case. Hence, the dissolution model closely follows the mean rate model in Eq. 4. In contrast to the multirate component model, the single rate constant model shows a homogeneous surface retreat of the calcite face as shown in Fig. 13b, showing the same region as Fig. 13a. The dissolution rate is color-coded from blue to red and is mostly determined by the estimated surface area of the corresponding voxels. Hence, it shows a much smaller variance than for the multi-rate model.

Model Comparison
Finally, a quantitative comparison is done between the two models in the sandstone sample. The resulting remaining mass and dissolution rate over time are shown in Fig. 15a and b. As the simulation time in the reactive transport model was much shorter than for the multi-rate model, the time-scale for the RTM was scaled by a factor of 300 for the comparison. The remaining mass plot shows that the models deviate due to the nonhomogeneous dissolution rate of the multi-rate model. The rate plot shows a higher dissolution rate for the first 3000 h of the simulation; however, for the rest of the dissolution process, the dissolution rates drop below the average rate in the mean rate model. This is due to the stronger dominance of the slow face dissolution rates after the first 3000 hours of the process, as is shown in Fig. 14, which can be attributed to the convex pore shapes, where the edge and corner rates are eliminated over time. Ultimately, this is reflected in the permeability over time plot (see Fig. 16), where the permeability increases faster for the multi-rate model in the beginning of the dissolution process, but later, the face rate modes become dominant, leading to a smaller increase in permeability than for the RTM model. Eventually, for both models the permeability converges to the value for the geometry without the calcite in 3 ( K = 3.5 * 10 −11 m 2 ). However, temporally, the dissolution process is qualitatively different for the models and cannot be brought into agreement, by scaling in time.

Discussion and Conclusion
In this paper, a new model for far-from-equilibrium crystal dissolution has been proposed, which assigns dissolution rates, based on the local geometry of the dissolving crystal. The model is parameterized using empirical rate-modes obtained from statistical analysis of single rate as well as rate map data. Simulation results of the new model for the dissolution of single crystals under far-from-equilibrium conditions have been presented, and their temporal evolution has been analyzed and compared to a representative substitute singlerate model. It has been shown that depending on the geometric setup, large contrasts in overall dissolution rate can be explained by the varying impact of different rate modes in specific geometries over time. In all test cases, the temporal evolution of the dissolution rate could not be reproduced using a substitute mean rate model. Further, it has been demonstrated that the dissolution process is fundamentally different for crystals grown into a convex pore compared to a fully exposed single-crystal. In addition, it has been shown that this effect is purely geometrical-controlled and rate mode related without any impact of transport control. Finally, the effect of a variable rate model on the dissolution simulation has been tested by using a micro-CT dataset of a sandstone, modified by adding calcite crystals to the segmented geometry. Strong deviations have been found in the temporal evolution of both porosity and permeability of the single vs. multi rate models. Major reason for the contrasting behavior is the dominating initial effect of elevated rates at crystal edges and corners, not present in the single rate model.
Even though this study was conducted using synthetic datasets, the results imply critical relevance of the suggested approach for a multitude of natural and technical materials. Hence, we suggest to extend this approach for the simulation of reactive transport processes at larger field-of-views, including more complex spatial distributions of the reacting crystals.
As stated in Sect. 2.1.2, we conclude that our present modeling can be extended within a class of models where the dissolution rate depends on local descriptors of the geometry or volume of solid function. An extension of the number of rate modes depending on a statistical analysis of rate maps, including a differentiation between acute and obtuse steps, is planned. In addition, randomized dissolution rates could be used to model the formation of etch-pits from nano-scale defects, as well as variations in defect densities of the crystal.
A further improvement is expected by using machine learning techniques and directly learning distributions of rates depending on the local neighborhood of surface voxels on Fig. 16 Computed effective permeability for the artificially generated rock overtime micro-CT or VSI data. An artificial neural network could be trained on 3D-dissolution-rate data and local dissolution rates on different -CT-data could be predicted by using the local geometry as a parameter.
Overall, by using this data-driven technique, variations in crystal reactivity can be learned on a representative data set and then utilized to predict their effect in additional porous materials.
Author contributions TP contributed to conceptualization of the study, model development and numerical implementations, simulations, writing of the manuscript; CF contributed to parameterization based on experimental data, conceptualization of the study; PG contributed to software development; OI helped in conceptualization of the study and proofreading of the manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL. We gratefully acknowledge financial support by the BMBF, Grant 03G0871A to T.P., C.F. and O.L.

Data availability
The data generated in this study can be made available on reasonable request.
Code availability PoreChem is a proprietary software developed at Fraunhofer ITWM; therefore, the source code is not publicly available.

Conflict of interest The authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.