# Modeling and simulation of coupled biochemical and two-phase compositional flow in underground hydrogen storage

Elyes Ahmed, Brahim Amaziane, Salaheddine Chabab, Stéphanie Delage

Santacreu, Guillaume Galliéro, Olav Møyner, Xavier Raynaud

► To cite this version:

Elyes Ahmed, Brahim Amaziane, Salaheddine Chabab, Stéphanie Delage Santacreu, Guillaume Galliéro, et al.. Modeling and simulation of coupled biochemical and two-phase compositional flow in underground hydrogen storage. International Journal of Hydrogen Energy, 2025, 168, pp.150947. 10.1016/j.ijhydene.2025.150947. hal-05225645

HAL Id: hal-05225645

https://univ-pau.hal.science/hal-05225645v1

Submitted on 27 Aug 2025

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 Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

## Review Article

# Modeling and simulation of coupled biochemical and two-phase compositional flow in underground hydrogen storage

Elyes Ahmed<sup>a,●</sup>,*, Brahim Amaziane<sup>b,●</sup>, Salaheddine Chabab<sup>b,c,●</sup>, Stéphanie

Delage Santacreu<sup>b,●</sup>, Guillaume Galliéro<sup>b,●</sup>, Olav Møyner<sup>a,●</sup>, Xavier Raynaud<sup>a,●</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

### Keywords:

Søreide-Whitson model

Equation of state

Methanogenic reactions

Archaea

Underground hydrogen storage

Hydrogen solubility

Vapor-liquid equilibrium

Overall-composition formulation

MRST

Bio-clogging

## ABSTRACT

Integrating microbial activity into underground hydrogen storage models is crucial for simulating long-term reservoir behavior. In this work, we present a coupled framework that incorporates bio-geochemical reactions and compositional flow models within the Matlab Reservoir Simulation Toolbox (MRST). Microbial growth and decay are modeled using a double Monod formulation, with populations influenced by hydrogen and carbon dioxide availability. First, a refined Equation of State (EoS) is employed to accurately capture hydrogen dissolution, thereby improving phase behavior and modeling of microbial activity. The model is then discretized using a cell-centered finite-volume method with implicit Euler time discretization. A fully coupled fully implicit strategy is considered. Our implementation builds upon MRST's compositional module by incorporating the Søreide-Whitson EoS, microbial reaction kinetics, and specific effects such as bio-clogging and molecular diffusion. Through a series of 1D, 2D and 3D simulations, we analyze the effects of microbial-induced bio-geochemical transformations on underground hydrogen storage in porous media. These results highlight that accounting for bio-geochemical effects can substantially impact hydrogen loss, purity, and overall storage performance.

## Contents

<table><tr><td>1.</td><td>Introduction</td><td>2</td></tr><tr><td>2.</td><td>Models and methodology</td><td>3</td></tr><tr><td>2.1.</td><td>Two-phase multi-component transport with microbial activity model</td><td>3</td></tr><tr><td>2.1.1.</td><td>Methanogenic archaeal model</td><td>3</td></tr><tr><td>2.1.2.</td><td>Overall composition formulation</td><td>3</td></tr><tr><td>2.1.3.</td><td>Impact of microbial activity on rock properties</td><td>4</td></tr><tr><td>2.2.</td><td>Thermodynamic models</td><td>4</td></tr><tr><td>2.2.1.</td><td>EoS models adapted to capture gas solubility</td><td>4</td></tr><tr><td>2.2.2.</td><td>Hydrogen and carbon-dioxide solubility results: PR vs SW EoS</td><td>6</td></tr><tr><td>2.2.3.</td><td>Compositional validation</td><td>6</td></tr><tr><td>3.</td><td>Numerical simulations</td><td>6</td></tr><tr><td>3.1.</td><td>Test 1: 1D compositional simulation</td><td>6</td></tr><tr><td>3.1.1.</td><td>Simulation setup</td><td>6</td></tr><tr><td>3.1.2.</td><td>Results and discussion</td><td>7</td></tr><tr><td>3.2.</td><td>Test 2: 2D dome-shaped aquifer with different rock types</td><td>8</td></tr><tr><td>3.2.1.</td><td>Simulation setup</td><td>8</td></tr><tr><td>3.2.2.</td><td>Results and discussion</td><td>9</td></tr><tr><td>3.3.</td><td>Test 3: UHS in a 3D regular reservoir</td><td>11</td></tr></table>

<footnote><sup>a</sup> SINTEF Digital, Forskningsveien 1, 0373 Oslo, Norway</footnote>

<footnote><sup>b</sup> Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LMAP, UMR5142, Pau, France</footnote>

<footnote><sup>c</sup> Université de Pau et des Pays de l'Adour, E2S UPPA, LaTEP, Pau, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: Elyes.Ahmed@sintef.no (E. Ahmed), brahim.amaziane@univ-pau.fr (B. Amaziane), salaheddine.chabab@univ-pau.fr (S. Chabab), stephanie.delage@univ-pau.fr (S.D. Santacreu), guillaume.galliero@univ-pau.fr (G. Galliéro), Olav.Møyner@sintef.no (O. Møyner), Xavier.Raynaud@sintef.no (X. Raynaud).</footnote>
activity. In this work, we adopt a fully compositional approach, enabling a more accurate thermodynamic description and the inclusion of micro-organisms effects.

The rest of the paper is organized as follows. In Section 2.1, we present an overall compositional model for the two-phase multicomponent transport equations coupled to a methanogenic archaea model and estimate its effect on rock properties. In Section 2.2, we describe the EoS models that can capture gas solubility and, more particularly the Søreide-Whitson model adapted to MRST. The simulation results are compared with experiments when they are available. Then, in Section 3, we present a series of simulation cases designed to evaluate the model's ability to predict the solubility of H₂ and CO₂, capture microbial dynamics and their impact on multicomponent liquid-gas phase transport, and assess potential losses of H₂.

## 2. Models and methodology

### 2.1. Two-phase multi-component transport with microbial activity model

This section presents a two-phase overall-composition flow model incorporating microbial dynamics, particularly methanogenic archaea, for simulating subsurface interactions in underground hydrogen storage. Unlike the molar fraction, natural variable formulation in Hagemann [18], which uses mole conservation, our model is based on mass conservation and employs pressure and overall mole fractions for flow, along with normalized micro-organisms population for microbial dynamics, as primary variables.

#### 2.1.1. Methanogenic archaea model

The methanogenic archaea are anaerobic microorganisms present in water and able to grow under difficult thermodynamical conditions. They can live in underground reservoirs or aquifers where they play an essential ecological role by consuming substrates such as H₂ and CO₂ and by producing CH₄. This methanogenic reaction can occur at low temperatures and high pressures with optimum conditions around T = 308.15 K and P = 90 bar [4,9,18] and can be expressed as:

$$4\text{H}_2 + \text{CO}_2 \rightarrow \text{CH}_4 + 2\text{H}_2\text{O} \quad (1)$$

The microbial dynamics are governed by a growth rate function, $\psi^{growth}$, and a decay rate function, $\psi^{decay}$, leading to the following conservation equation [4,18,19]:

$$\frac{\partial(\Phi\rho_l S_l \mathcal{N})}{\partial t} + \nabla \cdot (\Phi\rho_l S_l \mathbf{D}_{\mathcal{N},l} \nabla \mathcal{N}) = (\psi^{growth} - \psi^{decay})\Phi\rho_l S_l \mathcal{N}, \quad (2)$$

where $\Phi$ is the medium porosity, $D_{\mathcal{N},l}$ the microbial diffusion coefficient (m² s⁻¹), $\rho_l$ the liquid phase density (kg/m³), and $S_l$ the liquid saturation. The microbial population is normalized as $\mathcal{N} = n/n_0$, where n (respectively $n_0$) denotes the current (respectively initial) microbial density (m⁻³). A double Monod model [20] is employed for the growth rate function meanwhile the decay rate is formulated as a linear rate [10]:

$$\psi^{growth} = \psi_{\max}^{growth} \left( \frac{X_{l,H_2}}{l_{H_2} + X_{l,H_2}} \right) \left( \frac{X_{l,CO_2}}{l_{CO_2} + X_{l,CO_2}} \right), \quad (3)$$

$$\psi^{decay} = b_{decay} \mathcal{N} = b_{decay} \frac{n}{n_0}, \quad (4)$$

Here, $X_{l,H_2}$ and $X_{l,CO_2}$ denote the molar fractions of H₂ and CO₂, respectively, in the liquid phase; $l_{H_2}$ and $l_{CO_2}$ are the corresponding half-saturation constants; $\psi_{\max}^{growth}$ is the maximum growth rate (s⁻¹) and $b_{decay}$ is the decay rate (s⁻¹). Several experimental values obtained essentially in laboratory are available in the literature (see Table 1). Note that our results and implementations can be easily extended to other reactions involving H₂ consumption, such as sulfide reduction [13]. Several experiments showed that the dimensionless form of the Monod model mitigates numerical instability caused by extreme microbial density values, ensuring better conditioning and

scaling with the overall-composition mass-fraction formulation [21]. Although methanogenic archaea can displace in the liquid phase due to chemotaxis or thanks to dispersion or advection phenomena, [8, 22,23], the transport of the microbial population is not studied in the present work. Furthermore, we only focus on the microbial activity of the methanogenic archaea in the present paper.

#### 2.1.2. Overall composition formulation

The microbial model is coupled to the two-phase compositional model for fluid flow through a bioreactive source term in the mass conservation equations. This model consists of mass conservation equations for each component based on the overall composition formulation [34] and the multiphase extension of Darcy's law, which governs fluid phase velocity:

$$\partial_t(\Phi\rho_l z_l) + \nabla \cdot \mathbf{u}_l + q_l = \Phi\gamma_l^{H_2} \frac{\psi^{growth} \mathcal{N} S_l}{Y_{H_2}}, \quad i \in \{\text{H}_2, \text{H}_2\text{O}, \dots\}. \quad (5)$$

Here, $\rho_l$ denotes the total mass density, $z_l$ the overall mole fraction of the component i, $\mathbf{u}_l$ the flux, and $q_l$ represents external sources, such as well sources. Note that the pressure and the overall molar composition are chosen as the primary variables. The microbial contribution is modeled through the source term $\Phi\gamma_l^{H_2} \frac{\psi^{growth} \mathcal{N} S_l}{Y_{H_2}}$ [4,19], where $Y_{H_2}$ is the yield coefficient and $\gamma_l^{H_2} = n_0\gamma_l M_l/\gamma_{H_2}$ accounts for the stoichiometric conversion of hydrogen into other species i, with the stoichiometric coefficients in reaction (1) given by $[\gamma_{H_2}, \gamma_{CO_2}, \gamma_{CH_4}, \gamma_{H_2\text{O}}] = [-4, -1, 1, 2]$. The total density $\rho_l$ and the overall molar fraction $z_l$ are defined as:

$$\rho_l = \sum_{\zeta \in \{l,g\}} \rho_{\zeta} S_{\zeta}, \quad (6a)$$

$$z_l = \sum_{\zeta \in \{l,g\}} \frac{\rho_{\zeta} S_{\zeta}}{\rho_l} X_{\zeta,l} = \sum_{\zeta \in \{l,g\}} A_{\zeta} X_{\zeta,l}. \quad (6b)$$

where $A_{\zeta} = \frac{\rho_{\zeta} S_{\zeta}}{\rho_l}$ represents the contribution of each phase ($\zeta \in \{l, g\}$) to the total mass density $\rho_l$, and $\mathbf{u}_l$ is the species flux made up of two contributions: the standard multiphase extension of Darcy's law, $v_{\zeta}$, and the diffusive flux, $j_{\zeta,l}$, of component i in phase $\zeta$. These are given by:

$$\mathbf{u}_l = \sum_{\zeta \in \{l,g\}} \rho_{\zeta} X_{\zeta,l} v_{\zeta} - j_{\zeta,l}, \quad (7a)$$

$$v_{\zeta} = -\frac{k_{r\zeta}}{\mu_{\zeta}} \mathbf{K}(\nabla \rho_{\zeta} - \rho_{\zeta} \mathbf{g}), \quad \zeta \in \{l, g\}, \quad (7b)$$

$$j_{\zeta,l} = -\Phi \rho_{\zeta} S_{\zeta} \mathbf{D}_{\zeta,l} \nabla X_{\zeta,l}, \quad i \in \{\text{H}_2, \text{H}_2\text{O}, \dots\}. \quad (7c)$$

Here, $v_{\zeta}$ is the phase velocity from the generalized Darcy's law, where $k_{r\zeta}$ is the relative permeability, $\mu_{\zeta}$ is the viscosity, $\mathbf{K}$ is the absolute permeability tensor, and $\rho_{\zeta}$ is the phase pressure. The term $j_{\zeta,l}$ represents molecular diffusion flux, driven by the gradient of the mole fraction $X_{\zeta,l}$, with $\mathbf{D}_{\zeta,l}$ being the diffusion coefficient. Besides $\mathbf{D}_{\zeta,l}$ depends on the tortuosity $\tau_{\zeta}$ and binary molecular diffusion coefficients $\mathbf{D}_{\zeta,l}^{\psi}$ (see Table 2) such as $\mathbf{D}_{\zeta,l} = \mathbf{D}_{\zeta,l}^{\psi} \tau_{\zeta}$. The Millington and Quirk model is used for the tortuosity [35] and writes as:

$$\tau_{\zeta} = \frac{1}{\Phi^2} (\Phi S_{\zeta})^{7/3}. \quad (8)$$

To fully determine the unknowns in the coupled system (2)–(8), we must compute the phase saturations and the mole fractions of all components in each phase. Moreover, since both phase densities and flow potentials may exhibit pressure dependence, knowledge of the individual phase pressures is required. Crucially, evaluating phase behavior and flow properties demands an accurate description of component distribution between the fluid phases. This is achieved by solving the following set of nonlinear closure relations:

$$\sum_{\zeta \in \{l,g\}} S_{\zeta} = \sum_{\zeta \in \{l,g\}} A_{\zeta} = 1, \quad (9a)$$
<table><caption>**Table 1**
Range of parameters for simulating compositional flow with bacteria.</caption><thead><tr><th>Parameter</th><th>Range</th><th>Reference</th></tr></thead><tbody><tr><td>ψ<sub>max</sub><sup>growth</sup> (s<sup>-1</sup>)</td><td>1.805 × 10<sup>-6</sup>–3.009 × 10<sup>-5</sup></td><td>[24–26]</td></tr><tr><td rowspan="3">l<sub>H₂</sub> (mol/mol)</td><td rowspan="3">9 × 10<sup>-12</sup>–3.240 × 10<sup>-7</sup></td><td>[27,28]</td></tr><tr><td>[29–31]</td></tr><tr><td>[24–26]</td></tr><tr><td rowspan="5">l<sub>CO₂</sub> (mol/mol)<br>Y<sub>H₂</sub> (mol<sup>-1</sup>)</td><td>2.34 × 10<sup>-8</sup>–5.4 × 10<sup>-6</sup></td><td>[27,28]</td></tr><tr><td>7.697 × 10<sup>10</sup>–1.093 × 10<sup>13</sup></td><td>[29,30]</td></tr><tr><td>[24–26]</td></tr><tr><td>[27,29,31]</td></tr><tr><td>6.944 × 10<sup>-7</sup>–1.019 × 10<sup>-6</sup></td><td>[26,29]</td></tr><tr><td>b<sub>decay</sub> (s<sup>-1</sup>)<br>T (K)</td><td>298.15–313.15</td><td>[19,31,33]</td></tr><tr><td>P (bar)</td><td>80–100</td><td>[4,19]</td></tr></tbody></table>

<table><caption>**Table 2**
Binary molecular diffusion coefficients of components H₂, C₁, H₂O, CO₂, N₂, C₂, C₃, and NC₄ in both liquid and gas phases [41,42].</caption><thead><tr><th></th><th>H₂</th><th>C₁</th><th>H₂O</th><th>CO₂</th><th>N₂</th><th>C₂</th><th>C₃</th><th>NC₄</th></tr></thead><tbody><tr><td>D<sub>l,i</sub><sup>+</sup>(m² s<sup>-1</sup>)</td><td>4.5 × 10<sup>-9</sup></td><td>2.6 × 10<sup>-9</sup></td><td>2.3 × 10<sup>-9</sup></td><td>1.9 × 10<sup>-9</sup></td><td>2.1 × 10<sup>-9</sup></td><td>1.2 × 10<sup>-9</sup></td><td>0.9 × 10<sup>-9</sup></td><td>0.7 × 10<sup>-9</sup></td></tr><tr><td>D<sub>g,i</sub><sup>+</sup>(m² s<sup>-1</sup>)</td><td>6.1 × 10<sup>-5</sup></td><td>2.5 × 10<sup>-5</sup></td><td>1.5 × 10<sup>-5</sup></td><td>1.4 × 10<sup>-5</sup></td><td>1.8 × 10<sup>-5</sup></td><td>1.6 × 10<sup>-5</sup></td><td>1.2 × 10<sup>-5</sup></td><td>0.9 × 10<sup>-5</sup></td></tr></tbody></table>

$$ \sum_{l}^{N_c} z_l = \sum_{l}^{N_c} X_{\zeta,l} = 1, \quad \forall \zeta \in \{l, g\}, \qquad (9b) $$

$$ p_l - p_g - p_{c,l,g} = 0, \quad f_{l,i} - f_{g,i} = 0, \quad i \in \{\text{H}_2, \text{H}_2\text{O}, \dots\}, \qquad (9c) $$

Given that the equation of state is properly calibrated against experimental data, the compositional model can accurately capture phase behavior and interphase mass transfer processes. This is particularly important in the presence of microbial activity, where chemical transformations such as hydrogen consumption, carbon dioxide production, and biomass accumulation can significantly alter phase compositions, fluid properties, and ultimately reservoir performance. Note that these equations are coupled with well equations and specific initial and boundary conditions. For more details on the model, see [36]. The compositional fluid flow model was solved using the MRST object-oriented automatic differentiation framework [37,38]. By default, MRST employs a fully implicit two-point flow approximation (TPFA) with a single-point upwinding to solve the flow equations. Stability and flash calculations are performed at each Newton iteration for the compositional model. Note that [39] presented coupled bio-geochemical reactions through iterative integration of PHREEQC [40] with MRST's compositional module. Here, a versatile all-in-MRST framework with fully-implicit coupling is implemented with focus on accurate EoS modeling, bio-clogging effects, and flexible bacterial dynamics.

### 2.1.3. Impact of microbial activity on rock properties

When methanogenic archaea grow, they obstruct the pore space affecting the fluid flow through the porous medium [9]. Hence, the porosity and the permeability are impacted. The porosity $\Phi$ and permeability $K$ vary with bacteria density $N$ as:

$$ \Phi(\mathcal{N}) = \frac{\Phi_0}{1 + c_p(\mathcal{N}/\mathcal{N}_c)^2}, \qquad (10) $$

$$ \mathbf{K}(\mathcal{N}) = \mathbf{K}_0 \left( \frac{1 - \Phi_0}{1 - \Phi} \right)^2 \left( \frac{\Phi}{\Phi_0} \right)^3, \qquad (11) $$

where $N_c$ is the critical density for maximal clogging, $\Phi_0$ and $\mathbf{K}_0$ are reference values, and $c_p$ controls porosity reduction [9].

## 2.2. Thermodynamic models

In this section, we aim at employing a thermodynamic model able to take into account the solubility of gases such as hydrogen and carbon dioxide in two-phase and multi-component problems. An accurate estimation of the solubility of gases reveals to be essential to predict the Vapor-Liquid Equilibrium (VLE) and the phase compositions of the multi-component system.

### 2.2.1. EoS models adapted to capture gas solubility

In an isothermal and isobaric system, the coexistence of several phases results in the equality of chemical potentials of each component in all phases with respect to the Gibbs-Duhem equation [43]. For VLE, the equality of chemical potentials is expressed as the equality of fugacities [44].

In the literature, two approaches are available to calculate VLE: the symmetric and asymmetric approach. In the symmetric method, the fugacity coefficients of each component on both phases are calculated from a cubic (e.g. Peng-Robinson (PR), Soave-Redlich-Kwong (SRK), etc.) or a non-cubic (e.g. SAFT, CPA, multiparameter Helmholtz energy, ..., -type) equation of state (EoS). In the asymmetric approach, the fugacity coefficients of each component in the vapor phase are calculated as in the symmetric approach, while those of the liquid phase are determined using the model G-excess [45].

Furthermore, in their work on the VLE of H₂-H₂O mixtures, Rahbari et al. [46] have proved that using cubic EoS, such as PR or SRK, with classical mixing rules are unable to accurately predict the coexistence composition of the non aqueous and aqueous phases, even if adapted Binary Interactions Parameters (BIPs) were introduced.

However, Søreide and Whitson [47] have elaborated a model specifically tailored for brine solutions which consists in adding salt molality and temperature dependencies in the alpha-function of the PR-EoS and the BIPs between water and gases in the aqueous phase. Then, the model has been improved by Chabab et al. for CO₂ and H₂ solubilities in water and NaCl brine [48,49].

To fully describe the thermodynamic behavior of the multicomponent fluid, our compositional model implemented in MRST must be complemented with an appropriate equation of state. The fluid is assumed to consist of brine, treated as a pseudo-component representing H₂O and NaCl, along with n-1 additional light components, such as CO₂, N₂, and CH₄. In this work, NaCl is chosen as the representative salt.

To model the phase equilibrium of this multi-component system, we have implemented the enhanced version of the SW EoS found in the literature [3,49]. The PR EoS is expressed as [3]:

$$ P(T, V) = \frac{RT}{V-b} - \frac{a}{V(V+b) + b(V-b)}, \qquad (12) $$

where R is the universal gas constant, T is the temperature and V is the molar volume, a and b are the EoS parameters.

In a mixture, a and b are the EoS parameters of the mixture and depend on the temperature, the molar fraction and the phase. They
<table><caption>**Table 3**
Coefficients for generalized BIPs equations in liquid/gas phases [49–51].</caption><thead><tr><th colspan="5">Liquid-phase</th></tr><tr><th></th><th>CO<sub>2</sub></th><th>H<sub>2</sub></th><th>C<sub>1</sub></th><th>H<sub>2</sub>S</th></tr></thead><tbody><tr><td>a<sub>1</sub></td><td>-0.43575155</td><td>-2.11917</td><td>-1.625685</td><td>-0.42619</td></tr><tr><td>a<sub>2</sub></td><td>-0.05766906744</td><td>0.14888</td><td>1.114873</td><td>0.673586</td></tr><tr><td>a<sub>3</sub></td><td>0.00826464849</td><td>-</td><td>-</td><td>-</td></tr><tr><td>a<sub>4</sub></td><td>-</td><td>-13.01835</td><td>-</td><td>-</td></tr><tr><td>b<sub>1</sub></td><td>-0.00129539193</td><td>-0.0226322</td><td>8.590105 × 10<sup>-21</sup></td><td>-0.0575</td></tr><tr><td>b<sub>2</sub></td><td>-0.0016698848</td><td>-0.0044736</td><td>0.001812763</td><td>-0.078823343</td></tr><tr><td>b<sub>3</sub></td><td>-</td><td>-</td><td>-0.169968</td><td>-0.21625</td></tr><tr><td>b<sub>4</sub></td><td>-</td><td>-0.43946</td><td>-0.04198569</td><td>-0.160085087</td></tr><tr><td>c</td><td>-0.47866096</td><td>1</td><td>1</td><td>1</td></tr></tbody><tbody><tr><th colspan="5">Liquid-phase (Cont.)</th></tr><tr><th></th><th>N<sub>2</sub></th><th>C<sub>2</sub></th><th>C<sub>3</sub></th><th>NC<sub>4</sub></th></tr><tr><td>a<sub>1</sub></td><td>-1.702359096</td><td>1.1120 - 1.7369ω<sup>0.1</sup></td><td>1.1120 - 1.7369ω<sup>0.1</sup></td><td>1.1120 - 1.7369ω<sup>0.1</sup></td></tr><tr><td>a<sub>2</sub></td><td>0.450487</td><td>1.1001 + 0.8360ω</td><td>1.1001 + 0.8360ω</td><td>1.1001 + 0.8360ω</td></tr><tr><td>a<sub>3</sub></td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>a<sub>4</sub></td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>b<sub>1</sub></td><td>0.01792130</td><td>0.017407</td><td>0.017407</td><td>0.017407</td></tr><tr><td>b<sub>2</sub></td><td>0.066426</td><td>0.033516</td><td>0.033516</td><td>0.033516</td></tr><tr><td>b<sub>3</sub></td><td>-</td><td>-0.15742 - 1.0988ω</td><td>-0.15742 - 1.0988ω</td><td>-0.15742 - 1.0988ω</td></tr><tr><td>b<sub>4</sub></td><td>-</td><td>0.011478</td><td>0.011478</td><td>0.011478</td></tr><tr><td>c</td><td>0.8</td><td>1</td><td>1</td><td>1</td></tr></tbody><tbody><tr><th colspan="5">Gas-phase</th></tr><tr><th></th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th></th><th></th></tr><tr><td>CO<sub>2</sub></td><td>0.207440935</td><td>-2.066623464504 × 10<sup>-2</sup></td><td></td><td></td></tr><tr><td>H<sub>2</sub></td><td>0.042834</td><td>0.01993</td><td></td><td></td></tr><tr><td>C<sub>1</sub></td><td>-</td><td>0.494435</td><td></td><td></td></tr><tr><td>H<sub>2</sub>S</td><td>-0.05965</td><td>0.19031</td><td></td><td></td></tr><tr><td>N<sub>2</sub></td><td>-</td><td>0.385438</td><td></td><td></td></tr><tr><td>C<sub>2</sub></td><td>-</td><td>0.4920</td><td></td><td></td></tr><tr><td>C<sub>3</sub></td><td>-</td><td>0.5525</td><td></td><td></td></tr><tr><td>NC<sub>4</sub></td><td>-</td><td>0.5091</td><td></td><td></td></tr></tbody></table>

are noted $a_ζ(T, X)$ and $b(X)$ and are determined using standard mixing rules, namely:

$$a_{\zeta}(T, X) = \sum_{i,j=1}^{N_c} \sqrt{a_i(T)a_j(T)}(1 - k_{\zeta,ij})X_{\zeta,i}X_{\zeta,j}, \quad b(X) = \sum_{i=1}^{N_c} b_i X_{\zeta,i}, \quad (13)$$

where

$$a_i(T) = \Omega_a \frac{R^2 T_{c,i}}{P_{c,i}} \alpha_i(T), \quad b_i = \Omega_b \frac{RT_{c,i}}{P_{c,i}}, \quad (14)$$

with $\Omega_a = 0.45724$, $\Omega_b = 0.0778$ constants. $P_{c,i}$ and $T_{c,i}$ represent the critical pressure and temperature for component $i$, respectively.

In both gaseous and liquid phases, we use the generalized alpha function, $\alpha_i$, proposed by Soave [52] and formulated as:

$$\sqrt{\alpha_i} = 1 + (0.37464 + 1.54226\omega_i - 0.26992\omega_i^2)(1 - \sqrt{T_{r,j}}). \quad (15)$$

Here $T_{r,j} = T/T_{c,i}$ is reduced temperature corresponding to the critical temperature $T_{c,i}$ and $\omega_i$ is the acentric factor for component $i$.

The first modification applied in the SW model [47] consists in including a salt molality, $m_{salt}$ (mol kgw$^{-1}$), dependence for the component H$_2$O in the liquid phase, namely:

$$\sqrt{\alpha_{H_2O}} = 1 + 0.453(1 - T_{r,H_2O}(1 - 0.0103m_{salt}^{1.1})) + 0.0034(T_{r,H_2O}^{-3} - 1). \quad (16)$$

The second modification in the SW model consists in adding particular (phase-specific) BIPs, $k_{\zeta,ij}$, in the mixing rules in order to represent the cross interactions between the $i$th and $j$th components and more particularly between component H$_2$O and the others in each phase. The BIPs are symmetric, that is $k_{\zeta,ij} = k_{\zeta,ji}$. Given that the thermodynamical conditions in UHS reveal to be far from the H$_2$O critical point, the use of phase-specific BIPs remains suitable [49].

We compile a comprehensive dataset from the literature for BIPs across a wide range of mixtures. For isothermal compositional models,

this approach adds no computational cost to flash calculations since binary interaction parameters are computed once. However, it may significantly impact thermal simulation performance. Additionally, BIP calibration may be required if experimental data are available in the literature. As an illustrative example, the BIPs formulations can be unified in the following generalized forms:

$$\begin{align}
k_{r,CO_2-H_2O} &= T_{r,CO_2}(a_1 + a_2T_r + a_3T_{r,CO_2}m_{salt}) + m_{salt}^2(b_1 + b_2T_{r,CO_2}) + b_3, \nonumber \\
k_{r,j-H_2O} &= a_1(1 + b_1m_{salt}^c) + a_2(1 + b_2m_{salt}^c)T_{r,j} + a_3(1 + b_3m_{salt}^c)T_{r,j}^2 + a_4 \exp(b_4T_{r,j}), \nonumber \\
k_{g,j-H_2O} &= a_1T_{r,j} + a_2. \tag{17}
\end{align}$$

As an example, for the mixture [H$_2$, H$_2$S, N$_2$, C$_1$, C$_2$, C$_3$, NC$_4$], the corresponding coefficients are provided in Table 3.

Like PR EoS, the input parameters for each fluid component are critical properties ($P_c$, $T_c$), acentric factor ($\omega$), and molar mass ($M$). These properties are essential for accurately modeling the phase behavior of the multicomponent system. The critical properties and acentric factors for the components considered in this study come from the Coolprop library and are listed in Table 4.

Concerning the implementation of SW EoS in MRST, for reasons of simplicity, salinity in SW EoS has been considered as a constant input parameter and not as a composition that can vary (e.g. through water evaporation) as part of the stability and flash calculations. For a two-phase simulation and at conditions far from the critical point of the gas/water mixture, this assumption is fairly acceptable [53]. It is worth noting that recently, methods have been developed which consider salinity in SW EoS as a composition (and not a parameter) and include it in equilibrium calculations [53,54]. These modifications, which make the SW model thermodynamically consistent, are especially necessary for applications with more than two phases, at very high temperatures (close to the critical point of water), and when trapping by dissolution
<table><caption>**Table 4**
Components properties from the Coolprop library.</caption><thead><tr><th>Component</th><th>P<sub>c</sub>(Pa)</th><th>T<sub>c</sub>(K)</th><th>V<sub>c</sub>(m<sup>3</sup>)</th><th>ω</th><th>M(kgmol<sup>-1</sup>)</th></tr></thead><tbody><tr><td>H<sub>2</sub>O</td><td>2.21 × 10<sup>7</sup></td><td>647.1</td><td>5.595 × 10<sup>-5</sup></td><td>0.344</td><td>0.0180153</td></tr><tr><td>H<sub>2</sub></td><td>1.30 × 10<sup>6</sup></td><td>33.1</td><td>6.448 × 10<sup>-5</sup></td><td>-0.219</td><td>0.0020159</td></tr><tr><td>CO<sub>2</sub></td><td>7.38 × 10<sup>6</sup></td><td>304.1</td><td>9.412 × 10<sup>-5</sup></td><td>0.224</td><td>0.0440098</td></tr><tr><td>CH<sub>4</sub>/C<sub>1</sub></td><td>4.60 × 10<sup>6</sup></td><td>190.6</td><td>9.863 × 10<sup>-5</sup></td><td>0.011</td><td>0.0160428</td></tr><tr><td>C<sub>2</sub></td><td>4.87 × 10<sup>6</sup></td><td>305.3</td><td>1.458 × 10<sup>-4</sup></td><td>0.099</td><td>0.0300690</td></tr><tr><td>C<sub>3</sub></td><td>4.26 × 10<sup>6</sup></td><td>369.8</td><td>9.418 × 10<sup>-5</sup></td><td>0.152</td><td>0.0440962</td></tr><tr><td>NC<sub>4</sub></td><td>3.80 × 10<sup>6</sup></td><td>425.1</td><td>2.549 × 10<sup>-4</sup></td><td>0.201</td><td>0.0581222</td></tr><tr><td>N<sub>2</sub></td><td>3.40 × 10<sup>6</sup></td><td>126.2</td><td>8.941 × 10<sup>-5</sup></td><td>0.037</td><td>0.0280135</td></tr></tbody></table>

<table><caption>**Table 5**
Calculated errors between the EoS models and the experiment data.</caption><thead><tr><th>EoS model</th><th>Component solubility</th><th>Salinity (mol kgw<sup>-1</sup>)</th><th>error<sub>min</sub></th><th>error<sub>mean</sub></th><th>error<sub>max</sub></th></tr></thead><tbody><tr><td>PR</td><td>H<sub>2</sub></td><td>0</td><td>0.7710</td><td>0.9054</td><td>0.9805</td></tr><tr><td>SW</td><td>H<sub>2</sub></td><td>0</td><td>0.0002</td><td>0.0138</td><td>0.0284</td></tr><tr><td>SW</td><td>H<sub>2</sub></td><td>1</td><td>0.0013</td><td>0.0227</td><td>0.0622</td></tr><tr><td>SW</td><td>H<sub>2</sub></td><td>4</td><td>0.0046</td><td>0.0301</td><td>0.0678</td></tr><tr><td>PR</td><td>CO<sub>2</sub></td><td>0</td><td>0.5209</td><td>0.6401</td><td>0.7385</td></tr><tr><td>SW</td><td>CO<sub>2</sub></td><td>0</td><td>0.0208</td><td>0.0566</td><td>0.1074</td></tr><tr><td>SW</td><td>CO<sub>2</sub></td><td>1</td><td>0.0166</td><td>0.0239</td><td>0.0306</td></tr><tr><td>SW</td><td>CO<sub>2</sub></td><td>3</td><td>0.0048</td><td>0.0463</td><td>0.0969</td></tr></tbody></table>

is the dominant mechanism (e.g. CO₂ sequestration). Note that for accurate density calculations, volume translation techniques can be employed to modify the EoS model's volume output. While this approach improves mixture density predictions, we have not implemented it in the current work.

### 2.2.2. *Hydrogen and carbon-dioxide solubility results: PR vs SW EoS*

In MRST the VLE problem is solved using a symmetric approach with a PR or SRK Equation of State. Then, the obtained Flash equations are solved as a nested non linear system [38]. As previously mentioned, in this work, the Søreide-Whitson model has been implemented in MRST to account for the solubility of H₂ and CO₂ in water and brine.

Then, several numerical simulations have been performed on two-component H₂ – H₂O and CO₂ – H₂O and two-phase systems to verify the efficiency of the SW model developed in MRST over a wide range of thermodynamical conditions, i.e.: 100 bar ≤ P ≤ 200.11 bar and 298.15 K ≤ T ≤ 373.85 K for the H₂ – H₂O two-phase system and 17 bar ≤ P ≤ 141.1 bar and 323.2 K ≤ T ≤ 373.41 K for the CO₂ – H₂O two-phase system. Then, they are compared with the experimental data provided by [3,48].

As demonstrated in Rahbari et al. [46], Fig. 1 (left) shows the limits of the PR model to capture the H₂ solubility in water whereas the SW model matches very well with the experimental data in water and brine, see Fig. 1 (right). Relative errors are computed pointwise as |X<sub>i,l</sub> - X<sub>i,l</sub><sup>exp</sup>|/|X<sub>i,l</sub><sup>exp</sup>|, where i indexes the component i and l denotes the liquid phase. We report the minimum, mean, and maximum values over all points, see Table 5. Therein, the PR model exhibits a mean error of 90.5% (min error 77.1%) with the experiment data whereas the SW model does not exceed an error of 2.8%. Furthermore, it has to be noted that the solubility of H₂ decreases when the salinity increases. Similar results have been obtained for solubility of CO₂ in water and brine, see Fig. 2.

### 2.2.3. *Compositional validation*

The compositional MRST model was validated against commercial solvers using the Peng-Robinson EoS. The test case is a 1D problem from [34], later verified in [38], injects pure CO₂ (500 steps) into a 1000-cell reservoir initially containing a CO₂-CH₄-C₁₀H₂₂ mixture (liquid-dominated, with supercritical CO₂ at 423.15 K, 75 bar). A similar validation using the SW-EoS (not shown) confirmed agreement with MRST-PR and other commercial simulators.

## 3. Numerical simulations

In this work, we implemented the SW equation of state (see Section 2.2) in MRST to accurately capture gas dissolution effects. This thermodynamic model allows for a consistent coupling with microbial dynamics, enabling the analysis of their impact on hydrogen loss, transport behavior, and bio-clogging mechanisms. Although the biochemical model primarily uses literature-derived (or averaged) parameters, which provide a valid scientific foundation for simulations, calibration with site-specific laboratory observations could further enhance prediction accuracy.

These developments are built on the compositional module of MRST [16], which provides advanced simulation capabilities for multicomponent systems. The module supports integration with external, pre-compiled linear solvers such as AMGCL [16,36]. Through the use of memory-efficient backends for linearization and assembly, combined with high-performance solvers written in compiled languages, MRST can efficiently simulate both synthetic and field-scale reservoir models involving millions of unknowns.

The numerical results are structured as follows. In Section 3.1, we demonstrate the effects of microbial activity on H₂ loss and methanation using a simple 1D model without wells, focusing on reservoir pressure, saturation, and H₂ dissolution. In Section 3.2, we apply our fluid model to a 2D dome-shaped model with different rock types and heterogeneous entry pressures. This section showcases the capability of our module to handle a highly nonlinear model due to bio-clogging effects, varying capillary pressure, and different residual saturations, as well as the impact of bacteria on H₂ impurity. Section 3.3 demonstrates a 3D test case, highlighting the interaction between salinity effects, the molecular diffusion and microbial activity, and their influence on H₂ loss and production efficiency. Finally, in Section 3.4, we test our model on a large UHS 3D test case with a complex mixture and realistic injection scenarios.

### 3.1. Test 1: 1D compositional simulation

#### 3.1.1. Simulation setup

This is a simple 1D reservoir domain, inspired by Møyner and Tchelepi [55], consisting of 1000 cells with a total length of 100 m. This setup represents a simplified aquifer system where H₂ and CO₂ are present to study microbial activity and its impact on hydrogen loss. No wells are included in the simulation. We considered two initial compositions:
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** H₂ solubility as a function of pressure at different temperatures. Solubility in pure water is modeled using PR and SW equations of state (left), while solubility in saline water is shown using the SW model (right).</figcaption></figure>

Source: Data collected from [49].

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** CO₂ solubility along with pressure for different temperatures. CO₂ solubility in pure captured with the PR and SW EoS models (left). CO₂ solubility in salt water captured with the SW model (right).</figcaption></figure>

Source: Data collected from [3].

<table><caption>**Table 6**

Microbial model parameters used in the MRST simulations.</caption><thead><tr><th>Parameter</th><th>ψ<sub>max</sub><sup>growth</sup> (s<sup>-1</sup>)</th><th>b<sub>decay</sub> (s<sup>-1</sup>)</th><th>Y<sub>H₂</sub> (mol<sup>-1</sup>)</th><th>l<sub>H₂</sub> (-)</th><th>l<sub>CO₂</sub> (-)</th><th>n<sub>0</sub> (m<sup>-3</sup>)</th></tr></thead><tbody><tr><td>Value</td><td>1.338 × 10<sup>-4</sup></td><td>1.35 × 10<sup>-6</sup></td><td>3.9 × 10<sup>11</sup></td><td>3.6 × 10<sup>-7</sup></td><td>1.98 × 10<sup>-6</sup></td><td>1.0 × 10<sup>9</sup></td></tr></tbody></table>

* **Gas-rich case:** 40% water, 44.5% H₂, 5.5% CH₄, and 10% CO₂.

* **Liquid-rich case:** 95% water, 4% H₂, 0.5% CH₄, and 0.5% CO₂.

We analyze both low- and high-H₂ concentrations to address key phases of UHS development. The liquid-rich case (low H₂) represents initial field projects (e.g., [2]), while the gas-rich case (high H₂) anticipates future pure-H₂ systems, enabling evaluation of both early feasibility and long-term microbial-dominant behavior.

The temperature is set at 313.75 K, and the initial pressure is 82 bar. The reservoir has a uniform initial permeability of 10 mD and an initial porosity of 0.25, except in the final cells, where it is set to 0.9 to emulate source/sink behavior. Quadratic relative permeability functions are used for water and gas, with residual saturations $S_{lr} = 0.2$ and $S_{gr} = 0.1$. Effective permeabilities are obtained by shifting saturations and reinterpolating over the intervals [$S_{lr}$, 0.9] for water and [$S_{gr}$, 0.8] for gas. It should be noted that molecular and microbial diffusion are neglected. For the porosity reduction model, we set $c_p = 1$

and $N_c = n_c/n_0 = 180$. The densities of liquid and gas are set to their values under the surface conditions ($P = 1.013$ bar and $T = 15.56$ °C), that is, $\rho_l = 999.014$ kg/m³ and $\rho_g = 0.084$ kg/m³.

The first simulation is abiotic; the second includes microbial effects but no clogging; and the third includes microbial clogging effects. All simulations are performed with the SW EoS. The parameters for the microbial model are given in Table 6.

### 3.1.2. Results and discussion

In Fig. 3 (top), we plot the dissolved H₂ concentration in the liquid phase. The water-rich case shows values close to literature reports (1.22%), while the gas-rich initialization yields significantly lower dissolution (0.072%). Although bioclogging reduces dissolution in both cases, it has minimal impact on the total H₂ consumption at the simulation endpoint. The primary effects of bioclogging appear in the pressure and saturation behavior (Fig. 3, bottom) rather than in H₂ consumption, as expected for a well-free system. Mean gas
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** 1D case: Evolution of dissolved H₂ in the aqueous phase (top) and comparison of mean pressure and saturation (bottom) for three scenarios: (i) abiotic, (ii) with microbial activity, and (iii) with microbial activity and clogging effects.</figcaption></figure>

saturation decreases with H₂ and CO₂ consumption but increases under bioclogging due to gas trapping caused by reduced pore space and restricted fluid flow. Pressure declines are significant: the gas-rich case shows a 5% decrease without clogging versus 17% with clogging, while the water-rich case shows 50% and 62% decreases, respectively.

In Fig. 4 (top), we plot the changes in the mixture composition due to microbial activity. Microbial effects begin immediately, as H₂ and CO₂ are initially present in both gas- and water-rich cases. For the gas-rich case, H₂ consumption reaches 6.2% after 5 years, whereas the water-rich case achieves complete H₂ consumption within 3.5 years. Fig. 4 (bottom) reveals the evolution of mole fractions: CH₄ production increases as H₂ and CO₂ decrease, with more pronounced effects in the high-liquid-content case.

In the last experiment, we compare the SW EoS with the PR EoS. The results for the gas-rich case, shown in Fig. 5, indicate a slower consumption rates of H₂ and CO₂ when using the PR EoS, resulting in lower CH₄ production. This behavior aligns with theoretical expectations, as the PR EoS is known to less accurately capture dissolution effects, particularly the thermodynamic coupling between the gas and aqueous phases, which consequently reduces the extent of gas dissolution in the aqueous phase. These findings emphasize the importance of EoS selection for accurate modeling of microbial conversion processes in subsurface systems.

## 3.2. Test 2: 2D dome-shaped aquifer with different rock types

### 3.2.1. Simulation setup

This example, adapted from Ahmed et al. [17], investigates hydrogen storage in a saline aquifer using a compositional model instead of

<figure><img src="image_7.png" /><figcaption>**Table 7**
2D aquifer: Formation and petrophysical properties [17].</figcaption></figure>

<table><thead><tr><th>Symbol</th><th>Description</th><th>Caprock</th><th>Storage zone</th><th>Bedrock</th></tr></thead><tbody><tr><td>K (mD)</td><td>Permeability</td><td>1.0 × 10<sup>-4</sup></td><td>10</td><td>1.0 × 10<sup>-2</sup></td></tr><tr><td>Φ<sub>ref</sub> (-)</td><td>Porosity</td><td>0.25</td><td>0.1</td><td>0.1</td></tr><tr><td>S<sub>lg</sub> (-)</td><td>Gas residual saturation</td><td>0.1</td><td>0.05</td><td>0.1</td></tr><tr><td>S<sub>lg</sub> (-)</td><td>Water residual saturation</td><td>0.1</td><td>0.1</td><td>0.1</td></tr><tr><td>P<sub>e</sub> (bar)</td><td>Entry pressure</td><td>10</td><td>0.4</td><td>5</td></tr><tr><td>n<sub>f</sub> (-)</td><td>Water exponent</td><td>4</td><td>4</td><td>4</td></tr><tr><td>n<sub>r</sub> (-)</td><td>Gas exponent</td><td>3.5</td><td>3.5</td><td>3.5</td></tr><tr><td>λ (-)</td><td>Pore size index</td><td>1.2</td><td>1.2</td><td>1.2</td></tr></tbody></table>

the black-oil model used in the original setup. The 2D domain (Fig. 6) consists of a dome-shaped aquifer bounded by caprock and bedrock layers, extending 50 m horizontally and vertically. The upper boundary follows a sinusoidal profile $F(x) = \sigma + r \sin(\pi x)$ with $\sigma = 25$ m and $r = 5$ m, which enhances structural trapping of buoyant fluids. We analyze how microbial activity affects hydrogen plume evolution, mass loss, and H₂ recoverability during storage operations.

Hydrogen retention is ensured by the combined effects of high entry capillary pressure (Fig. 6) and low permeability in both caprock and bedrock. Lateral boundaries impose hydrostatic pressure conditions with a reference pressure $p_r = 40$ bar at $z = 0$, defined by $p_w = p_r + \rho g \Delta z$, while the top and bottom boundaries are no-flux. The domain is initially filled with brine in hydrostatic equilibrium. Fluid properties are listed in Table 7, with Brooks–Corey relative permeability and capillary pressure models governing multiphase flow. We simulate H₂ injection into pure water, with an initial composition of 84.8% water and 15.2% methane, resulting in an equilibrium saturation of [0.1522, 0.8478].
<figure><img src="image_8.png" /><figcaption>Fig. 4. 1D case: H₂ and CO₂ loss vs. CH₄ production (%) (top) and dissolved H₂ in the liquid phase (bottom), given by the ratio of liquid phase H₂ mass to the total H₂ mass.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 5. 1D case: Comparing PR and SW.</figcaption></figure>

A single vertical well, positioned at the aquifer center, is perforated exclusively below the caprock to ensure containment. The injection scenario consists of three distinct phases:

*   **Initial injection:** A gas mixture (90% CO₂, 10% H₂) is injected at 10 kg d⁻¹ for 120 days, followed by a 30-day shut-in period to allow for phase redistribution

*   **Cyclic storage operation:** 10 identical cycles (total duration: 24 months), each comprising:

- 30 days of pure H₂ injection at 18 kg d⁻¹

- 10-day idle phase for microbial activity

- 30-day production period at the same withdrawal rate

This benchmark case exhibits substantial physical complexity arising from the coupled interactions of multiphase flow, compositional effects, and geochemical reactions across heterogeneous rock types. While Ahmed et al. [17] originally used this configuration to examine molecular diffusion effects on plume migration, pressure evolution, and caprock integrity, our work extends the analysis by quantifying how microbial transformations affect both H₂ loss mechanisms and long-term recoverability. We expect this setup to pose significant numerical challenges due to strongly coupled nonlinearities from reaction-transport interactions, particularly through localized effects such as H₂ capillary leakage into low-permeability caprock zones and spatially variable bioclogging patterns that dynamically alter flow pathways.

### 3.2.2. Results and discussion

In Fig. 7, we present the spatio-temporal distribution of microbial density (after the build-up phase, first, and last injection cycles) and the corresponding mole fractions of CO₂, H₂, and CH₄ (after the first and last injection cycles). Note that bio-clogging effects are excluded in this first analysis. A distinct gas plume forms after the build-up phase, where the coexistence of H₂ and CO₂ supports microbial activity that progressively expands over time. During the first injection cycle, which uses a reduced CO₂ injection rate, we observe a noticeable decrease in CO₂ near the wellbore (with the maximum CO₂ mole fraction decreasing by 50%). The component plumes exhibit highly mixed gas saturation. By the final injection cycle, the maximum CH₄ mole fraction increases from 0.16 to 0.2, signaling enhanced microbial conversion. Additionally, dissolved H₂ infiltrates the caprock, though its concentration diminishes due to microbial activity (see the top-left figure for microbial distribution). The results indicate that in methanogenic UHS reservoirs, pure H₂ injection (rather than H₂-CO₂ mixtures)
<figure><img src="image_10.png" /><figcaption>Fig. 6. 2D aquifer:Ahmed et al. [17]: (left) 2D aquifer geometry and (right) capillary pressure curves used in the numerical experiments.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 7. 2D aquifer: Distribution of microbial (normalized) density and mole fractions after the build-up phase, as well as the first and last injection cycles.</figcaption></figure>

may be advantageous during build-up or conversion phases to control microbial activity.

Fig. 8 illustrates the impact of microbial activity on (i) component mass variation (top), (ii) chemical conversion rates (middle), and (iii) well efficiency (bottom). Methanogenesis commences immediately, resulting in an 18% reduction in H₂ mass by the simulation's end. While the total H₂ loss minimally affects production efficiency, the resulting CH₄ production compromises H₂ purity. The cyclic production patterns of H₂, CO₂, and CH₄ align with the scenario outlined earlier: H₂ and

CO₂ quantities increase during injection periods until reaching threshold values, then decline during production phases. Conversely, CH₄ levels (initially present in the aquifer) show a continuous rise during operations, as methane is not injected but produced solely through microbial activity. This leads to progressively higher H₂ impurity levels at the production stage.

Notably, bio-clogging amplifies chemical transformation (see Fig. 8, middle left) by extending residence times, thereby promoting additional reactions and restricting H₂ migration. Table 8 quantifies the final mixture composition changes. Although molecular diffusion
<figure><img src="image_12.png" /><figcaption>Fig. 8. 2D aquifer: Impact of biochemical transformation on cyclic H₂ injection.</figcaption></figure>

<table><caption>**Table 8**

2D aquifer: H₂ loss, CO₂ consumption, and CH₄ production with vs. without bio-clogging.</caption><thead><tr><th></th><th>H₂ loss (%)</th><th>CO₂ consumption (%)</th><th>CH₄ production (%)</th></tr></thead><tbody><tr><td>No clogging</td><td>18.21</td><td>8.44</td><td>1.20</td></tr><tr><td>Bio clogging</td><td>21.68</td><td>9.36</td><td>1.22</td></tr></tbody></table>

is neglected in this case, dissolved H₂ becomes effectively trapped in the aqueous phase through both limited mobility and microbial consumption, significantly reducing its migration potential toward the caprock.

Across all cases, the H₂ loss due to solubility remains approximately 1.5%, with potential for further reduction in higher salinity conditions. Microbial activity maintains significant transformation rates even under elevated salinity conditions, indicating primary control by microbial yield coefficients, decay rates, and gas phase composition. While pressure and temperature variations influence dissolution patterns similarly to salinity effects — remaining within ranges favorable for methanogens — they demonstrate negligible impact on overall H₂ loss, in agreement with findings from Safari et al. [21]. However, total losses from bacterial consumption strongly depend on reservoir conditions, particularly the initial microbial density. For example, losses can be reduced by half when $n_0 = 1.0 \times 10^8$ cells/m³ compared to the current higher seeding densities (see Fig. 9).

### 3.3. Test 3: UHS in a 3D regular reservoir

This case study is based on the geometry and scenario originally proposed by Khoshnevis et al. [10]. It examines the transport and evolution of H₂, CO₂, CH₄, and H₂O components and microbial populations in

<table><caption>**Table 9**

3D reservoir: Fluid properties.</caption><thead><tr><th>Phase</th><th>Density (kg m<sup>-3</sup>)</th><th>Viscosity (Pa s)</th><th>Compressibility (Pa<sup>-1</sup>)</th></tr></thead><tbody><tr><td>liquid</td><td>999.7</td><td>1.3059 &times; 10<sup>-3</sup></td><td>5.0015 &times; 10<sup>-10</sup></td></tr><tr><td>gas</td><td>1.2243</td><td>1.763 &times; 10<sup>-5</sup></td><td>1.0009 &times; 10<sup>-5</sup></td></tr></tbody></table>

a three-dimensional reservoir where H₂ and CO₂ (as cushion gas) are injected. In the current work, the reservoir is modeled as an aquifer containing either pure water or saline water.

#### 3.3.1. Simulation setup

The top of the reservoir is located at a depth of 1000 m below surface, with horizontal dimensions of 1525 m × 1525 m and a vertical thickness of 50 m. The reservoir rock has a homogeneous porosity $\Phi = 0.2$ and a heterogeneous absolute permeability tensor $K = [100 \text{ mD}, 100 \text{ mD}, 10 \text{ mD}]$. A vertical well is placed at the center of the reservoir's top surface.

The system is initialized to enable methanogenic archaea activity, with initial conditions of $P_0 = 100 \text{ bar}$, $T_0 = 313.15 \text{ K}$, and $S_{f0} = 0.2$. The initial composition comprises a liquid phase (70% H₂O) and a gas phase (2% CO₂ and 28% CH₄).

Fluid properties (Table 9) use standard water properties for the liquid phase and approximate gas properties at surface conditions ($P_{atm} = 10^5 \text{ Pa}$, $T_{atm} = 283.15 \text{ K}$). Note that gas is approximated as a binary mixture of 80% N₂ and 20% O₂, with data sourced from the NIST webbook.

Relative permeabilities are determined using the Brooks–Corey model [56], with exponent parameters ($\lambda_l, \lambda_g$) = (2, 2), entry pressure
<figure><img src="image_13.png" /><figcaption>Fig. 9. 2D aquifer: Porosity and permeability reduction due to bio-clogging from biochemical transformation.</figcaption></figure>

$P_{en} = 10^4$ Pa and residual saturations ($S_{lr}$, $S_{gr}$) = (0.2, 0.05) for each phase. The microbial diffusion effect is neglected in this setup.

Homogeneous Neumann boundary conditions are imposed. A Cartesian mesh is used for the reservoir with 31 × 31 × 8 cells along each spatial direction. The linearized systems are solved using either the default LU factorization or the ILU(0)-preconditioned GMRES method [38], with computational acceleration provided by MRST's AD backend. All simulations use a fixed time step of 1 d.

The injection scenario begins with a gas mixture (60% H₂ and 40% CO₂) injected at $10^6$ m³/d for a 60-day build-up phase, followed by a 20-day rest period and then a 30-day injection of a 95% H₂/5% CO₂ mixture. Subsequent phases include a 20-day idle period, a 30-day production phase at $10^6$ m³/d, and a final 20-day idle period. This scenario, adapted from Khoshnevis et al. [10], is extended to six cycles where subsequent cycles omit the build-up and rest periods.

Six simulation cases examine methanogenic archaea activity: (1) pure water aquifer without microbes; (2) pure water with microbial population $N_0 = 1$ ($n_0 = 10^8$ m⁻³); (3) saltwater without microbes; (4) saltwater with microbes; (5) pure water with higher microbial density ($n_0 = 10^9$ m⁻³); and (6) saltwater with higher density and molecular diffusion effects. The performance of both linear solvers (LU and ILU(0)-GMRES) is evaluated across all cases.

### 3.3.2. Results and discussion

We first examine the impact of microbial activity on UHS in a pure water aquifer ($m_{salt} = 0$). Fig. 10 shows microbial activity concentrated near the H₂-rich well, with population dynamics tracking H₂ and CO₂ availability. The results demonstrate a strong coupling between methanogenic archaea activity and gas availability: injection-phase gas increases promote microbial growth, while production-phase decreases lead to population decay (Fig. 11, right). Near the well, H₂ solubility reaches $1.3 \times 10^{-3}$ (compared to CO₂ at $9.8 \times 10^{-3}$), consistent with Section 2.2.2 results.

Moreover, Fig. 11 shows that microbial activity causes a decrease in pressure due to the reaction process. This pressure reduction correlates with a decrease in production efficiency, which is calculated as the total mass of H₂ produced during the production period divided by the total mass of H₂ injected into the well (see Table 10). Fig. 12 shows the production of H₂, CO₂ and C₁ in the well following the cyclic scenario described. At the end of the six cycles, Fig. 12 shows only a weak impact of microbial activity on production efficiency (approximately 0.4%).

To quantify the total loss of H₂, the consumption of CO₂, and the production of C₁ due to microbial activity, we compare the total mass of H₂ (resp. CO₂, C₁) over time in scenarios without and with microorganisms. After six cycles, microbial activity results in a 9.78% loss of H₂, a 0.26% consumption of CO₂, and less than 0.06% production of C₁ (see Fig. 13 and Table 10).

It has to be noted that these simulations have been performed with two different linear solvers: LU factorization and ILU(0)-preconditioned GMRES. As expected, the latter one is more efficient in terms of CPU time (see Table 11). For this 3D case with approximately

$38,440$ (without bacteria) and $46,128$ (with bacteria) degrees of freedom (DOFs), the introduction of microbial bacterial activity increases the computational time by 47% for LU factorization and by 30% when using ILU(0)-preconditioned GMRES. Overall, the results suggest that ILU(0)-preconditioned GMRES is more efficient for medium systems incorporating bacterial effects.

Secondly, we examine how salinity affects microbial activity in UHS by conducting the simulations in a saline aquifer ($m_{salt} = 3$). As expected, the salinity induces a decrease in the solubility of H₂ and CO₂ independently of the presence of microorganisms. The resulting solubilities near the well reaches $7.6 \times 10^{-4}$ for H₂ and $6.2 \times 10^{-3}$ for CO₂, matching Section 2.2.2 results for saltwater aquifers.

This salinity increase causes average reductions of 42% in H₂ solubility and 36% in CO₂ solubility, leading to a 1.8% decrease in archaea population (Fig. 11, right). Unlike microbial effects, salinity impacts neither reservoir pressure nor production efficiency, as shown by comparable results in pure and saltwater cases (Table 10).

Similarly, salinity has minimal influence (< 1%) on H₂ loss, CO₂ consumption, and CH₄ production (see Fig. 13 and Table 10) for more details.

In the third study, we analyze the impact of initial microbial density (we increase it to $n_0 = 10^9$ m⁻³) on production efficiency and gas losses. Fig. 11 reveals two key effects: a significant pressure decrease and microbial population dynamics synchronized with H₂ injection-production cycles. At this high density, the archaea population eventually declines as it exceeds the available H₂ in the liquid phase. Conversely, with $n_0 = 10^8$ m⁻³, the population increases throughout all 6 cycles but with diminishing growth rates, indicating a potential future threshold and decline.

With $n_0 = 10^9$ m⁻³, we observe H₂ loss reaching 11.01% (Fig. 13(a)) and CO₂ consumption (Fig. 13(b)) at 0.45% after the first cycle, while production efficiency (Fig. 13(c)) reaches 61.76% without microbes versus 60.98% with microbes. These measurements agree with Khoshnevis et al. [10]'s single-cycle findings.

The final analysis demonstrates that molecular diffusion in pure water aquifers reduces production efficiency by 0.73% while having minimal impact on gas losses (Table 10).

Overall, the results demonstrate that microbial activity in UHS significantly impacts H₂ storage efficiency, causing up to 9.8% H₂ loss and dynamic population shifts tied to gas availability. While salinity reduces gas solubilities, it does not significantly alter microbial pressure effects. Higher initial microbial densities can really amplify these losses substantially (up to 100 %). For computational performance, preconditioners (e.g., ILU(0)-GMRES or AMGCL-CPR [16,57]) are strongly recommended as both linear and nonlinear solvers — particularly flash calculations — incur significantly higher computational costs for complex mixtures. Finally, operational strategies including injection scheduling and gas composition control (CO₂ co-injection, staggered cycles, increasing salinity) can further improve storage efficiency while mitigating microbial impacts.
<figure><img src="image_14.png" /><figcaption>Fig. 10. 3D reservoir: Microbial activity (left) and H₂ solubility (right) in the vicinity of the well at the end of the first injection period.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 11. 3D reservoir: (Left) Mean pressure evolution; (right) total microbial population dynamics in pure and saline aquifers over six injection-production cycles.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12. 3D reservoir: H₂, CO₂ and C₁ well rates in the pure water aquifer.</figcaption></figure>

<caption>Table 10</caption>

3D reservoir: Gas loss and production from microbial activity in pure and saline aquifers, with resulting production efficiencies.

<table><thead><tr><th>Molecular diffusion (m s<sup>-2</sup>)</th><th>m<sub>salt</sub> (mol kgw<sup>-1</sup>)</th><th>n<sub>0</sub> (m<sup>-3</sup>)</th><th>H<sub>2</sub> loss (%)</th><th>CO<sub>2</sub> loss (%)</th><th>C<sub>1</sub> prod. (%)</th><th>Effic., no archae (%)</th><th>Effic., archae (%)</th></tr></thead><tbody><tr><td>No</td><td>0</td><td>10<sup>8</sup></td><td>9.78</td><td>0.26</td><td>0.006</td><td>94.28</td><td>93.88</td></tr><tr><td>No</td><td>0</td><td>10<sup>9</sup></td><td>47.63</td><td>1.38</td><td>0.12</td><td>94.28</td><td>90.13</td></tr><tr><td>No</td><td>3</td><td>10<sup>8</sup></td><td>9.02</td><td>0.24</td><td>0.01</td><td>94.27</td><td>93.01</td></tr><tr><td>No</td><td>3</td><td>10<sup>9</sup></td><td>46.63</td><td>1.35</td><td>0.12</td><td>94.27</td><td>90.42</td></tr><tr><td>Yes</td><td>0</td><td>10<sup>8</sup></td><td>10.42</td><td>0.3</td><td>0.0033</td><td>93.55</td><td>93.15</td></tr></tbody></table>
<figure><img src="image_17.png" /><figcaption>Fig. 13. 3D reservoir: H₂ loss, CO₂ consumption, and production efficiency in pure and saline aquifers over six injection-production cycles.</figcaption></figure>

<table><caption>**Table 11**
3D reservoir: CPU times using LU factorization and ILU(0) preconditioned GMRES.</caption><thead><tr><th>Linear solver</th><th>Archae</th><th>CpuTime (s)</th></tr></thead><tbody><tr><td>LU factorization</td><td>no archae</td><td>973</td></tr><tr><td>LU factorization</td><td>10<sup>8</sup></td><td>1263</td></tr><tr><td>ILU(0) preconditioned GMRES</td><td>no archae</td><td>320</td></tr><tr><td>ILU(0) preconditioned GMRES</td><td>10<sup>8</sup></td><td>471</td></tr></tbody></table>

<table><caption>**Table 12**
UHS benchmark: Initial and injected gas compositions in mole fractions.</caption><thead><tr><th>Component</th><th>Initial (%)</th><th>Injected (%)</th></tr></thead><tbody><tr><td>Methane (CH<sub>4</sub>)</td><td>7.015</td><td>17.522</td></tr><tr><td>Water (H<sub>2</sub>O)</td><td>90.230</td><td>0.000</td></tr><tr><td>Ethane (C<sub>2</sub>H<sub>6</sub>)</td><td>0.270</td><td>0.144</td></tr><tr><td>n-Butane (C<sub>4</sub>H<sub>10</sub>)</td><td>0.040</td><td>0.000</td></tr><tr><td>Pseudo-component (C<sub>3</sub><sup>+</sup>)</td><td>0.020</td><td>0.012</td></tr><tr><td>Hydrogen (H<sub>2</sub>)</td><td>0.000</td><td>80.000</td></tr><tr><td>Carbon Dioxide (CO<sub>2</sub>)</td><td>2.021</td><td>0.504</td></tr><tr><td>Nitrogen (N<sub>2</sub>)</td><td>0.405</td><td>1.818</td></tr><tr><td>Sum</td><td>100.000</td><td>100.000</td></tr></tbody></table>

## 3.4. Test 4: 3D underground hydrogen storage with complex mixture

### 3.4.1. Simulation setup

This test case is based on a large-scale 3D benchmark developed by the Institute of Subsurface Energy Systems (TU Clausthal) and studied in detail by Hogeweg et al. [58]. The model represents a typical sandstone reservoir located at a depth of 1210m, with horizontal dimensions of 3050m × 3050m and variable thickness. The domain is discretized using a corner-point grid consisting of 61 × 61 × 12 cells (a total of 44,652 cells). The reservoir has an average porosity of $\Phi = 0.15$ and an absolute permeability tensor given by $K = [143mD, 143mD, 3mD]$.

A vertical injection well is placed at the center of the top surface, injecting a gas mixture rich in H₂ (see Table 12 for the full composition). The system is initialized with a pressure of $P_0 = 81.6$ bar, and temperature of $T_0 = 333.15$ K. Unlike the original setup, which uses depth-dependent initial compositions and phase distributions, our MRST implementation assumes spatially uniform initial conditions.

The operational schedule comprises two phases. During the initial build-up phase, bottom-hole pressure (BHP) controlled injection gradually increases from 90 bar to 102 bar in 4 bar increments, enhancing reservoir pressure and hydrogen content. This phase consists of four 60 d cycles, each followed by 30 d idle periods. Subsequently, the regular storage phase commences, featuring alternating injection and production at a constant rate of 293.13 $mol s^{-1}$ ($6 \times 10^5 Sm^3/d$). Both injection and withdrawal periods last 90 d to maintain volume balance, with idle intervals separating cycles as in the conversion phase. The injected gas maintains a constant composition (Table 12) with 80% H₂ content as the primary component, balanced by 20% of cushion gas. All remaining model parameters are taken from Hogeweg et al. [58].

<table><caption>**Table 13**
UHS benchmark: Comparison of simulation performance for abiotic and microbial cases.</caption><thead><tr><th>Metric</th><th>Abiotic</th><th>n<sub>0</sub> = 10<sup>8</sup></th><th>n<sub>0</sub> = 10<sup>9</sup></th></tr></thead><tbody><tr><td>Total Time [s]</td><td>1132</td><td>2382</td><td>5174</td></tr><tr><td>Time per Step [s]</td><td>3.7</td><td>7.8</td><td>16.9</td></tr><tr><td>Linear Its (Total)</td><td>574</td><td>3714</td><td>7366</td></tr><tr><td>Linear Its/step (Mean)</td><td>1.8</td><td>12.2</td><td>23.2</td></tr><tr><td>Nonlinear Its (Total)</td><td>494</td><td>1238</td><td>2135</td></tr><tr><td>Nonlinear Its/step (Mean)</td><td>1.9</td><td>4.2</td><td>7.4</td></tr></tbody></table>

### 3.4.2. Results and discussion

This test focuses on the performance and behavior of simulations with and without microbial activity. The formulation of the overall composition is extended with microbial activity, resulting in approximately 446,520 DOFs (401,000 for abiotic). To enhance computational performance, we employ the accelerated AD backend in MRST, together with the AMGCL_CPR [57,59] preconditioner and a quasi-IMPES decoupling (block size = components + 1 for bacteria). We consider two different initial microbial densities, $n_0 \in \{10^8, 10^9\}$. Higher microbial concentrations intensify the coupling and increase the nonlinearity of the system as seen in the previous case.

Fig. 14 depicts the bottom-hole pressure (left) and gas injection/ production (right) rates throughout the simulation schedule. The key observation is that the simulation successfully maintains all schedule constraints, with no control switching interruptions observed during the operational cycles. Note that no differences in bottom-hole pressure (BHP) or gas well rates are observed in the abiotic simulation case. Fig. 15 depicts the spatial distribution of microbial population (left) and H₂ mole fraction (right) after the build-up phase. As expected, both are concentrated near the well region where the H₂ plume forms. In general, as in the other cases, the microbial activity expands and shrinks following the injection schedule of H₂. Interestingly, the microbial plume extends slightly beyond the H₂ plume, driven by the presence of even trace amounts of H₂ and CO₂ outside the core plume region. As illustrated in Fig. 16 (top left), the total loss of H₂ is approximately 7.5% for $n_0 = 10^8$ and 19.1% for $n_0 = 10^9$. In the high-density case ($n_0 = 10^9$), the produced gas shows increasing impurity from methane (bottom-left figure) production (partially offset by decreasing CO₂ content, see bottom-right figure). By contrast, for $n_0 = 10^8$, the produced gas composition remains largely intact, though the trend suggests a gradual impurity increase over time. Notably, microbial activity does not affect H₂ production in either case, as demonstrated in the top-right figure.

Table 13 reveals the important performance impact of including microbial activity. For $n_0 = 10^8$, the total simulation time increased from 1132 to 2382 (s) (+110%), and further to 5174 (s) (+357%) for $n_0 = 10^9$. This increase is attributed to the stronger nonlinear coupling introduced by microbial reactions, which led to more frequent time step cutting and a significant rise in solver effort due to more demanding
<figure><img src="image_18.png" /><figcaption>Fig. 14. UHS benchmark: BHP versus time (left) and gas rate versus time (right).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 15. UHS benchmark: Microbial population distribution (left) and H₂ mole fraction at the end of build-up phase (right).</figcaption></figure>

flash calculations, as the evolving composition required frequent and
more complex updates at each step.

Overall, the compositional-biochemistry simulator proves reliable for simulating fluid dynamics in UHS for small or medium size models. However, for large models with multiple components (with millions of degrees of freedom), it becomes necessary to consider analytical model reduction approaches [60], or treating the microbial dynamics as a fictitious pseudo-component in the aqueous phase [10]. Note also that the current study focuses on simulation infrastructure, with biochemical models primarily based on literature observations. For realistic scenarios, these models generally require calibration against actual laboratory investigations to match field conditions [13].

4. Conclusion

In this work, we have investigated the complex interplay of com-
positional, chemical, and biological processes involved in underground
hydrogen storage (UHS). To this end, we have implemented in MRST:
a SW EoS model to capture the solubility of gases, a double Monod
model to describe microbial dynamics, a molecular diffusion model,
and a bio-clogging model to assess the impact of microbial activity on
rock properties.

We conducted numerous simulations involving liquid-–gas phases
and multicomponent systems with microbial activity. The simulation
results are consistent with findings from the literature, particularly in
showing the influence of salinity on the solubility of H₂ and CO₂, which
in turn affects microbial dynamics and hydrogen loss due to microbial
consumption. This microbial activity subsequently alters rock proper-
ties, contributing further to H₂ loss in the aquifer. Simulations show
that hydrogen loss can reach up to 10% for initial microbial densities

reported in the literature, increasing to 47% at higher densities. This
is accompanied by H₂ impurities from CH₄ production. These results
highlight the significant impact of microbial activity on hydrogen stor-
age and underscore the need to incorporate bio-geochemical processes
into UHS simulations.

Moreover, molecular diffusion effects were found to be relatively weak under the simulated conditions. However, given the fact dispersive transport play a more significant role, incorporating dispersion mechanisms in future models is essential for a more realistic representation of multicomponent transport in porous media. In addition, microbial dynamics in UHS environments can involve multiple archaeal and bacterial populations coexisting and interacting. These populations are not only capable of substrate-specific metabolism but can also exhibit mobility driven by chemotaxis toward favorable gradients. Such behaviors can lead to spatial heterogeneities and competitive interactions that significantly affect microbial growth, gas composition, and overall system dynamics. Therefore, in a future work, we aim at integrating both hydrodynamic dispersion and microbial chemotaxis, along with competitive and cooperative microbial kinetics, to better capture the complexity of subsurface bio-reactive transport phenomena.

Overall, these expansions retain value for small-to-medium scale
models and investigative studies, as previously demonstrated. For large-
scale compositional simulations of complex mixtures, achieving a good
balance between physical fidelity and computational efficiency be-
comes critical for field projects. In such cases, practitioners may con-
sider reduced-order models or the integration of microbial populations
as fictitious pseudo-components in the liquid phase. Methods based
on the Arrhenius equation can be used for the approximate calcu-
lation of reaction rates—an approach particularly advantageous for
implementation in commercial simulators with limited support for
<figure><img src="image_20.png" /><figcaption>Fig. 16. UHS benchmark: Impact of biochemical transformation on cyclic H₂ injection. Total H₂ loss amounts for 7.5% for n₀ = 1.0e8 and 19.1% for n₀ = 1.0e9.</figcaption></figure>

custom biochemical model integration. The Søreide–Whitson equation of state, while introducing minimal computational overhead compared to Peng–Robinson in the studied cases, necessitates calibration against experimental data where available. This requirement extends equally to microbial dynamics modeling, where laboratory and field-scale validation remains essential for ensuring biologically realistic predictions of subsurface hydrogen storage systems.

### CRediT authorship contribution statement

**Elyes Ahmed:** Writing – review & editing, Writing – original draft, Validation, Software, Project administration, Methodology, Investigation, Conceptualization.
**Brahim Amaziane:** Writing – review & editing, Supervision, Resources, Project administration, Methodology, Conceptualization.
**Salaheddine Chabab:** Writing – review & editing, Validation, Methodology.
**Stéphanie Delage Santacreu:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation.
**Guillaume Galliéro:** Writing – review & editing, Validation, Supervision, Formal analysis, Conceptualization.
**Olav Møyner:** Writing – review & editing, Validation, Software, Methodology, Conceptualization.
**Xavier Raynaud:** Writing – review & editing, Project administration, Methodology, Funding acquisition.

### Funding

The authors affiliated with Sintef Digital acknowledge funding from the Norwegian Research and Innovation Centre for Hydrogen and Ammonia HYDROGENi (Grant No. 333118). Elyes Ahmed acknowledges support from KS Hydrogen: Konsernsatsinger i SINTEF for hydrogen.

The work of Brahim Amaziane and Stéphanie Delage Santacreu is partially supported by the ANR (French National Research Agency) within the HyStorEn project n° 22-CE05-0024, the Carnot Institute ISIFoR, and CNRS. The work of Stéphanie Delage Santacreu, Brahim Amaziane, Guillaume Galliéro and Salaheddine Chabab is carried out within the framework of the GdR HydroGEMM :Hydrogène du sous-sol: étude intégrée de la Genèse... à la Modélisation Mathématique. Their support is gratefully acknowledged.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Elyes Ahmed reports financial support was provided by Research Council of Norway. Brahim Amaziane reports financial support was provided by French National Research Agency. Stephanie Delage Santacreu reports financial support was provided by French National Research Agency. Xavier Raynaud reports financial support was provided by Research Council of Norway. Olav Moyner reports financial support was provided by Research Council of Norway. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### References

[1] Wulf C, Linssen J, Zapp P. Chapter 9 - power-to-gas—Concepts, demonstration, and prospects. In: Azzaro-Pantel C, editor. Hydrogen supply chains. Academic Press; 2018, p. 309–45. http://dx.doi.org/10.1016/B978-0-12-811197-0.00009-9, URL https://www.sciencedirect.com/science/article/pii/B9780128111970000099.
[2] Zivar D, Kumar S, Foroozesh J. Underground hydrogen storage: A comprehensive review. Int J Hydrogen Energy 2021;46(45):23436–62. http://dx.doi.org/10.1016/j.ijhydene.2020.08.138, Hydrogen Separation, Production and Storage.

[3] Chabab S, Théveneau P, Coquelet C, Corvisier J, Paricaud P. Measurements and predictive models of high-pressure H2 solubility in brine (H2O+NaCl) for underground hydrogen storage application. Int J Hydrog Energy 2020;45(56):32206–20. http://dx.doi.org/10.1016/j.ijhydene.2020.08.192, URL https://www.sciencedirect.com/science/article/pii/S0360319920332195.

[4] Panfilov M. Physicochemical fluid dynamics in porous media: applications in geosciences and petroleum engineering. Wiley-VCH; 2018.

[5] Hemme C, Van Berk W. Hydrogeochemical modeling to identify potential risks of underground hydrogen storage in depleted gas fields. Appl Sci 2018;8(11). http://dx.doi.org/10.3390/app8112282, URL https://www.mdpi.com/2076-3417/8/11/2282.

[6] Siegert M, Cichocka D, Herrmann S, Gründger F, Feisthauer S, Richnow H-H, et al. Accelerated methanogenesis from aliphatic and aromatic hydrocarbons under iron- and sulfate-reducing conditions. FEMS Microbiol Lett 2011;315(1):6–16. http://dx.doi.org/10.1111/j.1574-6968.2010.02165.x.

[7] Ginn TR, Wood BD, Nelson KE, Scheibe TD, Murphy EM, Clement TP. Processes in microbial transport in the natural subsurface. Adv Water Resour 2002;25(8):1017–42.

[8] Maier R, Pepper IL, Gerba CP. Environmental microbiology. Maier and pepper set, Elsevier Science; 2009, URL https://books.google.fr/books?id=A2zl8YBXQfoC.

[9] Eddaoui N, Panfilov M, Ganzer L, Hagemann B. Impact of pore clogging by bacteria on underground hydrogen storage. Transp Porous Media 2021;139:89–108. http://dx.doi.org/10.1007/s11242-021-01647-6.

[10] Khoshnevis N, Hogeweg S, Goncalves Machado C, Hagemann B. Numerical modeling of bio-reactive transport during underground hydrogen storage – a benchmark study. Eur Assoc Geosci & Eng 2023;1:1–5. http://dx.doi.org/10.3997/2214-4609.202321087.

[11] Ali M, Arif M, Sedev R, Sánchez-Román M, Keshavarz A, Iglauer S. Underground hydrogen storage: The microbiotic influence on rock wettability. J Energy Storage 2023;72:108405. http://dx.doi.org/10.1016/j.est.2023.108405, URL https://www.sciencedirect.com/science/article/pii/S2352152X23018029.

[12] Thaysen EM, McMahon S, Strobel GJ, Butler IB, Ngwenya BT, Heinemann N, et al. Estimating microbial growth and hydrogen consumption in hydrogen storage in porous media. Renew Sustain Energy Rev 2021;151:111481. http://dx.doi.org/10.1016/j.rser.2021.111481, URL https://www.sciencedirect.com/science/article/pii/S1364032121007620.

[13] Hogeweg S, Hagemann B, Bobrov V, Ganzer L. Development and calibration of a bio-geo-reactive transport model for UHS. Front Energy Res 2024;12:1385273.

[14] Dopffel N, Mayers K, Kedir A, Alagic E, An-Stepec BA, Djurhuus K, et al. Microbial hydrogen consumption leads to a significant pH increase under high-saline-conditions: implications for hydrogen storage in salt caverns. Sci Rep 2023;13(1):10564.

[15] Xiaofei Z, Yuxin T, Xinwei T, Maozheng H, Xuefang L, Yun P, et al. Improved mathematical model simulates methane production during anaerobic digestion of food waste under stable and unstable acidification conditions. J Environ Chem Eng 2025;13(2):115833.

[16] Lie K-A. An introduction to reservoir simulation using MATLAB/GNU Octave: user guide for the matlab reservoir simulation toolbox (MRST). Cambridge University Press; 2019.

[17] Ahmed E, Møyner O, Raynaud X, Nilsen HM. Phase behavior and black-oil simulations of hydrogen storage in saline aquifers. Adv Water Resour 2014;191:104772. http://dx.doi.org/10.1016/j.adwatres.2014.104772, URL https://www.sciencedirect.com/science/article/pii/S0309170824001593.

[18] Hagemann B. numerical and analytical modeling of gas mixing and bio-reactive transport during underground hydrogen storage (Ph.D. thesis), (2017LORR0328). Université de Lorraine; 2017, URL https://theses.hal.science/tel-01735019.

[19] Eddaoui N. Patterns in bioreactive transport in underground storage of hydrogen: impact of natural and induced heterogeneity (Ph.D. thesis), Université de Lorraine en cotutelle avec Université Hassan II (Casablanca, Maroc); 2021, URL http://www.theses.fr/2021LORR0244/document, Thèse de doctorat dirigée par Panfilov, Mikhail et Aniss, Saïd Mécanique - Génie civil Université de Lorraine 2021.

[20] Monod J. The growth of bacterial cultures. Annu Rev Microbiol 1949;3(1):371–94.

[21] Safari Z, Fatehi R, Azin R. Developing a numerical model for microbial methanation in a depleted hydrocarbon reservoir. Renew Energy 2024;227:120426. http://dx.doi.org/10.1016/j.renene.2024.120426, URL https://www.sciencedirect.com/science/article/pii/S0960148124004919.

[22] Becker S, Kuznetsov A. Heat transfer and fluid flow in biological processes. Academic Press; 1948.

[23] Toleukhanov A, Panfilov M, Kaltayev A. Storage of hydrogenous gas mixture in geological formations: Self-organisation in presence of chemotaxis. Int J Hydrog Energy 2015;40(46):15952–62. http://dx.doi.org/10.1016/j.ijhydene.2015.10.033.

[24] Elferink SO, Visser A, Pol LWH, Stams AJ. Sulfate reduction in methanogenic bioreactors. FEMS Microbiol Rev 1994;15(2-3):119–36.

[25] Elferink SO, Maas RN, Harmsen H, Stams AJ. Desulforhab-dus amnigenus gen. nov. sp. nov. a sulfate reducer isolated from anaerobic granular sludge. Arch Microbiol 1995;164(2):119–24.

[26] Karadagli F, Rittmann BE. Kinetic characterization of methanobacterium bryanti moh. Environ Sci Technol 2005;39(13):4900–5.

[27] Robinson JA, Tiedje JM. Competition between sulfate-reducing and methanogenic bacteria for H2 under resting and growing conditions. Arch Microbiol 1984;137(1):26–32.

[28] Kristjansson JK, Schönheit P, Thauer RK. Different ks values for hydrogen of methanogenic bacteria and sulfate reducing bacteria: an explanation for the apparent inhibition of methanogenesis by sulfate. Arch Microbiol 1982;131(3):278–82.

[29] Ebigbo A, Golfier F, Quintard M. A coupled, pore-scale model for methanogenic microbial activity in underground hydrogen storage. Adv Water Resour 2013;61:74–85.

[30] Stams AJM, Elferink SO, Westermann P. Metabolic interactions between methanogenic consortia and anaerobic respiring bacteria. In: Biomethanation I. Springer; 2003, p. 31–56.

[31] Vavilin VA, Lokshina LY, Rytov SV, Kotsyurbenko OR, Nozhevnikova AN. Description of two-step kinetics in methane formation during psychrophilic H2/CO2 and mesophilic glucose conversions. Bioresour Technol 2000;71(3):195–209.

[32] Dorseiffer P, Meyer B, Heinzle E. Modeling of anaerobic formate kinetics in mixed biofilm culture using dynamic membrane mass spectrometric measurement. Biotechnol Bioeng 1995;45(3):219–28.

[33] Kotsyurbenko OR, Glagolev MV, Nozhevnikova AN, Conrad R. Competition between homoaceticogenic bacteria and methanogenic archaea for hydrogen at low temperature. FEMS Microbiol Ecol 2001;38(2-3):153–9. http://dx.doi.org/10.1111/j.1574-6941.2001.tb00893.x.

[34] Voskov DV, Tchelepi HA. Comparison of nonlinear formulations for two-phase multi-component EoS based simulation. J Pet Sci Eng 2012;82–83:101–11. http://dx.doi.org/10.1016/j.petrol.2011.10.012, URL https://www.sciencedirect.com/science/article/pii/S0920410511002683.

[35] Millington RJ, Quirk JP. Permeability of porous solids. Trans Faraday Soc 1961;57:1200–7. http://dx.doi.org/10.1039/TF9615701200.

[36] Møyner O. Compositional simulation with the AD-OO framework. Adv Model the MATLAB Reserv Simul Toolbox 2021;324–74.

[37] Krogstad S, Lie K-A, Møyner O, Nilsen HM, Raynaud X, Skaflestad B. MRST-AD-an open-source framework for rapid prototyping and evaluation of reservoir simulation problems. In: 2015 SPE reservoir simulation conference. Society of Petroleum Engineers; 2015.

[38] Møyner O. Better AD simulators with flexible state functions and accurate discretizations. In: Lie K-A, Møyner O, editors. Advanced modelling with the MATLAB reservoir simulation toolbox. Cambridge University Press; 2021, p. 151–91.

[39] Shojaee A, Ghanbari S, Wang G, Gregory S, Dopffel N, Mackay E. New flow simulation framework for underground hydrogen storage modelling considering microbial and geochemical reactions. Int J Hydrog Energy 2025;159:150453.

[40] Parkhurst David L, Appelo CAJ. User's guide to PHREEQC (Version 2): A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. Technical report, US Geological Survey; 1999.

[41] Wilke CR, Chang Pin. Correlation of diffusion coefficients in dilute solutions. AIChE J 1955;1(2):264–70. http://dx.doi.org/10.1002/aic.690010222, arXiv:https://aiche.onlinelibrary.wiley.com/doi/pdf/10.1002/aic.690010222, URL https://aiche.onlinelibrary.wiley.com/doi/abs/10.1002/aic.690010222.

[42] Cusler EL. Diffusion: Mass transfer in fluid systems. Cambridge series in chemical engineering, third ed. Cambridge University Press; 2009.

[43] Gibbs JW. The collected works of J. Willard Gibbs. NewHaven, CT: Yale University Press; 1948.

[44] Helmig R. Multiphase flow and transport processes in the subsurface: a contribution to the modeling of hydrosystems. 1997.

[45] Lee LL. Molecular thermodynamics of electrolyte solutions. 2nd ed. World Scientific; 2021, http://dx.doi.org/10.1142/12181, arXiv:https://www.worldscientific.com/doi/pdf/10.1142/12181, URL https://www.worldscientific.com/doi/abs/10.1142/12181.

[46] Rahbari A, Brenkman J, Hens R, Ramdin M, van den Broeke LJ, Schoon R, et al. Solubility of water in hydrogen at high pressures: A molecular simulation study. J Chem Eng Data 2019;64(9):4103–15. http://dx.doi.org/10.1021/acs.jced.9b00513.

[47] Soreide I, Whitson CH. Peng-robinson predictions for hydrocarbons, CO2, N2, and H2 s with pure water and NaCl brine. Fluid Phase Equilib 1992;77:217–40. http://dx.doi.org/10.1016/0378-3812(92)85105-H, URL https://www.sciencedirect.com/science/article/pii/037838129285105H.

[48] Chabab S, Théveneau P, Corvisier J, Coquelet C, Paricaud P, Houriez C, et al. Thermodynamic study of the CO2 – H2O – NaCl system: Measurements of CO2 solubility and modeling of phase equilibria using Soreide and Whitson, electrolyte CPA and SIT models. Int J Greenh Gas Control 2019;91:102825. http://dx.doi.org/10.1016/j.ijggc.2019.102825, URL https://www.sciencedirect.com/science/article/pii/S1750583619302622.
[49] Chabab S, Kerkache H, Bouchkira I, Poulain M, Baudouin O, Moine E, et al. Solubility of H2 in water and NaCl brine under subsurface storage conditions: Measurements and thermodynamic modeling. Int J Hydrog Energy 2024;50:648–58. http://dx.doi.org/10.1016/j.ijhydene.2023.10.290, URL https://www.sciencedirect.com/science/article/pii/S0360319923054976.

[50] Afanasyev A, Vedeneeva E. Compositional modeling of multicomponent gas injection into saline aquifers with the MUFITS simulator. J Nat Gas Sci Eng 2021;94:103988. http://dx.doi.org/10.1016/j.jngse.2021.103988.

[51] Chabab S, Cruz JL, Poulain M, Ducoussou M, Contamine F, Serin JP, et al. Thermodynamic modeling of mutual solubilities in gas-laden brines systems containing CO2, CH4, N2, O2, H2, H2o, NaCl, CaCl2, and KCl: Application to degassing in geothermal processes. Energies 2021;14(17). http://dx.doi.org/10.3390/en14175239, URL https://www.mdpi.com/1996-1073/14/17/5239.

[52] Soave G. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem Eng Sci 1972;27:1197–203.

[53] Petitfrere M, Patacchini L, de Loubens R. Three-phase eos-based reservoir simulation with salinity dependent phase-equilibrium calculations. In: ECMOR XV-15th European conference on the mathematics of oil recovery. European Association of Geoscientists & Engineers; 2016, p. cp–494.

[54] Panfili P, Patacchini L, Ferrari A, Garipov T, Esler K, Cominelli A. Implementation of soreide and whitson EoS in a GPU-based reservoir simulator. Comput Geosci 2024;28(2):341–54.

[55] Møyner O, Tchelepi HA. A mass-conservative sequential implicit multi-scale method for isothermal equation-of-state compositional problems. SPE J 2018;23(06):2376–93.

[56] Brooks RH, Corey AT. Hydraulic properties of porous media. In: Colorado state university hydrology paper. Colorado: Colorado State University; 1964, p. 1–37.

[57] Nilsen HM, Ahmed E, Rasmussen AF, Bao K, Skille T. Constrained pressure residual preconditioner including wells for reservoir simulation. In: SPE reservoir simulation conference. 2023, http://dx.doi.org/10.2118/212172-MS, arXiv: https://onepetro.org/spersc/proceedings-pdf/23RSC/23RSC/D021S009R001/3091886/spe-212172-ms.pdf, D021S009R001.

[58] Hogeweg S, Strobel G, Hagemann B. Benchmark study for the simulation of underground hydrogen storage operations. Comput Geosci 2022;26(6):1367–78. http://dx.doi.org/10.1007/s10596-022-10163-5.

[59] Demidov Denis. AMGCL: An efficient, flexible, and extensible algebraic multigrid implementation. Lobachevskii J Math 2019;40:535–46.

[60] Hagemann Birger, Hogeweg Sebastian, Strobel Gion. Application of analytical solutions of the reactive transport equation for underground methanation reactors. Transp Porous Media 2024;151(14):2601–23.