Simulation of mineral dissolution at the pore scale with evolving fluid-solid interfaces: review of approaches and benchmark problem set

This manuscript presents a benchmark problem for the simulation of single-phase flow, reactive transport, and solid geometry evolution at the pore scale. The problem is organized in three parts that focus on specific aspects: flow and reactive transport (part I), dissolution-driven geometry evolution in two dimensions (part II), and an experimental validation of three-dimensional dissolution-driven geometry evolution (part III). Five codes are used to obtain the solution to this benchmark problem, including Chombo-Crunch, OpenFOAM-DBS, a lattice Boltzman code, Vortex, and dissolFoam. These codes cover a good portion of the wide range of approaches typically employed for solving pore-scale problems in the literature, including discretization methods, characterization of the fluid-solid interfaces, and methods to move these interfaces as a result of fluid-solid reactions. A short review of these approaches is given in relation to selected published studies. Results from the simulations performed by the five codes show remarkable agreement both quantitatively—based on upscaled parameters such as surface area, solid volume, and effective reaction rate—and qualitatively—based on comparisons of shape evolution. This outcome is especially notable given the disparity of approaches used by the codes. Therefore, these results establish a strong benchmark for the validation and testing of pore-scale codes developed for the simulation of flow and reactive transport with evolving geometries. They also underscore the significant advances seen in the last decade in tools and approaches for simulating this type of problem.


Introduction
Study of flow, transport, and reactions in geological materials has historically relied on treating the porous medium as a continuum. The underlying assumption is that a representative volume or REV can be defined where at each point in space all phases are assumed to exist simultaneously [11]. Bulk parameters such as porosity, permeability, or reactive surface area are then used to characterize this medium in what can be referred to as the Darcy scale.
The pore scale, on the other hand, is defined as the scale at which each point of space is occupied by a specific phase, whether fluid or solid. Thus, the porescale approach requires knowledge and correct representation of the spatial distribution of the fluid and solid phases and its evolution with time. At the pore scale, the continuum assumption applies to each phase separately. Hence, medium properties such as viscosity, diffusivity, and to an extent, reaction rate constants are measurable directly and independently. In contrast, Darcy scale parameters such permeability, effective diffusivity, dispersion, or reactive surface area must be fitted to aggregate experimental data. The functional form for these relations is usually based on empirical or theoretical relationships, e.g., Kozeny-Carman for permeability [16,48], Archie's law for effective diffusivity [6], or the two-thirds expression for reactive surface area [54]. As a result, the pore-scale approach provides insights not available with the porous medium treatment.
The last decade has seen an explosion of pore scale studies: these include imaging and characterization techniques, laboratory experiments, numerical simulations, and combinations of these approaches. Many of the applications of these studies have focused on the relations between porescale processes in fully resolved geometries and Darcy scale parameters such as permeability, effective diffusion, dispersion, and reactive surface area [12, 17, 30, 33, 36, 49, 56-58, 61, 62, 70, 95, 114, 119]. Numerical modeling plays an important role in the investigation of pore-scale processes as it provides a mechanistic understanding of the relevant coupled processes. Further, simulation results resolve variables that are not easily available from experiments, for example concentration gradients within the pore space. Such data enables additional insights to be drawn, especially when comparison to experimental data is possible [15,66,75,96,97,99].
A central theme of pore-scale numerical investigations has been the study of the effects of pore structure and mineral heterogeneity on reaction rates [30,51,65,75]. A significant number of pore-scale modeling contributions have focused on the interactions of CO 2 -rich fluids with carbonate-rich porous media in the context of CO 2 sequestration, including the development of the reactive infiltration instabilities [35,66,75,96,99,108,119]. Microbially driven evolution of porous media, including biofilm formation, has also received broad attention [14,93,109].
These applications involve a set of coupled processes including multiphase fluid flow, multi-component transport, biochemical or geochemical reactions, and evolving pore geometries. However, the fully coupled problem in complex image-based computational domains within a single model has been rarely considered, in part due to the complexity of code development and the computational costs involved in performing the simulations. Further, this complexity has also made it difficult to derive reference solutions that can be used to validate newly-developed codes. While benchmarks for different subsets of the full problem have been or are being developed, e.g., single-phase flow in micro-CT images [91], wettability controls on multiphase flow [121], or flow and conservative transport [72], there is still a lack of a benchmark problem that includes geochemical processes along with evolving pore geometries coupled to these reactive processes.
In this manuscript, we set out to develop and present a benchmark problem set for the simulation of single-phase flow and reactive transport processes at the pore scale with evolving pore-space geometries. One benchmark includes validation against experimental data. The manuscript is organized as follows. We begin by reviewing the equations that govern the processes outlined above (Section 2). Because the equations that describe the processes pose a coupled, multi-faceted problem, a number of numerical approaches are available to tackle the different aspects of the problem. Next, we review these approaches from different perspectives: (i) the characterization of the solid/fluid interfaces, (ii) the evolution of the interfaces, (iii) the approaches to couple different processes, and (iv) the discretization of the computational domain (Sections 3 and 4). This review allows us to classify pore-scale codes, specifically the codes that participate in the benchmark problem set (Section 5). The benchmark problem set is then described in Section 6. Results for the benchmark problems obtained with the participant codes are presented, compared, and discussed in Section 7. We conclude with a brief discussion of the value of the benchmark presented in the manuscript and an outlook on pore-scale reactive transport modeling (Section 8).

Pore-scale governing equations
The benchmark problem set, described in detail in Section 6, entails solving the flow equation for the motion of an aqueous solution, the advection-diffusion equation for transport of a single dissolved component in the aqueous solution, the reaction equation for the dissolution of a solid phase, and the equation that describes the evolution of the solid geometry as a result of dissolution.
Although these processes are coupled, the time scales associated with them are very different. For example, the velocity of the dissolving mineral surface can be assumed to be much slower than the fluid velocity [54]. Further, at typical flow rates in subsurface environments, inertial forces can be neglected and Stokes flow assumed. We will make use of these assumptions to present the governing equations strictly needed to solve the benchmark problem set. However, this does not preclude the use of a more general form of the equations to solve the problem. For example, one can use transient equations instead of steady-state equations as is the case of some codes (see Section 5). These different equation forms will be considered part of the individual methods and discussed in Sections 3, 5, and 7. Appendix A lists all the equations for completeness.
Fluid flow at the pore scale can be described by the Stokes equations, which express mass and momentum conservation in the absence of inertia: where ρ is the fluid density, ν is the kinematic viscosity, p is the fluid pressure, and u is the fluid velocity. Transport of aqueous species is described by the advection-diffusion equation: where c is the concentration of the dissolved component and D is the molecular diffusion coefficient. Mineral dissolution takes place at the fluid-solid interface ( ) and can be expressed as a Robin boundary condition on the transport (3), here n is the outward surface normal to the fluid region, r(c) is the mineral dissolution rate, and ξ is the stoichiometric coefficient of species in the dissolution reaction. In general, r is a non-linear function of concentration and the specific form used in the benchmark problems is given in Section 6.
The mineral surfaces evolve in the direction of the local normal, with a velocity u that follows from a mass balance across the interface u = −rV m n, (5) where V m is the molar volume of the dissolving mineral. Equation 5 implies that the tangential fluid velocity at the fluid-solid interface is zero but not the normal velocity. However, because the aqueous concentration is much smaller than the mineral concentration V −1 m , the normal velocity is usually negligible.

Computational methods
A number of methodologies have been used to solve the equations presented in Section 2 ( Table 1). They can be divided into those that seek to solve the equations from Section 2 directly and those that solve a different set of equations that under certain conditions recover (1)- (5). In this section, we will review these approaches, including two of the commonly used methods, which belong to the second group: the mesoscale lattice Boltzmann method (Section 3.3) and the micro-continuum Darcy-Brinkman-Stokes approach (Section 3.5).
We will start by reviewing the methods used to discretize the equations (e.g., Fig. 1), either the equations presented in Section 2 or the equations resulting from the mesoscale and micro-continuum approaches. Broadly, discretization methods have been classified in the literature in three groups. The first group includes those models that discretize the equations using finite volume (FV), finite element (FE), or finite differences (FD) methods in what are often regarded as traditional computational fluid dynamics (CFD) approaches (Section 3.1). A second group comprises particle methods that rely on a Lagrangian description of flow and transport (Section 3.2). The lattice Boltzmann method makes up the third group and is unique in that it can be viewed as both a solution approach and a discretization method (Section 3.3).

Traditional computational fluid dynamics methods
Computational fluid dynamics (CFD) methods involve the discretization of Eqs. 1-4 or Eqs. 13-17 -using finite volume, finite differences, or finite element methods. In reactive transport modeling, the finite volume method and the finite differences have been the most common choices ( Table 1, Fig. 1a, d, and e). These discretization methods have historically been considered computationally expensive. However, increasing computational power and, more recently, the availability of high-level, open-source libraries such as OpenFOAM ® (http://www.openfoam.org) or Chombo [3,25] have made it possible to adopt these methods to the simulation of pore-scale reactive transport processes [64-66, 96, 98-100]. Structured Cartesian meshes are more common than unstructured meshes, owing to their computational efficiency and ease of implementation. Cut-cell methods, such as the embedded boundary, provide a powerful approach to capture complex geometries with structured Cartesian grids [111,112]. Codes using unstructured meshes are better able to conform to the shape of the boundary, although the mesh must be able to adapt to the changing geometry. As points on the boundary move according to the computed normal displacements (5), the mesh as a whole must be smoothed to maintain mesh quality, while simultaneously maintaining the shape of the surface [99,100].

Lagrangian approach and particle methods
Particle methods, such as the smooth particle hydrodynamics (SPH), particle-in-cell (PIC), and moving particle semi-implicit (MPS) methods represent fluids by particles with intensive properties (e.g., mass m i ) that are tracked in time as they move in the pore space. Many of these methods have been successfully used in reactive transport at the pore scale (Table 1).
Continuous variables (e.g., density) are represented as a superposition of kernel functions centered on a set of discrete particle points r i . In continuous form, this superposition for variable A(r) at point r can be expressed in integral form as the following convolution: where W is an interpolation kernel, h is the kernel size that defines the domain of influence for each particle, and W h (r) = h −d W (r/h) is the rescaled kernel for the domain of influence (d = 2 or 3 is the dimension of the space). To reduce the computational costs in numerical calculations, h is chosen to be finite such that W h (r) = 0 only near the particle.
In particle methods, the kernel function satisfies W (r)dr = 1 and W (r)r q k dr = 0 ∀k = 1..d, (7) for all q strictly smaller than the order of the kernel. This polynomial reproducibility defines the order of the interpolation and thus is the key parameter for the accuracy of the method, together with the regularity of the function W [26]. A particle description means that A is a set of Dirac masses of weight A i at location r i and representing a volume v i . Thus, Eq. 6 may be written as where (v i remains constant in divergence-free fields) and can be solved by conventional methods for ordinary differential equations (ODEs). This removes the difficulty of solving a partial differential equation as it eliminates ∇ · (uA) and the associated CFL condition. Many applications take advantage of this fact [19-21, 76, 80, 90]. An alternative solution to move particles is the semi-analytical Pollock method [75].
Second, the formulation in Eq. 8 has given birth to several spatial discretizations. For the linear differential operators -e.g., the Laplace operator ∇ 2 -the convolution yields ∇ 2 A * = A * ∇ 2 W h . In smoothed particle hydrodynamics (SPH), ∇ 2 W h becomes a set of weights that has to satisfy the moment (7). Further, direct computation of ∇ 2 A * = ∇ 2 A * W h entails computing diffusion between the pointwise sources of A by means of particle strength exchange schemes (e.g., [77]) or by random walk particle tracking [13,89]). The convolution is used in this case only as an interpolation to extend the particles over the whole space. The diffusion of particles ∇ 2 A can be computed by finite differences on a grid with values obtained from the convolution A * = A * W h , and interpolated back on particles via the convolution by W h [28] in what is referred to as hybrid grid-particle methods (e.g., Fig. 1c).
In the special case that a distribution A is a collection of point-wise vorticity, the so-called method of singularities can be used. The related velocity is then computed via the convolution with a Green kernel. This makes up the original vortex method [4,9,10], which has the advantage of eliminating the pressure computation. Vortex methods can be purely Lagrangian by using Biot-Savart laws [4,24] or grid-particles hybrid methods [32]. Moreover, fluid-solid interaction can also be managed, beside by penalization, by immersed boundaries [79] in order to enforce continuity of the normal velocity field, and by vortical flux [78] in order to enforce continuity of its tangential components.

The mesoscale lattice Boltzmann method
The lattice Boltzmann method (LB) is a special discretization of the Boltzmann equation, with applications that range from microflows to turbulent flows, and from colloid transport to multiphase flows [104]. Originating from gas kinetic theory, the elementary variables are statistical functions (populations f i ), which describe the probability of finding a particle with a given velocity in a particular location in space. Every grid point of the discretized domain, a lattice (e.g., Fig. 1b), is populated by a set of discrete velocity vectors [86,105]. The Bhatnagar-Gross-Krook (BGK) relaxation to equilibrium is the simplest and most commonly used collision model, resulting in the Boltzmann BGK equation: where τ is the relaxation parameter and is correlated to the macroscopic dynamic viscosity of the fluid, c i are the discrete velocities of the populations, and f eq i are the populations at thermodynamic equilibrium. After time and space discretization, the lattice BGK equation can be formulated as (10) where δ t is the respective time step. Such a model, in the macroscopic limit, can be shown to recover the compressible Navier-Stokes equations.
ρ ∂u ∂t + u · ∇u = −∇p + μ∇ 2 u + 1 3 μ∇(∇ · u). (12) The lattice Boltzmann equations recover the incompressible Navier-Stokes, when the fluid velocity is small compared to speed of sound c s , in the same medium u/c s << 1, i.e., when the Mach number is small (Ma << 1). Local mass and momentum conservation, along with only nearestneighbor communication, allows for efficient parallelization of the method.
For chemically reacting flows, several approaches have been developed, based fully or partially on the lattice Boltzmann framework (Table 1). Some models have made use of LB for the flow solution only ( [107,108,120], Section 3.4). The most common approach however relies solely on the LB framework, solving a standard LB fluid flow equation coupled with an LB solver for the advection-diffusion Eq. 3 of the chemical species [45,47]. Such models have been used to study simultaneous dissolution and precipitation processes for applications that range from CO 2 sequestration [44], to nuclear waste repositories [29,84], and multiphase reacting flow in magma chambers [40,74].

Combining discretizations
Combination of discretization methods within a single model, e.g., using a different method for each component of the problem, has been a common approach to take advantage of existing, time-tested software packages, or because specific methods are seen as computationally advantageous. Different combinations of methods have been reported in a number of studies, for example, LB for flow and a stochastic solver for transport [107,108], LB for flow and a FV multicomponent solver for reactive transport [120], a FD solver for flow and a random walk method for reactive transport [88], or a FV solver for flow and a semi-analytical particle method for reactive transport [75]. However, it appears that as methods have become more mature in the field, and the use of high-level libraries has increased, reliance on a single discretization method is a common approach (Table 1).

The micro-continuum Darcy-Brinkman-Stokes approach
While the equations in Section 2 imply that the location of interfaces between the phases that make up the porous medium is known explicitly, approaches are also available that do not require an explicit description of the interface. Rather than treating the mineral phase as 100% solid and the pore space as 100% fluid, the medium is characterized as a porous continuum with a given porosity and permeability. In the pore-scale limit (i.e., when the porosity and permeability are very small), such a treatment recovers the description in Section 2 [35,96]. The equation that describes flow in a such a medium is termed the Darcy-Brinkman-Stokes (DBS) equation. Soulaine and co-workers referred to it as micro-continuum approach [96], an extension to the definition given by Steefel and co-workers [102,103] that also includes open pore space.
The volume fraction of fluid, ε, is used to map the solid geometry spatially with ε = 1 indicating fluid only and ε = 0 solid only (e.g., Fig. 1a and c). Then, the mathematical problem posed by Eqs. 1-5 can be reformulated in terms of locally averaged equations and immersed boundary conditions at the solid surface. A single equation is solved for to obtain the fluid flow in the entire computational domain (written here to recover the transient incompressible Navier-Stokes equations in the fluid domains, i.e., Eq. 12 in lieu of Eq. 2): (13) where k is the permeability and μ = νρ is the viscosity.
In the pore spaces, ε = 1 and Eq. 13 asymptotically tends towards Navier-Stokes, while in the porous domain it asymptotically tends to Darcy's law. When coupled with a low-porosity / low-permeability matrix, this has the effect of making a nearly no-slip boundary condition on the fluid-rock interface [5,98]. During the dissolution process, the solid morphology evolves with the reaction at the solid surface. The boundary condition for the receding velocity of the fluid/solid interface (Eq. 5) is replaced by the mass balance for the solid phase, where ρ s is the solid density (V m = ρ −1 s ) andṁ is the rate of dissolution per bulk volume. The dissolution rate (ṁ) is related to the reaction rate byṁ = a v r, where a v is the mineral surface per bulk volume. The surface area is evaluated from the porosity gradient (a v = |∇ε|) (see discussion in [96]).
As the solid mass dissolves into the fluid phase, the divergence-free velocity (Eq. 1) is replaced by With dissolution, the volume fraction of fluid (ε) evolves from the intial value of zero in the solid -effectively simulating the motion of the fluid/solid interface -and the local permeability field is updated according to the new value of ε. This requires a constitutive relation linking permeability to porosity, for example, the Kozeny-Carman equation, The transport of the aqueous species is described by a locally averaged equation, where the dissolution is here described as a source-sink term.
The vorticity formulation of the DBS equation can be obtained by taking the curl of Eq. 13 and assuming a constant density of the solvent: where ω = curl u is the vorticity. The term curl(ε∇p) can be neglected compared to curl(εk −1 u), while the viscosity μ is allowed to depend on ε by means of Archie's law.

Geometry generation and interface evolution
The pore-scale approach requires incorporating directly the spatial distribution of the fluid and solid phases and then simulate their evolution with time. For model validation and testing of fundamental concepts, computational domains have been synthetically generated using simple geometric shapes that sometimes may reproduce those etched in micromodels (e.g., [120]). In actual rock or sediment samples, the initial phase distributions are available from experimental characterization techniques. Significant advancements in imaging techniques and processing capabilities have made it possible to resolve the porous media from the nanometer to the micrometer scale and beyond. These techniques include X-ray computed microtomography (XCMT, e.g., [15,66]), scanning electron microscopy (SEM), and back-scattered SEM (e.g., [33,69]), SEM Quantitative Evaluation of Minerals by Scanning Electron Microscopy (QEMSCAN ® , e.g., [12,55]), focused ion beam SEM (FIB-SEM, e.g., [50,114]), and optical microscopy (e.g., [96]). Generally, experimental images are two-or threedimensional gray scale representations of the medium, with each pixel providing a measure of the phase or phases that occupy that point in space (e.g., x-ray attenuation in XCMT). To convert this image to a numerical domain, an assumption needs to be made regarding the phase distribution, i.e., how the gray scale values translate to values of the volume occupied by each solid and fluid phases.
Approaches that use an explicit representation of the interface (Section 2) require binarization of the gray scale image (or ternary and higher-order segmentations for multiple phases, e.g., [55]), although this can also be automated in the code. An advantage of the Darcy-Brinkman-Stokes approach (Section 3.5) is that the voxelbased image data can be correlated to porosity fields and readily incorporated into the model [95,98]. Different methods may be used for determining where the interfaces are in the computational domain and how they evolve with dissolution: level sets, phase fields, the volume of solid method, and conforming meshes.
The level set method describes fluid-solid interfaces, (x) with a contour of a function φ(x, t) such that where φ o is a constant. The level set function φ is greater than φ o for one phase, and less than φ o for the other phase. The evolution of the level set is governed by the following advection equation, which describes the motion of the fluid-solid interface [53]: where n is the normal vector of the level set function (n = ∇φ/|∇φ|) at φ = φ o ). The phase field method is based on the idea that the free energy depends on an order parameter (φ, the phase field variable) that acts as a function indicating what phase a point in space is in. The concentration field (diffusion-only version of Eq. 3 as in [116]) as a function of the phase field evolution may be captured with where α is proportional to the molar volume of the mineral. The method replaces the boundary conditions at the interface (e.g., Eq. 4) with a partial differential equation for the evolution of the phase field [116]: where τ is the phase field characteristic time parameter, is closely related to the interface thickness, and λ controls the strength of the coupling between the phase field and the concentration c.
In the volume of solid approach (named by analogy with the volume of fluid approach commonly used in two-phase flow simulations [38]), the porosity (ε, i.e., the volume fraction of void) field maps the distribution of the solid phase onto the computational grid ( Fig. 1a and c). The fluid-solid interface is then located in cells containing 0 < ε < 1. The exact location of the interface within a grid block is not known explicitly and therefore, an accurate description requires grid refinement in the vicinity of the interface. The evolution of the medium is then captured by performing a mass balance of the solid phase in each point in space using Eq. 14. The average normal to the interface is estimated using ∇ε/|∇ε|. This is the technique used in the micro-continuum approach (Section 5.2).
Unstructured methods can adapt to the geometry of the interfaces with conforming meshes (Fig. 1e). In a conforming mesh approach, a subset of the vertices lie on the fluid-solid boundaries, with a simple conservative calculation of the fluxes [99,100]. When the mineral dissolves, the boundary points move in accordance with Eq. 5 and the flow and transport Eqs. 2 and 3 are solved again with the new geometry. This is a simplification of the Arbitrary-Lagrangian-Eulerian (ALE) method [37] because mesh motion is decoupled from the solution of the field equations by the large time scale separation. In the ALE method, the field equations are solved in a Lagrangian frame, which automatically implies mesh motion. The mesh is then relaxed to prevent entanglements and the fields interpolated to the new mesh. Here, the field equations can be solved in an Eulerian frame and no interpolation is required; for each mesh, the steady-state field equations are solved from scratch. This is very efficient computationally.

Chombo-Crunch
Chombo-Crunch is a code suite for the solution of flow, reactive transport and geometry evolution at the pore scale developed since 2010 by Trebotich and co-workers [64][65][66]. Flow, transport, and geometry evolution processes are implemented using the Chombo software package [3,25], while geochemical reactions are implemented in the CrunchFlow code [101] which is coupled to Chombo via a custom interface. Geometric multigrid solvers from Chombo and algebraic multigrid solver from PETSc [7] are available at runtime.
The governing equations are discretized directly on a Cartesian grid using an embedded-boundary (EB) finitevolume method. Chombo-Crunch solves the transient incompressible Navier-Stokes equations (Eqs. 33 and 34) or the transient Stokes equations (as is the case in this benchmark). The code solves separately flow, reactive transport, and boundary displacement over a given same time step, assuming that flow and reactive transport solutions change instantaneously as the geometry evolves. This makes it possible to update the geometry based on a reaction rate that does not change over the duration of the time step. The size of this time step is chosen as a fraction of the Courant-Friedrichs-Levy (CFL)-constraint CFL-constrained time step for boundary displacement.
The solutions of flow and reactive transport are obtained via sub-time stepping, which may be viewed as a relaxation method, until steady state is achieved in each sub-problem. This is carried out sequentially: that is, flow is solved first using the velocity of the interface in the previous time step as boundary condition (5). This is distinct from other code participants that assume this velocity is small enough that its effect on the flow and transport equations can be neglected (Section 2). Then, reactive transport (Eqs. 3 and 4) is solved using the fluid velocities obtained from the solution of the flow problem. Finally, the fluid-solid interface is moved according to the local dissolution rate via (5) via a level set method (20) from which the implicit functions describing the new embedded boundary are obtained. Geochemical reactions are coupled to transport using an operator splitting approach.
Fluid-solid interfaces are represented with embedded boundaries that intersect the Cartesian grid cells. The resulting cut cells account for the partial volumes occupied by both fluid and solid, and for the interfacial area between fluid and solid. Conservation (34), (33), and (3) are solved using a predictor-corrector projection method. A higherorder upwind method with a van Leer flux limiter is to be applied to advection terms in a semi-implicit Crank-Nicolson approach to minimize numerical dispersion. Flow and transport time steps are constrained by the CFL criterion. Flow time steps may also be constrained by the viscous time scale when viscous forces dominate.

OpenFOAM ® Darcy-Brinkman-Stokes
The OpenFOAM-DBS method solves the micro-continuum model formed by Eqs. 13-17. The solver is implemented in the open-source simulation platform OpenFOAM ® 5.0 (http://www.openfoam.org). OpenFOAM ® is a C++ library which solves partial differential equations using the finite volume method on an unstructured grid. OpenFOAM-DBS benefits from all the features offered by the OpenFOAM ® library, including code parallelization, and an entire set of tools including discretization schemes and geometricalgebraic multigrid solvers.
For the benchmark problem (Section 6), the computational domain is discretized using a Cartesian grid. The code also works on all kinds of unstructured grids. Equations 13-17 are solved using a sequential coupling strategy. The solution algorithm is presented in detail in [98]. Within a given time step, the geometry of the solid phase is updated solving (17), along with the dissolution rate computed at the previous time. Then the advection-diffusion-reaction equation (Eq. 14) transports the reactive fluid in the domain. To maintain a sharp concentration front, the advection term is discretized using van Leer flux limiter schemes. The pressure-velocity coupling (13) and (15) are handled by a predictor-corrector strategy adapted from the Pressure-Implicit with Splitting of Operators (PISO) algorithm [41]. Finally the algorithm marches in time. The PISO algorithm is not unconditionally stable and a time step restriction based on a CFL number is necessary to handle numerical instabilities.

Lattice Boltzmann code
Lattice Boltzmann algorithms are relatively simple, and in principle they do not rely on external computational mathematical libraries and solvers. Due to their simplicity, they can be programmed in many different computer languages and can be executed entirely, for example, on general purpose graphical units (GP-GPU's). A family of such codes has been developed at Paul Scherrer Institut over the past few years with applications that range from catalytic reactive flows [43], to micro-flows [81], to threedimensional porous media and multiphase flows [83,87], and to geochemical flows [85]. The lattice Boltzmann solver used here follows the passive scalar approach [84]. For the underlying fluid flow, the guided equilibrium model of Ref. [82], with enhanced Galilean invariance is implemented, thus solving the compressible Navier-Stokes equations. The chemical species are advected and diffused following the motion of the fluid.
The computational domain is discretized using a regular Cartesian grid. Each grid point represents a discretized volume, which can be fully or partially filled with fluid. The pore geometry is mapped as a staircase. In order to obtain the partial surface of the fluid-solid interface during dissolution or precipitation, an approximation is made, based on the composition of neighboring solid and fluid nodes, similar to the volume of fluid method for capturing interfaces in two-phase flow solvers. When a surface crosses a lattice cell, the wetted area and volume of solid depend on the nearest neighbor configurations (solid or fluid). In the assumption that a solid node shall start dissolving only if at least 3 out of the 8 neighbor cells are in a fluid/quasifluid state, seven possible types of sub-lattice configurations exist (2D). This allows for a continuous estimate of the position of the interface, and an estimation of the reactive surface area and fluid-solid volumes. The time step, the diffusion coefficients, the viscosity, and the reaction rates are parametrized using the non-dimensional numbers that characterize flow (Appendix E). A half-way bounce back condition to represent the no-slip boundary condition at the fluid-solid interface.

Hybrid grid-particle vortex code
The vortex code solves the Darcy-Brinkman-Stokes reaction (13)- (14) and (16) using the vortex method (Section 3.2) applied to the vorticity formulation of the DBS Eq. 18.
The splitting between Lagrangian transport and diffusion-stretching (i.e., ∇ · (μ∇ω) + ω · ∇V) was performed in [28] in three dimensions, and the Kozeny-Carman term has been implemented in its vorticity formulation in [32]. The diffusion equation was solved by an improved Particle-Strength-Exchange method [77], especially well suited to variable diffusion, including diffusing depending on shear-rate [94].
This equation is coupled to the displacement of reactive fluids (acid and product) by means of the diffusion-reactiontransport (Eqs. 17). Equations 17 and 18 are transformed into a set of ordinary differential equations (Section 3), with vorticity considered as a set of point-wise vortices to be transported. The code uses Runge-Kutta methods, of order 2 or 4, for the solution of this large set of ODE's.
The velocity u is computed from the vorticity ω by taking the curl of the stream function ψ solution of the Poisson equation −∇ 2 ψ = ω. The Poisson solvers are gridbased solvers, either the fast Fourier transform (FFT) solver FISHPACK using cyclic reduction [2,106] or the AMG solver MUDPACK [1], while the curl and gradient operations are performed by 4th order finite differences (FD) schemes on Cartesian grids.
A time-splitting algorithm is used to split transport and reactions on the one hand (Lagrangian feature), and the differential operators and coupling terms on the other hand (Eulerian feature). The communications between grids and particles are performed by kernel interpolations defined by Eq. 8, using the 3rd order kernel M 4 for hydrodynamic quantities and M 3 for chemical quantities that need to be sign preserving (cf. [18,68] for the kernel definitions). Improvement of the integration over the kernel support has been introduced in [27,60] by performing a sequence of directional interpolations. The choice of these two kernels is justified in Appendix D.2. An adaptive time step is triggered when the coupled reaction-flow reaches a quasi-stationary state, allowing to use relatively large time steps despite the restrictive stability condition induced by the benchmark problems.

DissolFoam
The solver is based on the OpenFOAM ® toolkit, together with customized libraries for mesh motion and boundary conditions that are not supported by the distributed source code. The steady-state flow and transport (1) and (35) are solved by a second-order finite volume discretization of the fields using an unstructured mesh [42]. The computational domain is decomposed into polyhedral cells, and the governing equations are solved by operator splitting, under the assumption [54] that the velocity of the dissolving mineral surface is much slower than the fluid velocity.
Equations 1-3 are closed by boundary conditions on the domain and on the exposed surfaces of the porous material. Chemical reactions on the mineral surfaces are included by imposing a Robin boundary condition (4) on the surface of the solid. The rate of dissolution controls the (normal) motion of points on the pore surfaces, where the local velocity of the boundary is given by u (5).
In order to allow for significant dissolution, a Laplacian smoothing of the mesh was implemented in conjunction with the mesh motion from Eq. 23. First, the displacements of the face centers that make up the fluid-solid boundaries are obtained from the concentration gradient via (23). Next, the displacements of the cell centers, δ C , are obtained from a solution of the Laplace equation, using the prescribed motion of the dissolving surfaces as a boundary condition. The normal displacements on the faces of the surface polygons are combined with a slip condition, which imposes a zero gradient on the tangential displacements. The vertex displacements, δ P k , are then obtained from the cell-center values (δ C ) by interpolation. Finally, the vertex positions r are updated, using a projection operator P = (I − nn) to ensure that points remain on the boundary surface: Because the projection operator is non-linear, the outer solution must be iterated to convergence for stable and accurate mesh evolution. Once the new vertex positions are determined, OpenFOAM functions can be used to recreate the geometric information: cell centers and volumes, and lists of faces, edges, and neighbors. A description of the OpenFOAM implementation, along with source codes and test cases, can be found in the Supplemental Information attached to a recent publication [99]. At present, there is no means within dissolFoam to change the topology of the mesh dynamically, although such functionality exists within the more cutting-edge versions of OpenFOAM [113]. However, when necessary, the pore spaces can be re-meshed externally, after extracting the coordinates of the surfaces.

Benchmark problem set
The benchmark is organized as a problem set divided in 3 separate parts: one that seeks to establish solutions for the flow and concentration fields around a dissolving mineral grain, together with the associated reaction rates, but without consideration of the evolving geometry; a second one that adds the evolution of the grain geometry as a result of dissolution; and a third that provides validation against experimental data of grain evolution. Parts I and II are based on the same synthetic geometry, which may be generated computationally, to form a two-dimensional domain. Part III relies on image data obtained in an experiment to produce a three-dimensional domain that is used for the simulations. While the governing equations and process models are the same, the parameter values for part I/II and part III are different because the parameters for part III were calibrated to experimental results [96].

Part I: reactive transport and flow
This part entails the simulation of flow and transport in a two-dimensional rectangular domain, with chemical reactions on the surface of a circular grain placed in the center in the rectangle (Fig. 2). A single irreversible heterogeneous reaction is considered that simulates calcite dissolution according to the following stoichiometric relationship: The mineral reaction rate is expressed as a first-order dependence on the concentration of H + , assuming that far from equilibrium conditions are maintained: where r has units of mol cm −2 s −1 , and the activity (γ H + ) has units of cm 3 mol −1 . The activity coefficient is assumed to be one in this benchmark problem, but the activity (γ H + = 1000 cm 3 mol −1 ) is retained to indicate that c H + must be in mol cm −3 for unit consistency (see Appendix C for a discussion on rate constant units). The dissolution rate only depends on [H + ], and thus, only the transport of H + needs to be considered. The benchmark problem does not include multi-component geochemistry or complexation reactions in the aqueous phase. The geometry of the grain does not evolve as a result of the reaction; thus, this part entails the solution of Eqs. 2-4 only. This part is intended to compare components of the problem without evolving geometry, as differences in these components may affect the moving boundary problem solution (Section 6.2). Due to feedback processes and potential instabilities, these differences may be amplified when geometry evolves. In particular, part I is aimed at validating the flow solution, and the concentration at the fluid-solid boundary that determines the reaction rate r in Eq. 26.
A single calcite crystal is placed in cross-flow conditions in a two-dimensional rectangular channel (Fig. 2). A uniform velocity boundary condition is applied across the inlet face, which determines the volumetric flux and flow rate throughout the domain: The grain has a cylindrical shape with a radius equal to 0.01 cm, or 1/5 of the width of the flow channel (Fig. 2). A solution, at pH 2, out of equilibrium with respect to calcite (Table 2), flows into the domain and drives dissolution of the grain. The mass flux at the inlet is assumed to take place by advection (i.e., c in · u in , with c in and u in being the inlet concentration and fluid velocity, see Eq. 27). This may be expressed as a Dirichlet boundary condition: Because of dissolution, the concentration of products, such as calcium ions, increases around and downstream of the mineral grain. The first-order reaction rate value obtained in [23] at 25 • C is used ( Table 2). The upper and  Table 2 The rate constant (k H + ) is from [23] lower walls of the channel are assumed to be no-flow, non-slip boundaries: where the non-slip condition does apply at the inlet face (x = 0) for consistency with Eq. 28. Solute mass flux across the upper and lower boundary walls is also zero: At the outlet boundary, a fixed pressure boundary condition is used. The solute is allowed to flow freely out of the domain by advection (no concentration gradients allowed): The simulation parameters and symbols are summarized in Table 2.
The solution for flow, transport, and reaction is expected to reach a steady state because the geometry is not allowed to evolve. The focus here is to compare the steady-state solutions. Therefore, matching the trajectory to steady state from an initial condition is not required and is provided in Section 7.1 for the purpose of illustration. Steady state may be identified by a plateau in the average effluent concentrations (38) and average reaction rates (37). The time to reach steady state from the initial conditions noted in Table 2 is of the order of 3 seconds; a simulation of 5 seconds is sufficient for this benchmark.
The resolution of the spatial discretization is not set as part of the problem definition as it is specific to each code and method. A component of part II (Section 6.2) addresses specifically the effect of the spatial resolution on results.

Part II: moving boundary problem
This section builds on part I by including the evolution of the grain shape resulting from the dissolution reaction. Equations 1-5 are solved, updating the position of the boundary at each time step. Conditions on the external boundaries are the same as in part I (27)- (32). The results in this section emphasize the moving boundary component of the simulations, with both quantitative (such as the average reaction rate) and qualitative (from the shape of the mineral grain) comparisons of the different codes. The simulations in part II do not reach a steady state because the grain surface and flow field evolve with time; results are therefore compared at specific times. The simulations are run for 45 min, period in which the size of the grain evolves sufficiently to make comparisons between code results but is still large enough that is captured by the grid resolution. Two additional components are considered. First, part II is used to explore the sensitivity of the results to the resolution of the spatial discretization, i.e., whether results converge to a solution as the discretization is refined. Second, part II is used to explore the sensitivity of the results to the dimensionless parameters, Péclet (P e = u in w/D) and Damköhler (Da I I = k H + γ H + SR 2 /D) number (Table 3, with each parameter defined in Table 2).

Part III: three-dimensional experimental validation
This part is intended to provide an experimental validation of the pore-scale simulators in three dimensions. For this purpose, experimental and simulation data from [96] are used. In the experiment, calcite dissolution was created by an acidic solution flowing past an hexagonal-shaped calcite column, which was placed within a microfluidic cell (Fig. 3). The experiment and methods are described in detail in [96].
The initial geometry is generated from an image of the experimental calcite column (Fig. 3a) provided in the Supplemental Information. This is a two-dimensional image that needs to be extruded in the third dimension, perpendicular to the plane, to obtain the three-dimensional column for the simulations. The height of the resulting column is 0.2 mm. The width of the channel is 1.496 mm and is defined as the distance between the two dark layers on top and bottom of the mask image. The top 72 and bottom 38 rows of pixels of the 1416 pixels are considered solid and do not need to be included in the simulation. The long dimension of the image (i.e., 2526 pixels) corresponds to 2.667 mm in the micromodel experiment. There is no strict requirement in simulating exactly the entirety of the length of the image as long as the chosen domain length is sufficient to capture the flow field and the diffusive boundary layer that form around the column. In other words, the solution should be unaffected by the boundary conditions. For example, if plug flow is used as a boundary condition at the inlet, the distance to the calcite column must be enough to develop a Poiseuilletype velocity profile before reaching the grain. Because the computational methods differ, the approach to incorporate the image data and construct the domain is part of the benchmark problem. As in parts I and II, the spatial resolution is not set specifically.
Top and bottom boundaries (defined by the edge of the black bands in the image) as well as front and back boundaries (resulting from the extrusion of the image in the direction perpendicular to the image) are no-flow, noslip boundaries (i.e., Eq. 29 is also applied at z = 0 and z = h). The left and right boundaries (the inlet and outlet) are the only ones open for flow. At the inlet, the flow rate is prescribed. When plug flow is assumed, the fluid velocity at the inlet face is 0.117 cm s −1 (27). A fixed pressure is used as boundary condition at the outlet (31). The external boundary conditions for the transport problem are the same as for parts I and II, with consideration of the front and back boundaries via (30) applied at z = 0 and z = h.
The parameters used for the simulations are those employed by [96] to obtain a match to experimental observations (Table 4). It is worth noting that the rate constant is an order of magnitude faster than that for parts I/II as it was obtained by fitting the evolution of the grain observed in the experiment to results from OpenFOAM-DBS. The simulation is run for 12000 s, as in the experiment.

Results and discussion
Results from the simulations are analyzed and compared in this section in two ways: (i) spatially, in terms of the grain shape at specified times and in terms of concentrations, and (ii) by means of aggregate measures of the simulated processes for the whole domain, namely, the fluid-solid interfacial area, volume of the solid grain and average reaction rate (Appendix B.1).

Part I
Results from part I show that concentrations of H + decreases towards the surface of the mineral as the reaction consumes it Eq. 25. The rate of dissolution at the surface is controlled by both the reaction kinetics (26) and the diffusion of reactant towards the surface. As a result, a diffusive boundary layer is established around the mineral grain as shown in Fig. 4. The boundary layer is wider on the downstream side of the grain than on the upstream side, and [H + ] values are slightly lower at the surface. This difference Fig. 3 a Masking image that defines the domain for part III (grain in black) and b 3-dimensional rendering of the grain as used in the Vortex simulations obtained by extruding the masking image in the z direction. An exactly octagonal perimeter was used to initialize the dissolFoam simulations rather than the bitmap mask is better observed when plotting pH values (Fig. 6), with pH = −log 10 (γ H + c H + ). These differences give rise to the evolution in shape investigated in part II.
The resolution of the mesh used for part I is 256 × 128 for all codes (Fig. 1), except for dissolFoam that used an unstructured mesh with an underlying structured resolution 200 × 100 and additional refinement as in Fig. 1e).
Average reaction rates calculated from the different codes agree quite well, with a steady-state value around 4.3 × 10 −8 mol cm −2 s −1 ( Table 5). The effective rate is much slower than under well-mixed conditions (i.e., under strict surface control), k H + γ H + c in = 8.9 × 10 −7 mol cm −2 s −1 ; this implies that the dissolution is strongly transport limited and insensitive to the actual reaction rate. Overall the reaction rates are consistent between all the codes (Table 5) While codes using an explicit representation of the interface produce an accurate surface area and grain volume, the Darcy-Brinkman-Stokes results are slightly different from the theoretical values for a circle (Table 5). This is due to the inherent uncertainty in the calculation of the surface area from the porosity gradient (Section 3.5). Although the lattice Boltzmann code represents the solid volume using a stepwise (or staircase) description of the boundary on a Cartesian grid (similar to the output of an X-ray tomogram), the solid mass balance is performed at each flow and transport time step. The dissolving voxels are then in a partially solid state (diffusion is allowed, fluid is stagnant) and the sub-lattice surfaces may be approximated (Section 5.3). The trade-off for robustness is the introduction of slight uncertainties in the flux calculation (dissolution rate is about 5% larger than the other methods).
The steady-state flux-averaged effluent concentration also shows a good agreement among the codes (Fig. 5b). This is also true for the Vortex code, which used a different The source for all parameter values is [96]. Conversion of the reaction rate constant from the units reported in [96] to the units used here is detailed in Appendix C initial condition (pH = 7 in lieu of pH = 2). As noted earlier, the time scale associated with transport is faster than that of the evolution of the fluid-solid interface. Consistently, this result shows that the same dissolution rates are obtained in about 1 s (Fig. 5a) even with initial conditions that vary 5 orders of magnitude, while the time scales of geometry evolution obtained in part II are in the order of several minutes (Section 7.2).
For convenience, the rate R is defined from the mass balance between outlet and inlet concentrations (Appendix B.1), and therefore, only at steady-state R matches the instantaneous average reactive flux dissolving the grain. Thus, after a short delay (while calcium ions are convected towards the outlet), the average reaction rate increases sharply with time (Fig. 5a). Nevertheless, the relaxation times for all simulations to achieve a steady state are similar. In calculating the grain surface area and volume, the height of the cylinder is taken as 1 cm In most simulations, the reaction rate approaches its steadystate value monotonically, but the lattice Boltzmann results show a small overshoot as well as a larger asymptotic value. The reason for the overshoot might be attributable to the initialization of the flow field and its subsequent relaxation to steady state, or to the sudden change of state of the dissolving voxels from solid to quasi-solid. Moreover, the LB code is not implementing a time scale separation between fluid motion and the transport of ions, but is solving a fully coupled problem. However, the steady-state flux-averaged effluent concentration shows good agreement among the codes (Fig. 5b), although the initial condition for the vortex method was different from the other methods. Figure 6 shows the concentration (in the form of pH) along different axes: horizontal, vertical, and at 45 • . The agreement between the codes is good, in particular in regard to the thickness of the diffusive layer and the gradient of pH values. As will be discussed in Section 7.2, in relation to the sensitivity of results to P e and Da, rates in the simulation in part I are to a large extent controlled by transport and therefore capturing the concentration gradients is critical to obtaining an accurate solution. Differences exist between results from the different codes, especially for the concentration values in cells at or near the surface of the solid, which is represented differently in the different codes (Section 4, Fig. 1). Relatively high pH values are observed for the lattice Boltzmann code, which represents the interface as a staircase (Fig. 1). Artifacts related to extracting the values and plotting them on a line that does not align with the meshes employed also contribute to differences, especially for the Vortex code that also stores concentration values in the solid.

Part II
In contrast to part I, here the geometry is allowed to evolve from the local surface reaction. Initially, the results correspond to the long-time limit of Fig. 5. The time scales in Fig. 7 are about 1000 times longer, reflecting the difference in molar volume between the solution (10 5 cm 3 mol −1 ) and mineral (36.9 cm 3 mol −1 ). As a result of dissolution, the surface area exposed to the fluid and the size of the grain decrease with time ( Fig. 7b and c). Although the shrinking of the grain reduces the effluent concentration, the effective reaction rate increases because of the area of the grain decreases more rapidly than the effluent concentration (Fig. 7). This is because, as the grain becomes smaller, the diffusive control on the overall reaction rate decreases. However, the average rate in Fig. 7 is still significantly smaller than the reaction limit 8.9 × 10 −7 mol cm −2 s −1 . Results from the different codes agree quite well with one another, with differences of up to 10% in the average dissolution rate and somewhat smaller deviations in the grain area and volume.
The lack of fore-aft symmetry in the concentration field (Fig. 4) leads to non-uniform shrinking of the disk, with most of the dissolution occurring on its leading (upstream) edge while the most downstream point moves much less. Figure 8 shows the evolving grain shape as calculated by the different codes. The overall grain shape is quite similar in all Results from OpenFOAM-DBS, lattice Boltzmann, vortex, Chombo-Crunch, and dissolFoam simulations are compared. Here x represents the distance along the line the codes, but some noticeable differences develop at longer times. The dissolution in the vortex method is more uniform over the surface of the disk than in the other codes. The dissolution on the trailing (downstream) edge of the disk predicted by Chombo-Crunch and dissolFoam are similar and slower than results from the other codes.
As noted in Section 2, the time scales associated with different simulated processes are very different. In particular, the displacement of the fluid-solid interface in response to dissolution is much slower than all other processes such that flow and transport may be solved at steady state with Eqs. 1, 2, and 35. Depending on the specific approach and particular development history, each of the participant codes use different strategies to step through time in each component in order to solve each the overall problem and capture the coupling between processes (Section 5). The good agreement between code results thus also serves as validation that the different approaches can simulate the benchmark problem accurately with their specific strategy for time stepping and process coupling. The steady-state assumption for flow and transport is taken advantage directly by dissolFoam (Section 5.5), which solves the equations in the steady-state form. The methods to update the position of the interface discretely, whether based on a threshold, as in the LB code (Section 3.3) or limited by the boundary CFL (Section 5.1), yield similar results. Although the grain geometry is updated at nonuniform intervals, average interval times are of the same Fig. 7 Evolution of a the average dissolution rate, b the surface area (perimeter in two dimensions), and c the volume (area in two dimensions) as a function of time during the dissolution of a two-dimensional disk in part II. Panel d shows the discrete times used to draw for the solid lines shown in a-c, which correspond to the discrete times the geometry was updated for dissolFoam and Chombo-Crunch order of magnitude for Chombo-Crunch (13.3 s) and dissolFoam (30.4 s) (Fig. 7d). While the LB code updates the geometry at discrete times, the mass balance of the solid is performed at each time step (Section 5.3). Therefore, the flow and transport solutions are obtained every 2.54 · 10 −6 s, and no assumption is made regarding the time scales of the processes. The Darcy-Brinkman-Stokes codes (OpenFOAM-DBS and Vortex) update the geometry every time step and solve for flow and transport at each step. OpenFOAM-DBS used a constant step of 10 −3 s. Figure 9 illustrates the convergence of the various codes for different resolutions. The average reaction rate is plotted as a function of time for a number of different mesh resolutions. All the codes are showing convergent solutions with each doubling of the resolution making a smaller difference in the predicted reaction rate. This behavior can be partially attributed to the extent of the convection and concentration boundary layers that are formed around the disk. A larger number of grid points (grid refinement) allow for a more accurate resolution of the physics in the vicinity of the disk. In the case of dissolFoam, the results are essentially independent of mesh resolution, since the use of unstructured mesh with local refinement on the disk surface allows for a smoother and more accurate representation of the processes occurring near the surface. The converged solutions from each code are similar. Figure 10 reports results at varying P e and Da I I but at the same resolution as reported earlier. In all these simulations, the Reynolds number is kept constant Re = 0.6, with a fixed inlet velocity u in = 0.12 cm s −1 . The codes Here, dissolution is almost entirely limited by the reaction rate at the surface of the disk (R = 8.9 × 10 −4 mol cm −3 s −1 ), which dissolves almost uniformly. The dissolution rate predicted by Chombo-Crunch is noticeably smaller for this case, but still within 10% of the other codes. OpenFOAM-DBS and vortex were not able to complete this test because the reactant penetrated in excess the fictitious porous solid and the fluid/solid interface likely spread over too many cells.
In the third case (P e = 600, Da I I = 17800), the reaction rate has been increased by a factor of 100, so the solution is essentially in chemical equilibrium all the time (transport limit). The dissolution is only slightly enhanced over the base state, which suggests the base state itself is deep in the transport-limited regime. The dissolution rate predicted by Chombo-Crunch is very similar to the lattice Boltzmann and dissolFoam results but the predicted surface area is significantly smaller. Chombo-Crunch simulation was affected by instabilities and the area could not be captured accurately but the rate was due to being fully transport limited. Results from the codes using the DBS approach (Vortex and OpenFOAM-DBS) are similar to those of the base case.
Finally we examine a reduced Péclet number (P e = 6), while keeping the Damköhler number (Da I I = 178) constant. Here the dissolution time scale is reduced by about an order of magnitude, due to the much higher reaction rate (R = 8.9 10 −1 mol cm −3 s −1 ). Result from all codes compare well but the Chombo-Crunch simulation could not be run to completion.

Part III
The agreement between code results for the two-dimensional simulations on the synthetic geometry in Section 7.2 is remarkably good. Part III challenges the participant codes in two additional aspects of pore-scale reactive transport simulations, namely in the consideration of three dimensions and in incorporating image data to construct the numerical domain. Experimental data of the evolution of the solid geometry are also available for comparison.  As noted in Section 4, different approaches are used to incorporate image data into pore-scale models. These approaches are often tightly connected to the solution approach. In this manuscript, the two codes based on the Darcy-Brinkman-Stokes formulation incorporate the solid phase by assigning an initial value of zero to volume fraction of fluid (ε = 0) at the solid cells, and 1 otherwise (ε = 1). In contrast, Chombo-Crunch and dissolFoam must have an explicit description of the grain surface. For Chombo-Crunch, the mask image was saved as a binary file (0 for pore and 255 for solid) and processed with a dilate Fig. 11 Evolution of the a grain volume and b surface area obtained from part III. Experimental results (dots) can be compared with simulations from OpenFOAM-DBS (red), vortex (blue), Chombo-Crunch (cyan), and dissolFoam (orange) codes. Panel (d) shows the discrete times used for the solid lines shown in (a-c), which correspond to the discrete times the geometry was updated for dissolFoam and Chombo-Crunch filter followed by an erode cycle from the imageJ software [92]. The resulting grayscale image-with a smoother grain surface-was read by the code, which then used a level set algorithm (19) with a threshold value of 128 to generate the implicit functions that define the embedded boundary. The lattice Boltzmann code did not participate in part III. DissolFoam took advantage of the octagonal shape of the grain to initialize the domain with an exact octagonal perimeter rather than using the image directly. Further, dissolFoam also took advantage of the ability to remesh pore spaces externally (Section 5.5), as the octagonal shape of the object caused more skewness in the mesh than the circular shape in part II. The results were insensitive to the frequency of re-meshing. Overall, the solution approach is very efficient computationally. The simulations required about 100 time steps to determine the boundary motion, with less than 100 iterations of the linear solvers per time step; a typical mesh contains of the order of 10 5 cells.

Evolution results
The evolution of the grain in terms of both volume and surface area is qualitatively very similar to that of the disk in part II (Section 7.2, Fig. 11). The grain volume initially  Fig. 11, [96]). There is qualitative agreement between codes with OpenFOAM-DBS results being the closest to the experimental data. However, the results from the other codes are similar to one another and indicate a more rapid dissolution than OpenFOAM-DBS. The discrepancy in time scale between vortex, Chombo-Crunch, and dissolFoam simulations, and the experiment is approximately 20%.
As in part II, although the grain geometry is updated at non-uniform intervals, average interval times for all codes can be calculated. A wider range of interval values are observed, however: 7.3 seconds for Chombo-Crunch and 86 s for dissolFoam. For Chombo-Crunch, geometry update steps are limited by the cell with the fastest rate. In this rough grain geometry where embedded-boundary cells have relatively large variations in volume fractions, this led to very strict limitations on the update intervals. The mesh relaxation employed by dissolFoam allowed for larger time steps (Fig. 11c).
Differences in the initial shape and size of the grain are evident as the different codes used different ways to incorporate the geometry information from the image (Figs. 12, 13). For this reason, initial grain volume and surface area were slightly different and results in Fig. 11 were presented normalized to initial values of grain volume and grain surface area for each code. However, the staircase outline for OpenFOAM-DBS is an artifact of the visualization software after a threshold is applied to the output porosity (ε) field. These initial differences likely have an effect on the different evolution of the grain for the different codes. In spite of this, however, the agreement is in general very good and the simulated evolution of the shape is qualitatively similar between codes (Figs. 12). As in part II, the dissolution of the trailing edge of the grain predicted by Chombo-Crunch is slower than for other codes. In this case, it also it dissolves faster in the leading edge. There is however good agreement among all codes in capturing the width of the evolving grain. The vortex code predicts a slightly wider grain than the other codes.
In the direction perpendicular to the plane shown in Fig. 12), the top and bottom boundaries were considered non-slip boundaries. Because flow is faster in this midplane and slower towards the boundaries, dissolution is slightly faster here. As a result, the grain surface recedes faster in this area, especially in the leading edge of the grain (Fig. 13). The predictions for the trailing edge by the different codes diverge more significantly. While the vortex code predicts a similar pattern as in the leading edge, dissolFoam and OpenFOAM-DBS predict instead that trailing edge remains a vertical wall (only slightly concave for OpenFOAM-DBS), and Chombo-Crunch predicts a convex evolution. Although there is no experimental data to ascertain the evolution, additional simulations with Chombo-Crunch did not change this prediction. This result however explains the longer trailing edge for this code discussed earlier as Fig. 12 was obtained for the mid-plane, i.e., the section where the trailing edge extends the farthest for this code.

Conclusion and outlook
The overall good agreement between codes in the results of Section 7 establishes a strong benchmark that can be used for testing and validation of new and existing codes for the simulation of reactive transport at the pore scale with consideration of fluid-solid geometry evolution, in a variety of transport and reactive conditions.
The problem set was organized with increasing levels of complexity. Part I addresses only pore-scale reactive transport without geometry evolution and can be solved as a first step to part II. Part II includes geometry evolution for the same domain geometry and conditions. Additional components within part II such as analysis of grid convergence and sensitivity to Pe-Da conditions are helpful in validating and testing new codes. By varying the parameters (Fig. 10), they also investigate regimes where rates are closer to surface control, although carbonate dissolution is often closer to the transport limit. Part III can be addressed separately from parts I and II but new layers of difficulty are added; namely, the need to read a relatively rough fluid-solid interface from image data to generate the initial domain geometry and consideration of the third dimension.
The five codes that participated in this manuscript use different approaches including in the form of the equations they solve to recover the governing equations, the discretizations they employ, the characterization of the fluid-solid interfaces, and how they simulate the motion of these interfaces. These approaches and methods cover many of the approaches reported in the literature to solve reactive transport at the pore scale. This gives added confidence in that the results establish a well founded benchmark.
The benchmark problem presented focused exclusively on single-phase processes. Simulation of reactive processes in multiphase systems at the pore scale is a relatively new development, e.g., [59,97], but one that will likely receive increasing attention due to its relevance, e.g., in carbon sequestration or gas release from shales in hydraulic fracturing problems.
As a first effort, the benchmark considered a relatively simple geochemical problem: a single component reacts with a single mineral in a dissolution reaction, which is simulated with a first-order rate. Natural subsurface environments are however characterized by heterogeneous multi-mineral media involving multi-component aqueous solutions (e.g., [31]). Future benchmark efforts can build on this manuscript to consider more complex geochemical problem where multiple minerals are present, and may dissolve and precipitate in different areas of the domain.
The domain considered here is also relatively simple with single grain geometries, which included an experimental validation. Well-characterized natural media with relatively large porosity and relatively homogeneous mineralogy (e.g., [75]) offer a good opportunity for follow-up benchmarking efforts.
Simulation of multi-mineral systems at the pore scale has been performed but focused on media such as sandstones where uniform resolution was sufficient to represent [55,56]. It is often however the multiscale nature of mineral heterogeneity that challenges the pore-scale approach, where geochemical evolution leads to altered porous layers adjacent to open pore space [31,63]. Specifically, extreme resolution may be needed to capture the different mineral phases at multiple spatial resolutions.
Multiscale models are emerging as a powerful to tool to simulate systems that are heterogeneous at multiple scales as they offer a reasonable compromise between the fidelity of medium and process representation and feasibility of numerical simulation. In this sense, micro-continuum approaches (e.g., Section 5.2) or hybrid multiscale methods (e.g., [8,63]) offer appealing options to expand the range of applications of pore-scale modeling.
The focus on dissolution in the benchmark reflects the topic of relevance in most applications. Precipitation processes however have also been investigated in pore-scale models, e.g., [22,29,40,45,84,120], particularly as they can have a strong feedback on flow and transport properties by blocking pore spaces. Treatment of nucleation however varies significantly in these applications and warrants a dedicated benchmarking effort. effort. Chombo-Crunch was developed by D.T. and S.M., CrunchFlow by C.I.S., OpenFOAM-DBS by C.S., the lattice Boltzmann code by N.I.P., the vortex code by P.P., and dissolFoam by A.L. and V.S. Part III experiments were conducted by S.R.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A: Additional/alternative equations
Flow at the pore scale may be described by the incompressible Navier-Stokes (33) and (34): as well as the Stokes (1) and (2). In these benchmarks, the Reynolds number is sufficiently small that fluid inertia can be neglected; thus, these two approaches are equivalent.
In the dissolution benchmarks (parts II and III), codes may take advantage of the large time scale separation between boundary motion and transport processes to solve the steady-state transport equation directly, Time-dependent solutions of transport and reaction (part I) are more tightly coupled than dissolution (parts II and III), because t A and t R are often of the same order, especially for relatively fast reacting minerals such as carbonates. Both global implicit and operator splitting approaches have been used for time-dependent transport, with the time stepping in the operator splitting constrained by the Courant-Friedrichs-Levy (CFL) criterion

B.1. Upscaled parameters
Simulation results are compared in terms of the evolution with time of upscaled parameters. These upscaled parameters include the volume (V ) and surface area of all reacting reacting mineral surfaces (A) and the average reaction rate (R). The average rate is calculated as follows: where ξ is the stoichiometric coefficient, c in is the (uniform) concentration at the inlet, given by the boundary conditions, and c out is the flux-weighted-average outlet concentration, The volumetric flux at the outlet Q is found by integrating over the outlet area In addition to these upscaled parameters, simulation results are compared on the basis of the geometry of the grain at different time points and the concentration contours are prescribed times.

B.2. Grid convergence
As methods for simulating of moving boundary problems vary greatly, we want to investigate the impact of grid resolution on results for each method separately. For this purpose, the simulations were run at different resolutions (Figs. 14, 15, and 9) in the main text. Results for the grain volume and surface area were analyzed to ensure grid convergence of the methods, and choose a resolution for which results will be assumed to have converged within a reasonable accuracy.  concentration c H + must be in mol cm −3 so that the product k H + γ H + has units of cm s −1 . The conversion from the rate constant used in [96] (k H + γ H + = 0.5cm s −1 ) is k H + = 0.5 cm s −1 1000 cm 3 mol −1 = 5 × 10 −4 mol cm −2 s −1 . (42) However, in [96], this rate is applicable to the rate of HCl consumption when reacting with calcite according to the following stoichiometry CaCO3 (s) + 2HCl− > CaCl 2 + H 2 CO 3 .
To maintain consistency with the rate expressed for calcite, one must multiply by the stoichiometric coefficient of HCl in Eq. 43,

D.1. Lattice Boltzmann dimensionalization
Dimensionalization of the LB computations is a process that needs special care. Lattice Boltzmann unit conversion to physical units can be done after matching the characteristic non-dimensional Reynolds, Péclet, and Damköhler numbers. For a 256 × 128 discretization grid L y = 128 (in lattice units), each lattice space unit in parts I and II corresponds to w/128 = 3.91 × 10 −4 cm. For the current setup, viscosity is defined as ν = τ f ρT . The relaxation parameter for the fluid phase, τ f , is set to τ f = 0.5 in lattice units. By equating Re=Re LB =0.6, using the aforementioned viscosity, the inlet velocity can be calculated as u inLB = 0.00078125 (in lattice units), which corresponds to u in = 0.12cm s −1 .
Once the lattice velocity is set, the duration of the time step δt can be calculated by equating the inlet velocities: δt = 2.54 × 10 −6 s. Note that the time step is dictated by the slow advective flow, and by choosing to keep the same time step for all processes. This leads to a fully coupled advection-diffusion-reaction scheme applicable to all flow and chemical conditions. Separation of time scales is possible by solving for steady-state flow, then steady-state reactive transport, and finally the solid geometry update. Such an approach would be sufficient for these benchmarks and would greatly reduce the number of time steps to reach the solution. Diffusivity is defined as D = τ g T . By equating the Péclet numbers Pe=Pe LB =600, the relaxation parameter τ g , which corresponds to the diffusive time scale, is set to τ D = 0.0005, for the species that follow the advectiondiffusion equation. Finally, by equating the Damköler numbers Da II =Da II-LB =178.15, the rate constant k H + LB = 10 −3.2364 . For this dimensionalization Ma<Kn<<1, thus recovering the incompresible Navier-Stokes equations.

D.2. Discussion on interpolation kernel for Lagrangian methods
The choice of the kernel used for re-meshing the particle is crucial for the accuracy of vortex and particle methods. Indeed, in order to avoid holes and accumulation of particles that would ruin the convergence, particle information F p (including vorticity, concentration, ...) in volumes v p located at positions ξ p is remeshed on to a new structured mesh (with cell volumesṽ q ). This mesh defines a new set of particlesF q at locationsξ q by means of the following convolution: since the set of particles is mathematically defined by the generalized function F = p F p δ x p v p , based of Dirac functions at x p . In practice, when is the "hat" (or "tent") function, the reaction stays confined on the fluid/solid interface, but exhibits a pH over-estimation close to the stagnation points, thus over-estimating the reaction rate. When this kernel is smoother but positive in order to be sign preserving, such as the first-order cubic spline M 4 , the fluid/solid boundary becomes fuzzy and requires us to force the reaction on the interface by means of the function ∇ , as in [96]. When using the second-order kernel M 4 from [67], which is non sign preserving since the integral of x 2 M 4 (x) is zero, no negative concentration appears despite the jump of acid concentration at the body but it leads underestimation of reaction rate. However, the hydrodynamic flow is computed with better accuracy using M 4 , as expected [28]. Consequently, the short-supported function M 3 , smoother than the hat function with a support smaller than M 4 , has been chosen for interpolating and remeshing the chemical concentrations, while M 4 has been chosen for the interpolation hydrodynamic values (velocity and vorticity).
In practice for the present benchmark, for which the reaction properties (bounds and positivity) have to be strictly satisfied, the choice of the remeshing kernel is mainly driven by the following arguments: -The hat function, is good for the estimation of reaction rate but does respect the pH bounds (pH overshoots below 2 can occur), • The kernel M 4 is smooth but M 4 (0) = 2/3 = 1; thus, it is diffusive: pH bounds are good but reaction rate is under-estimated (see formula A.4 of [18] for definition), -M 4 (formula 4.5 of [20]) is algebraically massconservative, smooth, and second order, but its negative values induce oscillations at concentration jumps and over-estimate the reaction rate. Furthermore, it is not mathematically sign preserving, although negative concentrations were never been observed in this benchmark, -M 3 (formula A.3 of [18]) is smoother than hat, first order and sign preserving, with short support. It is the best choice for reactive flows like the one considered in the present study; the reaction rate is well estimated (a bit higher than the hat function and closer to other curves) and does not go lower than the initial pH=2 bound, consistent with this purely dissolution process, -M 6 and M 6 supports are too large for this geometry, and cannot handle correctly the final state of the dissolution.
Consequently, the kernel M 4 is the best choice for hydrodynamic computations (for particle remeshing and interpolation of velocity and vorticity from and to grids), while M 3 is the best choice for interpolation and transfer of concentrations.

D.3. Darcy-Brinkman-Stokes parameter values
Parameters specific to Darcy-Brinkman-Stokes code simulations are presented in Table 6.