# Improved Volume-of-Solid formulations for micro-continuum simulation of
mineral dissolution at the pore-scale

Julien Maes  Institute of GeoEnergy Engineering, Heriot-Watt University,
Edinburgh, U.K.  Cyprien Soulaine  CNRS, University of Orléans, Orléans,
France  Hannah P. Menke  Institute of GeoEnergy Engineering, Heriot-Watt
University, Edinburgh, U.K.

###### Abstract

We present two novel Volume-of-Solid (VoS) formulations for micro-continuum
simulation of mineral dissolution at the pore-scale. The traditional VoS
formulation (VoS-$\psi$) uses a diffuse interface localization function $\psi$
to ensure stability and limit diffusion of the reactive surface. The main
limitation of this formulation is that accuracy is strongly dependent on the
choice of the localization function. Our first novel improved formulation
(iVoS) uses the divergence of a reactive flux to localize the reaction at the
fluid-solid interface, so no localization function is required. Our second
novel formulation (VoS-$\psi$’) uses a localization function with a parameter
that is fitted to ensure that the reactive surface area is conserved globally.
Both novel methods are validated by comparison with experiments, numerical
simulations using an interface tracking method based on the Arbitrary Eulerian
Lagrangian (ALE) framework, and numerical simulations using the VoS-$\psi$.
All numerical methods are implemented in GeoChemFoam, our reactive transport
toolbox and three benchmark test cases in both synthetic and real pore
geometries are considered: (1) dissolution of a calcite post by acid injection
in a microchannel and experimental comparison, (2) dissolution in a 2D
polydisperse disc micromodel at different dissolution regimes and (3)
dissolution in a Ketton carbonate rock sample and comparison to in-situ micro-
CT experiments. We find that the iVoS results match accurately experimental
results and simulation results obtained with the ALE method, while the
VoS-$\psi$ method leads to inaccuracies that are mostly corrected by the
VoS-$\psi$’ formulation. In addition, the VoS methods are significantly faster
than the ALE method, with a speed-up factor of between 2 and 12.

##  1 Introduction

Prediction of solid mineral dissolution during reactive flow in porous media
is vital for a wide range of subsurface applications, including CO2
sequestration [1], geothermal systems [2] and enhanced oil recovery [3]. CO2
storage in underground reservoirs has the potential to significantly mitigate
the environmental impact of many industrial processes. However, mineral
dissolution is a potential barrier to the long-term storage of CO2 in the
subsurface, as CO2 reacts with water to make carbonic acid that can dissolve
solid minerals and threaten the structural integrity of a reservoir [4]. In
addition, most subsurface applications involve the injection of fluids with
chemical properties that are incompatible with existing reservoir fluids and
can lead to mineral precipitation or scaling in the pore structure. Scaling is
especially prevalent near well-bores, and can significantly reduce the
permeability, and thus productivity of a porous formation. Acid injection is
then often used to improve the flow in clogged wells [5]. Thus, accurate and
efficient modelling of mineral dissolution in porous media is crucial to
improve and optimise these engineering processes.

Modelling reactive transport at the field-scale relies on the assumption that
a representative elementary volume can be defined such that flow, transport
and reaction can be described in terms of bulk properties like porosity,
permeability and macro-scale reactive constant, in what is usually referred as
the Darcy scale [6, 7, 8, 9]. Mineral dissolution modifies the pore structure
and results in a change in these Darcy-scale properties. Pore-scale numerical
experiments can be used to predict the change in these Darcy-scale properties
during dissolution. At the pore-scale, these reactions are applied directly on
the solid surface while resolving flow and transport in a representative
elementary volume of pore space directly. The effects of these dissolution-
induced structural changes on the flow and transport properties of the bulk
medium can then be estimated for use in Darcy-scale simulations.

The last decade has seen an explosion in the study of flow and transport
behaviour at the pore-scale [10, 11, 12, 13, 14, 15, 16, 17]. Recent advances
in X-ray imaging techniques have enabled direct observation and quantification
of dissolution-induced changes in pore structures [10, 18, 19, 20, 21, 11,
22]. Numerical modeling has played an important role in these investigation of
pore-scale physics, as it provides a mechanistic understanding of the relevant
coupled processes. Furthermore, simulation results resolve variables that are
not easily available from experiments such as concentration gradients within
the pore space.

Numerical modelling of mineral dissolution at the pore-scale can be performed
using Pore-Network Modelling [12, 23, 13]. However, the evolution of the pore-
space can only be predicted using the finite range of geometrical parameters
of the network. Alternatively, computational microfluidics [24] has been
applied using a range of numerical methods [25]. Interface tracking models
explicitly deform and move the solid surface, either using solid balance with
a threshold on a lattice [26, 27, 28], a conforming mesh based on the
Arbitrary-Lagrangian-Eulerian (ALE) framework [15, 29] or smoothed particle
hydrodynamics [30]. Alternatively, the interface can be captured using a
level-set function [31, 14]. For all these methods, the boundary conditions on
the solid surface can be be applied directly, or using an immersed boundary
condition. However, they require additional treatment for interface
displacement, topological changes or remeshing, which usually lead to an
increase in their computational cost [32].

The micro-continuum approach [33, 16, 34, 35, 32] based on the Volume-of-Solid
(VoS) method offers an attractive substitute for interface tracking models.
Within this approach, the fluid-solid interface is captured using an indicator
function equal to the volume fraction of void space in each cell, and flow and
transport are solved using the Darcy-Brinkman-Stokes (DBS) equation. The VoS
method is computationally efficient as it does not require remeshing or any
special treatment for topological changes.

In the standard VoS approach, the surface area of the fluid-solid interface in
a control volume is computed through the gradient of a volume fraction. In
practice, a diffuse interface may emerge that spreads across a large number of
layers in the computational grid. To enforce the localization of the reactive
boundary condition at the fluid/solid interface, a diffuse interface
localization function $\psi$ is generally introduced [16] and this formulation
is labelled VoS-$\psi$. The main advantage of the VoS-$\psi$ method is that
standard Reactive Transport Modelling dedicated to Darcy-scale can be easily
applied to simulate geochemical processes at the pore-scale by simply changing
the way the fluid-rock interfacial area within control volume is estimated
[36]. While the surface area for Darcy-scale simulations is an input parameter
that is either constant or depends on complex function of porosity and flow
rates [10, 37], the surface area for pore-scale simulations using VoS is
directly calculated from the mapping of the solid volume fraction. The main
limitation of the VoS$-\psi$ is that the accuracy of the model depends
strongly on the choice of the localization function [38], and the optimal
choice depends on a large number of parameters, such as the geometry, the flow
rate, the reactive constant, the computational mesh and the discretization
method used for the computation of gradients.

In this paper we propose two novel VoS formulations. The first formulation
(iVoS) removes the need for a localization function by computing the reaction
rate using the divergence of flux. The second formulation (VoS-$\psi$’) uses a
localization function with a parameter that is fitted to ensure that the
reactive surface area is conserved globally. The numerical models are
presented in Section 2. The iVoS and VoS-$\psi^{\prime}$ methods are then
compared with the VoS-$\psi$ method and with an interface tracking method
based on the ALE framework on three benchmark test cases in Section 3. In each
case, we show that our new VoS methods match accurately experimental results
and/or simulation results using the ALE method while being significantly
faster.

##  2 Mathematical models

In this section, the governing equations and the micro-continuum approach are
first presented. Then, we describe the iVoS, VoS-$\psi$ and VoS-$\psi$’
models, which differ only in the way the reaction rate is computed. Further,
we show how dimensionless analysis lead to the quasi-static assumption that
reduces the computational time.

###  2.1 Governing equations

The governing equations consider flow, transport and reaction at the fluid-
solid interface. The domain $\Omega$ is partitioned into fluid $\Omega_{f}$
and solid $\Omega_{s}$. Under isothermal conditions and in the absence of
gravitational effects, fluid motion in $\Omega_{f}$ is governed by the
incompressible Navier-Stokes equations

$\nabla\cdot\mathbf{u}=0,$ (1) $\frac{\partial\mathbf{u}}{\partial
t}+\nabla\cdot\left(\mathbf{u}\otimes\mathbf{u}\right)=-\nabla
p+\nu\nabla^{2}\mathbf{u},$ (2)

with the continuity condition at the fluid-solid interface $\Gamma$,

$\rho\left(\mathbf{u}-\mathbf{w}_{s}\right)\cdot\mathbf{n}_{s}=-\rho_{s}\mathbf{w}_{s}\cdot\mathbf{n}_{s}\hskip
14.22636pt\text{at $\Gamma$},$ (3)

where $\mathbf{u}$ (m/s) is the velocity, $p$ (m2/s2) is the kinematic
pressure, $\nu$ (m2/s) is the kinematic viscosity, $\rho$ (kg/m3) is the fluid
density, $\rho_{s}$ (kg/m3) is the solid density, $\mathbf{n_{s}}$ is the
normal vector to the fluid-solid interface pointing toward the solid phase,
and $\mathbf{w_{s}}$ (m/s) is the velocity of the fluid-solid interface, which
is controlled by the surface reaction rate $R$ (kmol/m2/s) such that

$\mathbf{w}_{s}=\frac{M_{ws}}{\rho_{s}}R\mathbf{n}_{s},$ (4)

where $M_{ws}$ is the molecular weight of the solid. The concentration $c$
(kmol/m3) of a species in the system satisfies an advection-diffusion equation

$\frac{\partial c}{\partial
t}+\nabla\cdot\left(c\mathbf{u}\right)=\nabla\cdot\left(D\nabla c\right),$ (5)

where $D$ (m2/s) is the diffusion coefficient. The chemical reaction occurs at
the fluid-solid interface $\Gamma$, such that

$\left(c\left(\mathbf{u}-\mathbf{w}_{s}\right)-D\nabla
c\right)\cdot\mathbf{n}_{s}=\zeta R\hskip 14.22636pt\text{at $\Gamma$},$ (6)

where $\zeta$ is the stoichiometric coefficient of the species in the
reaction. In this work, we assume that the surface reaction rate depends only
on the concentration of one reactant species, following

$R=k_{c}c,$ (7)

where $k_{c}$ (m/s) is the reaction constant.

###  2.2 Micro-continuum approach with Volume-Of-Solid

In the micro-continuum approach, the entire domain $\Omega$ is considered, i.e
fluid $\Omega_{f}$ and solid $\Omega_{s}$, and the fluid-solid interface is
tracked in terms of $V_{f}$ and $V_{s}$, the volume of fluid and solid phase
in each control volume $V$, and their volume fraction $\varepsilon=V_{f}/V$
and $\varepsilon_{s}=1-\varepsilon$. The flow, transport and chemical reaction
are solved in term of the volume-averaged velocity

$\overline{\mathbf{u}}=\frac{1}{V}\int_{V_{f}}\mathbf{u}dV,$ (8)

and the phase-averaged pressure and reactant concentration

$\displaystyle\overline{p}=\frac{1}{V_{f}}\int_{V_{f}}pdV,$ (9)
$\displaystyle\overline{c}=\frac{1}{V_{f}}\int_{V_{f}}cdV.$ (10)

The averaging process results in an extension of the Darcy-Brinkman-Stokes
equation [16]

$\frac{1}{\varepsilon}\left(\frac{\partial\mathbf{u}}{\partial
t}+\nabla\cdot\left(\frac{\mathbf{u}\otimes\mathbf{u}}{\varepsilon}\right)\right)=-\nabla\overline{p}+\frac{\nu}{\varepsilon}\nabla^{2}\overline{\mathbf{u}}-\nu
k^{-1}\overline{\mathbf{u}},$ (11)

where $k$ (m2) is the permeability of the cell. $\nu
k^{-1}\overline{\mathbf{u}}$ represents the momentum exchange between the
fluid and the solid phase, i.e. the Darcy resistance. This term is dominant in
the solid phase and vanishes in the fluid phase. To model this, the local
permeability field $k$ is assumed to be a function of the local porosity
$\varepsilon$, following a Kozeny-Carman relationship

$k=k_{0}\frac{\varepsilon^{3}}{\left(1-\varepsilon\right)^{2}},$ (12)

where $k_{0}$ (m2) is the Kozeny-Carman constant. For the acid transport, the
mass-balance equation averaged over the control volume gives

$\frac{\partial\varepsilon_{f}\overline{c}}{\partial
t}+\nabla\cdot\left(\overline{c}\overline{\mathbf{u}}\right)-\nabla\cdot\left(\varepsilon
D^{*}\nabla\overline{c}\right)+\overline{R}_{f}=0,$ (13)

where $\varepsilon D^{*}$ (m2/s) is the effective diffusion coefficient and
$\overline{R}$ (kmol/m3/s) is the volume-averaged surface reaction rate. The
effective diffusion coefficient takes into account the reduction of the total
diffusion due to the presence of solid phase. In this paper, we take
$D^{*}=D$. The volume-averaged surface reaction rate is defined as

$\overline{R}=\frac{1}{V}\int_{A}k_{c}cdS,$ (14)

where $A=V\cap\Gamma$ is the reactive surface area in the control volume. The
specific surface area $a_{s}$ (m-1) in a control volume is defined as

$a_{s}=\frac{1}{V}\int_{A}dS.$ (15)

Finally, the mass balance equation for the solid phase writes

$\frac{\partial\varepsilon}{\partial t}=\overline{R}\frac{M_{ws}}{\rho_{s}}.$
(16)

and the volume averaged velocity satisfies

$\nabla\cdot\overline{\mathbf{u}}=\overline{R}M_{ws}\left(\frac{1}{\rho_{s}}-\frac{1}{\rho}\right).$
(17)

###  2.3 Improved Volume-of-Solid

The method presented here is analogue to the calculation of the mass transfer
across a multiphase interface presented in [39], for which the mass transfer
is calculated as the scalar product between a diffusive flux and the gradient
of the phase indicator function. To calculate the volume-averaged surface
reaction rate, the improved Volume-of-Solid (iVoS) therefore introduces the
reactive flux, $\mathbf{\Phi}_{R}$ (kmol/m2/s), defined as

$\mathbf{\Phi}_{R}=k_{c}c\mathbf{n}_{s},$ (18)

and the volume-averaged surface reaction rate can be rewritten as

$\overline{R}=\frac{1}{V}\int_{A}\mathbf{\Phi}_{R}\cdot\mathbf{n}_{s}dS.$ (19)

Assuming that the concentration of the reactant on the reactive surface can be
approximated by its volume-averaged on the control volume, and that the normal
vector to the interface can be approximated by

$\overline{\mathbf{n}}_{s}=-\frac{\nabla\varepsilon}{\|\nabla\varepsilon\|},$
(20)

the reactive flux can be approximated by

$\overline{\mathbf{\Phi}}_{R}=k_{c}\overline{c}\overline{\mathbf{n}}_{s}.$
(21)

Moreover, the average surface normal in a control volume can be calculated as
[6]

$\frac{1}{V}\int_{A}\mathbf{n}_{s}dS=-\nabla\varepsilon.$ (22)

Therefore, the volume-averaged surface reaction rate can be calculated as

$\overline{R}=-\overline{\mathbf{\Phi}}_{R}\cdot\nabla\varepsilon.$ (23)

To avoid problems related to the calculation of the gradient of $\varepsilon$
(see section 2.4), the divergence theorem is used to recast $\overline{R}$ as

$\overline{R}=\varepsilon\nabla\cdot\overline{\mathbf{\Phi}}_{R}-\nabla\cdot\left(\varepsilon\overline{\mathbf{\Phi}}_{R}\right),$
(24)

With this formulation, the reactive rate is the sum of two terms, an overall
mass transfer term
($\sum\overline{R}_{f}=\sum\left(\varepsilon\nabla\cdot\overline{\mathbf{\Phi}}_{R}\right)$)
and a conservative term
($-\sum\left(\nabla\cdot\left(\varepsilon\overline{\mathbf{\Phi}}_{R}\right)\right)=0$),
which balances the local reaction rate between two adjacent control volumes.
This means that the reaction can consume reactant in one cell and use it to
dissolve solid in a neighbor cell. By using a second-order high resolution
difference scheme [40], the reaction rate is balanced toward the reactive
surface and the diffusion of the solid interface is limited. This is an
accurate representation of a reaction at an interface between two cells with
$\varepsilon\approx 0$ and $\varepsilon=1$, where the reaction rate is
calculated using the concentration in the fluid where $\varepsilon=1$ and the
reaction dissolves the solid where $\varepsilon\approx 0$. This formulation is
labelled iVoS.

###  2.4 Volume-of-Solid with localization function

As an alternative to Equ. (24), the volume-averaged surface reaction rate can
be calculated as

$\overline{R}=k_{c}\overline{c}a_{s}.$ (25)

The specific surface area can be directly calculated as
$a_{s}=~{}\|\nabla\varepsilon\|$. However, this can lead to a diffuse
interface that spreads across a large number of layers in the computational
grid. To enforce localization of the dissolution front on the fluid-solid
interface, the VoS-$psi$ method introduces a diffuse interface localization
function $\psi$ [16] so that

$\overline{R}=k_{c}\overline{c}\psi\|\nabla\varepsilon\|.$ (26)

The main advantage of VoS-$\psi$ compared to iVoS is that it provides a direct
calculation of the reactive surface area in a control volume. Therefore, the
reaction rate can be calculated by a dedicated geochemical solver, such as
Phreeqc [41] or Reaktoro [42], as it is often done for standard Reactive
Transport Modelling dedicated to multi-scale applications [36].

While the VoS-$\psi$ method has proved to be a fast and flexible method to
match experimental results [25], this formulation has one main limitation. It
is strongly dependent on the choice of the localization function $\psi$.
Several functions have been proposed by [38], but their accuracy depends on
the case considered and on the discretization scheme used for the gradient.
For example, using a centered difference scheme requires $\psi(1.0)=0$ for
stability. However, using a decentered scheme (in the direction of
$\varepsilon=0$ to avoid instabilities) will result in a higher and more
diffuse reaction rate, due to a higher reactant concentration away from the
surface and a larger specific surface area when $\varepsilon=0$. Centered
difference schemes are less diffuse, but they result in incomplete
dissolution, since a cell with $\varepsilon<1$ but $\varepsilon=1$ for all its
neighbors will have a zero reaction rate. Currently, there is no consensus on
the ideal combination of localization function and discretization scheme to
use for every scenario, as this will depend on the geometry and flow
conditions.

In this work, we use a centered difference scheme for the gradient and
$\psi=\lambda\varepsilon(1-\varepsilon)$, which is the most accurate
combination of discretization scheme and localization function proposed by
[38] for the case of dissolution of a calcite post by acid injection [16].
Typically, $\lambda=4$, but since $\psi\leq 1.0$, this will lead to a
reduction in interfacial area. For this reason, the VoS-$\psi$ method gives an
overall lower dissolution rate that the iVoS method. Instead, $\lambda$ can be
calculated as a function of $\varepsilon$ so that

$\lambda\left(\varepsilon\right)=\frac{\int_{\Omega}\|\nabla\varepsilon\|}{\int_{\Omega}\|\nabla\varepsilon\|\varepsilon\left(1-\varepsilon\right)}.$
(27)

In this work, we label VoS-$\psi$ the formulation using
$\psi=4\varepsilon\left(1-\varepsilon\right)$ and VoS-$\psi$’ the formulation
using
$\psi=\lambda\left(\varepsilon\right)\varepsilon\left(1-\varepsilon\right)$.
Using the VoS-$\psi$’ formulation, the total surface area is conserved
globally. However, it is not conserved locally. At an interface between two
cells with $\varepsilon\approx 0$ and $\varepsilon=1$, the interface area will
remain close to 0 in both cells even after the correction.

###  2.5 Upscaling to the Darcy scale

In order to investigate the capabilities of our numerical model to calculate
upscaled properties for macro-scale simulations, the flow and reaction in the
whole domain are characterised by the the total porosity $\phi$, the Darcy
velocity $U_{D}$ (m/s),

$U_{D}=\frac{Q_{i}}{A_{i}},$ (28)

and the permeability $K$ (m2),

$K=\frac{\nu U_{D}L_{D}}{\Delta P}.$ (29)

where $Q_{i}$ (m3/s) is the inlet flow rate, $A_{i}$ (m2) is the inlet area,
$L_{D}$ (m) is the distance between the inlet and outlet, $\Delta P$ (m2/s2)
is the kinematic pressure drop. In addition, the chemical reaction is
characterised by the total specific surface area $a$ (m-1)

$a=\frac{1}{V_{\Omega}}\int_{\Gamma}dS,$ (30)

and the correction factor $\alpha$, that represents the reduction of the
reactive surface area accessible to reactant, and is defined as

$\alpha=\frac{1}{a}\frac{\int_{\Gamma}cdS}{\int_{\Omega_{f}}cdV}.$ (31)

Change in effective upscaling parameters $K$, $a$ and $\alpha$ are often
modelled as a function of porosity with power law functions [37, 43, 10]. The
parameters of these power-law functions depend strongly on the flow, transport
and reaction conditions, which are characterized by the Péclet number

$\displaystyle Pe=\frac{UL}{D},$ (32)

which quantifies the relative importance of advective and diffusive transport,
and the Damköhler number

$\displaystyle Da=\frac{k_{c}}{U},$ (33)

which quantifies the relative importance of chemical reaction and advective
transport. Here $U$ and $L$ are the reference velocity and length. The product
of the Damköhler and Péclet numbers is also a relevant quantity called the
Kinetic number, defined as

$\displaystyle Ki=DaPe=\frac{k_{c}L}{D}.$ (34)

In addition, the reactant strength is characterized by

$\displaystyle\beta=\frac{c_{i}M_{ws}}{\zeta\rho_{s}},$ (35)

where $c_{i}$ is the concentration of reactant at the inlet.

##  3 Benchmark cases

In this section, the numerical models are benchmarked based on experimental
results and simulation results using an interface tracking method based on the
ALE framework. All numerical methods are implemented in GeoChemFoam, our
reactive transport toolbox, and their implementation is presented in
appendices. In order to save on computational time, the equations are
implemented using the quasi-static assumption presented in A, using that, for
all our test cases, $\beta Da<<1$ and $\beta Ki<<1$. The solution procedures
are presented in B. The domains are meshed using Adaptive Mesh Refinement
(AMR) for the VoS methods and Local Mesh Refinement (LMR) for the ALE methods
presented in C, and adaptive time-stepping strategies presented in D. Three
benchmark test cases are considered. In the first benchmark case, the
numerical models are used to simulate the dissolution of a 3D calcite post in
a straight microchannel, and the results are compared with experimental
results [16]. Convergence, accuracy and efficiency of the methods are
compared. In the second benchmark, the methods are used to simulate
dissolution in a 2D micromodel at various dissolution regimes. The accuracy
and efficiency of the methods are compared using the ALE method as a
reference. The capability of each model to calculate upscaled coefficients is
then studied. In the third test case, the iVoS method is used to simulate
dissolution in a 3D micro-CT image of Ketton carbonate. The ALE method could
not be used in this case due to its computational cost. The accuracy of the
methods is compared based on experimental results [20] and simulation results
[44] using solid balance with a threshold. Two additional simulations at two
different dissolution regimes are then run and the accuracy of upscaling laws
[10, 11, 43] are explored.

###  3.1 Benchmark 1: calcite post dissolution - comparison with experiment

In benchmark 1, we dissolve a calcite post with all four numerical methods
(ALE, VoS-$\psi$, VoS-$\psi^{\prime}$ and iVoS) and compare the results with
the experimental data from [16]. In the experiment, an octagonal-shaped
calcite post is placed at the center of a straight microchannel and is
dissolved by an acidic solution that is flowing past. The experiment and
methods are described in detail in [16].

Simulations are performed using the VoS methods, and using the ALE method for
reference. The domain is a straight microchannel of size 2.67 mm $\times$ 1.5
mm $\times$ 0.2 mm with a calcite post of height 0.2 mm in its center. Images
of the post (in stl and h5 formats) are given in the supplementary material.
In each case, a cartesian mesh of resolution $\Delta x$=20 $\mu$m is
generated, which is snapped onto the calcite solid surfaces for the ALE method
(C). Increased resolutions of $\Delta x=10$ $\mu$m and $\Delta x=5$ $\mu$m at
the interface are obtained using LMR or AMR (C). The initial meshes have
respectively 95,900, 109,476 and 160,390 cells for the ALE simulations and
105,000, 114,430 and 174,280 for both of the micro-continuum simulations.

At t=0, a solution of hydrochloride acid is injected from the left boundary at
constant flow rate, extrapolated from a zero-gradient pressure velocity field
[45], and constant concentration $c_{i}=0.0126$ kmol/m3. The acid reacts with
the calcite surface to produce $Ca^{2+}$ and $HCO_{3}^{-}$. However, due to
the very low pH, $HCO_{3}^{-}$ reacts instantaneously with $H^{+}$ to give
$H_{2}CO_{3}^{*}$. The two reactions can be added and modelled as the single
reaction described below.

Parameter | Symbol | Value | Unit
---|---|---|---
Kinematic viscosity | $\nu$ | 2.61$\times 10^{-6}$ | m2/s
Diffusion coefficient | $D$ | $5\times 10^{-9}$ | m2/s
Inlet flow rate | $Q_{i}$ | 3.5$\times 10^{-10}$ | m3/s
Inlet acid concentration | $c_{i}$ | 0.0126 | kmol/m3
Reaction constant | $k_{c}$ | 8.9125$\times 10^{-4}$ | m/s
Stoichiometric coefficient | $\zeta$ | 2 | (-)
Calcite molecular weight | $M_{ws}$ | 100 | kg/kmol
Calcite density | $\rho_{s}$ | 2710 | kg/m3
Kozeny-Carman constant | $k_{0}$ | $10^{-12}$ | m2
Table 1: Simulation parameters for Benchmark 1.

$\displaystyle CaCO_{3}+H^{+}$ $\displaystyle\leftrightharpoons
Ca^{2+}+HCO_{3}^{-}$ (36) $\displaystyle HCO_{3}^{-}+H^{+}$
$\displaystyle\leftrightharpoons H_{2}CO_{3}$ $\displaystyle CaCO_{3}+2H^{+}$
$\displaystyle\leftrightharpoons Ca^{2+}+H_{2}CO_{3}^{*}$

The simulation parameters are summarized in Table 1. Each simulation is run
until $t=12000$ s or until all of the solid has been dissolved, whichever
happens first.

Figure 1: Concentration map at various times during dissolution of a calcite
post by acid injection using the ALE, VoS-$\psi$, VoS-$\psi$’ and iVoS
methods.

Fig. 1 shows the concentration map at different times for the different
methods, with $c=0$ in the solid phase. We see that for the ALE and iVoS
methods, the concentration maps are very similar and the calcite post has
fully disappeared at t=12000 s. However, for the VoS-$\psi$ method, the
dissolution is delayed and there is still a small but significant volume of
calcite at t=12000 s. The VoS-$\psi$’ method corrects most of this error and
the concentration maps are similar to the one obtained with the ALE and iVoS
methods.

Figure 2: Evolution of (a) the grain volume and (b) the grain surface area
obtained experimentally (black) and with simulation using the ALE (red),
VoS-$\psi$ (blue), , VoS-$\psi$’ (purple) and iVoS (green) methods at various
mesh resolutions during dissolution of a calcite post by acid injection.

The results of numerical simulations at different mesh resolutions are
compared with experimental results from [16] in Fig. 2. We observe that the
ALE method converges toward a solution close to the experimental results. The
difference between the experimental results and the numerical simulation at
mesh resolution $\Delta x=5$ $\mu$m in grain volume and grain area are less
than 1% of the initial values. Although the iVoS method gives a significantly
lower grain volume and surface area than the experiment at a resolution
$\Delta x=20$ $\mu$m, the results are very similar to the ALE results for
$\Delta x=10$ $\mu$m and $\Delta x=5$ $\mu$m. In addition, with the iVoS
method as with the ALE method, the difference between the experimental results
and the numerical simulation at mesh resolution $\Delta x=5$ $\mu$m in grain
volume and grain area are less than 1% of the initial values.

Although the VoS-$\psi$ method matches the trend of the experiment, it
overestimates the grain volume and surface area for all resolutions. The
results do not improve as the mesh resolution increases and converge toward a
solution with an error of 5% in the grain volume and 10% in the grain area.
This error is corrected by the surface area correction provided by the
VoS-$\psi$’ method. At all resolution, the VoS-$\psi$’ method gives errors in
grain volume and grain area that are less than 1% of the initial values.
During the simulation, $\lambda\left(\varepsilon\right)$ is a value between 8
and 16 that changes as $\varepsilon$ changes. This shows that the VoS-$\psi$
underestimates the overall surface area by a factor between 2 and 4.

| 20$\mu$m | 10$\mu$m | 5 $\mu$m
---|---|---|---
ALE | 39 | 78 | 147
VoS-$\psi$ | 11 | 22 | 62
VoS-$\psi^{\prime}$ | 15 | 30 | 64
iVoS | 16 | 31 | 65
Table 2: CPU time [min] obtained with the ALE, VoS-$\psi$, VoS-$\psi$’ and
iVoS methods at various mesh resolutions during dissolution of a calcite post
by acid injection.

Table 2 shows the CPU times for all simulations. The simulations are
significantly faster using the VoS methods, with the iVoS being approximately
2.5x faster than the ALE for each simulation. The VoS-$\psi$ method is
slightly faster than the iVoS and VoS-$\psi$’ methods, but this is mostly due
to the larger time steps resulting from the lower dissolution rates.

We conclude that all VoS methods are significantly faster than the ALE method.
The iVoS provides a result that converges toward a solution with an error less
than 1% when the mesh resolution reaches 5 $\mu$m. The VoS-$\psi$ results in a
small error which is due to the reduction of the overall surface area and is
corrected in the VoS-$\psi$’ method by fitting the constant $\lambda$ for each
time-step. The VoS-$\psi$’ method gives accurate results for all mesh
resolution.

###  3.2 Benchmark 2: dissolution regimes in a 2D model

Figure 3: Micromodel geometry for Benchmark 2.

In benchmark 2, we use the three VoS methods to model the various dissolution
regimes (i.e compact, wormholes, and uniform) that occur during mineral
dissolution in a 2D porous media model. ALE simulations are also run for
reference. The model is constructed from a homogeneous domain with discs
radius 270 $\mu$m by adding a random deviation of magnitude 270 $\mu$m in disc
radius and center position. The micromodel generation code is available open
source (https://github.com/hannahmenke/DrawMicromodels) and the method is
described in [46]. The geometry is presented in Fig. 3 and a high resolution
image can be found in the supplementary material.

The domain is meshed with a cartesian mesh with uniform resolution $\Delta
x=3$ $\mu$m, which is snapped on the solid surfaces for the ALE method. A band
of two cells width is added on each side of the model to avoid dissolution at
the boundaries. The final meshes include 1,008,016 cells for the VoS methods
and 457,455 for the ALE method. The porosity $\phi$ and the permeability $K$
of the full domain can be numerically calculated as $\phi=0.45$ and
$K=5.6\times 10^{-10}$ m2.

Parameter | Symbol | Value | Unit
---|---|---|---
Kinematic viscosity | $\nu$ | $10^{-6}$ | m2/s
Diffusion coefficient | $D$ | $1\times 10^{-9}$ | m2/s
Inlet acid concentration | $c_{i}$ | 0.01 | kmol/m3
Stoichiometric coefficient | $\zeta$ | 1 | (-)
Calcite molecular weight | $M_{ws}$ | 100 | kg/kmol
Calcite density | $\rho_{s}$ | 2710 | kg/m3
Kozeny-Carman constant | $k_{0}$ | 1.79$\times 10^{-11}$ | m2
Table 3: Simulation parameters for Benchmark 2.

At t=0, acid is injected from the left boundary at constant flow rate,
extrapolated from a zero-gradient pressure velocity field [45], and flows out
of the domain from the right boundary at constant pressure. The top and bottom
boundaries are no-flow, no-slip conditions. The fluid and solid properties are
summarized in Table 3. The Kozeny-Carman constant is fitted to obtain the same
permeability as in the direct method at $t=0$. For each simulation, the inlet
flow rate $Q_{i}$ and the chemical reaction constant $k_{c}$ are adapted to
obtain the correct $Pe$ and $Ki$, using the pore-scale length
$L=\sqrt{12K/\phi}$ as the reference length and the average pore velocity
$U=U_{D}/\phi$ as the reference velocity. The factor 12 is added so that the
pore-scale length corresponds to the channel size for an homogeneous bundle of
straight channels [47]. Four cases are considered that characterize four
different dissolution regimes: $Pe=0.01$ and $Ki=0.1$ (compact dissolution),
$Pe=0.3$ and $Ki=10$ (conical wormhole), $Pe=1$ and $Ki=1$ (dominant
wormhole), and $Pe=10$ and $Ki=0.01$ (uniform dissolution). The simulations
are performed until 20% of the solid has been dissolved.

Figure 4: Acid concentration map after 20% dissolution in a micromodel for
four different regimes (Benchmark 2) obtained by numerical simulations. Figure
5: Evolution of the solid volume during numerical simulation of dissolution in
a 2D micromodel at four different regimes.  Figure 6: Evolution of the solid
surface area during numerical simulation of dissolution in a 2D micromodel at
four different regimes.

Fig. 4 shows the concentration map at the end of the simulations. We observe
that all methods are able to model all regimes qualitatively. In particular,
the iVoS method reproduces the ALE results almost exactly. However, there are
several inaccuracies in the VoS-$\psi$ method. Due to the use of a centered
scheme for the computation of the gradient of $\varepsilon$, the VoS-$\psi$
method results in several grains with incomplete dissolution in the compact
and wormhole regimes. This problem can be resolved by applying a decentered
gradient, but this is done at the expense of accuracy, as a decentered
gradient generates more numerical diffusion. In addition, the wormhole
obtained with VoS-$\psi$ for $Pe=1$ and $Ki=1$ is more diffused and ramified
that the ones obtained with ALE and iVoS, due to the reduction of reaction
rate induced by $\psi$. For the uniform regime, the patterns are almost
identical, but the acid concentration obtained with the VoS-$\psi$ method is
higher than with ALE and iVoS, suggesting that the reaction rate is lower for
VoS-$\psi$. Most of these errors are corrected in the VoS-$\psi$’ method,
although there are still grains with incomplete dissolution in the compact and
wormhole regime, and the wormhole is still slightly more diffuse and ramified
for $Pe=1$ and $Ki=1$.

Fig. 5 and 6 show the evolution of the solid volume and solid surface area for
all regimes obtained with all numerical methods. We observe that the iVoS
method reproduces the ALE results with good accuracy for all regimes. For all
cases, the solid volume and surface area evolutions obtained with the ALE and
the iVoS methods are similar. The VoS-$\psi$ method underpredicts the amount
of dissolution occurring in the systems compared to the ALE and the iVoS
methods. This is particularly true for the cases with $Ki\leq 0.1$, i.e. in
the compact and uniform regime. In these cases, the concentration gradient on
the reactive surface is small and the reaction rate is mostly dependent on the
surface area, which is significantly reduced by the localization function
$\psi$. For the dominant wormhole regime, although the solid volume evolution
obtained with the VoS-$\psi$ method is similar to the ones obtained with the
ALE and iVoS methods, the evolution of the solid area is significantly
different. This suggests that, although the total amount of dissolution is
correctly predicted, it occurs in a slightly different pattern than with the
ALE and iVoS methods. The dissolution front is less sharp, and the reaction
occurs in the vicinity of the wormhole rather than at the tip of the wormhole
as predicted by the ALE and iVoS methods. This leads to a more ramified
dissolution pattern, as observed in Fig. 4. The VoS-$\psi$’ method corrects
some of these errors and the overall amount of dissolution is similar to the
ALE and iVoS results for all cases. However, the surface area is slightly
higher for all cases. This is because, since the dissolution is high when
$\varepsilon\approx 0.5$ and low when $\varepsilon\approx 0$ and
$\varepsilon\approx 1$, the interface becomes artificially sharp and leads to
a larger interfacial area.

Figure 7: Evolution of the permeability and macro-scale reaction constant
during numerical simulation of dissolution in a 2D micromodel at four
different regimes obtained with the ALE, the VoS-$\psi$ and the iVoS methods.

Fig. 7 shows the evolution of the permeability $K$ and the correction factor
$\alpha$ as a function of porosity for all regimes obtained with the four
methods. In all cases, the order of the permeability evolution is similar for
the three methods, except for the wormholing regimes, for which the order is
19 for ALE and 17 for iVoS, typical of the dominant wormhole regime, but only
8 for VoS-$\psi$ and 12 for VoS-$\psi$’, which is more typical of the ramified
wormhole regime. The macro-scale reaction constant is consistently lower with
VoS-$\psi$ than with ALE and iVoS, especially for the uniform regime where it
is more than twice as small.

| $Pe=0.01$ | $Pe=0.3$ | $Pe=1$ | $Pe=10$
---|---|---|---|---
| $Ki=0.1$ | $Ki=10$ | $Ki=1$ | $Ki=0.01$
ALE | 69 | 71 | 46 | 18
iVoS | 5.8 | 6.2 | 4.6 | 4.3
VoS-$\psi$ | 3.8 | 3.8 | 3.3 | 2.7
VoS-$\psi$’ | 7.1 | 7.3 | 6.2 | 6.3
Table 4: CPU time (in hours) obtained for all three methods for all four
regimes during dissolution in a 2D micromodel.

Table 4 shows the CPU time for all cases for all three methods. We observe
that the VoS methods are between 3 and 12 time faster than the ALE method. The
VoS methods are particularly efficient compared to the ALE method for the
cases with localized dissolution front, i.e. for the compact and the
wormholing regime, for which the ALE method performs a large number of
remeshing steps. The VoS-$\psi$ method is slightly faster than the iVoS
method, but this is mostly due to larger time steps resulting from the lower
dissolution rates. The VoS-$\psi$’ method is slightly slower than the iVoS
method, which we attribute to a slower convergence of the transport equation
due to a more localized dissolution.

We conclude that the iVoS method is capable of modelling all regimes during
dissolution in a 2D micromodel and calculating macro-scale coefficients with
similar results as the ones obtained with ALE. In addition, the iVoS method is
significantly faster than the ALE method. The VoS-$\psi$ method, using
$\psi=4\varepsilon_{f}(1-\varepsilon_{f})$ underpredicts the dissolution in
all cases due to a reduction of overall surface area, but this error is mostly
corrected by using the VoS-$\psi$’ method.

###  3.3 Benchmark 3: dissolution in a 3D micro-CT image and comparison with
experiment

In benchmark 3 we simulate the dissolution of a 3D micro-CT image of Ketton
limestone and compare the numerical results to in-situ dissolution
experiments. The main advantage of the micro-continuum approach is its
computational efficiency compared to the ALE method. In this part, we take
advantage of this to perform simulations with a VoS method, while the ALE
method is too computationally expensive to perform on a 3D image of this size.
The iVoS method was selected since it gave the most accurate results in
benchmark 2. The experiment is the one conducted in [20], where CO2-saturated
brine is injected in a Ketton carbonate core. A core of 4 mm diameter was
flooded with a brine solution representative of a typical saline aquifer
consisting of 1% KCl and 5% NaCl by weight, pre-equilibrated with
supercritical CO2 at 10 MPa and 50 °C. Micro-CT images were acquired after 17,
33, 50 and 67 min.

Figure 8: Initial condition for simulation of dissolution in a 3D micro-CT
image of carbonate: (a) initial local porosity field;(b) Initial calculated
velocity field; (c) Initial calculated concentration field ($Pe=190$,
$Da=4.4\times 10^{-2}$).

The simulated conditions are similar to the ones presented in [44]. The sample
is a cylinder of radius 1.7mm and length 3.5mm. The image has $911\times
902\times 922$ voxels with resolution 3.8 $\mu$m, however, we ran it on two-
level adpative mesh with maximum resolution of 7.6 $\mu$m to decrease
computation time. The initial grid includes 17 million cells. Fig. 8a shows
the initial porosity field. The velocity field is then initialised by solving
Equ. (17) and (11) in the domain using $k_{0}=2\times 10^{-13}$ m2 which was
fitted to obtain a permeability of 16D. The streamlines and corresponding
magnitude of velocity are shown on Fig. 8b. The reaction rate is described in
terms of the concentration of calcium cations Ca2+, following
$R=k_{eff}\left(C_{eq}-C\right)$, where $C_{eq}$ and $k_{eff}$ depend on the
concentrations of $CO_{3}^{2-}$, $H^{+}$ and $CO_{2}$. Because our model only
tracks one component, the reaction constant and equilibrium concentration are
fitted to match the reaction-limited constant obtained for a flat pure crystal
of calcite [48] $R=8.1\times 10^{-7}$ kmol/m2/s and the initial reaction rate
in the experiment $R=8.8\times 10^{-8}$ kmol/m2/s. Such rates are obtained
using a reactant concentration of 0.0035 kmol/m3 at the inlet and a reaction
constant $k_{c}=2.314\times 10^{-4}$ m/s. Fig. 8c shows the concentration
field in the pores at $T=0$ s obtained by solving Equ. (43).

Parameter | Symbol | Value | Unit
---|---|---|---
Kinematic viscosity | $\nu$ | $2.61\times 10^{-6}$ | m2/s
Diffusion coefficient | $D$ | $7.5\times 10^{-10}$ | m2/s
Inlet flow rate | $Q_{i}$ | 8.33$\times 10^{-9}$ | m3/s
Inlet concentration | $c$ | 0.0035 | kmol/m3
Reaction constant | $k_{c}$ | 2.314$\times 10^{-4}$ | m/s
Stoichiometric coefficient | $\zeta$ | 1 | (-)
Calcite molecular weight | $M_{ws}$ | 100 | kg/kmol
Calcite density | $\rho_{s}$ | 2390 | kg/m3
Kozeny-Carman constant | $k_{0}$ | 2$\times 10^{-13}$ | m2
Table 5: Simulation parameters for Benchmark 3.

The simulation parameters are summarized in Table 5. We use for reference
velocity the average pore velocity defined as

$U=\frac{U_{D}}{\phi}.$ (37)

The reference length is associated to the initial permeability so that

$L=\sqrt{\left(\frac{8K}{\phi}\right)},$ (38)

where the constant 8 is added in order to obtain the throat radius for a
capillary bundle of uniform size. With these definitions, we obtain
$U=5.3\times 10^{-3}$ m/s and $L=27$ $\mu$m. Using the diffusion coefficient
and reaction constant in Table 5, we obtain $Pe=190$ and $Da=4.4\times
10^{-2}$.

Figure 9: Experimental and simulated results showing the evolution of the
dissolution pattern in a 3D micro-CT image of Ketton Carbonate. The colors
represent different time intervals.

The simulation is run with an adaptive time-step until $t=4000$ s on 128 CPUs
using Oracle cloud computing. The total CPU time was 65 hours. Fig. 9 shows a
comparison of the evolution of the dissolution pattern between experiment and
simulation. The colors show the part of the rock that is dissolved for each
time interval. We observe that the simulation captures the correct patterns,
with pores being enlarged in the direction of the flow. Some expected
differences are observed, primarily due to uncertainties in the experimental
conditions, segmentation error due to reaction occurring during image
acquisition, and the fact that the simulation is performed on a sub-image of
the full sample used in the experiment. Furthermore, some of the differences
in the local porosity at different times in the experiment are due to
imperfect alignment of the images. Table 6 shows the evolution of the porosity
during dissolution for both experiment and simulation. We observe that the
porosity is accurately predicted by the simulation until $T>3000$ s, at which
time the simulation starts diverging from the experiment slightly. This can be
explained by the fact that the simulation is done on a sub-image of the full
sample, starting at 2 mm away from the inlet [20], and thus the concentration
of acid will have varied as dissolution occured in the unimaged portion of the
core.

| Simulation | Experiment
---|---|---
Time | $\phi$ |  $K$ (D) | $\phi$ |  $K$ (D)
0 | 0.173 | 16 | 0.173 | 16
1000 | 0.204 | 33 | 0.204 | 37
2000 | 0.228 | 57 | 0.227 | 62
3000 | 0.251 | 89 | 0.248 | 93
4000 | 0.273 | 129 | 0.267 | 129

Table 6: Experimental and simulated porosity and permeability at different
time during dissolution in a 3D micro-CT image of Ketton Carbonate.

Figure 10: Permeability as a function of evolving porosity during dissolution
in a 3D micro-CT image of Ketton Carbonate.

| Simulation | Experiment
---|---|---
Time (s) |  $R_{\Omega}$ (mol/m3/s) |  $C$ (mol/m3) | $\alpha$ |  $R_{\Omega}$ (mol/m3/s) | $\alpha^{*}$
0-1000 | 0.74 | 1.2 | 0.35 | 0.74 | 0.35
1000-2000 | 0.57 | 1.3 | 0.23 | 0.55 | 0.23
2000-3000 | 0.55 | 1.4 | 0.21 | 0.50 | 0.19
3000-4000 | 0.52 | 1.6 | 0.18 | 0.46 | 0.16
Table 7: Evolution of experimental and simulated reaction rate, reactant
concentration and macro-scale reaction constant. The experimental macro-scale
constant $\alpha^{*}$ is calculated using the simulated concentration.

In addition, permeability can be calculated during the simulation. Table 6
shows the porosity and permeability at different times for the simulation and
the experiment and we observe a good correspondence. The permeability is then
plotted as a function of the porosity in Fig. 10. We see that the permeability
in the simulation follows the same trend that of the experiment. Finally, the
average macro-scale reaction rate for each time interval can be calculated as

$R_{\Omega}=\frac{\rho_{s}\Delta\phi}{M_{w}\Delta t},$ (39)

for both simulation and experiment. The correction factor $\alpha$ can then be
calculated by dividing by the average concentration in the pore space,
calculated at the end of the time interval from the simulation results. The
values are summarized in Table 7. We observe a good correspondence between
simulation and experiment. The overall trend of decreasing reaction rates and
correction factors observed in the experiment is reproduced in the simulation,
albeit with a slightly slower rate. This decrease is not a typical
characteristic of the uniform regime, and is an indication that the
dissolution might be occuring in the channeling regime, identified in [11].

Two additional simulations are run in different regimes by dividing the
reaction rate by 100 and 10,000, which gives Péclet number of 1.9 and 0.019.
The simulations are run on 128 CPUs using Oracle cloud computing until the
porosity reaches an approximate value of 0.28. The CPU time was 75 hours for
$Pe=1.9$ and 80 hours for $Pe=0.019$. Figure 11 shows the reactant
concentration along reconstructed streamlines at the end of the simulation. At
$Pe=0.019$, the dissolution is in the compact regime. The reactant is consumed
and dissolves the rock close to the inlet face. At $Pe=1.9$, the dissolution
is in the wormholing regime. The reactant penetrates further in the domain and
flow instabilities kick in, leading to a preferential dissolution pathway. The
streamlines and reactant concentration for $Pe=190$ are also shown in Figure
11 and the dissolution appears to be in the uniform regime.

Figure 11: Reactant concentration along streamlines during Ketton dissolution
in three different regimes, shown when the porosity has reached
approximatively 0.28. Figure 12: Evolution of permeability and correction
factor as a function of porosity and corresponding correlations during Ketton
dissolution in three different regimes.

Figure 12a shows the evolution of permeability as a function of porosity
during the simulations for the three regimes. The permeability of the compact
and uniform can be fitted with power laws, and we obtain $K=30\phi^{0.36}$ for
the compact regime and $K=5.5\times 10^{4}\phi^{4.7}$ for the uniform regime.
For the wormhole regime, the permeability curve has a strong inflection point
around $\phi=0.25$, which corresponds to porosity at reactant breakthrough
time. The permeability curve can then be fitted with two power law curves,
$K=74\phi^{2.1}$ before breakthrough and $K=3.4\times 10^{10}\phi^{15}$ after
breakthrough. The orders of the uniform and wormhole regimes correspond those
observed in the literature [11, 49].

Similarly, Figure 12b shows the evolution of the correction factor $\alpha$ as
a function of porosity during the simulations for the three regimes. For each
regime, the evolution of $\alpha$ can be fitted with a power law [43] and we
obtain $\alpha=1.7\times 10^{-6}\phi^{-5.4}$ for the compact regime,
$\alpha=3.1\times 10^{-4}\phi^{-2.8}$ for the wormholing regime and
$\alpha=7.8\times 10^{-3}\phi^{-2.4}$ for the uniform regime. We observe that
the correction factor for the uniform regime is not constant but decreases
with a power-law. This indicates that the dissolution might not be in the
uniform regime, but could be in the channeling regime identified in [11].

We conclude that the accuracy and computational efficiency of the iVoS method
enables simulation of dissolution in a 3D micro-CT image of a real carbonate
sample, reproduce experimental results with good precision and can be used to
investigate upscaling parameters.

##  4 Conclusion

We have presented two novel numerical methods, iVoS and VoS-$\psi$’, to
simulate mineral dissolution in real pore-scale geometries that when compared
to existing methods has two main advantages. First, they are based on the
micro-continuum approach, and therefore do not require a complex algorithm for
interface tracking, nor any special treatment if topological changes occur.
Second, the iVoS method calculates a reaction rate based on the divergence of
a reactive flux, and thus does not require an interface localization function,
which greatly improved its accuracy. The VoS-$\psi$’ uses a localization
function with a constant that is fitted to ensure that the reactive surface
area is conserved globally and therefore avoids most of the numerical errors
present with VoS-$\psi$. The advantages of these methods were demonstrated in
three benchmark cases.

In benchmark 1, the iVoS method was used to simulate the dissolution of a 3D
calcite post in a straight microchannel, and the results were compared with
experimental results and with simulations obtained with the ALE method and
with the standard micro-continuum approach based on the VoS-$\psi$ method. We
observed that the iVoS method was capable of reproducing the experimental
results with similar accuracy but significantly less computational time than
the ALE method, while the VoS-$\psi$ method using
$\psi=4\varepsilon_{f}(1-\varepsilon_{f})$ showed significantly more error,
but this error was corrected by the VoS-$\psi$’ approach. All VoS methods were
significantly faster than the ALE method.

In benchmark 2, we simulated dissolution in a 2D micromodel and calculated
macro-scale coefficients using all numerical methods. Four cases in three
different dissolution regimes (i.e. compact, wormholing, and uniform) were
considered and we observed a qualitative match between the ALE and iVoS
simulations, and a good correspondence between the macro-scale coefficients
calculated. The iVoS method was between 4 and 12 times faster than the ALE
method. The VoS-$\psi$ method using $\psi=4\varepsilon_{f}(1-\varepsilon_{f})$
was much faster than the ALE method, but underpredicted the dissolution and
resulted in inaccurate macro-scale coefficients. Most of these errors were
corrected by the VoS-$\psi$’ approach, but the results were still slightly
less accurate than with iVoS.

In benchmark 3, the computational efficiency of the iVoS method was used to
perform a simulation in a 3D micro-CT image of a real carbonate rock (i.e.
Ketton) and the simulation results were compared to experimental results. We
observe a good correspondence between the experimental and simulated results
for the evolution of the total porosity change, total reaction rate and
permeability with time. Two additional simulations were performed in different
dissolution regimes and the accuracy of correlations for macro-scale
coefficients evaluated. We observed that the simulated results could be
matched by power-law correlations, and that the coefficients obtained for the
permeability in the wormhole and uniform regimes correspond to what has been
observed in experiments [20, 11].

In future work, the advantages of our novel approach will be used to perform a
large number of simulations in 2D micromodels and 3D micro-CT images that will
form a large database for data-driven research. This database will then be
used to identify precisely the boundaries between the regimes and decipher the
impact of pore-scale heterogeneities. Machine-learning algorithms will also be
used to estimate macro-scale coefficients in a similar way that for single-
phase flow and reactive transport [50, 51] but extended to predict their
evolution with the porosity change. Because the micro-continuum approach can
also be used to simulate dissolution at the Darcy-scale [52, 53, 54], our
method can be extended to simulate flow, transport and dissolution in multi-
scale porous media [55]. Further, the applicabilities of the methods to model
precipitation [32] will be investigated, Finally, our approach is compatible
with the Volume-Of-Fluid and the Continuous Species Transfer methods [56] for
simulation of multiphase flow and multiphase transport with interfacial
transfer, so it has the potential to be extended to simulate mineral
dissolution during multiphase processes, which is relevant to a number of
clean-energy applications, including CO2 storage and geothermal systems [57].

##  Appendix A Quasi-static assumption

Dissolution of a solid grain is typically orders of magnitude slower than
reactant transport [58]. Over typical dimensions of pores (10 $\mu$m to 1 mm),
the diffusion time scale of the reactant ranges from 1 ms to 10 h. By
contrast, to dissolve a layer of calcite 1 mm thick takes at least 1 year
[59]. In our dimensionless analysis, this is characterised by the fact that
$\beta Da<<1$ and $\beta Ki<<1$. In this case, flow (Equ. (2)) and transport
(Equ. 5)) are assumed to be in a quasi-static state

$\nabla\cdot\left(\mathbf{u}\otimes\mathbf{u}\right)=-\nabla
p+\nu\nabla^{2}\mathbf{u},$ (40)
$\nabla\cdot\left(c\mathbf{u}\right)=\nabla\cdot\left(D\nabla c\right).$ (41)

In this case, the micro-continuum apparoach gives

$\displaystyle\nabla\cdot\left(\frac{\mathbf{u}\otimes\mathbf{u}}{\varepsilon_{f}}\right)=-\nabla\overline{p}_{f}+\frac{\nu}{\varepsilon_{f}}\nabla^{2}\overline{\mathbf{u}}-\nu
K^{-1}\overline{\mathbf{u}},$ (42)
$\displaystyle\nabla\cdot\left(\overline{c}_{f}\overline{\mathbf{u}}\right)-\nabla\cdot\left(\varepsilon
D^{*}\nabla\overline{c}_{f}\right)+\overline{R}_{f}=0$ (43)

and the volume-averaged velocity satisfies the continuity equation

$\nabla\cdot\overline{\mathbf{u}}=0.$ (44)

The quasi-static assumption allows the models to run with a large time-step
controlled only by the velocity of the solid interface (see appendix) to save
on computational time.

##  Appendix B Implementation

The iVoS and VoS-$\psi$ methods have been implemented in GeoChemFoam [60, 61,
62, 63]. GeoChemFoam is an open-source OpenFOAM®-based [45] toolbox that
contains several additional packages for solving various flow processes
including multiphase transport with interface transfer, single-phase flow in
multiscale porous media, and reactive transport with mineral dissolution. The
full code can be downloaded from
www.github.com/geochemfoam.

The full solution procedure is presented in Fig. 13a. For each time-step, the
new volume of fluid fraction is solved using Equ. (16), and the permeability
is updated (Equ. (12)). The velocity (Equ. (42)), pressure (Equ. (44)) and
concentration (Equ. (43)) equations are solved using the Semi-Implicit Method
for Pressure Linked Equation (SIMPLE) [64]. The pressure equation is obtained
by combining the continuity (Eq. (44)) and momentum (Eq. (42)) equations. The
equations are solved sequentially with under-relaxation factors of 0.3 and
iterated until convergence. At the end of each SIMPLE loop iteration, the new
volume-averaged surface reaction rate is calculated, using Equ. (24) for the
iVoS method and Equ. (26) for the VoS-$\psi$ method. The SIMPLE loop is
iterated until steady-state is reached.

Figure 13: Solution procedures for reactiveTransportDBSFoam and
reactiveTransportALEFoam

In order to compare the iVoS and VoS-$\psi$ methods with an accurate interface
tracking method, the ALE method has also been implemented in GeoChemFoam in a
separate solver called reactiveTransportALEFoam and the full solution
procedure is presented in Fig. 13b. For each time-step, the mesh points are
moved with velocity $\mathbf{w}$, which satisfies the Laplace equations with
boundary condition (Equ. (4))

$\displaystyle\nabla\cdot D_{m}\nabla w_{j}=0\hskip 14.22636pt\text{j=x,y,z}$
(45) $\displaystyle w_{j}=\mathbf{w}_{s}\cdot\mathbf{e}_{j}\hskip
14.22636pt\text{at $\Gamma$},$ (46)

where $D_{m}$ is the diffusivity of the mesh motion, $w_{j}$ is the
j-directional component and $\mathbf{e}_{j}$ is the j-directional standard
basis vector. These equations will allow the mesh points to track the fluid-
solid interface, and the mesh motion is diffused to avoid large volume ratio
between neighbor cells. However, the skewness of the mesh can increase and
become unacceptably large, which can lead to failure of the transport solver.
To avoid this, the mesh quality is checked at the end of each time-step, and
if it fails, the domain is fully remeshed and the fields are mapped to the new
mesh. To avoid topological errors that can appear when two faces of the same
mineral grain overlap, leading to failure of the flow or transport solver, the
faces which are fully located in a topological error are eliminated. These
collapsing faces are identified by the following condition: a face defined as
faceI collapsed if a ray leading from its center following its normal vector
pointing toward the solid phase meets another face defined as faceJ at a
distance lower than the grid size, and faceI and faceJ do not intersect. Then,
a new mesh can be constructed. Following this remeshing algorithm, our
numerical simulations are stable and topological errors are eliminated. After
the mesh is updated. The velocity, pressure and concentration equations (Equ.
(1), (40) and (41)) are then solved using the SIMPLE algorithm with under-
relaxation factors of 0.3 and iterated until convergence.

The equations are discretized on a collocated Eulerian grid. The space
discretization of the convection terms is performed using the second-order
vanLeer scheme [40] while the diffusion term is discretized using the Gauss
linear limited corrected scheme, which is second order and conservative. For
the iVoS method, the discretization of the reactive flux term (Equ. (24)) is
also done with the vanLeer scheme.

##  Appendix C Meshing

Figure 14: Example of adaptive mesh refinement for the micro-continuum
approach (a) and local mesh refinement for the ALE method (b) at the surface
of an octagonal grain.

For the micro-continuum approach, the solid phase is described using the
grayscale value of a raw image, encapsulated in h5 format. A cartesian mesh is
generated and the volume fraction of fluid $\varepsilon_{f}$ in each grid cell
is calculated from the image. The volume fraction of fluid in the solid phase
is set to a small value $\varepsilon_{f}=10^{-4}$ to avoid division by zero.
An Adaptive Mesh Refinement (AMR) strategy can be used to refine the mesh
automatically near the tracked interface (Fig. 14a). This method splits the
computational cell when the volume fraction of fluid satisfies
$\varepsilon_{f,min}<\varepsilon_{f}<\varepsilon_{f,max}$ [65]. The volume
fraction of fluid $\varepsilon_{f}$ is then recalculated for the refined cells
with higher resolution before the start of the simulation. During the
simulation, we refine cells for which the volume fraction of fluid becomes
$\varepsilon_{f,min}<\varepsilon_{f}<\varepsilon_{f,max}$ and merge them when
it becomes $\varepsilon_{f}\leq\varepsilon_{f,min}$ or
$\varepsilon_{f}\geq\varepsilon_{f,max}$. To save on computational time, the
mesh is only modified every $n_{ref}$ time-steps. In our simulations, we use
$\varepsilon_{f,min}=0.01$, $\varepsilon_{f,max}=0.99$ and $n_{ref}=200$.

For the ALE method, the solid surface is described using an stl image. First,
a cartesian mesh is generated. Local Mesh Refinement (LMR) can be used to
refine the mesh near the solid boundaries (Fig. 14b). The mesh is then snapped
onto the solid surface using the snappyHexMesh utility [45], i.e. cell
containing solid are then removed and replaced by hexahedral or tetrahedral
cells that match the solid boundaries. An additional layer of cells of the
same resolution $\Delta x$ is then added around the solid boundary to improve
the representation of the solid surface.

##  Appendix D Time-stepping strategy

The simulations are performed using an adaptive time-stepping strategy. For
the ALE method, this is done using the mesh Courant-Friedrich-Lewy (CFL)
number defined as

$mCFL=\frac{\mathbf{w}\Delta t}{\Delta x},$ (47)

where $\Delta t$ is the time-step and $\Delta x$ is the mesh resolution. For
the micro-continuum approach, the mCFL number corresponds a maximum change of
$\varepsilon_{f}$ during a time-step $\Delta\varepsilon_{f,max}$. Our
experience indicates that maximum $mCFL$ number and
$\Delta\varepsilon_{f,max}$ of 0.005 offers a good compromise between
accuracy, robustness and efficiency, and that is what we are using in this
paper.

## References

  * [1] J.R. Black, S.A. Carroll, R.R. Haese, Chemical Geology 399, 134 (2015)
  * [2] S.N. Pandey, A. Chaudhuri, H. Rajaram, S. Kelkar, Geothermics 57, 111 (2015)
  * [3] M.U. Shafiq, H. Ben Mahmud, Journal of petroleum exploration and production technology 7(4), 1205 (2017)
  * [4] J.M. Nordbotten, M.A. Celia, _Geological storage of CO 2: Modeling approaches for large-scale simulation_ (John Wiley & Sons, 2011)
  * [5] B.B. Williams, J.L. Gidley, R.S. Schechter, _Acidizing Fundamentals_ (Society of Petroleum Enginners, 1979)
  * [6] M. Quintard, S. Whitaker, Adv. Water Resour. 17(4), 116 (1994)
  * [7] M. Shapiro, H. Brenner, Chem. Engng Sci. 43 3, 551–571 (1988)
  * [8] P.C. Lichtner, Geochimica et Cosmochimica Acta 52(1), 143 (1988)
  * [9] C.I. Steefel, C.A.J. Appelo, B. Arora, D. Jacques, T. Kalbacher, O. Kolditz, V. Lagneau, P.C. Lichtner, K.U. Mayer, J.C.L. Meeussen, S. Molins, D. Moulton, H. Shao, J. Simunek, N. Spycher, S.B. Yabusaki, G.T. Yeh, Computational Geosciences 19(3, SI), 445 (2015)
  * [10] C. Noiriel, L. L., B. Madé, L. Raimbault, P. Gouze, J. Van Der Lee, Chemical Geology 265((1-2)), 160 (2009)
  * [11] H.P. Menke, B. Bijeljic, M.J. Blunt, Geochimica Cosmochimica Acta 204, 267–285 (2017)
  * [12] A. Raoof, H. Nick, S. Hassanizadeh, C. Spiers, Computers & Geosciences 61, 160 (2013)
  * [13] C. Varloteaux, S. Bekri, P.M. Adler, Advance in Water Resources 53, 87–100 (2013)
  * [14] S. Molins, D. Trebotich, G.H. Miller, C.I. Steefel, Water Resource Research 53(5), 3645 (2017)
  * [15] V. Starchenko, C.J. Marra, L.A.J. C., Journal of Geophysical Research: Solid Earth 121, 6421–6444 (2016)
  * [16] C. Soulaine, S. Roman, A. Kovscek, H. Tchelepi, J. Fluid Mech. 827(457-483) (2017)
  * [17] Y. Yang, Y. Li, J. Yao, S. Iglauer, L. Luquot, K. Zhang, H. Sun, L. Zhang, W. Song, Z. Wang, Water Resources Research 56(4), e2019WR026112 (2020)
  * [18] Y. Hao, M. Smith, Y. Sholokhova, S. Carroll, Advance in Water Resources 62((C)), 388–408 (2013)
  * [19] L. Luquot, O. Rodriguez, P. Gouze, Transport in Porous Media 101(3), 507 (2014)
  * [20] H.P. Menke, B. Bijeljic, M.G. Andrew, M.J. Blunt, Environmental Scence and Technology 49(7), 4407 (2015)
  * [21] H.P. Menke, M.G. Andrew, M.J. Blunt, B. Bijeljic, Chemical Geology 428, 15 (2016)
  * [22] H.P. Menke, C.A. Reynolds, M.G. Andrew, J.P. Pereira Nunes, B. Bijeljic, M.J. Blunt, Chemical Geology 481, 27–37 (2018)
  * [23] J. Nogues, J. Fitts, M. A. Celia, C. Peters, Water Resources Research 49 (2013)
  * [24] C. Soulaine, J. Maes, S. Roman, Frontiers in Water 3(643714) (2021)
  * [25] S. Molins, C. Soulaine, N.I. Prasianakis, A.J.C. Ladd, V. Starchenko, S. Roman, D. Trebotich, H.A. Tchelepi, C.I. Steefel, Computational Geoscience (2020).  DOI 10.1007/s10596-019-09903-x
  * [26] Q. Kang, D. Zhang, S. Chen, X. He, Phys. Rev. E 65, 036318 (2002).  DOI 10.1103/PhysRevE.65.036318.  URL https://link.aps.org/doi/10.1103/PhysRevE.65.036318
  * [27] P. Szymczak, A.J.C. Ladd, Journal of Geophysical Research: Solid Earth 114(B6) (2009)
  * [28] N. Prasianakis, S. Ansumali, Communications in Computational Physics 9(5), 1128–1136 (2011)
  * [29] V. Starchenko, A.J.C. Ladd, Water Resource Research 54(10), 7946 (2018)
  * [30] A.M. Tartakovsky, P. Meakin, T.D. Scheibe, R.M. Eichler West, Journal of Computational Physics 222(2), 654 (2007)
  * [31] X. Li, H. Huang, P. Meakin, International Journal of Heat and Mass Transfer 53(13), 2908 (2010)
  * [32] F. Yang, G.A. Stack, V. Starchenko, Scientific reports 11(3495) (2021)
  * [33] C. Soulaine, H.A. Tchelepi, Transport In Porous Media 113, 431 (2016).  DOI 10.1007/s11242-016-0701-3
  * [34] C. Soulaine, P. Creux, H.A. Tchelepi, Transport in Porous Media (2019)
  * [35] Chatelin, Robin, Sanchez, David, Poncet, Philippe, ESAIM: M2AN 50(2), 565 (2016)
  * [36] C. Soulaine, S. Pavuluri, F. Claret, C. Tournassat, Environmental Modelling and Software 145(105199) (2021)
  * [37] H. Wen, L. Li, Geochimica et Cosmochimica Acta 210, 289 (2017)
  * [38] H. Luo, M. Quintard, G. Debenest, F. Laoufat, Comput. Geosci. 16, 913–932 (2012)
  * [39] J. Maes, C. Soulaine, Journal of Computational Physics 402, 109024 (2020)
  * [40] B. van Leer, Journal of Computational Physics 14(4), 361 (1974)
  * [41] D.L. Parkhurst, L. Wissmeier, Advances in Water Resources 83, 176 (2015)
  * [42] S.D.C. Walsh, N. Garapati, A.M.M. Leal, M.O. Saar, Geothermics 70, 146
  * [43] N. Seigneur, K.U. Mayer, C.I. Steefel, Reviews in Mineralogy & Geochemistry 85, 197 (2019)
  * [44] J.P. Pereira Nunes, M.J. Blunt, B. Bijeljic, Journal of Geophysical Research: Solid Earth 121, 558 (2016)
  * [45] OpenCFD, _OpenFOAM, the open source cfd toolbox, User Guide_.  OpenCFD Ltd (2016)
  * [46] A. Patsoukis-Dimou, H.P. Menke, J. Maes, Transport in Porous Media (2021).  DOI 10.1007/s11242-021-01718-8
  * [47] S. Pavuluri, J. Maes, F. Doster, microfluidics and nanofluidics 22(8) (2018).  DOI –10.1007/s10404-018-2106-9˝
  * [48] C. Peng, J.P. Crawshaw, G.C. Maitland, J.P.M. Trusler, Chemical geology 403, 74 (2015)
  * [49] H. Menke, J. Maes, S. Geiger, In preparation (2022)
  * [50] H.P. Menke, J. Maes, S. Geiger, Scientific Reports 11(2065) (2021).  DOI 10.1038/s41598-021-82029-2
  * [51] M. Liu, B. Kwon, P.K. Kang, (2022).  DOI 10.21203/rs.3.rs-1284059/v1
  * [52] X. Liu, A. Ormond, K. Bartko, P. Ortoleva, Journal of Petroleum Science and Engineering 17, 181 (1997)
  * [53] A. Ormond, P. Ortoleva, Journal of Geophysical Research 105, 16,737 (2000)
  * [54] F. Golfier, C. Zarcone, B. Bazin, R. Lenormand, D. Lasseux, M. Quintard, Journal of Fluid Mechanics 457, 213 (2002)
  * [55] A. Patsoukis-Dimou, A. Suzuki, H.P. Menke, S. Geiger, J. Maes, in _Proceedings of the 18th International Conference on Flow Dynamics (ICFD2021)_ (2020)
  * [56] C. Soulaine, S. Roman, A. Kovscek, H.A. Tchelepi, Journal of Fluid Mechanics 855, 616–645 (2018).  DOI 10.1017/jfm.2018.655
  * [57] P. Li, H. Deng, S. Molins, Frontiers in water (2022).  DOI 10.3389/frwa.2021.734518
  * [58] P. Szymczak, A.J.C. Ladd, J. Fluid Mech. 702, 239–264 (2012)
  * [59] A.N. Palmer, Geol. Soc. Am. Bull. 103, 1–21 (1991)
  * [60] J. Maes, H.P. Menke, in _Proceedings of the 17th International Conference on Flow Dynamics (ICFD2020)_ (2020)
  * [61] J. Maes, H.P. Menke, arXiv:2110.03311 (2021).  DOI arXiv:2110.03311
  * [62] J. Maes, H.P. Menke, arXiv:2105.10576 (2021).  DOI arXiv:2105.10576
  * [63] J. Maes, H.P. Menke, Transport in Porous Media 139, 271–299 (2021).  DOI 10.1007/s11242-021-01661-8
  * [64] S.V. Patankar, _Numerical Heat and Mass Transfer_ (Hemisphere Publ. Corp., Washington, 1980)
  * [65] J. Cooke, L. Armstrong, K. Luo, S. Gu, Computers Chemical Engineering 60, 297 (2014)