# porousMedia4Foam: Multi-scale open-source platform for hydro-geochemical simulations with OpenFOAM®

Cyprien Soulaine, Saideep Pavuluri, Francis Claret, Christophe Tournassat

► To cite this version:

Cyprien Soulaine, Saideep Pavuluri, Francis Claret, Christophe Tournassat. porousMedia4Foam: Multi-scale open-source platform for hydro-geochemical simulations with OpenFOAM®. Environmental Modelling and Software, 2021, 145, pp.105199. 10.1016/j.envsoft.2021.105199. insu-03346087v2

HAL Id: insu-03346087

https://insu.hal.science/insu-03346087v2

Submitted on 21 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# porousMedia4Foam: Multi-scale open-source platform for hydro-geochemical simulations with OpenFOAM®

Cyprien Soulaine<sup>a, *,</sup> Saideep Pavuluri<sup>a, b</sup>, Francis Claret<sup>b</sup>, Christophe Tournassat<sup>a, c</sup>

## ARTICLE INFO

### Keywords:

- Reactive transport modelling
- Multi-scale simulations
- Micro-continuum
- OpenFOAM
- PHREEQC

## ABSTRACT

porousMedia4Foam is a package for solving flow and transport in porous media using OpenFOAM® - a popular open-source numerical toolbox. We introduce and highlight the features of a new generation open-source hydro-geochemical module implemented within porousMedia4Foam, which relies on micro-continuum concept and which makes it possible to investigate hydro-geochemical processes occurring at multiple scales i.e. at the pore-scale, reservoir-scale and at the hybrid-scale. Geochemistry is handled by a third party package (e.g. PHREEQC) that is coupled to the flow and transport solver of OpenFOAM®. We conducted benchmarks across different scales to validate the accuracy of our simulator. We further looked at the evolution of mineral dissolution/precipitation in a fractured porous system. Application fields of this new package include the investigation of hydro-bio-geochemical processes in the critical zone, the modelling of contaminant transport in aquifers, as well as and the assessment of confinement performance for geological barriers.

## 1. Introduction

Over the last decades, reactive transport modelling (RTM) has become an essential tool for the study of subsurface processes involving flow, transport and geochemical reactions (Steefel et al., 2015a). This discipline is at the junction of two scientific communities, namely geochemistry and transport in porous media. RTM consists of computational models that describe the coupled physical, chemical, and biological processes that interact with each other over a broad range of spatial and temporal scales. RTM modelling tools enable the prediction of contaminant migration in polluted aquifers, and are used to design enhanced remediation techniques. Integration of physical and biogeochemical processes makes RTM also an ideal research instrument for elucidating the complex and non-linear interactions between roots, micro-organisms, water composition and minerals in the Critical Zone (Li et al., 2017). Other applications include the assessment of the long-term integrity of reservoirs for storing carbon dioxide, hydrogen or nuclear wastes in deep geological formations (DePaolo and Cole, 2013; Claret et al., 2018). In practice, three different kinds of models are used to describe reactive transport in porous media as illustrated in Fig. 1: (i) continuum models, (ii) pore-scale models, and (iii) hybrid models that combine both former approaches.

Continuum models (see Fig. 1c) are representative of the historical and standard approach for solving reactive transport in large-scale natural porous systems (Lichtner, 1985). Continuum-scale RTM codes include, among others, MIN3P (Mayer et al., 2002), CrunchFlow (Steefel et al., 2015a), TOUGHREACT (Xu et al., 2006), PFlotran (Lichtner et al., 2015) or HP1 (Jacques et al., 2008). In such codes, flow and transport equations are formulated in terms of volume-averaged equations with respect to a Representative Elementary Volume (REV) of the porous structure (Bear, 1972) and are coupled with geochemical reactions (see Steefel et al., 2015a for a comprehensive description of the coupling). The topology of the rock micro-structures is described using effective properties including porosity, tortuosity, permeability – or hydraulic conductivity – and specific surface area. Flow is usually modelled using Darcy’s law and multi-component species transport relies on a set of advection-dispersion-reaction equations. In addition to the classic challenges related to transport in porous media, including the medium heterogeneity awareness and the description of hydrodynamic dispersion, RTM has to consider the variation of rock properties in response to chemical reactions. Indeed, by enlarging or clogging pore throats and fractures, geochemical processes such as minerals dissolution and precipitation can alter the local flowfield and subsequently modify the rock properties, e.g permeability, tortuosity, accessible reactive surface area

<footnote><sup>a</sup> Institut des Sciences de la Terre d'Orléans, CNRS, Université d'Orléans, BRGM, Orléans, France</footnote>

<footnote><sup>b</sup> BRGM, French Geological Survey, Orléans, France</footnote>

<footnote><sup>c</sup> Energy Geosciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA, USA</footnote>

<footnote>* Corresponding author.
E-mail addresses: cyprien.soulaine@cnrs-orleans.fr, cyprien.soulaine@gmail.com (C. Soulaine).</footnote>
(Poonoosamy et al., 2020; Seigneur et al., 2019).

Changes in rock properties with chemical reaction is usually described in continuum-scale RTM as heuristic functions of the porosity. For example, in most of reactive transport codes, tortuosity is described using Archies's law, permeability change is modelled using Kozeny-Carman relationship, and mineral surface areas evolve as a two-third power law of the porosity (Xie et al., 2015). However, the complex interplay between reactions, advection, and diffusion can lead to highly nonlinear porosity feedback that is poorly captured using this kind of relationships that were not built on a strong theoretical background (Daccord and Lenormand, 1987; Garing et al., 2015). The limiting factor of the continuum-scale models is therefore the determination of empirical parameters and their evolution as a function of the progress of geochemical processes. To circumvent these challenges, recent efforts have focused on the numerical modelling of coupled hydro-geochemical processes at the pore-scale (Békri et al., 1997; Chen et al., 2013; Tar-takovsky et al., 2007; Molins et al., 2014; Molins et al., 2017).

In pore-scale models (see Fig. 1a), the pore network is fully resolved, i.e. each point of space is occupied by either a fluid or solid phase. As the exact knowledge of the phase distribution is known, continuum-scale concepts such as porosity, permeability, and reactive surface area do not apply at the pore-scale. They can be obtained, however, by aver- aging pore-scale simulation results if the computational domain is large enough to reach the size of a REV (Whitaker, 1999; Soulaine et al., 2013). The strategy that consists in simulating flow and transport in a three-dimensional image of a porous sample to characterize its continuum-scale properties has become an independent scientific discipline sometimes referred to as Digital Rock Physics (Blunt et al., 2013; Andrä et al., 2013a; Andrä et al., 2013b; Soulaine et al., 2021). Most of the efforts so far have been devoted to solve the Navier-Stokes equations under single (Spanne et al., 1994; Bijeljic et al., 2013; Gui- bert et al., 2015; Soulaine et al., 2016) and two-phase flow conditions (Horgue et al., 2013; Raeini et al., 2014; Graveleau et al., 2017; Maes and Soulaine, 2018; Pavuluri et al., 2020) to compute absolute and relative permeabilities. Despite the growing investment in the devel- opment of RTM at the pore-scale – pioneer simulators date back to the late 90s (Békri et al., 1997) – the pore-scale RTM field is still emerging. One of the main challenge of this approach consists in moving the flu- id/solid boundary with respect to chemical reactions at the mineral boundaries (Noiriel and Soulaine, 2021). A comprehensive review of the

different approaches used to solve this problem can be found in Molins et al. (Molins et al., 2020). It is only very recently that pore-scale sim- ulators have been proved mature for reproducing accurately and without any adjusting parameters the dissolution of a calcite crystal (Soulaine et al., 2017; Molins et al., 2020), or of a gypsum grain (Dutka et al., 2020). Actually, the benchmark presented in Molins et al. (Molins et al., 2020) is a first effort based on a relatively simple geochemical reaction (a single component that reacts with a single mineral using a first-order kinetics) to demonstrate the ability of current codes to accurately simulate mineral dissolution at the pore scale in a repro- ducible manner with several codes. Further developments and verifi- cations still need to be done for simulating multi-component aqueous solutions interacting with heterogeneous multi-mineral media using comprehensive reaction networks.

Naturally occurring porous media involve a wide range of spatial scales. For example, the important contrast in pore-size distributions in fractured porous rocks comes from much larger characteristic lengths for the fractures than for the surrounding porous matrix. Therefore, the domain size required to reach a REV limits the use of pure pore-scale modelling. Hybrid-scale models have been proposed to describe sys- tems that include multiple characteristic length-scales, for which some regions are described using pore-scale modelling while others are modelled with continuum approaches (see Fig. 1b) (Liu and Ortoleva, 1996; Liu et al., 1997). Two different approaches have been developed to solve hybrid-scale problems. On the one hand, the domain decom- position technique solves different physics on separate domains – one for Darcy flow, another for Stokes flow – linked together through appro- priate boundary conditions (Molins et al., 2019) including Beaver-Joseph and Ochoa-Tapia-Whitaker conditions (Beavers and Jo- seph, 1967; Ochoa-Tapia and Whitaker, 1995). On the other hand, micro-continuum models use a single set of partial differential equations throughout the computational domain regardless of the content of a grid block (Steefel et al., 2015b; Soulaine and Tchelepi, 2016). The latter approach is particularly well-suited to capture the dynamic displace- ment of the interface between the porous and solid-free regions without involving complex re-meshing strategies. For example, micro-continuum models have been used successfully to simulate the formation and growth of wormholes in acidic environments (Ormond and Ortoleva, 2000; Golfier et al., 2002; Soulaine and Tchelepi, 2016). Hybrid-scale modelling is also a powerful tool in image-based

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Porosity distribution considered in: a) a pure pore-scale approach (Navier-Stokes) for which the porosity is fully resolved, b) a micro-continuum approach (DBS) that handles region free of solid and porous region in the same framework, c) a pure continuum-scale approach (Darcy) for which all the control volumes contain an aggregate of fluid and solid.</figcaption></figure>
simulations to account for microscale features that are not visible in the images because they are smaller than the imaging instrument resolution (Arns et al., 2005; Apourvari and Arns, 2004; Scheibe et al., 2015; Soulaine et al., 2016; Soulaine et al., 2019).

In this study, we developed a comprehensive open-source simulator to model coupled hydro-geochemical processes at continuum-, pore- and hybrid-scales. This unique multi-scale framework relies on the micro-continuum model and its ability to tend asymptotically toward continuum-scale models if grid block contains solid content ($0 < \varphi < 1$) or towards pore-scale models otherwise ($\varphi = 1$) (Soulaine and Tchelepi, 2016). The resulting advanced RTM allows the treatment of complex reactions network as a function of flow conditions, water composition and minerals distribution within the rock including the complex porosity feedback between flow and chemistry. It is part of porousMedia4Foam, an open-source package developed by the authors to solve flow and transport in porous media within the popular simulation platform OpenFOAM®. Because of its versatility and advanced features such as three-dimensional unstructured grids, dynamic meshes, and high-performance computing, there is a growing interest in the community to develop mathematical models for solving flow and transport in porous media within OpenFOAM® (Horgue et al., 2015; Maes and Geiger, 2018; Orgogozo et al., 2014; Soulaine et al., 2017). We developed a generic interface to combine flow and transport models with existing geochemical packages. We illustrate the versatility of our coupling interface by combining flow models with geochemical models using PHREEQC (Parkhurst and Wissmeier, 2015), an open-source and popular geochemistry package that is used in many continuum-scale RTM (Rolle et al., 2018; Healy et al., 2018; Muniruzzaman and Rolle, 2019; Muniruzzaman et al., 2020; Moortgat et al., 2020).

In Section 2, we present the mathematical models implemented in porousMedia4Foam including a multi-scale and a continuum-scale flow solver, a wide range of porous properties models and the packages used to perform geochemical calculations. In Section 3, we verify the robustness of the coupled hydro-geochemical platform by simulating cases for which reference solutions exist both at the continuum-scale and at the pore-scale. Simulation results are compared with the results obtained with state-of-the-art reactive transport codes. Then in Section 4, we use the simulation framework to illustrate the potential of porousMedia4Foam to model hybrid-scale cases.

## 2. The porousMedia4Foam package

The multi-scale solver for simulating hydro-geochemical problems is part of porousMedia4Foam (https://github.com/csoulain/porousMedia4Foam), a generic platform for solving flow and transport in porous media at various scales of interest. porousMedia4Foam is an open-source platform developed by the authors using the C++ library OpenFOAM (http://www.openfoam.org). Hence, the package benefits from all the features of OpenFOAM including the solution of partial differential equations using the finite-volume method on three-dimensional unstructured grids as well as High Performance Computing. Although porousMedia4Foam has capabilities for solving two-phase flow (liquid-liquid and liquid-gas) in porous systems, the geochemistry coupling that is introduced in this paper only considers single-phase flow.

The code is organized in three interacting parts: a class that describes porous media properties (Section 2.4), the flow solvers (Section 2.2) and the geochemical packages (Section 2.3). As porousMedia4Foam intends to be a versatile platform, it is designed in a such way that other porous media models, other flow solvers or other geochemical packages can be easily implemented using the C++ code architecture. In this section, we introduce the models and their numerical implementation in the code. It is organized as a multiple-entry document in which each model is presented in a self-contained manner so that users can directly refer to the subsection associated with the said model.

### 2.1. Mineral distribution and porosity

A geological medium is made of an assembly of $N_s$ minerals whose porous properties are defined in Section 2.4. The distribution of each mineral $i$ on the computational grid is determined by the volume fraction,

$$Y_{s,i}(x, y, z, t) \text{ with } i \in [1, N_s], \quad (1)$$

in each grid block. The $N_s$ solid volume fraction fields are dimensionless. They can be initialized with uniform or distributed values. The evolution of $Y_{s,i}$ due to geochemical reactions is dictated by the geochemical packages that are described in Section 2.3. In some simulations, it is relevant to define an inert mineral, $Y_{s,inert}$, that is not part of the geochemical calculation.

The porosity field is computed by,

$$\varphi = 1 - \sum_{i}^{N_f-1} Y_{s,i} - Y_{s,inert}, \quad (2)$$

and can be updated at any moment following dissolution or precipitation processes. The porosity update at every time-step is optional.

### 2.2. Flow solvers

porousMedia4Foam for hydro-geochemical simulations includes three flow models: a multi-scale flow solver based on the micro-continuum approach, a continuum-scale Darcy solver and a constant velocity solver (see Table 1).

#### 2.2.1. dbsFoam: Multi-scale micro-continuum flow model

dbsFoam is a multi-scale flow solver based on the micro-continuum modelling approach developed in Soulaine and Tchelepi (Soulaine and Tchelepi, 2016). Micro-continuum approaches are intermediate between a pure Navier-Stokes description of the transport for which all the porosity is fully resolved (see Fig. 1a), and a pure continuum-scale modelling for which the flow is governed by Darcy's law (see Fig. 1c). This hybrid-scale approach relies on the Darcy-Brinkman-Stokes (DBS) equation (Brinkman, 1947) that allows for the modelling of flow and transport in regions free of solid and porous regions in a single framework (Neale and Nader, 1974; Soulaine and Tchelepi, 2016). DBS equation arises from the integration of Navier-Stokes equations over a control volume (Vafai and Tien, 1981; Hsu and Cheng, 1990; Bouquet-Melou et al., 2002; Goyeau et al., 2003). The momentum equation reads,

$$\frac{1}{\varphi} \left( \frac{\partial \rho_f}{\partial t} + \nabla \cdot \left( \frac{\rho_f}{\varphi} \mathbf{v}_f \mathbf{v}_f \right) \right) = -\nabla p_f + \rho_f \mathbf{g} + \nabla \cdot \left( \frac{\mu_f}{\varphi} (\nabla \mathbf{v}_f + i \nabla \mathbf{v}_f) \right) - \mu_f k^{-1} \mathbf{v}_f, \quad (3)$$

where $\varphi$ is the porosity, $\mathbf{v}_f$ is the seepage velocity, $p_f$ is the fluid pressure, $\mathbf{g}$ is the gravity, $\rho_f$ is the fluid density, $\mu_f$ is the fluid viscosity and $k$ is the cell permeability. The porous media properties including porosity and

<table><caption>**Table 1**

Summary of the flow solvers implemented in porousMedia4Foam for simulating hydro-geochemical processes.</caption><thead><tr><th>Name</th><th>Model</th><th>Setion</th><th>Comments</th></tr></thead><tbody><tr><td>dbsFoam</td><td>Micro-continuum (Darcy-Brinkman-Stokes)</td><td>Sec. 2.2.1</td><td>pore-scale, hybrid-scale, continuum-scale, (Soulaine and Tchelepi, 2016).</td></tr><tr><td>darcyFoam</td><td>Darcy's law</td><td>Sec. 2.2.1</td><td>continuum-scale only.</td></tr><tr><td>constantVelocityFoam</td><td>constant velocity profile</td><td>Sec. 2.2.3</td><td>uniform or non-uniform velocity profiles.</td></tr></tbody></table>
permeability change dynamically with geochemical processes and are updated at every time steps.

Eq. (3) is valid throughout the computational domain regardless the content of a cell. In regions that contain fluid only, $\varphi = 1$, and the drag force $\mu_k k^{-1}v_f$ vanishes so that the momentum equation tends to the Navier-Stokes equation. In regions that contain an aggregate of fluid and solid, $0 < \varphi < 1$, and the drag force is dominant over the inertial and viscous forces so that Eq. (3) tends asymptotically to Darcy's law.

The momentum equation, Eq. (3), can be used to model pore-scale flows. Indeed, if a solid region is approximated by a low-permeability low-porosity matrix, the velocity in this region goes to near zero which forces a no-slip boundary condition at the fluid/solid interface. This feature of the DBS equation is particularly relevant to solve Navier-Stokes problems using Cartesian grids only (also called penalized approach) (Angot et al., 1999; Soulaine and Tchelepi, 2016) and to move the fluid/solid interface with respect to geochemical processes such as precipitation/dissolution (Soulaine et al., 2017; Molins et al., 2020) or swelling by using the local porosity field, $\varphi$, as a phase indicator function (Carrillo and Bourg, 2019).

The pressure-velocity coupling is achieved by solving the momentum equation along with the micro-continuum continuity equation for multiple minerals. For an incompressible Newtonian aqueous fluid, the latter reads,

$$ \nabla \cdot v_f = \sum_{i=1}^{N_i} \dot{m}_{s,i} \left( \frac{1}{\rho_f} - \frac{1}{\rho_{s,i}} \right), \quad (4) $$

where $\rho_{s,i}$ is the solid density of mineral i and $\dot{m}_{s,i}$ represents the rate of phase change of solid into fluid, or of a fluid into solid. For example, it can represent the rate of solid minerals that is dissolved into aqueous solution. Inversely, it can describe an amount of fluid that is removed of a control volume because it has precipitated. The right-hand side of Eq. (4) is provided by the geochemistry calculation (Section 2.3). Although this term if often neglected in continuum-scale models, it ensures the mass balance at the fluid/solid interface in pore-scale simulations (Soulaine et al., 2017), as well as in continuum-scale simulations (Seigneur et al., 2018).

The flow model formed by Eqs (3) and (4) is discretized using the finite volume method and solved sequentially. The pressure-velocity coupling is handled by a predictor-corrector strategy based on the PIMPLE algorithm implemented in OpenFOAM. It consists in a combination of PISO (Pressure Implicit with Splitting of Operator, Issa, 1985) and SIMPLE (Semi-Implicit Method for Pressure Linked Equations, Patankar, 1980). PIMPLE algorithm allows both transient and steady-state simulations. Moreover, PIMPLE enables larger time steps than PISO. Further information regarding the numerics is found in Soulaine and Tchelepi (2016).

### 2.2.2. darcyFoam: Darcy's law

darcyFoam is a standard continuum-scale solver that is based on Darcy's law,

$$ v_f = -\frac{k}{\mu_f} (\nabla p_f - \rho_f g), \quad (5) $$

for describing flow in porous media. Numerically, Eq. (5) is combined along with Eq. (4) to form a Laplace equation solving implicitly for the pressure field, $p_f$. Then, the velocity field is calculated point-wise using Eq. (5) and $p_f$. If activatePorosityFeedback is switched on, Darcy's law is recalculated at every time steps to update the velocity and pressure fields according to the new permeability value.

Boundary conditions can be described by imposing fixed pressure or fixed velocity values on the domain edges. However, as darcyFoam solves implicitly for the pressure field, the boundary conditions on the velocity are transformed into pressure gradient conditions using Darcy's law:

$$ \mathbf{n} \cdot \nabla p_f = - \mathbf{n} \cdot (\mu_f k^{-1} \mathbf{v}_f - \rho_f \mathbf{g}), \quad (6) $$

where $\mathbf{n}$ is the vector normal to the domain boundary. In the code, Eq. (6) is achieved using the boundary condition darcyGradPressure (Horgue et al., 2015).

### 2.2.3. constantVelocityFoam: Constant flow rate

The flow solver constantVelocityFoam is used to model cases in which the chemical species are transported using a steady-state velocity field, $v_f$-uniform or non-uniform- provided as an input data that can come from a separate flow simulation. constantVelocityFoam is particularly useful if the feedback between geochemical reactions and the flow is negligible. Indeed, in such a case, the characteristic timescale of flow changes is much longer than the characteristic time of species transport and the calculation of the velocity profile can be decoupled from the species transport.

## 2.3. Geochemical packages

In porousMedia4Foam, complex reaction networks are handled by geochemical packages. The aqueous components are transported using the velocity profile, $v_f$, computed by the flow solver (see section 2.2) and surface reactions rely on the reactive surface area, $A_e$, calculated with the porous media models (see section 2.4.2). The code architecture of porousMedia4Foam is generic so that a wide variety of third-party geochemical packages can be coupled with our platform for solving hydro-geochemical processes at the pore-scale and at the continuum-scale. In this paper, we illustrate the potential of the coupled simulation framework using the popular geochemistry package PHREEQC (Parkhurst and Appelo, 2013; Parkhurst and Wissmeier, 2015). Models currently implemented in porousMedia4Foam to account for geochemistry are summarized in Table 2.

The geochemical packages update the water composition, $C_j$, and the distribution of the solid minerals, $Y_{s,i}$, and return the rate of solid changes,

$$ \dot{m}_{s,i} = - \frac{\partial \rho_{s,i} Y_{s,i}}{\partial t}, \quad (7) $$

where $\rho_{s,i}$ is the density of solid mineral i.

### 2.3.1. phreeqcRM

The phreeqcRM class calls the general-purpose geochemical reaction model PHREEQC through the PhreeqcRM module. It carries out the transport of the aqueous solution composition, $C_j$ (in mol/kgwater), along with equilibrium and kinetic reactions with the solid minerals described by $Y_{s,i}$. phreeqcRM is set with components, i.e. SOLUTION_MASTER-species total concentration.

The geochemistry setup is carried out using an input file that follows PHREEQC format. Hence, the aqueous composition is defined in the

<caption>Table 2</caption>

Summary of the geochemical packages implemented in porousMedia4Foam for simulating hydro-geochemical processes.

<table><thead><tr><th>Name</th><th>Model</th><th>Section</th><th>Comments</th></tr></thead><tbody><tr><td>phreeqcRM</td><td>PHREEQC</td><td>Sec. 2.3.1</td><td>Parkhurst and Wissmeier (2015)</td></tr><tr><td>simpleFirstOrderKineticMole</td><td>first order kinetic, C<sub>i</sub> in mol/m<sup>3</sup></td><td>Sec. 2.3.2</td><td>Molins et al. (2020); Soulaine et al. (2017)</td></tr><tr><td>transportOnly</td><td>no geochemistry</td><td>Sec. 2.3.3</td><td>–</td></tr><tr><td>flowOnly</td><td>no transport, no geochemistry</td><td>Sec. 2.3.4</td><td>–</td></tr></tbody></table>
block SOLUTION (0 for the composition of the injected fluid at the inlet boundary, 1 for the initial aqueous composition in the bulk). The EQUILIBRIUM smart PHASES and KINETICS blocks are generated automatically within the code and the user only has to assign before the calculation which mode of reactions is used for each mineral. Moreover, porousMedia4Foam can load any customized database using PHREEQC format.

The coupling between transport and reactions relies on an operator-splitting approach based on the Strang's algorithm (Strang, 1968). First, all species concentration fields, C<sub>j</sub>, are transported sequentially using the advection-dispersion equations,

$$ \frac{\partial \varphi C_j}{\partial t} + \nabla \cdot (\mathbf{v}_j C_j) - \nabla \cdot (\varphi D_j^* \nabla C_j) = 0, \quad (8) $$

where v<sub>j</sub> is the fluid velocity computed with the flow solver (see Section 2.2) and D<sub>j</sub><sup>*</sup> is an effective diffusion tensor that accounts for tortuosity and hydrodynamic dispersion effects (see Section 2.4.3). The transport equation is discretized on the computational domain using the finite-volume method and solved implicitly using OpenFOAM's engines.

Then, the volume fractions of solid minerals, Y<sub>s,i</sub>, are updated according to phase equilibrium and/or kinetic reaction calculations provided by PHREEQC. Reaction kinetics use the surface area computed at every time steps using the surface area models in section 2.4.2. It corresponds to the surface area per volume and its units are m<sup>2</sup>/m<sup>3</sup> (or m<sup>-1</sup>). Hence, the RATE block provided in PHREEQC database to describe reaction rates has to be defined accordingly.

### 2.3.2. simpleFirstOrderKineticMole

simpleFirstOrderKineticMole is a simple geochemical engine for solving the transport of a single species labelled A that reacts with solid minerals using first order kinetic reactions. It is an extension to multiple minerals of the model used in the benchmark presented in Molins et al. (2020) in which pore-scale simulators were used to model the dissolution of a calcite crystal by hydrochloric acid.

The chemical reaction reads,

$$ \text{Mineral}_{i(s)} + A = B_i. \quad (9) $$

The mass balance equation for species A reads,

$$ \frac{\partial \varphi C_A}{\partial t} + \nabla \cdot (\mathbf{v}_f C_A) - \nabla \cdot (\varphi D_j^* \nabla C_A) = - \left( \sum_{j=1}^{N_g} A_{s,j} (k_{j,A} \gamma_A) \right) C_A, \quad (10) $$

where v<sub>f</sub> is the fluid velocity, D<sub>j</sub><sup>*</sup> is an effective diffusion tensor, A<sub>s,j</sub> is the reactive surface area of mineral j (in m<sup>-1</sup>), and (k<sub>j,A</sub>γ<sub>A</sub>) is the constant of reaction of the species A with the mineral j (in m/s) following the notations adopted in Molins et al. (2020). In simpleFirstOrderKineticMole, the concentration field, C<sub>j</sub>, is defined in mol/m<sup>3</sup>. The equation is discretized on the computational grid using the finite-volume method and solved implicitly.

The distribution of solid minerals evolves according to,

$$ \frac{\partial Y_{s,i}}{\partial t} = -A_{s,A} (k_{i,A} \gamma_A) V_{m,s,i} C_A, \quad (11) $$

where V<sub>m,s,i</sub> (in m<sup>3</sup>/mol) is the molar volume of the reacting mineral.

### 2.3.3. transportOnly

transportOnly solves the advection-dispersion equation,

$$ \frac{\partial \varphi C_j}{\partial t} + \nabla \cdot (\mathbf{v}_f C_j) - \nabla \cdot (\varphi D_j^* \nabla C_j) = 0, \quad (12) $$

without considering geochemistry. It allows the transport of species using the dispersion models implemented in porousMedia4Foam (see Table 5).

### 2.3.4. flowOnly

flowOnly is an empty class for computing velocity profiles without species transport nor geochemistry. For example, Poonoosamy et al. (2020) used the multi-scale flow solver of porousMedia4Foam to compute the steady-state velocity profile in absence of geochemical reactions within a two-scale domain, i.e. a domain that contains both solid-free regions and porous regions (see Fig. 1b). This option is particularly interesting in cases for which geochemistry and flow can be treated independently from each other.

## 2.4. Porous media models

The flow solvers and geochemistry modules rely on porous media properties including absolute permeability, specific surface area and dispersion tensor. These properties describe pore-scale effects related to the micro-structure geometry of the porous medium. Hence, they may change if the micro-structure evolves with geochemical reactions.

### 2.4.1. Absolute permeability models

The absolute permeability describes the ability of a porous medium to conduct the flow. This property is intrinsic to the medium micro-structure and therefore evolves with geochemical processes including precipitation and dissolution. porousMedia4Foam includes several porosity-permeability relationships summarized in Table 3.

### 2.4.2. Surface area models

The estimation of the reactive surface area is crucial to model geochemical processes described by kinetic reactions. Actually, reactive surface area is a difficult quantity to assess as only a portion of the geometric surface area is accessible to the reactants. For example, in an advection-dominated transport, only the surfaces at the vicinity of the faster flowlines react (Soulaine et al., 2017) leading to a reactive surface area smaller than the geometric surface area. Moreover, the evolution of the specific surface area as the mineral volume fractions change due to dissolution or precipitation is not necessarily monotonic (Noiriel et al., 2009). Table 4 summarizes the models implemented in porousMedia4-Foam to describe the surface area as a function of the mineral volume fraction.

Unlike continuum-scale simulations, the surface area in pore-scale modelling is not an input parameter but is a direct output of the simulation. Indeed at the pore-scale, the micro-structure of the porous medium is fully resolved in the computational grid and there is a sharp interface between the fluid and the solid mineral. The Volume-of-Solid

<caption>Table 3</caption>

Summary of the permeability-porosity models implemented in porousMedia4-Foam. In the table, subscripts 0 refers to variable data at initial time. Optionally, φ<sub>0</sub> and k<sub>0</sub> are updated at every time-steps.

<table><thead><tr><th>Name</th><th>Expression</th><th>Comments</th></tr></thead><tbody><tr><td>none</td><td>k = 0</td><td>-</td></tr><tr><td>constant</td><td>k = k<sub>0</sub></td><td>k<sub>0</sub> is uniform or non-uniform.</td></tr><tr><td>Power-law</td><td>k = k<sub>0</sub> (φ/φ<sub>0</sub>)<sup>n</sup></td><td>n is a user defined variable.</td></tr><tr><td>Kozeny-Carman</td><td>k =<br>k<sub>0</sub> (φ/φ<sub>0</sub>)<sup>n</sup> (1 - φ<sub>0</sub>/φ)</td><td>by default, n = 2 and m = 3 (Kozeny, 1927; Carman, 1937).</td></tr><tr><td>Verma-Pruess</td><td>k = k<sub>0</sub> (φ/φ<sub>0</sub>)<sup>n</sup></td><td>n is a model parameter. φ<sub>c</sub> refers to the critical porositywhere permeability reduces to 0 (Verma and Pruess, 1988).</td></tr><tr><td>Hele-Shaw</td><td>k = h<sup>2</sup>/12</td><td>for simulating 2D depth-averaged flow in micromodels (e.g. Poonoosamy et al., 2020; Roman et al., 2016).</td></tr></tbody></table>
<table><caption>**Table 4**
Summary of the specific surface area models implemented in porousMedia4-Foam. Units of specific surface area are m⁻¹.</caption><thead><tr><th>Name</th><th>Expression</th><th>Comments</th></tr></thead><tbody><tr><td>None</td><td>A<sub>s</sub> = 0</td><td>If specific surface area is not necessary, e.g. for phase equilibrium calculation.</td></tr><tr><td>Constant</td><td>A<sub>s</sub> = A<sub>0</sub></td><td>–</td></tr><tr><td>Volume of solid</td><td>A<sub>s</sub> = |∇Y<sub>s</sub>|ψ</td><td>for pore-scale simulations only. Compute thelocal surface area based on the mineral mapping. ψ is a diffuse interface function (Soulaine et al., 2017).</td></tr><tr><td>Power-law</td><td>A<sub>s</sub> = A<sub>0</sub>(Y<sub>s</sub>)<sup>n</sup></td><td>n is a user defined variable</td></tr><tr><td>Sugar-lump</td><td>A<sub>s</sub> = (A<sub>0</sub> + A<sub>m</sub>(1 - (Y<sub>s</sub>/Y<sub>0</sub>)<sup>n<sub>1</sub></sup>)<sup>n<sub>2</sub></sup>)(Y<sub>s</sub>/Y<sub>0</sub>)<sup>n<sub>3</sub></sup></td><td>Evolution of the surface area of an aggregateduring dissolution (Noiriel et al., 2009). A<sub>m</sub> is the maximum surface area given by the sum of the surface areas of all individual particles, n<sub>1</sub>, n<sub>2</sub> and n<sub>3</sub> are user-defined parameters. including surface reduction due to hydro- geochemical coupling ( Soulaine et al., 2017). n, p, q are user defined parameters.</td></tr><tr><td>Hydro- geochemical coupling</td><td>A<sub>s</sub> = A<sub>0</sub>(Y<sub>s</sub>/Y<sub>0</sub>)<sup>n</sup>(1 - exp(Pe<sup>-p</sup>Da<sup>-q</sup>))</td><td></td></tr></tbody></table>

model computes the surface area of an explicit fluid/solid interface using the gradient of the volume fraction of mineral (see Soulaine et al., 2017 for additional details on the technique).

### 2.4.3. Dispersion models

In porous media, the spreading of a solute is not governed only by molecular diffusion ($D_i$) but also by the micro-structure and the local velocity field. On the one hand, the tortuosity of the porous structure tends to slow down the spreading. On the other hand, hydrodynamic dispersion stretches a solute band in the flow direction during its transport. In porousMedia4Foam, a single effective diffusion tensor, $D_i$, is used to represent both mechanisms. The models implemented in the code are summarized in Table 5.

<table><caption>**Table 5**
Summary of the dispersion models implemented in porousMedia4Foam. I is the unit tensor.</caption><thead><tr><th>Name</th><th>Expression</th><th>Comments</th></tr></thead><tbody><tr><td>none</td><td>D<sub>i</sub><sup>+</sup> = 0</td><td>for modelling transport by advection only.</td></tr><tr><td>diffusionOnly</td><td>D<sub>i</sub><sup>-</sup> = D<sub>i</sub>|l</td><td>no tortuosity effects, no hydrodynamic dispersion.</td></tr><tr><td>archiesLaw</td><td>D<sub>i</sub><sup>+</sup> = φ<sup>n</sup>D<sub>i</sub>|l</td><td>tortuosity is represented by φ<sup>n</sup>. By default, n = 0.</td></tr><tr><td>linearDispersion</td><td>D<sub>i</sub><sup>+</sup> = φ<sup>n</sup>(D<sub>i</sub> + α<sub>T</sub>|v|)l + (α<sub>L</sub> - α<sub>T</sub>)v<sub>ν</sub></td><td>tortuosity is represented by φ<sup>n</sup>, α<sub>L</sub> and α<sub>T</sub> are model parameters describing lateral and longitudinal dispersion respectively.</td></tr></tbody></table>

## 3. Verification of the hydro-geochemical simulation platform

In this section, the multi-scale hydro-geochemical simulation package porousMedia4Foam introduced in Section 2 is used along with PHREEQC to investigate several scenarios for which reference solutions exist. The verification of the results is achieved by comparison against benchmarks published in literature both at the continuum-scale and at the pore-scale. Essential files required to run all the test cases presented in this section are available as examples within the package. All simulations were run on Intel Xeon with 2.60 GHz.

### 3.1. Verification at the continuum-scale

We verify the ability of porousMedia4Foam to simulate coupled hydro-geochemical processes that include porosity feedback on the transport properties at the continuum-scale. We also verify that the multi-scale solver dbsFoam tends asymptotically towards Darcy's law in porous domains. The case is based on the Benchmark 1 described in Xie et al. (2015). It consists of a 2 meters long 1D column initially filled with 35% of inert mineral and 30% of calcite. An acid (pH = 3) is continuously injected at the inlet to initiate the dissolution of calcite according to,

$$ \begin{aligned} & \mathrm{CaCO}_3 + \mathrm{H}^+ \rightleftharpoons \mathrm{Ca}^{2+} + \mathrm{HCO}_3\text{-}, \\ & \mathrm{HCO}_3\text{-} \rightleftharpoons \mathrm{H}^+ + \mathrm{CO}_3^{2-}. \end{aligned} \qquad (13) $$

Table 6 provides the initial and boundary conditions data specific to the primary components. A difference of 0.007 m in hydraulic head is applied between the inlet and outlet (Xie et al., 2015) by fixing the pressure at 70 Pa and 0 Pa respectively at the inlet and outlet boundaries throughout the simulation.

The calcite dissolution with porosity feedback is simulated using both the multi-scale solver dbsFoam and the continuum-scale solver darcyFoam. The aqueous species are transported by advection only. The calcite dissolution is modelled using a kinetic rate of reaction with $k_{calcite} = 5 \times 10^{-5}$ mol/m²/s and the initial specific area $A_0 = 1$ m²/m³. As the calcite dissolves, the surface area decreases according to a power-law function with $n = 2/3$ (see Table 4). The porosity-permeability relationship is described by the Kozeny-Carman equation (Table 3). The initial permeability of the column is set to $k_0 = 1.186 \times 10^{-11}$ m². The column is spatially discretized with $\Delta x = 25$ mm (80 cells). 150 years are simulated with $\Delta t = 21600$ s.

There is a perfect match between dbsFoam and darcyFoam confirming that the multi-scale solver converges well to the continuum-scale solution predicted by Darcy's law (Fig. 2). The analysis of results includes the evolution of porosity (Fig. 2a), calcite volume fraction (Fig. 2b), hydraulic head (Fig. 2c) along the column length and the outflux over time (Fig. 2d). To be consistent with the benchmark of Xie et al. (2015), we consider the cross-sectional area of the column to be 1 m². As the dissolution of calcite occurs, the calcite volume fraction decreases and the porosity increases over time. In Fig. 2c, we notice different slopes of hydraulic head at different times. The slope is minimal where porosity is large and vice-versa. The velocity – and therefore the outflux – increases over time as the permeability (and porosity) of the system increases. The evolution of porosity, calcite volume fraction, hydraulic head and outflux predicted by porousMedia4Foam solvers are in close agreement with

<table><caption>**Table 6**
Initial and boundary conditions of the primary components for calcite dissolution under kinetic conditions considering porosity feedback.</caption><thead><tr><th>Primary components</th><th>Units</th><th>Initial conditions</th><th>Boundary condition</th></tr></thead><tbody><tr><td>pH</td><td>–</td><td>9.38</td><td>3</td></tr><tr><td>–</td><td>–</td><td>9.38</td><td>3</td></tr><tr><td>Ca</td><td>mol/kg<sub>water</sub></td><td>1.56 × 10<sup>-4</sup></td><td>9.97 × 10<sup>-5</sup></td></tr><tr><td>C(4)</td><td>mol/kg<sub>water</sub></td><td>2.56 × 10<sup>-4</sup></td><td>9.97 × 10<sup>-3</sup></td></tr><tr><td>S(6)</td><td>mol/kg<sub>water</sub></td><td>9.97 × 10<sup>-11</sup></td><td>6.44 × 10<sup>-4</sup></td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Fig. 2. Calcite dissolution under kinetic conditions considering feedback of porous media properties. Evolution of (a) porosity, (b) calcite volume fraction, (c) hydraulic head along the channel and (d) evolution of outflux. MIN3P data is from Xie et al. (2015) for comparison.</figcaption></figure>

those of MIN3P which demonstrates the ability of our platform to simulate hydro-geochemical processes with porosity feedback.

## 3.2. Verification at the pore-scale

In this section, we highlight the capabilities of our OpenFOAM package to model hydro-geochemical interactions occurring at the pore-scale using PHREEQC.

In porousMedia4Foam, pore-scale simulations are run using the micro-continuum approach through the flow solver dbsFoam. At the pore-scale, the reaction rates are directly applied at the fluid-mineral interface that is described explicitly in the computational grid using the mineral volume fraction. The micro-continuum approach has been used to simulate the dissolution of a calcite crystal at the pore-scale and compared successfully with microfluidic experiments (Soulaine et al., 2017). In Molins et al. (2020), the approach is compared with state-of-the-art RTM at the pore-scale using various numerical techniques including Chombo-Cruch with Level-Set (Molins et al., 2017), Lattice Boltzmann Method (Prasianakis et al., 2018), dissolFoam moving grids with conformal mapping (Starchenko et al., 2016), and Vortex

methods (Sanchez et al., 2019). The benchmark consists of a 0.2 mm diameter calcite crystal posted in a 1 mm long 0.5 width channel (see Fig. 3). An acidic solution of pH = 2 is continuously injected from the inlet at a rate of $U_{inj} = 1.2 \times 10^{-3}$ m/s. The calcite crystal dissolution is described considering a kinetic rate,

$$r = A_{calcite} (k_{calcite} \gamma) c_{H^+}, \quad (14)$$

where r is the reaction rate in mol/m³/s, $A_{calcite}$ is the reactive surface area in m²/m³ computed using the volume-of-solid approach (see Table 4), $(k_{calcite} \gamma) = 0.89 \times 10^{-3}$ m/s is the reaction rate constant of calcite and, $c_{H^+}$ is the concentration of $H^+$ in mol/m³. This reaction rate may not be fully representative of the underlying geochemical processes. It has been chosen in Molins et al. (2020) to demonstrate the ability of various approaches to move the fluid-mineral interface according to geochemical processes. All the numerical methods were able to capture accurately the shape evolution of the calcite crystal, giving confidence in pore-scale simulators for moving fluid-mineral interfaces along with geochemical processes.

Actually, in Soulaine et al. (2017) and Molins et al. (2020), the micro-continuum approach, dbsFoam, was combined with the

<figure><img src="image_5.png" /><figcaption>Fig. 3. The considered model set-up along with initial and boundary conditions to investigate calcite grain dissolution in a microchannel.</figcaption></figure>
geochemical package simpleFirstOrderKineticsMole (see Table 2) that solves Eq. (14) using OpenFOAM's internal engines. This limits drastically the applicability of the approach to comprehensive reaction networks. In this section, we reproduce the two-dimensional case presented in Molins et al. (2020) using dbsFoam and phreeqcRM to demonstrate the robustness of our coupling between OpenFOAM and PHREEQC at the pore-scale. The kinetic rate in PHREEQC input file has been modified to match Eq. (14). The system is spatially discretized using a Cartesian mesh of 128 × 64 cells. The simulation is run for 45 minutes using a time step size Δt = 5 ms.

The shape evolution of the calcite grain determined by the two approaches matches perfectly (Fig. 4) which verifies, therefore, that in our modelling plateform, PHREEQC can be used to model hydrogeochemical interactions occurring at the pore-scale. This gives us confidence for further investigations that rely on more complex reactive transport phenomena occurring at the pore-scale.

## 4. Hybrid-scale simulation in fractured porous media

In most subsurface environments, fractures intercept porous media domains. These fractures act as free-flow zones transporting substantial quantities of fluids alongside chemical species compared to the flow and transport that occur within the porous medium (Noiriel et al., 2007; Ajo-Franklin et al., 2018). The complex interplay between advection, diffusion, and reaction can lead to very different dissolution and precipitation patterns. For example, Poonoosamy et al. (2020) uses micro-Raman spectroscopy to visualise the replacement of celestite with barite in a fractured porous media flooded with a solution that contains barium ions. They observe that the mineral replacement occurs either uniformly or at the vicinity of the fracture-matrix interface. This difference in the mineral distribution was attributed to the injection flow rates leading to advection or diffusion-dominated transport.

We investigate such a multiscale system, where a fracture is sandwiched in between a reactive porous matrix made of 50% celestite (SrSO₄) having specific reactive surface area of A₀ = 20000 m²/m³ as shown in Fig. 5. The fracture has a length of ℓ = 0.03 m and height h = 0.002 m. Transport phenomena in the fracture is fully resolved, i.e. the flow is governed by Navier-Stokes equations whereas the flow in the matrix is described by Darcy's law. This hybrid-scale case is modelled using the dbsFoam solver. The initial porosity and permeability of the porous medium are φ₀ = 0.5 and k₀ = 10⁻¹² m², respectively. A solution

containing 300 mol/m³ of barium (Ba²⁺) is continuously injected through the inlet at a constant velocity U_inj for 200 hours. The dispersivity of species within the porous matrix are taken into account considering linear dispersion law (Table 5) with molecular diffusion set to Dⱼ = 10⁻⁹ m²/s, hydrodynamic dispersion coefficient set to αₘ = 10⁻⁵ m and tortuosity exponent set to n = 2. Once the barium ions reach the porous matrix, celestite dissociates into strontium (Sr²⁺) and sulphate (SO₄²⁻) ions. The barium ions react with sulphate ions resulting in the precipitation of barite (BaSO₄) according to the following reaction (Poonoosamy et al., 2018),

$$Ba^{2+} + SrSO_4 \rightarrow BaSO_4 + Sr^{2+}. \quad (15)$$

Celestite dissolution is taken into account considering kinetics with k_celestite = 10⁻⁵.⁶⁶ mol/m²/s whereas, barite precipitation is accounted considering phase equilibrium. Celestite reactive surface area evolves linearly with its volume fraction. The matrix permeability is updated according to Kozeny-Carman.

We investigate the ongoing hydrogeochemistry within this system considering two different injection velocities, U_inj = 10⁻² m/s and U_inj = 10⁻⁶ m/s. The Péclet number, Pe = U_injℓ/Dⱼ (where the reference length-scale is the fracture aperture), characterizes the importance of advection with respect to diffusion within the fracture. The highest velocity corresponds to advection-dominated transport (Pe ≈ 10⁴) while the lowest corresponds to diffusion-dominated regime (Pe ≈ 1). For both cases, the second Damkhöler number that determines the timescale of reaction with respect to the timescale of species diffusion at the mineral surface is Daᵢᵢ = k_celestiteℓ/(cₛBa²·Dⱼ) ≈ 3.6 × 10⁻⁴ (where the reference length-scale is the inverse of the specific surface area, ℓ = A₀⁻¹, according to Soulaine et al. (2017)). We notice differences in the pattern of celestite dissolution and barite precipitation whether the transport in the fracture is dominated by advection or by diffusion in agreement with Poonoosamy et al. (2020) observations (Fig. 6). For advection dominated regime (Pe > 1), there are two characteristic timescales for the solute transport: first, the barium ions flow through the fracture by advection, then they diffuse laterally into the matrix. Because of the timescale contrast between the two processes, the concentration profile of barium ions is uniform along the porous matrix which leads to a uniform pattern of celestite dissolution and barite precipitation as seen in Figs. 6 and 7. For diffusion-dominated regimes (Pe ≤ 1), the characteristic transport timescales both in the fracture and in the matrix are of the same order of magnitude. Therefore, the front of barium ions in the matrix follows the

<figure><img src="image_6.png" /><figcaption>Fig. 4. Evolution of the shape of the calcite grain as a function of time ((a) t = 0 min, (b) t = 15 min, (c) t = 30 min and (d) t = 45 min) predicted by dbsFoam + first order kinetics geochemical module (simpleFirstOrderKineticsMole, red line) and dbsFoam + PHREEQC (phreeqcRM geochemical module of porousMedia4Foam, black line). The red and black lines represent cells having calcite volume fraction of 0.5.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 5. Numerical set-up for the hybrid-scale case study. A fracture is sandwiched in between two layers of reactive porous medium. The reactive porous medium comprises of celestite. Barium is injected at the inlet. We investigate this scenario considering two different injection velocities, $U_{inj} = 10^{-2}$ m/s and $U_{inj} = 10^{-6}$ m/s.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Evolution of mineral volume fractions - celestite on the left and barite on the right - at different injection velocities. (a) Initial ($t = 0$ s) mineral volume fractions. (b) Mineral volume fractions at 100 h and 200 h for $U_{inj} = 1 \times 10^{-2}$ m/s, and (c) mineral volume fractions at 100 h and 200 h for $U_{inj} = 1 \times 10^{-6}$ m/s.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Plots comparing the mineral volume fractions of celestite and barite at three different time intervals (a) $t = 0$ s, (b) $t = 100$ h and (c) $t = 200$ h for different injection velocities. The mineral volume fraction data is collected along the length of the domain at a distance of 0.001 m from the top wall boundary as highlighted by dashed black line in Fig. 6 (see initial mineral volume fraction).</figcaption></figure>
diffusive front within the fracture. Subsequently, we observe mineral dissolution (celestite) and precipitation (barite) fronts within the porous matrix (Figs. 6 and 7).

This illustration emphasizes the capabilities of *porousMedia4Foam* to model dual-porosity systems in reactive environments using hybrid-scale approach. Our platform is therefore a powerful tool to complement and augment reactive transport experiments including high-resolution imaging of the evolution the fracture aperture including the effect of the weathered zone (Noiriel et al., 2007; Noiriel et al., 2009; Ajo-Franklin et al., 2018; Deng et al., 2020) and two-scale reactive microfluidic experiments (Poonoosamy et al., 2020; Nissan et al., 2021; Osselin et al., 2016).

In a logic of cascade of scales nested within each other, fractured porous media can be modelled by: i) pore-scale approaches, ii) discrete fracture networks, iii) dual porosity models. The hybrid-scale approach that we propose in this paper is intermediate between a full pore-scale description in which the fracture and the pores in the matrix are fully resolved and a discrete fracture network in which the matrix is modelled as a porous medium and the fractures as discrete elements that exchange matter with the matrix. The hybrid-scale approach is therefore crucial to characterize and improve the effective parameters (e.g. transfer functions between the porous matrix and the fracture) used in discrete fracture networks and larger-scale models.

## 5. Conclusion

We developed an integrated open-source simulator to model hydro-geochemical processes at various scales of interest including pore-scale and reservoir-scale. The simulation platform is part of *porousMedia4-Foam*, a package that solves flow and transport in porous media using the open-source library OpenFOAM. The modelling framework handles complex reactions network as a function of flow conditions, water composition and minerals distribution within the rock including the complex porosity feedback between flow and chemistry. Moreover, *porousMedia4Foam* benefits from all features of OpenFOAM libraries. Hence, the code is fully parallel and handles structured as well as un-structured grids in one, two and three dimensions. The interface between the flow simulator and the geochemistry is generic and can be used to couple a large variety of geochemical packages. In this paper, we illustrated the hydro-geochemical capabilities of the coupled solver using PHREEQC.

Unlike other reactive transport simulators, *porousMedia4Foam* is multi-scale, i.e. a unique flow solver describes transport processes both at the continuum-scale and the pore-scale. Importantly, the two scales can be solved simultaneously in geological formations that feature large contrast of permeability and porosity. For example, in fractured rocks, *porousMedia4Foam* solves Stokes flow in the fracture network and Darcy's law in the porous matrix. This multi-scale model is achieved using the micro-continuum approach, hybrid-scale approach based on the Darcy-Brinkman-Stokes equation. Indeed, this approach is intermediate between a pure Navier-Stokes description of the transport for which all the porosity is fully resolved and pure continuum-scale modeling based on Darcy's law. Besides this hybrid-scale approach, *porousMedia4Foam* also includes a standard Darcy solver for continuum-scale simulations. Therefore, the same simulator can be used to simulate flow, transport, and geochemical reactions in an reservoir and in 3D micro-tomography images.

The coupled hydro-geochemical simulator was verified by running cases for which reference solutions exist. These solutions are well-established and used in the reactive transport community to benchmark state-of-the-art codes available both at the continuum-scale (Xie et al., 2015) and at the pore-scale (Molins et al., 2020). Finally, we demonstrated the ability of our advanced modelling framework to simulate dissolution and precipitation processes in fractured porous media at the pore-scale using the hybrid-scale approach. Here, the reactive medium consisted of celestite grains that reacted with a barium

chloride solution injected into the system, leading to the dissolution of celestite and the growth of barite. We observed differences in mineral precipitation - dissolution patterns by varying the injection rates.

Because *porousMedia4Foam* has already capabilities for modelling two-phase flow in porous media both at the pore- and Darcy's scales using two-phase micro-continuum technique (Soulaine et al., 2019; Soulaine et al., 2018; Carrillo et al., 2020), true multiscale and multi-phase RTM is envisioned to be implemented in *porousMedia4Foam* framework.

## Software and data availability

*porousMedia4Foam*, the software introduced in this paper is built using open-source libraries including OpenFOAM and PHREEQC. The source code and the cases presented in the paper are available on GitHub (https://github.com/csoulain/porousMedia4Foam).

## Contributions of the authors

CS is the *porousMedia4Foam* architect. CS and SP implemented new models in *porousMedia4Foam* and corrected bugs. SP and CT designed and setup the benchmark problems. SP run the simulations. CS, SP, CT and FC discussed, interpreted the results and wrote the paper. CS, CT and FC applied for funding.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

The research leading to these results has received funding from the French Agency for Research (Agence Nationale de la Recherche, ANR) through the labex Voltaire ANR-10-LABX-100-01, the grant CATCH ANR-18-CE05-0035, and through the FraMatI project under contract ANR-19-CE05-0002. It has also received financial support from the CNRS through the MITI interdisciplinary programs. This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreements No 847593 (WP DONUT) and No 850626 (REFLECT Project). SP postdoctoral fellowship was funded by BRGM through the PORE-REACTIF project from the Alliance Nationale de Coordination de la Recherche pour l'Energie (ANCRE). The authors benefitted from the use of the cluster at the Centre de Calcul Scientifique en région Centre-Val de Loire.

## References

* Ajo-Franklin, J., Voltolini, M., Molins, S., Yang, L., 2018. Coupled processes in a fractured reactive system: a dolomite dissolution study with relevance to GCS caprock integrity. *Caprock Integrity in Geological Storage: Hydrogeochemical and Hydrogeomechanical Processes and Their Impact on Storage Security*. Wiley Publishing, New York.

* Andrä, H., Combaret, N., Dvorkin, J., Glatt, E., Han, J., Kabel, M., Keehm, Y., Krzikalla, F., Lee, M., Madonna, C., Marsh, M., Mukerji, T., Saenger, E.H., Sain, R., Saxena, N., Ricker, S., Wiegmann, A., Zhan, X., 2013a. Digital rock physics benchmarks Part I: imaging and segmentation. Comput. Geosci. 50, 25–32. ISSN 0098-3004. https://doi.org/10.1016/j.cageo.2012.09.005. http://www.sciencedirect.com/science/article/pii/S0098300412003147, (benchmark problems, datasets and methodologies for the computational geosciences).

* Andrä, H., Combaret, N., Dvorkin, J., Glatt, E., Han, J., Kabel, M., Keehm, Y., Krzikalla, F., Lee, M., Madonna, C., Marsh, M., Mukerji, T., Saenger, E.H., Sain, R., Saxena, N., Ricker, S., Wiegmann, A., Zhan, X., 2013b. Digital rock physics benchmarks Part II: computing effective properties. Comput. Geosci. 50, 33–43. ISSN 0098-3004. https://doi.org/10.1016/j.cageo.2012.09.008. http://www.sciencedirect.com/science/article/pii/S0098300412003172, (benchmark problems, datasets and methodologies for the computational geosciences).

* Angot, P., Bruneau, C.-H., Fabrie, P., 1999. A penalization method to take into account obstacles in incompressible viscous flows. Numer. Math. 81 (4), 497–520.
Arns, C., Bauget, F., Limaye, A., Sakellariou, A., Senden, T., Sheppard, A., Sok, R., Pinczewski, W., Bakke, S., Berge, L., Oeren, P.-E., Knackstedt, M., 2005. Pore-scale characterization of carbonates using X-ray microtomography. *SPE J.* 10 (4), 475–484.

Bear, J., 1972. Dynamics of Fluids in Porous Media. Elsevier, New York. URL
*Dynamicsoffluidsinporousmedia*.

Beavers, G.S., Joseph, D.D., 1967. Boundary conditions at a naturally permeable wall.
*J. Fluid Mech.* 30, 197–207. https://doi.org/10.1017/S0022112067001375.

Békri, S., Thovert, J.-F., Adler, P., 1997. Dissolution and deposition in fractures. *Eng. Geol.* 48 (3-4), 283–308. https://doi.org/10.1016/S0013-7952(97)00044-6.

Bijeljic, B., Raeini, A., Mostaghimi, P., Blunt, M.J., 2013. Predictions of non-Fickian
solute transport in different classes of porous media using direct simulation on pore-
scale images. *Phys. Rev. E* 87, 013011. https://doi.org/10.1103/
*PhysRevE.87.013011*. URL.

Blunt, M.J., Bijeljic, B., Dong, H., Gharbi, O., Iglauer, S., Mostaghimi, P., Paluszny, A., Pentland, C., 2013. Pore-scale imaging and modelling. *Adv. Water Resour.* 51, 197–216.

Bousquet-Melou, P., Goyeau, B., Quintard, M., Fichot, F., Gobin, D., 2002. Average momentum equation for interdendritic flow in a solidifying columnar mushy zone. ISSN 0017-9310 Int. J. Heat Mass Tran. 45 (17), 3651–3665. https://doi.org/10.1016/S0017-9310(02)00077-7. URL. http://www.sciencedirect.com/science/article/pii/S0017931002000777.

Brinkman, H.C., 1947. A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles. *Appl. Sci. Res.* A1, 27–34.

Carman, P.C., 1937. Fluid flow through granular beds. Trans. Inst. Chem. Eng. 15, 150–166.

Carrillo, F.J., Bourg, I.C., 2019. A Darcy-Brinkman-Biot approach to modeling the
hydrology and mechanics of porous media containing macropores and deformable
microporous regions. *Water Resour. Res.* 55 (10), 8096–8121. https://doi.org/
10.1029/2019WR024712. URL.

Carrillo, F.J., Bourg, I.C., Soulaine, C., 2020. Multiphase flow modeling in multiscale porous media: an open-source micro-continuum approach. J. Comput. Phys. X 8, 100073. https://doi.org/10.1016/j.jcxp.2020.100073. ISSN 2590-5552.

Chen, L., Kang, Q., Robinson, B.A., He, Y.-L., Tao, W.-Q., 2013. Pore-scale modeling of multiphase reactive transport with phase transitions and dissolution-precipitation processes in closed systems. *Phys. Rev. E* 87, 043306. https://doi.org/10.1103/
*PhysRevE.87.043306*.

Claret, F., Marty, N., Tournassat, C., 2018. Modeling the Long-Term Stability of Multi-
Barrier Systems for Nuclear Waste Disposal in Geological Clay Formations. chap. 8.
John Wiley & Sons, Ltd. ISBN 9781119060031, 395–451. https://doi.org/10.1002/9781119060031.ch8.

Daccord, G., Lenormand, R., 1987. Fractal patterns from chemical dissolution. *Nature* 325 (6099), 41–43.

Deng, H., Fitts, J.P., Tappero, R.V., Kim, J.J., Peters, C.A., 2020. Acid erosion of
carbonate fractures and accessibility of arsenic-bearing minerals: in operando
synchrotron-based microfluidic experiment. *Environ. Sci. Technol.* 54 (19),
12502–12510. https://doi.org/10.1021/acs.est.0c03736 pMID: 32845141.

DePaolo, D.J., Cole, D.R., 2013. Geochemistry of geologic carbon sequestration: an
overview. *Rev. Mineral. Geochem.* 77 (1), 1–14. https://doi.org/10.2138/
rmg.2013.77.1.

Dutka, F., Starchenko, V., Osselin, F., Magni, S., Szymczak, P., Ladd, A.J., 2020. Time-
dependent shapes of a dissolving mineral grain: comparisons of simulations with
microfluidic experiments. *Chem. Geol.* 540, 119459.

Garing, C., Gouze, P., Kassab, M., Riva, M., Guadagnini, A., 2015. Anti-correlated
porosity–permeability changes during the dissolution of carbonate rocks:
experimental evidences and modeling. *Transport Porous Media* 107 (2), 595–621.

Golfier, F., Zarcone, C., Bazin, B., Lenormand, R., Lasseux, D., Quintard, M., 2002. On the
ability of a Darcy-scale model to capture wormhole formation during the dissolution
of a porous medium. J. Fluid Mech. 457, 213–254.

Goyeau, B., Lhuillier, D., Gobin, D., Velarde, M., 2003. Momentum transport at a
fluid-porous interface. Int. J. Heat Mass Tran. 46 (21), 4071–4081.

Graveleau, M., Soulaine, C., Tchelepi, H.A., 2017. Pore-scale simulation of interphase
multicomponent mass transfer for subsurface flow. *Transport Porous Media* 120 (2),
287–308.

Guibert, R., Nazarova, M., Horgue, P., Hamon, G., Creux, P., Debenest, G., 2015.
Computational permeability determination from pore-scale imaging: sample size,
mesh and method sensitivities. *Transport Porous Media* 107 (3), 641–656. https://
doi.org/10.1007/s11242-015-0458-0. ISSN 0169-3913.

Healy, R.W., Haile, S.S., Parkhurst, D.L., Charlton, S.R., 2018. VS2DRTI: simulating heat and reactive sulfate transport in variably saturated porous media. *Ground Water* 56 (5), 810–815. https://doi.org/10.1111/gwat.12640.

Horgue, P., Augier, F., Duru, P., Prat, M., Quintard, M., 2013. Experimental and
numerical study of two-phase flows in arrays of cylinders. *Chem. Eng. Sci.* 102,
335–345. https://doi.org/10.1016/j.ces.2013.08.031. ISSN 0009-2509. http://
www.sciencedirect.com/science/article/pii/S0009250913005745.

Horgue, P., Soulaine, C., Franc, J., Guibert, R., Debenest, G., 2015. An open-source toolbox for multiphase flow in porous media. ISSN 0010-4655 Comput. Phys. Commun. 187, 217–226. https://doi.org/10.1016/j.cpc.2014.10.005. URL. http://www.sciencedirect.com/science/article/pii/S0010465514003400.

Hsu, C., Cheng, P., 1990. Thermal dispersion in a porous medium. Int. J. Heat Mass Tran. 33 (8), 1587–1597.

Issa, R.I., 1985. Solution of the implicitly discretised fluid flow equations by operator-splitting. J. Comput. Phys. 62, 40–65.

Jacques, D., Simunek, J., Mallants, D., Van Genuchten, M.T., 2008. Modeling coupled hydrologic and chemical processes: long-term uranium transport following phosphorus fertilization. *Vadose Zone J.* 7 (2), 698–711.

Kozeny, J., 1927. Über kapillare leitung der wasser in boden, der Vienna Academy of Science, Vienna, Proc. Class. Irel. 136, 271–306.

Li, L., Maher, K., Navarri-Sitchler, A., Druhan, J., Meile, C., Lawrence, C., Moore, J.,
Perdrial, J., Sullivan, P., Thompson, A., et al., 2017. Expanding the role of reactive
transport models in critical zone processes. *Earth Sci. Rev.* 165, 280–301.

Lichtner, P.C., 1985. Continuum model for simultaneous chemical reactions and mass transport in hydrothermal systems. *Geochem. Cosmochim.* Acta 49 (3), 779–800.

Lichtner, P.C., Hammond, G.E., Lu, C., Karra, S., Bisht, G., Andre, B., Mills, R., Kumar, J., 2015. PFLOTRAN user manual: a massively parallel reactive flow and transport model for describing surface and subsurface processes. Tech. Rep., United States. URL. https://www.osti.gov/servlets/purl/1168703.

Liu, X., Orteleva, P., 1996. A general-purpose, geochemical reservoir simulator. *SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers.*

Liu, X., Ormond, A., Bartko, K., Ying, L., Orteleva, P., 1997. A geochemical reaction-transport simulator for matrix acidizing analysis and design. J. Petrol. Sci. Eng. 17 (1), 181–196.

Maes, J., Geiger, S., 2018. Direct pore-scale reactive transport modelling of dynamic wettability changes induced by surface complexation. *Adv. Water Resour.* 111, 6–19.

Maes, J., Soulaine, C., 2018. A new compressive scheme to simulate species transfer across fluid interfaces using the Volume-Of-Fluid method. *Chem. Eng. Sci.* 190, 405–418.

Mayer, K.U., Frind, E.O., Blowes, D.W., 2002. Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions. *Water Resour. Res.* 38 (9), 13–1.

Molins, S., Trebotich, D., Yang, L., Ajo-Franklin, J.B., Ligocki, T.J., Shen, C., Steefel, C.I., 2014. Pore-scale controls on calcite dissolution rates from flow-through laboratory and numerical experiments. *Environ. Sci. Technol.* 48 (13), 7453–7460.

Molins, S., Trebotich, D., Miller, G.H., Steefel, C.I., 2017. Mineralogical and transport controls on the evolution of porous media texture using direct numerical simulation. *Water Resour. Res.* 53 (5), 3645–3661.

Molins, S., Trebotich, D., Arora, B., Steefel, C.I., Deng, H., 2019. Multi-scale model of reactive transport in fractured media: diffusion limitations on rates. *Transport Porous Media* 128 (2), 701–721.

Molins, S., Soulaine, C., Prasianakis, N., Abbasi, A., Poncet, P., Ladd, A., Starchenko, V., Roman, S., Trebotich, D., Tchelepi, H., Steefel, C., 2020. Simulation of mineral dissolution at the pore scale with evolving fluid-solid interfaces: review of approaches and benchmark problem set. *Comput. Geosci.* 1–34ISSN 1573–1499.

Moortgat, J., Li, M., Amooie, M., Zhu, D., 2020. A higher-order finite element reactive transport model for unstructured and fractured grids. *Sci. Rep.* 10, 15572. https://doi.org/10.1038/s41598-020-72354-3.

Muniruzzaman, M., Rolle, M., 2019. Multicomponent ionic transport modeling in physically and electrostatically heterogeneous porous media with PhreeqCRM coupling for geochemical reactions. *Water Resour. Res.* 55 (12), 11121–11143. https://doi.org/10.1029/2019WR026373. URL.

Muniruzzaman, M., Karlsson, T., Ahnadi, N., Rolle, M., 2020. Multiphase and multicomponent simulation of acid mine drainage in unsaturated mine waste: modeling approach, benchmarks and application examples. *ISSN 0883-2927 Appl. Geochem.* 120, 104677. https://doi.org/10.1016/j.apgeochem.2020.104677. URL. https://www.sciencedirect.com/science/article/pii/S0883292720301682.

Neale, G., Nader, W., 1974. Practical significance of Brinkman's extension of Darcy's law: coupled parallel flows within a channel and a bounding porous medium. *Can. J. Eng.* 52 (4), 475–478.

Apourvari, S.N., Arns, C.H., 2004. An Assessment of the Influence of Micro-porosity for Effective Permeability Using Local Flux Analysis on Tomographic Images. In: *International Petroleum Technology Conference, 19-22 January, Doha, Qatar.*

Noiriel, C., Madé, B., Gouze, P., 2007. Impact of coating development on the hydraulic and transport properties in argillaceous limestone fracture. *Water Resour. Res.* 43 (9), 1–16.

Nissan, A., Alcolombri, U., de Schaetzen, F., Berkowitz, B., Jimenez-Martinez, J., 2021. Reactive transport with fluid-solid interactions in dual-porosity media. *ACS EST Water* https://doi.org/10.1021/acsestwater.0cc0043. ISSN 2690-0637.

Noiriel, C., Luquot, L., Madé, B., Raimbault, L., Gouze, P., van der Lee, J., 2009. Changes in reactive surface area during limestone dissolution: an experimental and modelling study. *Chem. Geol.* 265 (1–2), 160–170. ISSN 0009-2541. https://doi.org/10.1016/j.chemgeo.2009.01.032. http://www.sciencedirect.com/science/article/pii/S000925410900643. c02 geological storage: Integrating geochemical, hydrodynamical, mechanical and biological processes from the pore to the reservoir scale.

Noiriel, C., Soulaine, C., 2021. Pore-scale imaging and modelling of reactive flow in evolving porous media: tracking the dynamics of the fluid-rock interface. *Transport Porous Media* https://doi.org/10.1007/s11242-021-01613-2.

Ochoa-Tapia, J.A., Whitaker, S., 1995. Momentum transfer at the boundary between a porous medium and a homogeneous fluid: I. Theoretical development. *Int. J. Heat Mass Tran. 38* (14), 2635–2646.

Orgogozo, L., Renon, N., Soulaine, C., Hénon, F., Tomer, S., Labat, D., Pokrovsky, O., Sekhar, M., Ababou, R., Quintard, M., 2014. An open source massively parallel solver for Richards equation: mechanistic modelling of water fluxes at the watershed scale. *ISSN 0010-4655 Comput. Phys. Commun.* 185 (12), 3358–3371. https://doi.org/10.1016/j.cpc.2014.08.004. URL. http://www.sciencedirect.com/science/article/artic/le/pii/S0010465514002719.

Ormond, A., Orteleva, P., 2000. Numerical modeling of reaction-induced cavities in a porous rock. *J. Geophys. Res.*: Solid Earth 105 (B7), 16737–16747.

Osselin, F., Kondratiuk, P., Budek, A., Cybulski, O., Garstecki, P., Szymczak, P., 2016. Microfluidic observation of the onset of reactive-infiltration instability in an analog fracture. *Geophys. Res. Lett.* 43, 6907–6915. https://doi.org/10.1002/2016gl069261. ISSN 0094-8276.
Parkhurst, D., Appelo, C., 2013. Description of input and examples for PHREEQC Version 3—a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations, 6-A43. U.S. Department of the Interior, U.S. Geological Survey Techniques and Methods.

Parkhurst, D.L., Wissmeier, L., 2015. PhreeqCRM: a reaction module for transport simulators based on the geochemical model PHREEQC. Adv. Water Resour. 83, 176–189.

Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Taylor & Francis, Washington, DC.

Pavuluri, S., Maes, J., Yang, J., Regaieg, M., Moncorgé, A., Doster, F., 2020. Towards pore network modelling of spontaneous imbibition: contact angle dependent invasion patterns and the occurrence of dynamic capillary barriers. Comput. Geosci. 24, 951–969. https://doi.org/10.1007/s10596-019-09842-7.

Poonoosamy, J., Wanner, C., Alt Epping, P., Aguila, J.F., Samper, J., Montenegro, L., Xie, M., Su, D., Mayer, K.U., Mäder, U., Van Loon, L.R., Kosakowski, G., 2018. Benchmarking of reactive transport codes for 2D simulations with mineral dissolution-precipitation reactions and feedback on transport parameters. Comput. Geosci. 1–22. https://doi.org/10.1007/s10596-018-9793-x.

Poonoosamy, J., Soulaine, C., Burmeister, A., Deissmann, G., Bosbach, D., Roman, S., 2020. Microfluidic flow-through reactor and 3D Raman imaging for in situ assessment of mineral reactivity in porous and fractured porous media. Lab-on-a-Chip 20 (14), 2562–2571. https://doi.org/10.1039/d0lc00360c.

Prasianakis, N.I., Gatschet, M., Abbasi, A., Churakov, S.V., 2018. Upscaling strategies of porosity-permeability correlations in reacting environments from pore-scale simulations, 2018 Geofluids 1–8. https://doi.org/10.1155/2018/9260603. ISSN 1468-8115.

Raeini, A.Q., Blunt, M.J., Bijeljic, B., 2014. Direct simulations of two-phase flow on micro-CT images of porous media and upscaling of pore-scale forces. Adv. Water Resour. 74, 116–126. ISSN 0309-1708, URL http://www.sciencedirect.com/science/article/pii/S0309170814001730.

Rolle, M., Sprocati, R., Masi, M., Jin, B., Muniruzzaman, M., 2018. Nernst-Planck-based description of transport, coulombic interactions, and geochemical reactions in porous media: modeling approach and benchmark experiments. Water Resour. Res. 54 (4), 3176–3195. https://doi.org/10.1002/2017WR022344. URL.

Roman, S., Soulaine, C., AlSaud, M.A., Kovscek, A., Tchelepi, H., 2016. Particle velocimetry analysis of immiscible two-phase flow in micromodels. ISSN 0309-1708 Adv. Water Resour. 95, 199–211. https://doi.org/10.1016/j.adwatres.2015.08.015. URL http://www.sciencedirect.com/science/article/pii/S0309170815002018.

Sanchez, D., Hume, L., Chatelin, R., Poncet, P., 2019. Analysis of the 3D non-linear Stokes problem coupled to transport-diffusion for shear-thinning heterogeneous microscale flows, applications to digital rock physics and mucociliary clearance. Math. Model.Numer. Anal. 53 (1124), 1083. https://doi.org/10.1051/m2an/2019013. ISSN 0764-583X.

Scheibe, T.D., Perkins, W.A., Richmond, M.C., McKinley, M.I., Romero-Gomez, P.D.J., Oostrom, M., Wietsma, T.W., Serkowski, J.A., Zachara, J.M., 2015. Pore-scale and multiscale numerical simulation of flow and transport in a laboratory-scale column. ISSN 1944-7973, doi:10.1002/2014WR015959 Water Resour. Res. 51 (2), 1023–1035. https://doi.org/10.1002/2014WR015959. URL.

Seigneur, N., Lagneau, V., Corvisier, J., Dauzères, A., 2018. Recoupling flow and chemistry in variably saturated reactive transport modelling - an algorithm to accurately couple the feedback of chemistry on water consumption, variable porosity and flow. ISSN 0309-1708 Adv. Water Resour. 122, 355–366. https://doi.org/10.1016/j.adwatres.2018.10.025. URL https://www.sciencedirect.com/science/article/pii/S0309170818306304.

Seigneur, N., Mayer, K.U., Steefel, C.I., 2019. Reactive transport in evolving porous media. Rev. Mineral. Geochem. 85 (1), 197–238. https://doi.org/10.2138/rmg.2019.85.7.

Soulaine, C., Tchelepi, H.A., 2016. Micro-continuum approach for pore-scale simulation of subsurface processes. Transport Porous Media 113, 431–456. https://doi.org/10.1007/s11242-016-0701-3.

Soulaine, C., Davit, Y., Quintard, M., 2013. A two-pressure model for slightly compressible single phase flow in bi-structured porous media. ISSN 0009-2509 Chem. Eng. Sci. 96, 55–70. https://doi.org/10.1016/j.ces.2013.03.060. URL http://www.sciencedirect.com/science/article/pii/S0009250913002492.

Soulaine, C., Gjetvaj, F., Garing, C., Roman, S., Russian, A., Gouze, P., Tchelepi, H., 2016. The impact of sub-resolution porosity of X-ray microtomography images on the permeability. Transport Porous Media 113 (1), 227–243. https://doi.org/10.1007/s11242-016-0690-2.

Soulaine, C., Roman, S., Kovscek, A., Tchelepi, H.A., 2017. Mineral dissolution and wormholing from a pore-scale perspective. J. Fluid Mech. 827, 457–483. https://doi.org/10.1017/jfm.2017.499.

Soulaine, C., Roman, S., Kovscek, A., Tchelepi, H.A., 2018. Pore-scale modelling of multiphase reactive flow. Application to mineral dissolution with production of CO2. J. Fluid Mech. 855, 616–645. https://doi.org/10.1017/jfm.2018.655.

Soulaine, C., Creux, P., Tchelepi, H.A., 2019. Micro-continuum framework for pore-scale multiphase fluid transport in shale formations. Transport Porous Media 127, 85–112.

Soulaine, C., Maes, J., Roman, S., 2021. Computational microfluidics for geosciences. Frontiers in Water 3, 1–11. https://doi.org/10.3389/frwa.2021.643714.

Spanne, P., Thovert, J., Jacquin, C., Lindquist, W., Jones, K., Adler, P., 1994. Synchrotron computed microtomography of porous media: topology and transports. Phys. Rev. Lett. 73 (14), 2001.

Starchenko, V., Marra, C.J., Ladd, A.J., 2016. Three-dimensional simulations of fracture dissolution. J. Geophys. Res.: Solid Earth 121, 6421–6444. https://doi.org/10.1002/2016JB013321.

Steefel, C., Appelo, C., Arora, B., Jacques, D., Kalbacher, T., Kolditz, O., Lagneau, V., Lichtner, P., Mayer, K.U., Meeussen, J., et al., 2015a. Reactive transport codes for subsurface environmental simulation. Comput. Geosci. 19 (3), 445–478.

Steefel, C.I., Beckingham, L.E., Landrot, G., 2015b. Micro-continuum approaches for modeling pore-scale geochemical processes. Rev. Mineral. Geochem. 80, 217–246. https://doi.org/10.2138/rmg.2015.80.07.

Strang, G., 1968. On the construction and comparison of difference schemes. SIAM J. Numer. Anal. 5 (3), 506–517. ISSN 00361429, URL http://www.jstor.org/stable/2949700.

Tartakovsky, A.M., Meakin, P., Scheibe, T.D., West, R.M.E., 2007. Simulations of reactive transport and precipitation with smoothed particle hydrodynamics. J. Comput. Phys. 222 (2), 654–672.

Vafai, K., Tien, C., 1981. Boundary and inertia effects on flow and heat transfer in porous media. Int. J. Heat Mass Tran. 24 (2), 195–203.

Verma, A., Pruess, K., 1988. Thermohydrological conditions and silica redistribution near high-level nuclear wastes emplaced in saturated geological formations. J. Geophys. Res.: Solid Earth 93 (B2), 1159–1173. https://doi.org/10.1029/JB093iB02p01159. URL.

Whitaker, S., 1999. The Method of Volume Averaging, vol. 13. Kluwer Academic, Dordrecht of Theory and Applications of Transport in Porous Media.

Xie, M., Mayer, K.U., Claret, F., Alt-Epping, P., Jacques, D., Steefel, C., Chiaberge, C., Simunek, J., 2015. Implementation and evaluation of permeability-porosity and tortuosity-porosity relationships linked to mineral dissolution-precipitation. Comput. Geosci. 19 (3), 655–671.

Xu, T., Sonnenthal, E., Spycher, N., Pruess, K., 2006. TOUGHREACT—a simulation program for non-isothermal multiphase reactive geochemical transport in variably saturated geologic media: applications to geothermal injectivity and CO2 geological sequestration. ISSN 0098-3004 Comput. Geosci. 32 (2), 145–165. https://doi.org/10.1016/j.cageo.2005.06.014. URL http://www.sciencedirect.com/science/article/pii/S0098300405001500.