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

.


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 an comprehensive description of the coupling).The topology of the rock micro-structures is described using effective properties including porosity, tortuosity, permeabilityor hydraulic conductivityand 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 (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;Tartakovsky 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 averaging 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;Guibert 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 development of RTM at the pore-scalepioneer 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 fluid/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 simulators 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 reproducible manner with several codes.Further developments and verifications 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 systems 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 decomposition technique solves different physics on separate domainsone for Darcy flow, another for Stokes flowlinked together through appropriate boundary conditions (Molins et al., 2019) including Beaver-Joseph and Ochoa-Tapia-Whitaker conditions (Beavers and Joseph, 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 displacement 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 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.
C. Soulaine et al. 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 microcontinuum model and its ability to tend asymptotically toward continuum-scale models if grid block contains solid content (0 < φ < 1) or towards pore-scale models otherwise (φ = 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 porousMe-dia4Foam, 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 porous-Media4Foam to model hybrid-scale cases.

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 (liquidliquid 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.

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, 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, and can be updated at any moment following dissolution or precipitation processes.The porosity update at every time-step is optional.

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).

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;Bousquet-Melou et al., 2002;Goyeau et al., 2003).The momentum equation reads, where φ is the porosity, v f is the seepage velocity, p f is the fluid pressure, g is the gravity, ρ f is the fluid density, μ f is the fluid viscosity and k is the cell permeability.The porous media properties including porosity and 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, φ = 1, and the drag force μ f 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 < φ < 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, φ, 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, where ρ s,i is the solid density of mineral i and ṁ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).

darcyFoam: Darcy's law
darcyFoam is a standard continuum-scale solver that is based on Darcy's law, 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: where 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).

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.

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 continuumscale.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, where ρ s,i is the density of solid mineral i.

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/kg water ), along with equilibrium and kinetic reactions with the solid minerals described by Y s,i .phreeqcRM is set with components, i.e.SOL-UTION_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 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_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 operatorsplitting approach based on the Strang's algorithm (Strang, 1968).First, all species concentration fields, C j , are transported sequentially using the advection-dispersion equations, ∂φC j ∂t where v f is the fluid velocity computed with the flow solver (see Section 2.2) and D * j 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 finitevolume method and solved implicitly using OpenFOAM's engines.
Then, the volume fractions of solid minerals, Y s,i , 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 2 /m 3 (or m − 1 ).Hence, the RATE block provided in PHREEQC database to describe reaction rates has to be defined accordingly.

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, The mass balance equation for species A reads, where v f is the fluid velocity, D * j is an effective diffusion tensor, A s,j is the reactive surface area of mineral j (in m − 1 ), and 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 j , is defined in mol/m 3 .The equation is discretized on the computational grid using the finite-volume method and solved implicitly.
The distribution of solid minerals evolves according to,

∂Y s,i ∂t
where V ms,i (in m 3 /mol) is the molar volume of the reacting mineral.

transportOnly
transportOnly solves the advection-dispersion equation, ∂φC j ∂t without considering geochemistry.It allows the transport of species using the dispersion models implemented in porousMedia4Foam (see Table 5).

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.

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.

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

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

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

Name Expression Comments
) n n is a user defined variable.
) n n is a model parameter.φ c refers to the critical porositywhere permeability reduces to 0 (Verma and Pruess, 1988).
C. Soulaine et al. 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).

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.

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.

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, 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 × 10 − 5 mol/m 2 /s and the initial specific area A 0 = 1 m 2 /m 3 .As the calcite dissolves, the surface area decreases according to a powerlaw 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 × 10 − 11 m 2 .The column is spatially discretized with Δx = 25 mm (80 cells).150 years are simulated with Δ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 2 .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 velocityand therefore the outfluxincreases 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 4
Summary of the specific surface area models implemented in porousMedia4-Foam.Units of specific surface area are m − 1 .

Name Expression Comments
None A s = 0 If specific surface area is not necessary, e.g. for phase equilibrium calculation.

Constant
A s = A 0 -

Volume of solid
As = |∇Ys|ψ for pore-scale simulations only.
Compute thelocal surface area based on the mineral mapping.ψ is a diffuse interface function (Soulaine et al., 2017).

Power-law
As = A0(Ys) n n is a user defined variable Sugar-lump ) n3 Evolution of the surface area of an aggregateduring dissolution (Noiriel et al., 2009).A m is the maximum surface area given by the sum of the surface areas of all individual particles, n 1 , n 2 and n 3 are user-defined parameters.Hydrogeochemical coupling including surface reduction due to hydrogeochemical coupling ( Soulaine et al., 2017).n, p, q are user defined parameters.

Table 5
Summary of the dispersion models implemented in porousMedia4Foam.I is the unit tensor.those of MIN3P which demonstrates the ability of our platform to simulate hydro-geochemical processes with porosity feedback.

Verification at the pore-scale
In this section, we highlight the capabilities of our OpenFOAM package to model hydro-geochemical interactions occurring at the porescale 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 × 10 − 3 m/s.The calcite crystal dissolution is described considering a kinetic rate, where r is the reaction rate in mol/m 3 /s, A calcite is the reactive surface area in m 2 /m 3 computed using the volume-of-solid approach (see Table 4), (k calcite γ) = 0.89 × 10 − 3 m/s is the reaction rate constant of calcite and, c H + is the concentration of H + in mol/m 3 .This reaction rate may not be fully representative of the underlying geochemical processes.It has been chosen in Molins et al. (2020) to demonstraste 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  C. Soulaine et al. 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.

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 4 ) having specific reactive surface area of A 0 = 20000 m 2 /m 3 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 = 0.5 and k 0 = 10 − 12 m 2 , respectively.A solution containing 300 mol/m 3 of barium (Ba 2+ ) 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 j = 10 − 9 m 2 /s, hydrodynamic dispersion coefficient set to α L = 10 − 5 m and tortuosity exponent set to n = 2. Once the barium ions reach the porous matrix, celestite dissociates into strontium (Sr 2+ ) and sulphate (SO 2− 4 ) ions.The barium ions react with sulphate ions resulting in the precipitation of barite (BaSO 4 ) according to the following reaction (Poonoosamy et al., 2018), Celestite dissolution is taken into account considering kinetics with k celestite = 10 − 5.66 mol/m 2 /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 − 2 m/s and U inj = 10 − 6 m/s.The Péclet number, Pe = U inj ℓ/D j (where the reference lengthscale 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 4 ) 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 II = k celestite ℓ/(c Ba 2+ D j ) ≈ 3.6 × 10 − 4 (where the reference lengthscale is the inverse of the specific surface area, ℓ = A − 1 0 , 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   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 hybridscale approach.Our platform is therefore a powerful tool to complement and augment reactive transport experiments including highresolution 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.

Conclusion
We developed an integrated open-source simulator to model hydrogeochemical 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 unstructured 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 wellestablished 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 multiphase 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.

Fig. 2 .
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.

Fig. 3 .
Fig. 3.The considered model set-up along with initial and boundary conditions to investigate calcite grain dissolution in a microchannel.

Fig. 4 .
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.

Fig. 5 .
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.

Fig. 6 .Fig. 7 .
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 × 10 − 2 m/s, and (c) mineral volume fractions at 100 h and 200 h for U inj = 1 × 10 − 6 m/s.

Table 1
Summary of the flow solvers implemented in porousMedia4Foam for simulating hydro-geochemical processes.

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

Table 6
Initial and boundary conditions of the primary components for calcite dissolution under kinetic conditions considering porosity feedback.