# A benchmark study on reactive two-phase flow in porous media: Part I -model description

Stephan de Hoop, Denis Voskov, Etienne Ahusborde, Brahim Amaziane, Michel Kern

► To cite this version:

Stephan de Hoop, Denis Voskov, Etienne Ahusborde, Brahim Amaziane, Michel Kern. A bench-
mark study on reactive two-phase flow in porous media: Part I -model description. Computational
Geosciences, 2024, 10.1007/s10596-024-10268-z. http://hal-04237764

HAL Id: hal-04237764

https://hal.science/hal-04237764v1

Submitted on 11 Oct 2023

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>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A benchmark study on reactive two-phase flow in porous media:
## Part I - model description

Stephan de Hoop<fnref n="1" />, Denis Voskov<fnref n="1" />,<fnref n="2" />, Etienne Ahusborde<fnref n="3" />, Brahim Amaziane<fnref n="3" />, and
Michel Kern<fnref n="4" />,<fnref n="5" />

Emails: S.deHoop-1@tudelft.nl, D.V.Voskov@tudelft.nl, etienne.ahusborde@univ-pau.fr, brahim.amaziane@univ-pau.fr, Michel.Kern@inria.fr

October 11, 2023

## Abstract

This paper proposes a benchmark study for reactive multiphase multicomponent flow in porous media. Modeling such problem leads to a highly nonlinear coupled system of partial differential equations, ordinary differential equations and algebraic constraints which requires special numerical treatment. The benchmark consists of five test problems in total (both in 1D and in 2D), with varying degrees of difficulty, designed to verify the algorithms and the codes dedicated to simulating coupled isothermal Hydro-Chemical processes during injection and storage of CO₂ in the subsurface. It is intended to be used as a basis for comparing codes in order to better understand different couplings such as chemical reactions with two-phase flow, phase behavior with equilibrium reactions, dissolution and precipitation.

# 1 Introduction

Geological Carbon Sequestration (GCS) is viewed as one of the possible mitigation strategies for reducing worldwide carbon emissions [25]. After being captured from large carbon-emitting facilities, the gas is injected into adequate geologic formations, such as deep saline aquifers. A major concern about GCS is the prediction of plume movement to ensure subsurface retention. CO₂ injection projects should assure safe storage in the subsurface for thousands of years. Primary containment mechanisms include structural trapping, dissolution, mineralization, and capillary trapping. Simulations are required to make sure that the stored gas does not escape to the surface (see [45] for a full discussion of the issues involved).

Reactive multiphase flow models simulate the migration of CO₂ by solving a highly nonlinear system of degenerate PDEs coupled to algebraic and/or ordinary differential equations requiring special numerical treatment. From a physical point of view, the system to be described involves water (H₂O) and gas (CO₂), as well as several chemical components derived from the interaction between water, gas and the rock matrix. Both water and gas can exist in two forms: liquid and vapor. One thus has to model a two-phase system, with exchanges between the two phases (gas can dissolve in the liquid phase and water can evaporate into the gas phase). Additionally, the chemical components can react with the rock matrix to either dissolve the rocks or precipitate into one of several possible minerals (calcite being the most frequently found).

<footnote n="1">¹T.U. Delft, Department of Geoscience & Engineering, Delft, Netherlands</footnote>

<footnote n="2">²Department of Energy, Science and Engineering, School of Earth Sciences, Stanford University, CA 94305, USA</footnote>

<footnote n="3">³Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LMAP, Pau, France</footnote>

<footnote n="4">⁴Inria, Paris research Center, 2 rue Simone Iff, 75012, Paris, France</footnote>

<footnote n="5">⁵CERMICS, ENPC, 77455, Marne-la-Vallée, France</footnote>
Several methods have been developed for numerical modeling of such problems and we only give here a few references. The numerical modeling of the reactive single-phase flow model is quite well understood and many methods have been developed to solve it, see [10] for a recent review.

However, the situation is quite different for the reactive multiphase flows since such simulations pose a significant challenge due to the complexity of the coupled processes. While a majority of published works and reactive transport codes still deal with single-phase flow, studies dealing with two- or multi-phase reactive flow are now much more frequent. According to [51], “one of the most significant achievements in reactive transport analysis in recent years has been the application to multiphase systems”.

The book [64] and a more recent survey [49] describe several codes having multiphase capa-
bility: IPARS [58], PFLOTRAN [20, 33], NUFT [21], TOUGHREACT [61, 60], HYTEC [53, 46],
eSTOMP [59], MIN3P [35, 34], HYDROGEOCHEM [63] and OpenGeoSys [29, 30]. Without being
exhaustive, we may also mention some other codes dealing with reactive multiphase flow, such as
ADGPRS [16], CORE²D V4 [44], PROOST [19], DARTS [55].

Among these codes, one finds the usual two approaches for solving reactive transport problems, going back to the seminal review paper [62]: in the Sequential Approach, multiphase multi-component transport and chemistry are solved sequentially (either once per time-step or in an iterative loop, in which case one obtains the Sequential Iterative Approach (SIA)), while in the Global Implicit Approach (GIA), the chemical equations (mass action laws for equilibrium or rate laws for kinetics) and mass conservation laws are gathered and solved simultaneously in a single system of equations. Two recent surveys [47, 10] (the first one is focused on multiphase flow, while the second one mostly deals with one phase flow) include a detailed discussion of the pros and cons of SIA vs GIA, as well as a review of existing codes. They also classify the reviewed codes according to whether they follow the SIA or the GIA approach. Other works dealing with reactive transport for two-phase or multiphase flow also follow the SIA-GIA dichotomy: [43], as well as the recent work in [57], use the SIA, while [15, 27] follow a GIA.

Benchmarks (or sometimes more properly described code intercomparisons) have long been rec-
ognized as an important contribution to the community. They serve a dual purpose. First, they are
instrumental in raising awareness of the particular challenges one finds when attempting to simulate
these problems. Second, they offer well-documented test cases for researchers who enter the domain,
or who want to test new methods in their codes. One can find such benchmarks for different porous
media models, including reactive transport, see the book [30] for example. An early and influential
effort on the topic of CO₂ storage, though not focused on geochemistry, was [12]. Later, the MoMaS
Benchmark [9, 8] was concentrated on the numerical issues of reactive transport, in a single-phase
setting. The SeSBench family of benchmarks described in [50, 14] and the references therein, pro-
vides a wide-ranging set of examples for testing the codes in more realistic physical situations. A
very recent benchmark for CO₂ storage is the FluidFlower International Benchmark (see [39] for
a description of the benchmark, and [18] for a discussion of the results). That benchmark has a
different set of goals than the present proposal: it features a model for a physical lab experiment
and aims at a comparison between simulation and experimental results. On the other hand, the
physics is not fully specified, as participants are expected to provide the full set of equations of
states and other constitutive laws. The difference is between a *verification* (this proposal) and a
*validation* (FluidFlower [39]) benchmark (see for instance [40] for more on these two notions). Other
benchmarking efforts include [28, 48, 2, 42] and [22] for hydrogen storage simulations.

The aim of this paper is to present a benchmark that was designed to test the properties of
existing discretization schemes for reactive multiphase multi-component flow in porous media. The
benchmark provides all necessary data to ensure the reproducibility of the proposed tests. The
present paper is an expanded version better suited for publication of the original benchmark proposal
that was submitted to potential participants in 2021-2022, and is available at [24]. The results
obtained by the original participants are discussed in the companion paper [4].

The benchmark was written to address several challenges commonly met in applications:

1. Robust coupling of chemical reactions with multiphase flow in porous media,

2. Phase behavior coupling with equilibrium reactions,

3. Conservative treatment of solid phase dissolution and precipitation,

4. Effective coupling of equations in the case of multiple (concurrent) reactions.
The Benchmark is a set of five test cases of increasing difficulty. The first two tests are a 1D homogeneous model with two cases corresponding to whether the chemical reaction is modeled as kinetic or equilibrium. The next three test cases are based on a 2D heterogeneous model with two cases corresponding to a simple chemical model (with two sub-cases with or without gravity) and an extended chemical model composed of four reactions of which three are at equilibrium and one is kinetic.

We detail the model equations, including compositional flow, chemical reactions and how they are coupled, in Section 2. Section 3 includes a summary of the equations and a discussion of the main hypotheses. We emphasize that one major simplification made to formulate the model is that capillary pressure is not taken into account. We also provide some background information on how to deal with the issue of phase appearance and disappearance. Section 4 presents the actual benchmark test cases, their geometry, the chemical systems and specifies all the relevant data. The data are also provided in readable form on [24]. In Section 5, we recall which quantities of interest were required for the benchmark participants, so that scientists interested in running the cases with their software can easily validate their results. Some additional data are given in an Appendix.

## 2 Governing equations

This section briefly covers the governing equations of the multiphase multi-component reactive transport framework for the proposed benchmark study.

We consider a set of C fluid species (also called components), distributed over P fluid phases, and M mineral phases (each mineral phase contains only one mineral species). The species are involved in chemical reactions, of which $n_K$ are kinetic reactions, while $n_Q$ are modeled as equilibrium reactions. For fluid species, the total number of moles per unit bulk volume [mol/m³] is

$$n_c = \phi \sum_{\alpha=1}^{P} \rho_{\alpha} s_{\alpha} x_{c\alpha}, \quad c = 1, \dots, C, \quad (1)$$

where $\rho_{\alpha}$ [mol/m³] is the molar density of phase $\alpha$, $x_{c\alpha}$ [-] is the molar fraction of component $c$ in phase $\alpha$ and $\phi$ is the porosity [-].

For the solid mineral species, we denote the total number of moles per unit bulk volume in phase $m$ by $n_m$.

### 2.1 Conservation laws

The basic mass balance equations including the effect of chemical reactions as source/sink term are a natural extension of the equations for multiphase, multicomponent flow. They are by now well known, we refer the reader to any of the following references, among others [32, 38, 16, 27, 26, 3]

$$\begin{align}
\frac{\partial n_c}{\partial t} + l_c + q_c &= \sum_{k=1}^{n_K} v_{ck} r_k^K + \sum_{q=1}^{n_Q} v_{cq} r_q^Q, \tag{2} \\
c &= 1, \dots, C, \nonumber
\end{align}$$

$$\begin{align}
\frac{\partial n_m}{\partial t} &= \sum_{k=1}^{n_K} v_{mk} r_k^K + \sum_{q=1}^{n_Q} v_{mq} r_q^Q, \tag{3} \\
m &= C + 1, \dots, C + M. \nonumber
\end{align}$$

For each component $c$, $l_c$ is the total flux associated with that component, $q_c$ is the total well flow rate associated with that component, $v_{ik}$ is the stoichiometric coefficient associated with kinetic reaction $k$ for component $c$, $v_{iq}$ is the stoichiometric coefficient associated with equilibrium reaction $q$ for component $c$, $r_k^K$ [mol/m³/s] is the rate for kinetic reaction $k$ and $r_q^Q$ [mol/m³/s] is the equilibrium reaction rate for reaction $q$. We emphasize that for equilibrium reactions the rates are unknown. In practice, they are eliminated from the system by taking appropriate linear combinations of the conservation equations in (2)-(3). We come back to more explanations in Section 3.3.
The flux of component c, c = 1, ..., C is given by:

$$l_c = \nabla \cdot \sum_{\alpha=1}^{P} [\rho_{\alpha} x_{c\alpha} \mathbf{u}_{\alpha} - \phi s_{\alpha} d_{c\alpha} \nabla(\rho_{\alpha} x_{c\alpha})], \quad (4)$$

where the term $d_{c\alpha}$ [m$^2$/s] corresponds to the molecular diffusion of component c in phase $\alpha$. We note that dispersion could have been included (e.g., using Scheidegger's model [6]), but the simpler model has been retained for this benchmark. The term $\mathbf{u}_{\alpha}$ [m/s] is the Darcy velocity of phase $\alpha$ and is defined by Darcy's law:

$$\mathbf{u}_{\alpha} = -\mathbb{K} \frac{k_{r\alpha}}{\mu_{\alpha}} (\nabla p - \hat{\rho}_{\alpha} g \nabla h) \quad \alpha = 1, \dots, P, \quad (5)$$

where $\mathbb{K}$ [D]$^1$ is the rock permeability, $\hat{\rho}_{\alpha}$ [kg/m$^3$] is the mass density of phase $\alpha$, $\mu_{\alpha}$ [P]$^2$ its dynamic viscosity, $g$ [m/s$^2$] is the gravity constant, $h$ [m] is a vertical coordinate and $p$ [Pa] is the pressure. In this benchmark capillary pressure effects are ignored, so there is only a single pressure for both phases. Furthermore, $k_{r\alpha} = k_{r\alpha}(s_{\alpha})$ [-] is the relative permeability function for phase $\alpha$, and is discussed further in Section 2.4.

Equations (2)-(3) can be written in a vector form:

$$\frac{\partial \mathbf{n}}{\partial t} + \mathbf{l} + \mathbf{q} = \mathbf{V}^{\mathbf{Q}} \mathbf{r}^{\mathbf{Q}} + \mathbf{V}^{\mathbf{K}} \mathbf{r}^{\mathbf{K}}, \quad (6)$$

with vectors $\mathbf{n} = (n_1, \dots, n_{C+M})^T$, $\mathbf{l} = (l_1, \dots, l_C, 0, \dots, 0)^T$, and $\mathbf{q} = (q_1, \dots, q_C, 0, \dots, 0)^T$. The matrices $\mathbf{V}^{\mathbf{Q}}$ and $\mathbf{V}^{\mathbf{K}}$ are the stoichiometric matrices respectively for the equilibrium and kinetic reactions, while $\mathbf{r}^{\mathbf{Q}} = (r_1^Q, \dots, r_{n_Q}^Q)^T$ and $\mathbf{r}^{\mathbf{K}} = (r_1^K, \dots, r_{n_K}^K)^T$ are the equilibrium and kinetic reaction rate vectors.

The system of equations is closed with the following algebraic constraints:

$$\sum_{\alpha=1}^{P} s_{\alpha} = 1, \quad (7)$$

and

$$\sum_{c=1}^{C} x_{c\alpha} = 1, \quad \alpha = 1, \dots, P. \quad (8)$$

## 2.2 Phase equilibrium model

The following equations are used for thermodynamic equilibrium of the multicomponent system. A component is in thermodynamic equilibrium if the fugacities of the components in both phases are equal:

$$f_{c1} - f_{c\alpha} = 0, \quad c = 1, \dots, C, \quad \alpha = 2, \dots, P. \quad (9)$$

The fugacity of a component in a phase is given by

$$f_{c\alpha} = \phi_{c\alpha} x_{c\alpha} p, \quad c = 1, \dots, C, \quad \alpha = 1, \dots, P, \quad (10)$$

where $\phi_{c\alpha}$ [-] is the fugacity coefficient of an ideal mixture. Equation (9) can also be written in terms of the partition coefficients $K_{c\alpha} = \phi_{c1}/\phi_{c\alpha}$:

$$K_{c\alpha} x_{c,1} - x_{c\alpha} = 0, \quad c = 1, \dots, C, \quad \alpha = 2, \dots, P. \quad (11)$$

The phase equilibrium assumption introduces additional complexity since it poses auxiliary non-linear constraints to the solution of the system that can lead to a possible appearance or disappearance of the fluid phases. How this is treated is intimately linked to the way the nonlinear system is solved and is an important implementation choice. We come back to this question in Section 3.3.

<footnote n="1">¹the SI unit for permeability is [m$^2$] but the Darcy is commonly used, with 1 D = 9.87 × 10$^{-13}$ m$^2$</footnote>

<footnote n="2">²Poise, with 1 P=0.1 Pa s</footnote>
## 2.3 Equations for geochemistry

The other important part of the model is the chemical part, which is the specification of the kinetic and equilibrium rates. Before we describe it, we note that chemical quantities are usually given in terms of the activities of the components, which are defined in terms of the molality of the corresponding component, modified by an activity coefficient. In this work, we make the simplifying assumption of an ideal solution, so that we do not include complex activity coefficients, and the activity of the component is simply taken equal to its molality:

$$a_{cw} = M_w \left( \frac{x_{cw}}{x_{ww}} \right), \qquad (12)$$

where $M_w$ is the number of moles of H₂O per kilogram of aqueous phase (typically taken as 55.508) and $x_{ww}$ is the mole fraction of H₂O in the aqueous phase.

The kinetic rates $r_k^k$ are given functions of the activities. Since these functions are typically model-dependent, they will be specified when we describe the benchmark scenarios, see equation (21).

For the equilibrium reactions, we recall that the reaction rates themselves are unknown. However, for each equilibrium reaction, we need to add the law of mass action to either the global or the local system (depending on the preferred nonlinear formulation). The law of mass action for equilibrium reaction q is given as:

$$Q_q - K_q = \prod_{c=1}^{C} a_{cw}^{v_{cq}} - K_q = 0, \quad q = 1, \dots, n_Q. \qquad (13)$$

Here $Q_q$ is the reaction quotient, whereas $K_q$ [-] is the equilibrium reaction constant or equilibrium solubility limit for dissolution/precipitation of minerals and $a_{cw}$ is the activity of the component c in the aqueous phase.

One important consequence of the presence of reactions involving minerals is that these reactions may change the properties of the rock, both its porosity and its permeability. The most usual model [52] for letting porosity change as a result of a change in the mineral concentration is given by

$$\phi = 1 - \sum_{m=1}^{M} \frac{M_m n_m}{\hat{\rho}_m}, \qquad (14)$$

where $M_m$ [kg/mol] is the molar mass of species m, $n_t = \sum n_m$ and $\hat{\rho}_m$ [kg/m³] is its mass density. A more precise model, taking into account the presence of non-reactive minerals can be found in [26, 17].

The permeability dependence on porosity is approximated using the following power-law equation, following [23]

$$K = K_0 \left( \frac{\phi}{\phi_0} \right)^A, \qquad (15)$$

where $K_0$ and $\phi_0$ are initial porosity and permeability respectively, and A is a coefficient that takes into account the strength of the dependence.

## 2.4 Fluid and rock parameters

The relative permeability functions used in this benchmark follow the Brooks-Corey description, more precisely

$$k_{r\alpha} = k_{r,\alpha}^{e} \left( \frac{s_{\alpha} - s_{r\alpha}}{1 - \sum_{\beta \in P} s_{r\beta}} \right)^{n_{\alpha}},$$

where $k_{r\alpha}$ is the relative permeability, $k_{r\alpha}^e$ is the maximum relative permeability, $s_{r\alpha}$ is the residual saturation, and $n_{\alpha}$ is the Corey exponent of phase $\alpha$ respectively. In the absence of any residual saturation and taking $P = \{w, g\}$ (i.e., assuming liquid (water) and vapor (gas) are the fluid phases present in the system), this results in

$$k_{rw} = k_{rw}^{e} (s_w)^{n_w}, \qquad (16)$$

for the water and

$$k_{rg} = k_{rg}^{e} (1 - s_w)^{n_g}, \qquad (17)$$

for the gas relative permeability.
For the phase densities, a simple linear compressibility model is assumed, particularly

$$ \rho_{\alpha} = \rho_{\alpha,0} [1 + c_{\alpha}(p - p_0)]. \qquad (18) $$

Here $c_{\alpha}$ [1/Pa] is the compressibility of the phase $\alpha$ and $\rho_{\alpha,0}$ is its density at pressure $p_0$. This is assumed to hold for each of the three phases present in the system, water, gas, and solid. Additional physical complexity can be obtained by adopting a fully compressible model for the gas phase.

## 2.5 Initial and boundary conditions

The model obviously needs to be completed by specifying initial and boundary conditions. Since these depend both on the specific test case and on the chosen numerical formulation, as discussed in Section 3.3, they will be specified when we discuss the benchmark scenarios in Section 4

# 3 Discussion of the model

## 3.1 Summary of the main equations

The main unknowns of the system are the pressure $p$, the phase saturations $(s_{\alpha})_{\alpha=1,...,P}$, the number of moles of minerals $(n_{m})_{m=1,...,M}$ and the molar fractions $(x_{c\alpha})_{c=1,...,C, \alpha=1,...,P}$. In addition, the reaction rates for the equilibrium reactions are unknown and must be added to the system. This gives a number of unknowns equal to $1 + P + M + CP + n_Q$.

Conversely, the reactive system to be solved is composed of:

* the conservation equations (6);

* the closure relations (7) and (8);

* the phase equilibrium relations (11);

* the mass action laws for equilibrium reactions (13), with the activities given by (12).

This nonlinear system is completed by specifying:

* Darcy's law (5), together with the definition of the relative permeabilities, given by equations (16)-(17), as well as an equation of state relating the density to the pressure as in eq. (18);

* the evolution of porosity (14) and permeability (15) as a result of mineral precipitation or dissolution;

* the definition for the rates of the kinetic reactions, to be provided in Section 4.

It is worth noting that the number of equations is $C + M + 1 + P + C(P - 1) + n_Q$, which is indeed the same as the number of unknowns.

## 3.2 Main assumptions

In this section we review the main assumptions that were made with respect to a full compositional model in order to keep the model simple.

**No capillary pressure** The most important simplification is that capillary pressure is not taken into account. This choice was done for simplicity reason as capillary pressure effects add significant difficulties when one needs to solve the nonlinear problem at each iteration. If capillary pressure effects were to be added, Darcy's law (5) would need to include the phase pressure. The system would have an extra unknown, and the capillary pressure relation

$$ p_g - p_w = P_c(s_w), $$

where $P_c(s_w)$ is the capillary pressure function, can be added to close the system. For high capillary pressure (usually in tight media), the definitions of the fugacities would have to be changed accordingly. We emphasize again that the formal changes to the system are minor, but the consequences in terms of the numerical challenges would be significant;
**Equations of states and fugacities** The definition of all the functional relationships has been kept as simple as possible, while still keeping nonlinear dependence between the related quantities. This includes the relative permeability functions (16)-(17) as well as the equations of state for the density-pressure and fugacity relations (a cubic equation of state of Peng-Robinson [41] type could be used). We decided to use the simple equation of state stated above, as we felt it was more important to make sure all participants use the same equations than to let participants use their own favorite equation of state;

**Ideal chemical solution** This assumption was again made to "remove one layer of non-linearity". Equation (12) would have to be modified to include an *activity coefficient* depending on the molalities of *all* species (most models depend on the ionic strength, which is a global quantity). In chemical codes, such activity coefficients are usually evaluated in an extra outer loop, but they could also be included in the global non-linear system.

## 3.3 Different formulations and solution methods: phase appearance and disappearance

In this Section, we briefly address some common issues that, while not strictly part of the benchmark formulation, have a direct influence on the way the numerical problem will be formulated and solved.

### 3.3.1 Phase appearance and disappearance

A defining feature of compositional flows is that phases can appear or disappear, as a result of the thermodynamical exchanges between the phases. Indeed, it is quite clear that gas can dissolve into the liquid phase, while water can vaporize into the gas phase. The equations given in Section 2.1 were written as if both phases exist. If one of the phases disappears, the molar fractions for that phase are no longer meaningful. Of course, which phases exist at a given instant in time and point in space is not known in advance, and the set of present phases is itself part of the unknowns of the problem. We briefly review two of the main approaches to handle this issue. The discussion is mainly taken from the survey paper [56], to which the reader is referred for more complete information, as well as other possible formulations.

Before we describe the possible formulations, we note that the equations described in Section 2 can be split into two sets: the conservation equations (6), which are global (couple variables from different spatial locations) and the algebraic constraints, which are local at a single point in space. It is possible to use the local equations to eliminate some variables, leaving a system with just the conservation equations. The number of such equations is always equal to the number of components $C+M$, so one must identify a set of the same number of *primary* variables, while the other secondary variables can be computed based on the substitution of the secondary equations. One must be careful to select the primary equations so as to make sure the local set of equations for the secondary variables is actually solvable.

**Natural variables** This formulation is also often called the “Coats” formulation, after [13]. As long as both phases exist, the equations from Section 2.1 apply. But some mechanism must be added to decide which phases actually do exist, and to modify both the equations and the set of unknowns when one of the phases disappears. This mechanism must be applied *during* the nonlinear iterations, and it usually involves some switch of variables, as well as a stability test to check on which phases a represent at a given grid point. We refer to the above reference, as well as to [11] for further details. Alternatively, a formulation based on complementarity conditions [31, 7] or negative flash [1, 54] can be used.

As explained above, the primary unknowns for this formulation may need to change during the Newton iterations (the pressure will always be part of the primary unknowns), but we reiterate that the number of primary variables remains fixed during the whole simulation.

**Total molar variables** In this formulation, the main unknowns are pressure and the total molar fractions of fluid species, defined as

$$z_c = \sum_{\alpha=1}^{P} x_{c\alpha} \nu_{\alpha}, \quad c = 1, \dots, C, \qquad (19)$$
where the fluid phase fractions are defined as

$$ \nu_{\alpha} = \frac{\rho_{\alpha} s_{\alpha}}{\sum_{k=1}^{P} \rho_{k} s_{k}}, \quad \alpha = 1, \dots, P, $$

which obviously satisfy

$$ \sum_{\alpha} \nu_{\alpha} = 1. $$

The main advantage of this formulation as compared with the natural variables formulation is that the total molar fractions are always valid unknowns, so one does not have to use a switch of variables. On the other hand, a phase equilibrium computation (known as a “Flash” computation [36]) involving (a subset of) equations (8), (11) and (19) must be performed for each grid point, adding a notable cost to the overall solution procedure. We add that in the presence of equilibrium reactions, the flash must be modified, as explained in [26].

Note that, for the convenience of the potential participants, we have provided the initial and boundary conditions both in terms of natural variables and total molar variables. Note also that molar fraction for the C-th component can be obtained by $z_C = 1 - \sum_{\alpha=1}^{C-1} z_{\alpha}$ and is not a primary unknown (hence the initial and injection composition do not contain the composition of $z_C$), and the primary unknowns in this system are $X = [p, z_1, \dots, z_{C-1}]$.

### 3.3.2 Elimination of equilibrium reaction rates

As already noted, the reaction rates for the equilibrium equations are not specified and are in principle part of the unknowns of the problem. As their value is obviously of no interest, the usual practice has been to eliminate them by a linear combination of the conservation equations (6). This is done by introducing a kernel matrix $U$ (often called elimination matrix) such that $UV^Q = 0$, together with the auxiliary variables $\xi = \mathbf{Un}$. Given the values of $\xi$ from the flow equations, one recovers the original unknowns $\mathbf{n}$ by using the definition of $\xi$ together with the mass action laws (13). One may also inject directly the mass action laws into the conservation laws.

The procedure is by now a standard for one-phase reactive transport and has been described in various forms in numerous places (see [6] for a recent presentation in book form). It has recently been extended to deal with multi-phase flow, see [16, 27, 5, 3] for the natural variables formulation, as well as [26] for an extension to the formulation in total molar variables.

## 4 Benchmark scenarios

We describe here the specific scenarios that constitute the benchmark. The first two test cases involve a 1D geometry, while the remaining three cases are posed on a 2D heterogeneous geometry.

In all cases, there will be two fluid phases: liquid and gas phases, and one solid phase. The system consists of the following components: $[H_2O, CO_2, Ca^{+2}, CO_3^{-2}, CaCO_3]$. The components and their distribution among the phases are given in Table 1.

<table><caption>Table 1: Component-Phase distribution matrix.</caption><thead><tr><th>Component</th><th>Liquid (water)</th><th>Vapor (gas)</th><th>Solid</th></tr></thead><tbody><tr><td>H<sub>2</sub>O</td><td>✓</td><td>✓</td><td>✗</td></tr><tr><td>CO<sub>2</sub></td><td>✓</td><td>✓</td><td>✗</td></tr><tr><td>Ca<sup>+2</sup></td><td>✓</td><td>✗</td><td>✗</td></tr><tr><td>CO<sub>3</sub><sup>-2</sup></td><td>✓</td><td>✗</td><td>✗</td></tr><tr><td>CaCO<sub>3</sub></td><td>✗</td><td>✗</td><td>✓</td></tr></tbody></table>

For all but the last test case, a simple chemical system is used, with a single reaction approximating the precipitation or dissolution of calcite:

$$ CaCO_3 \rightleftharpoons Ca^{+2} + CO_3^{\cdot 2}. \qquad (20) $$
This reaction will be modeled as kinetic or equilibrium depending on the specific cases. The last test case features a more complex chemical model, with several reactions in the aqueous phase, in addition to the mineral reaction given above.

<table><caption>Table 2: Values for fluid properties.</caption><thead><tr><th>Property</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>Phase density at p<sub>0</sub>, ρ<sub>w,g,s</sub></td><td>[1000, 100, 2000]</td><td>[kg/m<sup>3</sup>]</td></tr><tr><td>Molar mass of CaCO<sub>3</sub>, M<sub>CaCO<sub>3</sub></sub></td><td>100.09</td><td>[g/mol]</td></tr><tr><td>Phase compressibility, c<sub>w,g,s</sub></td><td>[10<sup>-6</sup>, 10<sup>-4</sup>, 10<sup>-7</sup>]</td><td>[Pa]</td></tr><tr><td>Phase viscosity, μ<sub>w,g</sub></td><td>[1, 0.1]</td><td>[cP]</td></tr><tr><td>End-point relative permeability, k<sup>e</sup><sub>rw,rg</sub></td><td>[1, 1]</td><td>[-]</td></tr><tr><td>Corey exponents, n<sub>w,g</sub></td><td>[2, 2]</td><td>[-]</td></tr><tr><td>Residual saturation, s<sub>rw,rg</sub></td><td>[0, 0]</td><td>[-]</td></tr><tr><td>Phase partition coefficients, K<sub>H<sub>2</sub>O,CO<sub>2</sub></sub></td><td>[0.1, 10]</td><td>[-]</td></tr><tr><td>Diffusion coefficients, d<sub>cα</sub> = d</td><td>10<sup>-9</sup></td><td>[m<sup>2</sup>/s]</td></tr><tr><td>Activity coefficients, γ<sub>cw</sub> = γ</td><td>1</td><td>[-]</td></tr><tr><td>Porosity-permeability exponent (eq. (15)), A</td><td>3</td><td>[-]</td></tr></tbody></table>

Finally, Table 2 gives the values of the fluid parameters, which are used for all the test cases.

## 4.1 1D homogeneous domain

The first test case is a basic 1D model, shown in Figure 1. A pure CO₂ gas stream is injected from the left of a horizontal column containing pure water at a constant pressure. An outflow boundary condition, with constant pressure, is imposed at the right end of the column.

The properties describing the geometry of the domain, its discretization and the rock properties are shown in Table 3, while Table 4 summarizes the initial, injection and outflow conditions and simulation time.

<figure><img src="image_2.png" /><figcaption>Figure 1: One dimensional domain setup. Constant injection rate on the left boundary (Neumann) and constant pressure (Dirichlet) on the right. The lateral (y-direction) and vertical (z-direction) dimensions of the domain are 1[m].</figcaption></figure>

The initial and injection composition expressed in terms of the molar fraction of individual species and saturation of each phase are given in Appendix A. We distinguish two cases according to whether the reaction is modeled as kinetic or at equilibrium.
<table><caption>Table 3: Spatial data for 1D models</caption><thead><tr><th>Property</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>Initial permeability, K<sub>0</sub></td><td>100</td><td>[mD]</td></tr><tr><td>Initial porosity, φ<sub>0</sub></td><td>0.3</td><td>[-]</td></tr><tr><td>Control volume dims, Δx, Δy, Δz</td><td>[1, 1, 1]</td><td>[m]</td></tr><tr><td>Number of control volumes, N<sub>x</sub></td><td>1000</td><td>[-]</td></tr></tbody></table>

<table><caption>Table 4: Boundary conditions and other simulation parameters (1D models).</caption><thead><tr><th>Property</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>Injection rate, Q<sub>inj</sub></td><td>0.2</td><td>[m<sup>3</sup>/day]</td></tr><tr><td>Injection composition, z<sub>c,inj</sub>, c = 1, ..., C − 1,</td><td>[0, 1, 0, 0]</td><td>[-]</td></tr><tr><td>Injection composition in molar fractions x<sub>cw,ini</sub>, c = 1, ..., C − 1,</td><td>[0, 1, 0, 0]</td><td>[-]</td></tr><tr><td>Initial pressure, P<sub>ini</sub></td><td>95</td><td>[bar]</td></tr><tr><td>Initial composition, z<sub>c,ini</sub>, c = 1, ..., C − 1,</td><td>[0.15, 0, 0.075, 0.075]</td><td>[-]</td></tr><tr><td>Initial composition in molar fractions x<sub>cw,ini</sub>, c = 1, ..., C − 1,</td><td>[0.5, 0, 0.25, 0.25]</td><td>[-]</td></tr><tr><td>Outflow pressure, P<sub>out</sub></td><td>95</td><td>[bar]</td></tr><tr><td>Simulation time, T</td><td>1000</td><td>[day]</td></tr></tbody></table>

### 4.1.1 Kinetic chemistry

Here we assume that the chemical reaction (20) is kinetic. The kinetic rate (i.e., the right-hand side of equation (6)) is written as

$$r_k = A_s K_k \left( 1 - \frac{Q}{K_{sp}} \right), \qquad (21)$$

where $A_s$ [m²/m³] is the reactive surface area, which is a linear function of the solid concentration ($A_s = A_0 n_m = (1 - \phi_0)n_m$), $K_k$ [mol/m²/s] is the kinetic reaction constant, $Q$ is the activity product (to simplify $Q = x_{Ca^{2+},w} \times x_{CO_3^{2-},w}$) and $K_{sp}$ is the equilibrium constant.

The values of the reaction constants are given in Table 5. $K_{sp}$ is equal to $0.25 \times 0.25 = 0.0625$ to ensure that the initial state is in equilibrium and no dissolution occurs.

<table><caption>Table 5: Kinetic and equilibrium constants.</caption><thead><tr><th>Property</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>Kinetic constant, K<sub>k</sub></td><td>1</td><td>[kmol/m<sup>2</sup>/day]</td></tr><tr><td>Solubility constant, K<sub>sp</sub></td><td>0.0625</td><td>[-]</td></tr></tbody></table>

### 4.1.2 Equilibrium chemistry

The second test case is similar to the first one, except that now the reaction is treated as an equilibrium reaction. Mathematically, this adds an additional constraint equation of the form

$$Q - K_{sp} = 0, \qquad (22)$$

where $Q$ is the activity product of the equilibrium reaction as defined in equation (13) (which is taken here to have the same form as in Section 4.1.1) and $K_{sp}$ is the solubility constant, with the value given in Table 4. All the other parameters, fluid/rock/boundary condition/simulation parameters (as specified in Table 2 and 4), are the same as for the previous case (including the stoichiometry of the reaction).
## 4.2 2D Heterogeneous domain

The remaining test cases are set up in a two-dimensional heterogeneous domain. In the model, a zone of high porosity and permeability is embedded within a lower porosity and permeability zone. All the dimensions are stated in Figure 2. The boundary conditions are a constant injection rate on the left (bottom half of the domain pure CO₂, top half pure H₂O) and constant pressure on the right boundary (outflow) with no flow on top and bottom.

<figure><img src="image_3.png" /><figcaption>Figure 2: Configuration of the 2D test case (Section 4.2). Constant injection rate on the left boundary (Neumann) and constant pressure (Dirichlet) on the right. The lateral (y-direction) of the domain is 10m.</figcaption></figure>

Kinetic chemistry is used to model the dissolution of CaCO₃. See Tables 6 and 7 for all the parameters used in this model. Note that we have provided the initial values both in terms of the overall mole fractions and in terms of the individual mole fractions and also that the concentration of the calcite can be directly computed as a function of the porosity (see equation (14)). All fluid and chemical parameters (e.g., kinetic constants, reference mass density, etc.) are the same as in test case 4.1.1, and were given in Table 2.

<table><caption>Table 6: Spatial data for 2D model</caption><thead><tr><td>Property</td><td>Value</td><td>Units</td></tr></thead><tbody><tr><td>Porosity in Ω₁, φ</td><td>0.8</td><td>[-]</td></tr><tr><td>Permeability in Ω₁, k<sub>x,y,z</sub></td><td>[1896, 1896, 1896]</td><td>[mD]</td></tr><tr><td>Porosity in Ω₂, φ</td><td>0.3</td><td>[-]</td></tr><tr><td>Permeability in Ω₂, k<sub>x,y,z</sub></td><td>[100, 100, 100]</td><td>[mD]</td></tr><tr><td>Control volume dimension, Δx, y, z</td><td>[10, 10, 10]</td><td>[m]</td></tr><tr><td>Number of control volumes, N<sub>x</sub> × N<sub>y</sub> × N<sub>z</sub></td><td>120 × 1 × 48</td><td>[-]</td></tr><tr><td>Diffusion coefficients, d<sub>cα</sub> = d</td><td>10<sup>-9</sup></td><td>[m<sup>2</sup>/s]</td></tr><tr><td>Gravitational acceleration, g</td><td>9.8</td><td>[m/s<sup>2</sup>]</td></tr></tbody></table>

In addition to the two chemical systems described later, this test case can be executed with or without gravity (i.e., g = 0). We note that when gravity is included, it would have been more natural for the initial pressure to follow a hydrostatic law. However, the effect is quite small, and the simpler constant initial pressure was retained.
<table><caption>Table 7: Boundary conditions and other simulation parameters (2D model).</caption><thead><tr><th>Property</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>Gas injection rate, Q<sub>inj</sub></td><td>1000</td><td>[m<sup>3</sup>/day]</td></tr><tr><td>Water injection rate, Q<sub>inj</sub></td><td>200</td><td>[m<sup>3</sup>/day]</td></tr><tr><td>Gas injection composition, z<sub>c,inj</sub>, c = 1, ..., C − 1,</td><td>[0, 1, 0, 0]</td><td>[-]</td></tr><tr><td>Gas injection composition in molar fractions, x<sub>cg,inj</sub>, c = 1, ..., C − 1,</td><td>[0, 1, 0, 0]</td><td>[-]</td></tr><tr><td>Water injection composition, z<sub>c,inj</sub>, c = 1, ..., C − 1,</td><td>[1, 0, 0, 0]</td><td>[-]</td></tr><tr><td>Water injection composition in molar fractions, x<sub>cw,inj</sub>, c = 1, ..., C − 1,</td><td>[1, 0, 0, 0]</td><td>[-]</td></tr><tr><td>Initial pressure in Ω<sub>1</sub> ∪ Ω<sub>2</sub>, P<sub>ini</sub></td><td>95</td><td>[bar]</td></tr><tr><td>Initial composition in Ω<sub>1</sub>, z<sub>c,ini</sub>, c = 1, ..., C − 1,</td><td>[0.4, 0, 0.20, 0.20]</td><td>[-]</td></tr><tr><td>Initial composition in Ω<sub>2</sub>, z<sub>c,ini</sub>, c = 1, ..., C − 1,</td><td>[0.15, 0, 0.075, 0.075]</td><td>[-]</td></tr><tr><td>Initial fluid composition in Ω<sub>1</sub> ∪ Ω<sub>2</sub> in molar fractions, x<sub>cw,ini</sub>, c = 1, ..., C − 1,</td><td>[0.5, 0, 0.25, 0.25]</td><td>[-]</td></tr><tr><td>Outflow pressure, P<sub>out</sub></td><td>95</td><td>[bar]</td></tr><tr><td>Simulation time, T</td><td>1000</td><td>[days]</td></tr></tbody></table>

### 4.2.1 Simple chemical model

Here, only one chemical reaction is included, and the chemical model is the same as in Section 4.1.1.

One further simplification may be necessary: one may encounter convergence problems when running the system as described previously. If that is the case, it may be helpful to consider a single “meta-ion” Ca<sup>2+</sup> + CO<sub>3</sub><sup>2−</sup> in the liquid phase.

### 4.2.2 Extended chemical model

We consider a somewhat more realistic chemical system, including the dissociation of water and carbon dioxide, as this makes it possible to take into account the influence of pH.

The system is composed of 4 reactions:

$$ 
\begin{align} 
& H_2O \rightleftharpoons H^+ + OH^-, \nonumber \\ 
& CO_2 + H_2O \rightleftharpoons HCO_3^- + H^+, \nonumber \\ 
& HCO_3^- \rightleftharpoons H^+ + CO_3^{2-}, \nonumber \\ 
& CaCO_3 + H^+ \rightleftharpoons Ca^{2+} + HCO_3^-. \tag{23}
\end{align} 
$$

Note that this increases the total number of species by three, particularly to z<sub>c</sub> = [H<sub>2</sub>O, CO<sub>2</sub>, Ca<sup>2+</sup>, CO<sub>3</sub><sup>2−</sup>, H<sup>+</sup>, OH<sup>−</sup>, HCO<sub>3</sub><sup>−</sup>, CaCO<sub>3</sub>]. However, this makes it possible to represent the CaCO<sub>3</sub> reaction (last equation in (23)) with an explicit dependency on the pH of the solution.

The first three reactions are at equilibrium, the logarithms of the equilibrium constants are given in Table 8, while the fourth reaction is kinetic, and the rate for the last reaction is given by equation (21), with Q now defined as:

$$ Q = \frac{a_{Ca^{2+},w} a_{HCO_3^-,w}}{a_{H^+,w}}. \qquad (24) $$

<table><caption>Table 8: Log<sub>10</sub> of the equilibrium constants for extended chemical system (23).</caption><thead><tr><th>K<sub>1</sub></th><th>K<sub>2</sub></th><th>K<sub>3</sub></th><th>K<sub>sp</sub></th></tr></thead><tbody><tr><td>-13.95</td><td>-6.293</td><td>-10.279</td><td>-1.899</td></tr></tbody></table>
For this last test case, the mass actions laws are expressed in activities:

$$
\begin{align*}
K_1 \ a_{H_2O,w} &= a_{H^+,w} \ a_{OH^-,w}, \\
K_2 \ a_{CO_2,w}a_{H_2O,w} &= a_{HCO_3^-,w} \ a_{H^+,w}, \\
K_3 \ a_{HCO_3^-,w} &= a_{H^+,w} \ a_{CO_3^{2-},w}.
\end{align*}
\tag{25}
$$

Initial conditions (computed using Arxim [37]) are given in Table 9. All other physical parameters are the same as in test 4.2.1.

<table><caption>Table 9: Initial molar fractions for the extended chemical model</caption>
  <thead>
    <tr>
      <th>Property</th>
      <th>Value</th>
      <th>Units</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>x<sub>CO<sub>2</sub>,w</sub></td>
      <td>3.9624 &times; 10<sup>-10</sup></td>
      <td>[-]</td>
    </tr>
    <tr>
      <td>x<sub>Ca<sup>2+</sup>,w</sub></td>
      <td>2.1703 &times; 10<sup>-6</sup></td>
      <td>[-]</td>
    </tr>
    <tr>
      <td>x<sub>H<sup>+</sup>,w</sub></td>
      <td>2.3507 &times; 10<sup>-12</sup></td>
      <td>[-]</td>
    </tr>
    <tr>
      <td>x<sub>OH<sup>-</sup>,w</sub></td>
      <td>1.5475 &times; 10<sup>-6</sup></td>
      <td>[-]</td>
    </tr>
    <tr>
      <td>x<sub>HCO<sub>3</sub><sup>-</sup>,w</sub></td>
      <td>1.5467 &times; 10<sup>-6</sup></td>
      <td>[-]</td>
    </tr>
    <tr>
      <td>x<sub>CO<sub>3</sub><sup>2-</sup>,w</sub></td>
      <td>6.2315 &times; 10<sup>-7</sup></td>
      <td>[-]</td>
    </tr>
  </tbody>
</table>

5 Expected output

In this section, we gather the quantities that were part of the output that participating teams had to provide. The corresponding files are now available on a dedicated web site https://github.com/ eahusbor/Reactive-Multiphase-Benchmark. For 1D quantities (either as functions of space or as functions of time), the data files (in a CSV format specified on the web site) were to be provided. For 2D output, only contour plots were requested. Researchers who would like to compare their methods on the models presented in Section 4 will know which quantities of interest they need to output.

5.1 1D homogeneous domain

The following quantities were required: S<sub>g</sub>, p, φ, x<sub>H₂O,w</sub>, x<sub>CO₂,w</sub>, x<sub>Ca²⁺,w</sub> + x<sub>HCO₃⁻,w</sub>

• as a function of space at time t = 1000 days;

• as a function of time at x = 25 m.

5.2 2D heterogeneous domain

5.2.1 Simple chemical model

The participants were asked to provide the following quantities at time 1000 days:

• 2D plots of the CO₂ fraction x<sub>CO₂,w</sub>, the gas saturation S<sub>g</sub> and the porosity φ;

• 1D output of S<sub>g</sub>, p, φ, x<sub>H₂O,w</sub>, x<sub>CO₂,w</sub>, x<sub>Ca²⁺,w</sub> + x<sub>HCO₃⁻,w</sub>, along the vertical line x = 40 m and along the vertical line y = 50 m.

5.2.2 Extended chemical model

The participants were asked to provide the following quantities at time 1000 days:

• 2D plots of the CO₂ fraction x<sub>CO₂,w</sub>, the gas saturation S<sub>g</sub>, the porosity φ and the (logarithm) of the mole fraction for all ions;

• 1D plots of S<sub>g</sub>, p, x<sub>H₂O,w</sub>, x<sub>CO₂,w</sub>, x<sub>Ca²⁺,w</sub>, x<sub>H⁺,w</sub>, x<sub>HCO₃⁻,w</sub>, x<sub>CO₃^{2-},w</sub>, φ along the vertical line x = 40 m and along the vertical line.
## Acknowledgments

The work of E. Ahusborde and B. Amaziane has been partly supported by the Carnot ISIFOR Institute, and “la Région Nouvelle-Aquitaine”, France. These supports are gratefully acknowledged. The authors are grateful to E. Flauraud (IFPEN, France) for providing the initial values for the “Extended chemical model” in Section 4.2.2. We wish also to thank all of the teams who took part in the Benchmark presented herein for their very active participation, as well as for a careful reading of the paper: M. El Ossmani (UPPA, France), E. Flauraud (IFPEN, France), F. Hamon (TotalEnergies), A. Socié, K. U. Mayer, D. Su (University of British Columbia, Canada), M. Tóth (University of Heidelberg, Germany). The authors also thank the anonymous reviewers whose comments led to significant improvements to the paper.

## Data availability

The original statement of the benchmark, as well as a set of files containing the data given in Tables 2 to 9 and Tables 10 to 15 are available from doi:10.5281/zenodo.8319570

# Appendix A Molar composition for initial and injection conditions

Tables 10-15 give the initial and injection composition in terms of molar fraction of individual species and saturation of each phases.

<table><caption>Table 10: Properties based on 1D injection state: $[P, z_{H_2O}, z_{CO_2}, z_{Ca}, z_{CO_3}] =$ [165, 10⁻¹², 1, 10⁻¹², 10⁻¹²]</caption><thead><tr><th></th><th>H<sub>2</sub>O</th><th>CO<sub>2</sub></th><th>Ca<sup>2+</sup></th><th>CO<sub>3</sub><sup>2-</sup></th><th>CaCO<sub>3</sub></th></tr></thead><tbody><tr><th>Composition, z<sub>c</sub></th><td>10<sup>-12</sup></td><td>1</td><td>10<sup>-12</sup></td><td>10<sup>-12</sup></td><td>10<sup>-12</sup></td></tr><tr><th>Liquid molar fraction</th><td>10<sup>-11</sup></td><td>0.01</td><td>0.494</td><td>0.494</td><td>0</td></tr><tr><th>Vapor molar fraction</th><td>10<sup>-11</sup></td><td>1</td><td>4.94 &times; 10<sup>-13</sup></td><td>4.94 &times; 10<sup>-13</sup></td><td>0</td></tr><tr><th>Solid molar fraction</th><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td></tr></tbody></table>

<table><caption>Table 11: Properties based on 1D initial state: $[P, z_{H_2O}, z_{CO_2}, z_{Ca}, z_{CO_3}] =$ [95, 0.15, 10⁻¹², 0.075, 0.075]</caption><thead><tr><th></th><th>H<sub>2</sub>O</th><th>CO<sub>2</sub></th><th>Ca<sup>2+</sup></th><th>CO<sub>3</sub><sup>2-</sup></th><th>CaCO<sub>3</sub></th></tr></thead><tbody><tr><th>Composition, z<sub>c</sub></th><td>0.15</td><td>10<sup>-12</sup></td><td>0.075</td><td>0.075</td><td>0.7</td></tr><tr><th>Liquid molar fraction</th><td>0.5</td><td>3.33 &times; 10<sup>-12</sup></td><td>0.25</td><td>0.25</td><td>0</td></tr><tr><th>Vapor molar fraction</th><td>10<sup>-11</sup></td><td>1</td><td>4.94 &times; 10<sup>-13</sup></td><td>4.94 &times; 10<sup>-13</sup></td><td>0</td></tr><tr><th>Solid molar fraction</th><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td></tr></tbody></table>

<table><caption>Table 12: Properties based on 2D initial in $\Omega_2$ state: $[P, z_{H_2O}, z_{CO_2}, z_{Ions}] =$ [95, 0.15, 10⁻¹², 0.15]</caption><thead><tr><th></th><th>H<sub>2</sub>O</th><th>CO<sub>2</sub></th><th>Ions</th><th>CaCO<sub>3</sub></th></tr></thead><tbody><tr><th>Composition, z<sub>c</sub></th><td>0.15</td><td>10<sup>-12</sup></td><td>0.15</td><td>0.7</td></tr><tr><th>Liquid molar fraction</th><td>0.5</td><td>3.33 &times; 10<sup>-12</sup></td><td>0.5</td><td>0</td></tr><tr><th>Vapor molar fraction</th><td>0.045454</td><td>0.95454</td><td>4.5 &times; 10<sup>-13</sup></td><td>0</td></tr><tr><th>Solid molar fraction</th><td>0</td><td>0</td><td>0</td><td>1</td></tr></tbody></table>
<table><caption>Table 13: Properties based on 2D initial in $\Omega_1$ state: $[P, z_{H_2O}, z_{CO_2}, z_{Ions}] = [95, 0.4, 2.66 \times 10^{-12}, 0.4]$</caption><thead><tr><th></th><th>H<sub>2</sub>O</th><th>CO<sub>2</sub></th><th>Ions</th><th>CaCO<sub>3</sub></th></tr></thead><tbody><tr><th>Composition, z<sub>c</sub></th><td>0.4</td><td>2.66 × 10<sup>-12</sup></td><td>0.4</td><td>0.2</td></tr><tr><th>Liquid molar fraction</th><td>0.5</td><td>3.33 × 10<sup>-12</sup></td><td>0.5</td><td>0</td></tr><tr><th>Vapor molar fraction</th><td>0.045454</td><td>0.95454</td><td>4.5 × 10<sup>-13</sup></td><td>0</td></tr><tr><th>Solid molar fraction</th><td>0</td><td>0</td><td>0</td><td>1</td></tr></tbody></table>

<table><caption>Table 14: Properties based on 2D injection (water stream) state: $[P, z_{H_2O}, z_{CO_2}, z_{Ions}] = [95, 1, 2 \times 10^{-12}, 10^{-12}]$</caption><thead><tr><th></th><th>H<sub>2</sub>O</th><th>CO<sub>2</sub></th><th>Ions</th><th>CaCO<sub>3</sub></th></tr></thead><tbody><tr><th>Composition, z<sub>c</sub></th><td>1</td><td>2 × 10<sup>-12</sup></td><td>10<sup>-12</sup></td><td>10<sup>-14</sup></td></tr><tr><th>Liquid molar fraction</th><td>1</td><td>2 × 10<sup>-12</sup></td><td>0</td><td>0</td></tr><tr><th>Vapor molar fraction</th><td>0.909</td><td>0</td><td>0.909</td><td>0</td></tr><tr><th>Solid molar fraction</th><td>0</td><td>0</td><td>0</td><td>1</td></tr></tbody></table>

<table><caption>Table 15: Properties based on 2D injection (gas stream) state: $[P, z_{H_2O}, z_{CO_2}, z_{Ions}] = [95, 0, 1, 10^{-12}]$</caption><thead><tr><th></th><th>H<sub>2</sub>O</th><th>CO<sub>2</sub></th><th>Ions</th><th>CaCO<sub>3</sub></th></tr></thead><tbody><tr><th>Composition, z<sub>c</sub></th><td>0</td><td>1</td><td>10<sup>-12</sup></td><td>10<sup>-14</sup></td></tr><tr><th>Liquid molar fraction</th><td>0</td><td>0.1</td><td>0.9</td><td>0</td></tr><tr><th>Vapor molar fraction</th><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><th>Solid molar fraction</th><td>0</td><td>0</td><td>0</td><td>1</td></tr></tbody></table>

## References

[1] A. Abadpour and M. Panfilov. "Method of negative saturations for modeling two-phase compositional flow with oversaturated zones". *Transport in Porous Media* 79.2 (2009), pp. 197-214. DOI: 10.1007/s11242-008-9310-0.

[2] J. Águila, V. Montoya, J. Samper, et al. "Modeling cesium migration through Opalinus clay: a benchmark for single- and multi-species sorption-diffusion models". *Comput. Geosci.* 25 (2021), pp. 1405-1436. DOI: 10.1007/s10596-021-10050-5.

[3] E. Ahusborde, B. Amaziane, and M. Id Moulay. "High performance computing of 3D reactive multiphase flow in porous media: Application to geological storage of CO<sub>2</sub>". *Computational Geosciences* 25 (2021), pp. 2131-2147. DOI: 10.1007/s10596-021-10082-x.

[4] E. Ahusborde, B. Amaziane, S. deHoop, et al. "A benchmark study on reactive two-phase flow in porous media: Part II - results and description". *Comput. Geosci., under review* (2023).

[5] E. Ahusborde, M. Kern, and V. Vostrikov. "Numerical simulation of two-phase multicomponent flow with reactive transport in porous media: application to geological sequestration of CO<sub>2</sub>". *ESAIM: Proc.* 50 (2015), pp. 21-39. DOI: 10.1051/proc/201550002.

[6] J. Bear and A. H.-D. Cheng. *Modeling Groundwater Flow and Contaminant Transport*. Heidelberg, New-York: Springer, 2010. ISBN: 978-1-4020-6681-8. DOI: 10.1007/978-1-4020-6682-5.

[7] I. Ben Gharbia and E. Flauraud. "Study of compositional multiphase flow formulation using complementarity conditions". *Oil Gas Sci. Technol. - Rev. IFP Energies nouvelles* 74 (2019), p. 43. DOI: 10.2516/ogst/2019012.
[8] J. Carrayrou, J. Hoffmann, P. Knabner, et al. "Comparison of numerical methods for simulating strongly nonlinear and heterogeneous reactive transport problems-the Mo-MaS benchmark case". *Comput. Geosci.* 14 (2010), pp. 483–502. DOI: 10.1007/s10596-010-9178-2.

[9] J. Carrayrou, M. Kern, and P. Knabner. "Reactive transport benchmark of MoMaS". *Comput. Geosci.* 14 (2010), pp. 385–392. DOI: 10.1007/s10596-009-9157-7.

[10] J. Carrera, M. W. Saaltink, J. Soler-Sagarra, J. Wang, and C. Valhondo. "Reactive Transport: A Review of Basic Concepts with Emphasis on Biochemical Processes". *Energies* 15.3 (2022). ISSN: 1996-1073. DOI: 10.3390/en15030925.

[11] H. Class, R. Helmig, and P. Bastian. "Numerical simulation of non-isothermal multi-phase multicomponent processes in porous media.: 1. An efficient solution technique". *Adv. Water Res.* 25.5 (2002), pp. 533–550. ISSN: 0309-1708. DOI: 10.1016/S0309-1708(02)00014-3.

[12] H. Class, A. Ebigbo, R. Helmig, et al. "A benchmark study on problems related to CO₂ storage in geologic formations". *Comput. Geosci.* 13 (2009), pp. 409–434. DOI: 10.1007/s10596-009-9146-x.

[13] K. H. Coats. "An Equation of State Compositional Model". *Society of Petroleum Engineers Journal* 20.05 (Oct. 1980), pp. 363–376. ISSN: 0197-7520. DOI: 10.2118/8284-PA.

[14] D. Dwivedi, B. Arora, S. Molins, and C. I. Steefel. "Benchmarking reactive transport codes for subsurface environmental problems". *Groundwater Assessment, Modeling, and Management*. Ed. by M. Thangarajan and V. P. Singh. Boca-Raton, USA: CRC Press, 2016, pp. 299–316. DOI: 10.1201/9781315369044.

[15] Y. Fan, L. J. Durlofsky, and H. A. Tchelepi. "A fully-coupled flow-reactive-transport formulation based on element conservation, with application to CO₂ storage simulations". *Adv. Water Resour.* 42 (2012), pp. 47–61. DOI: 10.1016/j.advwatres.2012.03.012.

[16] S. F. Farshidi, Y. Fan, L. J. Durlofsky, and H. A. Tchelepi. "Chemical Reaction Modeling in a Compositional Reservoir-Simulation Framework". SPE Reservoir Simulation Conference. 2013. DOI: 10.2118/163677-MS.

[17] S. F. Farshidi. "Compositional reservoir simulation-based reactive-transport formulations, with application to CO₂ storage in sandstone and ultramafic formations". PhD thesis. Stanford University, 2016. URL: https://pangea.stanford.edu/ERE/db/pereports/record_detail.php?filename=Farshidi2016.pdf.

[18] B. Flemisch, J. M. Nordbotten, M. Fernø, et al. *The FluidFlower International Benchmark Study: Process, Modeling Results, and Comparison to Experimental Data*. 2023. DOI: 10.48550/ARXIV.2302.10986.

[19] P. Gamazo, L. J. Slooten, J. Carrera, et al. "PROOST: object-oriented approach to multiphase reactive transport modeling in porous media". *J. Hydroinf.* 18.2 (2015), pp. 310–328. DOI: 10.2166/hydro.2015.126.

[20] G. Hammond, P. Lichtner, C. Lu, and R. Mills. "PFLOTRAN: Reactive flow & transport code for use on laptops to leadership-class supercomputers". *Groundwater Reactive Transport Models*. Ed. by F. Zhang, G. Yeh, and J. Parker. United Arab Emirates: Bentham Publishers, 2012, pp. 141–159. DOI: 10.2174/978160805306311201010141.

[21] Y. Hao, Y. Sun, and J. Nitao. "Overview of NUFT: A versatile numerical model for simulating flow and reactive transport in porous media". *Groundwater Reactive Transport Models*. Ed. by F. Zhang, G. Yeh, and J. Parker. United Arab Emirates: Bentham Publishers, 2012, pp. 212–239. DOI: 10.2174/978160805306311201010212.

[22] S. Hogeweg, G. Strobel, and B. Hagemann. "Benchmark study for the simulation of Underground Hydrogen Storage operations". *Comput. Geosci.* 26 (2022), pp. 1367–1378. DOI: 10.1007/s10596-022-10163-5.
[23] J. Hommel, E. Coltman, and H. Class. "Porosity-Permeability Relations for Evolving Pore Space: A Review with a Focus on (Bio-)geochemically Altered Porous Media". *Transp Porous Med* 124 (2018), pp. 589-629. DOI: 10.1007/s11242-018-1086-2.

[24] S. de Hoop, D. Voskov, E. Ahusborde, B. Amaziane, and M. Kern. *Reactive Multiphase Flow in Porous Media at the Darcy Scale: a Benchmark proposal*. 2023. DOI: 10.5281/zenodo.8319570.

[25] IPCC. *Special report on carbon dioxide capture and storage*. Ed. by B. Metz, O. Davidson, H. de Coninck, M. Loos, and L. Meyer. Prepared by working group III of the Intergovernmental Panel on Climate Change. Cambridge, UK: Cambridge University Press, 2005, p. 431. URL: https://www.ipcc.ch/report/carbon-dioxide-capture-and-storage/.

[26] K. Kala and D. Voskov. "Element balance formulation in reactive compositional flow and transport with parameterization technique". *Comput.Geosci.* 24.2 (2020), pp. 609-624. DOI: 10.1007/s10596-019-9828-y.

[27] M. Knodel, S. Kräutle, and P. Knabner. "Global implicit solver for multiphase multi-component flow in porous media with multiple gas components and general reactions". *Comput. Geosci.* 26 (2022), pp. 697-724. DOI: 10.1007/s10596-022-10140-y.

[28] O. Kolditz, S. Bauer, C. Beyer, et al. "A systematic benchmarking approach for geologic CO2 injection and storage". *Environ. Earth Sci.* 67.2 (2012), pp. 613-632. DOI: 10.1007/s12665-012-1656-5.

[29] O. Kolditz, S. Bauer, L. Bilke, et al. "OpenGeoSys: An open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media". *Environ. Earth Sci.* 67 (2012), pp. 589-599. DOI: 10.1007/s12665-012-1546-x.

[30] O. Kolditz, U.-J. Görke, H. Shao, and W. Wang, eds. *Thermo-Hydro-Mechanical-Chemical Processes in Porous Media. Benchmarks and Examples*. Vol. 86. Lecture Notes in Computational Science and Engineering. Berlin, Heidelberg: Springer, 2012. DOI: 10.1007/978-3-642-27177-9.

[31] A. Lauser, C. Hager, R. Helmig, and B. Wohlmuth. "A new approach for phase transitions in miscible multi-phase flow in porous media". *Advances in Water Resources* 34.8 (2011), pp. 957-966. ISSN: 0309-1708. DOI: 10.1016/j.advwatres.2011.04.021.

[32] P. C. Lichtner. "Reactive Transport in Porous Media". Ed. by P. C. Lichtner, C. I. Steefel, and E. H. Oelkers. Berlin, Boston: De Gruyter, 1996. Chap. Continuum Formulation of Multicomponent Multiphase Reactive Transport, pp. 1-82. ISBN: 9781501509797. DOI: 10.1515/9781501509797-004.

[33] P. C. Lichtner, G. E. Hammond, C. Lu, et al. *PFLOTRAN Web page*. 2020. URL: http://www.pflotran.org.

[34] K. U. Mayer, R. T. Amos, S. Molins, and F. Gérard. *Groundwater Reactive Transport Models*. Ed. by F. Zhang, G. Yeh, and J. Parker. United Arab Emirates: Bentham Publishers, 2012. Chap. Reactive Transport Modeling in Variably Saturated Media with MIN3P: Basic Model Formulation and Model Enhancements, pp. 186-211. DOI: 10.2174/978160805306311201010186.

[35] K. U. Mayer, E. O. Frind, and D. W. Blowes. "Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions". *Water Resour. Res.* (2002), 38 (9), 1-21. DOI: 10.1029/2001WR000862.

[36] M. L. Michelsen. "The isothermal flash problem. Part II. Phase-split calculation". *Fluid Phase Equilibria* 9.1 (1982), pp. 21-40. ISSN: 0378-3812. DOI: 10.1016/0378-3812(82)85002-4.
[37] J. Moutte, A. Michel, G. Battaia, et al. "Arxim, a library for thermodynamic modeling of reactive heterogeneous systems, with applications to the simulation of fluid-rock system". *21st Congress of IUPAC. Conference on Chemical Thermodynamic*. Tsukuba, Japan, 2010.

[38] L. Nghiem, P. Sammon, J. Grabenstetter, and H. Ohkuma. "Modeling CO2 Storage in Aquifers with a Fully-Coupled Geochemical EOS Compositional Simulator". SPE Improved Oil Recovery Conference. Apr. 2004, paper SPE-89474-MS. DOI: 10.2118/89474-MS.

[39] J. M. Nordbotten, M. Fernø, B. Flemisch, R. Juanes, and M. Jørgensen. *Final Benchmark Description: FluidFlower International Benchmark Study*. 2022. DOI: 10.5281/zenodo.6807102.

[40] W. L. Oberkampf, T. G. Trucano, and C. Hirsch. "Verification, validation, and predictive capability in computational engineering and physics". *Appl. Mech. Rev.* 57.5 (Dec. 2004), pp. 345-384. ISSN: 0003-6900. DOI: 10.1115/1.1767847.

[41] D.-Y. Peng and D. B. Robinson. "A New Two-Constant Equation of State". *Industrial & Engineering Chemistry Fundamentals* 15.1 (1976), pp. 59-64. DOI: 10.1021/i160057a011.

[42] J. Poonoosamy, C. Wanner, P. Alt Epping, et al. "Benchmarking of reactive transport codes for 2D simulations with mineral dissolution-precipitation reactions and feedback on transport parameters." *Comput. Geosci.* 25 (2021), pp. 1337-1358. DOI: 10.1007/s10596-018-9793-x.

[43] M. Saaltink, V. Vilarrasa, F. De Gaspari, et al. "A method for incorporating equilibrium chemical reactions into multiphase flow models for CO₂ storage". *Adv. Water Resour.* 62 (2013), pp. 431-441. DOI: 10.1016/j.advwatres.2013.09.013.

[44] J. Samper, C. Yang, L. Zheng, et al. "CORE2D V4: A code for water flow, heat and solute transport, geochemical reactions, and microbial processes". *Groundwater Reactive Transport Models*. Ed. by F. Zhang, G. Yeh, and J. Parker. United Arab Emirates: Bentham Publishers, 2012, pp. 160-185. DOI: 10.2174/978160805306311201010160.

[45] D. Scheer, H. Class, and B. Flemisch. "Geologic Carbon Sequestration". *Subsurface Environmental Modelling Between Science and Policy*. Switzerland: Springer Cham, 2021, pp. 109-152. DOI: 10.1007/978-3-030-51178-4.

[46] I. Sin, V. Lagneau, and J. Corvisier. "Integrating a compressible multicomponent two-phase flow into an existing reactive transport simulator." *Adv. Water Resour.* 100 (2017), pp. 62-77. DOI: 10.1016/j.advwatres.2016.11.014.

[47] I. Sin and J. Corvisier. "Multiphase Multicomponent Reactive Transport and Flow Modeling". *Rev. Mineral. Geochem.* 85.1 (2019), pp. 143-195. DOI: 10.2138/rmg.2019.85.6.

[48] I. Sin, V. Lagneau, L. De Windt, and J. Corvisier. "2D simulation of natural gas reservoir by two-phase multicomponent reactive flow and transport—Description of a benchmarking exercise". *Math. Comput. Simul.* 137 (2017), pp. 431-447. DOI: 10.1016/j.matcom.2016.12.003.

[49] C. I. Steefel, C. A. J. Appelo, B. Arora, et al. "Reactive transport codes for subsurface environmental simulation". *Comput. Geosci.* 19.3 (2015), pp. 445-478. DOI: 10.1007/s10596-014-9443-x.

[50] C. Steefel, S. Yabusaki, and K. Mayer. "Reactive transport benchmarks for subsurface environmental simulation". *Comput. Geosci.* 19 (2015), pp. 439-443. DOI: 10.1007/s10596-015-9499-2.
[51] C. I. Steefel. "Reactive Transport at the Crossroads". *Rev. Mineral. Geochem.* 85.1 (2019), pp. 1-26. DOI: 10.2138/rmg.2019.85.1.

[52] C. I. Steefel and A. C. Lasaga. "A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow in single phase hydrothermal systems". *American Journal of Science* 294.5 (May 1, 1994), pp. 529-592. DOI: 10.2475/ajs.294.5.529.

[53] J. van der Lee, L. De Windt, V. Lagneau, and P. Goblet. "Module-oriented modeling of reactive transport with HYTEC". *Computers & Geosciences* 29.3 (2003), pp. 265-275. DOI: 10.1016/S0098-3004(03)00004-9.

[54] D. Voskov. "An Extended Natural Variable Formulation for Compositional Simulation Based on Tie-Line Parameterization". *Transport in Porous Media* 92.3 (2012), pp. 541-557. DOI: 10.1007/s11242-011-9919-2.

[55] D. Voskov, I. Saifullin, M. Wapperom, et al. *open Delft Advanced Research Terra Simulator (open-DARTS)*. Version 1.0.2. July 2023. DOI: 10.5281/zenodo.8116928.

[56] D. V. Voskov and H. A. Tchelepi. "Comparison of nonlinear formulations for two-phase multi-component EoS based simulation". *Journal of Petroleum Science and Engineering* 82-83 (2012), pp. 101-111. ISSN: 0920-4105. DOI: 10.1016/j.petrol.2011.10.012.

[57] Y. Wang, D. Fernández-Garcia, and M. W. Saaltink. "Modeling reactive multi-component multi-phase flow for Geological Carbon Sequestration (GCS) with Matlab". *Computers & Geosciences* 172 (2023), p. 105300. DOI: 10.1016/j.cageo.2023.105300.

[58] M. Wheeler, S. Sun, and S. Thomas. "Modeling of flow and reactive transport in IPARS". *Groundwater Reactive Transport Models*. Ed. by F. Zhang, G. Yeh, and J. Parker. United Arab Emirates: Bentham Publishers, 2012, pp. 42-73. DOI: 10.2174/978160805306311201010042.

[59] M. White and M. Oostrom. *STOMP subsurface transport over multiple phases version 4.0 user's guide*. Tech. rep. Pacific Northwest National Laboratory, 2006. URL: https://stomp-userguide.pnnl.gov/.

[60] T. Xu, E. Sonnenthal, N. Spycher, et al. "TOUGHREACT: A simulation program for subsurface reactive chemical transport under non-isothermal multiphase flow conditions". *Groundwater Reactive Transport Models*. Ed. by F. Zhang, G. Yeh, and J. Parker. United Arab Emirates: Bentham Publishers, 2012, pp. 74-95. DOI: 10.2174/978160805306311201010074.

[61] T. Xu, E. Sonnenthal, N. Spycher, and K. Pruess. "TOUGHREACT – A simulation program for non-isothermal multiphase reactive geochemical transport in variably saturated geologic media: Applications to geothermal injectivity and CO₂ geological sequestration". *Computers and Geosciences* 32.2 (2006), pp. 145-165. DOI: 10.1016/j.cageo.2005.06.014.

[62] G. T. Yeh and V. S. Tripathi. "A critical evaluation of recent developments in hydrogeochemical transport models of reactive multichemical components". *Water Resour. Res.* 25 (1989), pp. 93-108. DOI: 10.1029/WR025i001p00093.

[63] G. Yeh, V. Tripathi, J. Gwo, et al. "HYDROGEOCHEM: A coupled model of variably saturated flow, thermal transport, and reactive biogeochemical transport." *Groundwater Reactive Transport Models*. Ed. by F. Zhang, G. Yeh, and J. Parker. United Arab Emirates: Bentham Publishers, 2012, pp. 3-41. DOI: 10.2174/978160805306311201010003.

[64] F. Zhang, G. Yeh, and J. Parker, eds. *Groundwater Reactive Transport Models*. United Arab Emirates: Bentham Publishers, 2012. DOI: 10.2174/97816080530631120101.