# A compressible pore-scale numerical simulation of hydrogen flow into brine: Application to underground hydrogen storage

Bourama Diarra, Antonio Rodríguez de Castro, Sakir Amiroudine

► To cite this version:

Bourama Diarra, Antonio Rodríguez de Castro, Sakir Amiroudine. A compressible pore-scale numerical simulation of hydrogen flow into brine: Application to underground hydrogen storage. Chemical Engineering Science, 2026, 320, pp.122583. 10.1016/j.ces.2025.122583. hal-05379972

HAL Id: hal-05379972

https://artsetmetiers.hal.science/hal-05379972v1

Submitted on 24 Nov 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>

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

# A compressible pore-scale numerical simulation of hydrogen flow into brine: Application to underground hydrogen storage

Bourama Diarra<sup>a</sup>, Antonio Rodríguez de Castro<sup>b,*</sup>, Sakir Amiroudine<sup>a,†</sup>

## ARTICLE INFO

**Keywords:**

Multiphase flow
Porous Media
Compressible fluid
Hydrogen Storage
OpenFOAM

## ABSTRACT

This study investigates the dynamics of hydrogen-brine flow in idealized pore geometries, with the aim of improving underground hydrogen storage in saline aquifers. Direct numerical simulations (DNS) using Open-FOAM are carried out to assess the influence of key parameters including injection flow rate, pore geometry, dynamic contact angle, and fluid compressibility on immiscible displacement at the pore scale. Results demonstrate marked differences between compressible and incompressible models in terms of brine sweep efficiency, interface displacement patterns, and pressure drop for a given flow rate. Incompressible simulations fail to capture critical phenomena such as hydrogen bubble formation and associated pressure fluctuations. Variations in flow rate and geometric constriction significantly impact residual brine saturation and inlet pressure dynamics; lower rates reduce pressure buildup and leakage risk while increasing storage efficiency. Incorporating a dynamic contact angle reduces the capillary resistance, accelerates the flow, and mitigates pressure peaks compared to static angle models. Overall, this work demonstrates the applicability of OpenFOAM to multiphase compressible hydrogen-brine flows at the pore scale, providing application-specific validation and novel insights to guide the design and optimization of underground hydrogen storage systems.

## 1. Introduction

The urgent need to reduce greenhouse gas emissions has driven the increasing adoption of renewable energy sources such as solar and wind power. However, their intermittent nature presents significant challenges in balancing energy supply and demand. A promising solution involves storing excess renewable energy in the form of hydrogen within underground geological formations, such as saline aquifers and depleted hydrocarbon reservoirs (Shi et al., 2020; Tarkowski, 2019; Yan et al., 2025). Underground engineering, in its various energy and environmental applications, relies on a deep understanding of multiphase flows in soils composed of layers of unconsolidated granular media and consolidated rocks with pores of different sizes. This knowledge is crucial for sectors such as CO₂ storage, soil remediation, hydrocarbon recovery, geothermal energy production, and natural gas storage in geological formations (Lake et al., 2014; Li and Benson, 2015; Rocco et al., 2017; Yin et al., 2018). A precise understanding of these flows enables the development of sustainable solutions to current energy and environmental challenges.

Numerous studies have explored multiphase flows in porous media

using numerical simulations. Mowla and Ioannidis (2020) analyzed the dynamic effect of wettability (dynamic contact angle) on multiphase flow dynamics as a function of porous medium wettability, using a laboratory model of a gas-liquid-solid fluidized bed system. Hashemi et al. (2021) employed pore network modeling to estimate hydrogen-brine transport properties, specifically capillary pressure and relative permeability. Their results demonstrated that these properties are particularly sensitive to the contact angle (wettability), pore structure, and clay content. In another study, Dai et al. (2017) simulated crude oil invasion into a water-filled rough fracture using the level set method. They found that roughness and wettability affect invasion pressure and water retention. Higher pressure and wettability reduce breakthrough time and residual water, with roughness having the strongest influence. Furthermore, Abdelsamie et al. (2016) developed a DNS model for reactive and turbulent diphasic flows at low Mach numbers, assuming an incompressible fluid. Also, Rodríguez de Castro and Goyeau (2022) used circle matrices to model solid grains, analyzing the impact of the contact angle, capillary number, fluid rheology, and dynamic viscosity ratio on displacement efficiency. Their pore-scale study, based on incompressible fluids with different rheological behaviors, modeled interface

<footnote><sup>a</sup> Univ. Bordeaux, CNRS, Bordeaux INP, Arts et Métiers Paris Tech, I2M UMR 5295, Talence F-33400, France</footnote>

<footnote><sup>b</sup> Arts et Métiers Institute of Technology, CNRS, Univ. Bordeaux, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: antonio.rodriguezdecastro@ensam.eu (A. Rodríguez de Castro).</footnote>
dynamics using the Level Set method in COMSOL software. Wang et al. (2023) analyzed the effect of pressure on trapping mechanisms during drainage and imbibition cycles using Direct Numerical Simulations (DNS). However, their study focused solely on the influence of pressure on the contact angle and the interfacial tension between hydrogen and water. The effect of hydrogen compressibility on trapping mechanisms in pore-scale multiphase flow was not considered in their analysis.

Derijani et al. (2023) evaluated the interFoam solver from Open-FOAM for simulating immiscible two-phase flows at the pore scale using a pore doublet geometry. Their work, based on incompressible fluids, was validated against experimental data for both drainage and imbibition processes, successfully capturing key mechanisms such as preferential flow through wider pores and early invasion of the wetting phase in narrower constrictions. A detailed sensitivity analysis was also conducted, revealing that mesh resolution, contact angle dynamics, and wall slip conditions significantly affect interface stability and wetting film development. Similarly, Rabbani et al. (2017) employed high-resolution DNS with OpenFOAM to study diphasic flows in porous media under intermediate-wet conditions. Their results, also based on incompressible fluids, demonstrated that pore geometry critically shapes interface dynamics, with the coexistence of concave and convex menisci. The complexity of interfacial motions observed under these wettability conditions led to counter-intuitive macroscopic behaviors, such as reduced recovery of the defending phase compared to classical drainage or imbibition regimes.

Most pore-scale models for hydrogen-brine systems in porous media have historically assumed incompressible fluids, neglecting hydrogen compressibility (Tarkowski, 2019; Huang et al., 2023; Wang et al., 2023). Outside the context of porous media, a few authors (Caiden et al., 2001; Hao and Prosperetti, 2004), coupled compressibility with interface tracking methods by omitting thermal and viscous dissipation effects. Also, Dellacherie (2005) developed a Diphasic Low Mach Number framework for simulating two-phase flows without phase change, applicable to general equations of state. While this model captures natural properties related to bubble dilation and compression in low Mach number flows, it does not specifically address hydrogen transport in porous media. More recently, several authors (Bagheri et al., 2023; Nazari et al., 2025) introduced pore-scale studies that explicitly account for hydrogen compressibility under reservoir conditions. Bagheri et al. (2023) used DNS to investigate the influence of capillary number, compressibility, and hysteresis on hydrogen displacement and trapping efficiency. Their simulations identified that misinterpreting hydrogen compressibility can lead to incorrect predictions of breakthrough dynamics and storage potential, especially in cyclic injection production scenarios. Nazari et al. (2025) on the other hand, incorporated wettability heterogeneity and compressibility effects in pore-scale DNS, focusing on how spatial distributions of wettability affect displacement paths and storage recovery.

While these recent contributions have highlighted the importance of compressibility in multiphase hydrogen transport, they primarily focused on drainage-imbibition cycles and heterogeneous geological formations. However, they do not provide a detailed analysis of local field variables such as pressure, velocity, and density variations and they do not offer a clear comparison between static and dynamic contact angles.

In addition to pore-scale modeling in saline aquifers, other storage strategies have been investigated. For example, Yu et al. (2023) reported that PAN-based activated carbon nanofibers can store hydrogen under ambient pressure, highlighting the potential of adsorption-based methods. Hosseini et al. (2024) analyzed secondary imbibition processes in carbonate and basalt reservoirs for H₂ and CO₂, showing the critical role of lithology on wettability and storage efficiency. These complementary studies underline the diversity of hydrogen storage approaches and motivate the need for pore-scale mechanistic investigations in geological systems.

In the present work, we focus on a pore-scale study on a single

idealized constricted pore geometry, with specific emphasis on:

* the influence of hydrogen compressibility on interface dynamics and storage efficiency, quantified by hydrogen pore occupancy;

* the onset of hydrogen bubble formation in the constrictions of the pores;

* the effect of injection conditions and geometric variations on pressure evolution, hydrogen breakthrough times, and fluid distributions at breakthrough;

* the role of static and dynamic contact angles on gas drainage efficiency and interface motion, with dynamic contact angle effects and wall slip conditions incorporated to improve the physical realism of the model beyond previous studies.

Furthermore, a comparison is made between compressible and incompressible hydrogen flows under identical conditions, providing a quantitative estimate of the error introduced by neglecting compressibility, an aspect that remains underexplored in previous studies. This nuanced contribution offers new insights into the transient pore-scale dynamics relevant to hydrogen storage, particularly for designing and validating larger-scale models such as pore network models. To achieve this, a series of DNS are performed using the open-source software OpenFOAM under varying injection boundary conditions and with different pore geometries.

## 2. Mathematical modelling

In the considered multiphase flow, which involves both compressible and incompressible fluids, the mass conservation equation is expressed as follows:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0 \qquad (1) $$

where $\rho(\text{kg}\cdot\text{m}^{-3})$ and $\mathbf{u}(\text{m}\cdot\text{s}^{-1})$ are the average density of the mixture and the velocity vector. The overall density is defined as:

$$ \rho = \alpha\rho_1 + (1 - \alpha)\rho_2 \qquad (2) $$

with $\alpha$, $\rho_1$ and $\rho_2$ being the volume fraction of hydrogen, the density of hydrogen and the density of brine, respectively. For the compressible gas phase (hydrogen), the density is related to the pressure and temperature through the ideal gas law:

$$ \rho_1 = \frac{P_1}{rT_1} \qquad (3) $$

where $r = 4124 \text{ J/(kg}\cdot\text{K)}$ is the specific gas constant for hydrogen, $P_1(\text{Pa})$ and $T_1(\text{K})$ are the pressure and temperature of the hydrogen.

The equation for the conservation of momentum considers the interaction between inertia, pressure, viscous stress and surface tension forces:

$$ \frac{\partial(\rho\mathbf{u})}{\partial t} + \nabla \cdot (\rho\mathbf{u}\mathbf{u}) = -\nabla P + \nabla \cdot (\mu\nabla\mathbf{u}) + F_{\sigma} \qquad (4) $$

where $\mu(\text{Pa}\cdot\text{s})$ and $F_{\sigma}(\text{N}/\text{m}^3)$ are respectively the dynamic viscosity and the surface tension force which is defined as follows (Brackbill et al., 1992):

$$ F_{\sigma} = \sigma\kappa\nabla\alpha \qquad (5) $$

where $\sigma(\text{N}\cdot\text{m}^{-1})$ and $\kappa(\text{m}^{-1})$ are the interfacial tension between both fluids and the curvature of the fluid-fluid interface, defined as:

$$ \kappa = \nabla \left( \frac{\nabla \alpha}{|\nabla \alpha|} \right) \qquad (6) $$

The average dynamic viscosity is given by the sum of the contributions from each phase:
$$ \mu = \alpha\mu_1 + (1 - \alpha)\mu_2 \qquad (7) $$

where $\mu_1$ and $\mu_2$ are the dynamic viscosities of the gas phase (hydrogen) and liquid phase (brine), respectively.

In this study, the Volume of Fluid (VOF) method is used to capture the evolution of the interface between the liquid and gas phases. The simulation domain consists of two subdomains (Fig. 1): one corresponding to the gas (hydrogen) and the other to the liquid (brine). The interface between these phases evolves according to the transport equation for the volume fraction of hydrogen:

$$ \frac{\partial \alpha}{\partial t} + \nabla \cdot (\alpha \mathbf{u}) + \nabla \cdot (\mathbf{u}_c \alpha (1 - \alpha)) = 0 \qquad (8) $$

where $\mathbf{u}_c(\mathbf{m} \cdot \mathbf{s}^{-1})$ is an artificial velocity, named compression velocity, introduced to refine the representation of the interface between phases. The term $\nabla \cdot (\mathbf{u}_c \alpha (1 - \alpha))$ is only active in the fluid–fluid interface region. This velocity has no direct physical significance but improves interface sharpness (Weller, HG, 2008), $\mathbf{u}_c$ is computed as:

$$ \mathbf{u}_c = \min(C_\alpha, \max(|\mathbf{u}|)) \frac{\nabla \alpha}{|\nabla \alpha|} \qquad (9) $$

where $C_\alpha$ is the compression factor, a key parameter that controls interface sharpness in the VOF method. While $C_\alpha$ can theoretically take any value greater than or equal to 0, it is often limited to a maximum value of 1 (Wardle and Weller, 2013). When $C_\alpha$ is set to 0, the interface compression is disabled, leading to a diffuse interface representation, which may cause numerical diffusion of the phases. In the solver implementation, $C_\alpha$ is defined and applied independently for each phase pair, ensuring a sharp interface for specific phases. However, this method can introduce numerical artifacts, particularly spurious currents at the interface, which are especially problematic for small-scale flows influenced by surface tension, such as capillary rise. These spurious currents can be minimized by maintaining $C_\alpha = 1$ (Cerne et al., 2002), as adopted in this study.

The equation for the conservation of energy accounts for the thermal exchanges between phases, the effects of compression, and the dissipation of viscous energy:

$$ \frac{\partial(\rho C_v T)}{\partial t} + \nabla \cdot (\rho \mathbf{u} C_v T) = -p \nabla \cdot \mathbf{u} + \nabla (\lambda \nabla T) + \Phi \qquad (10) $$

where $C_v(J \cdot kg^{-1} \cdot K^{-1})$ is the specific heat at constant volume, $\lambda(J \cdot m^{-1} \cdot K^{-1})$ is the thermal conductivity and $\Phi(J \cdot m^{-3})$ is the viscous dissipation term. The average thermal conductivity and specific heat at constant volume are given by the sum of the contributions from each phase:

$$ \lambda = \alpha \lambda_1 + (1 - \alpha) \lambda_2 \qquad (11) $$

$$ C_v = \alpha C_{v1} + (1 - \alpha) C_{v2} \qquad (12) $$

where $\lambda_1$ and $\lambda_2$ are the thermal conductivities of the hydrogen and brine, respectively and $C_{v1}$ and $C_{v2}$ are the heat capacities at constant volume of the hydrogen and brine, respectively.

The above equations are solved with the open-source code OpenFOAM which is based on Finite-Volume method. The governing equations (mass, momentum and energy, along with the volume fraction equations) are discretized using the combined PISO and SIMPLE algorithms, known as the PIMPLE algorithm, to ensure a robust and accurate solution of the conservation laws. For incompressible multiphase flow simulations, the interFoam (InterFoam, 2024) solver is employed, whereas for compressible multiphase flow scenarios, the compressibleInterFoam version of InterFoam solver is used to accurately account for variations in density, temperature, and pressure. The Courant-Friedrichs-Lewy (CFL) stability condition is crucial for ensuring that numerical simulations converge to the exact solution (Courant et al., 1928). The CFL number is defined as:

$$ \text{CFL} = \frac{|\mathbf{u}| \Delta t}{\Delta h} \qquad (13) $$

where $|\mathbf{u}|(\mathbf{m} \cdot \mathbf{s}^{-1})$, $\Delta t(s)$, $\Delta h(m)$ represent velocity magnitude, time step, and mesh size, respectively. For the time discretization, the first-order Euler discretization scheme is used. A strict limit on the time step is imposed by the maximum CFL condition $CFL_{max} = 1$ (Chen et al., 2013; Da Riva and Del Col, 2009) recommended setting $CFL_{max} = 0.25$ to ensure numerical stability. This recommendation is followed in the current simulations.

## 3. Computational setup and test cases

### 3.1. Geometry

In this study, the considered geometries are cylindrical pores with a diameter of $D_p = 20 \ \mu\text{m}$ and a length of 110 $\ \mu\text{m}$ (see Fig. 1). These pores may include one or more constrictions, representing, in a simplified manner, the variations in flow cross-section found in granular media, such as the interstices of a sphere packing. The presence of these constrictions allows the study of the influence of section variations on the flow and thermal transfers. At the initial state, the first 10 $\ \mu\text{m}$ of the domain, in the axial direction, are filled with hydrogen ($\alpha = 1$). Three geometric configurations are examined, as shown in Fig. 1:

* **Geometry 1:** A cylindrical pore with a constant cross-section.

* **Geometry 2:** A cylindrical pore with a single constriction.

* **Geometry 3:** A cylindrical pore with two separated constrictions.

For Geometry 2, three different pore-to-constriction ratios are also investigated, leading to a total of five defined variants, as detailed in Table 1: a variant with constant cross-section (V1), a variant with a constriction of diameter 0.6 $D_p$ (V2), a variant with a constriction of diameter 0.4 $D_p$ (V3), a variant with a constriction of diameter 0.3 $D_p$ (V4), and a variant with two constrictions of diameter 0.6 $D_p$ (V5).

This study uses a simplified axisymmetric pore-throat geometry with one or two constrictions to model a porous medium element. While many previous pore-scale studies have used 2D models (e.g., Derijani et al., 2023) to isolate displacement mechanisms, our axisymmetric

<figure><img src="image_4.png" /><figcaption>Fig. 1. Different geometries investigated in the current simulations: cylindrical pore with constant cross-section (geometry 1), cylindrical pore with one constriction (geometry 2) and cylindrical pore with two constrictions (geometry 3). Red color represents hydrogen (gas) and blue color represents brine (liquid).</figcaption></figure>
<table><caption>**Table 1**

Geometry variants used for this study.</caption><thead><tr><th>Variants</th><th>L1</th><th>L2</th><th>L3</th><th>C1</th><th>C2</th></tr></thead><tbody><tr><td>V1 : Geometry1</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>V2 : Geometry2</td><td>50 μm</td><td>—</td><td>—</td><td>10 μm</td><td>4 μm</td></tr><tr><td>V3 : Geometry2</td><td>50 μm</td><td>—</td><td>—</td><td>10 μm</td><td>6 μm</td></tr><tr><td>V4 : Geometry2</td><td>50 μm</td><td>—</td><td>—</td><td>10 μm</td><td>7 μm</td></tr><tr><td>V5 : Geometry3</td><td>—</td><td>20 μm</td><td>30 μm</td><td>10 μm</td><td>4 μm</td></tr></tbody></table>

setup better represents radial injection conditions while maintaining numerical stability and interpretability. The solid matrix is assumed rigid, focusing the analysis on hydrogen compressibility, pore geometry, and wettability. The effect of solid deformation on flow dynamics is not considered here. This aspect has, however, been examined in other contexts, for example in a numerical study of single-phase flow through 2D individual pores, where a detailed mathematical model was developed (Rodríguez de Castro et al., 2023a).

This work aims to investigate the influence of compressibility and contact angle, under a defined slip velocity, in a model porous medium representative of reservoir length scales and pore-to-throat size ratios. The pore-to-throat size ratios considered in the investigated geometries range from 1 to 10/3, encompassing typical values observed in many sandstone and carbonate reservoirs (Zhao et al., 2024), which also applies to the pore diameters used in this work. Geometry 1 represents a capillary with a constant cross section, as employed in traditional bundle-of-capillaries models for porous media, whereas Geometry 2 represents the shape of individual pores in pore-network modeling approaches (Gostick et al., 2019; Rodríguez de Castro et al., 2021, 2023b). Previous work has examined the effect of pore shape on the accuracy of pore-network modeling predictions, in comparison with experiments in model porous media (Rodríguez de Castro et al., 2021), as well as the applicability of ball-and-stick models in pore-network simulations for different pore-to-throat size ratios (Dai et al., 2025). Furthermore, pore-network modeling approaches have been applied to study brine-H₂ flows in geological porous media (Hashemi et al., 2021; Zhao et al., 2024). The effects of pore interconnectivity have not been specifically addressed in the present study. Geometry 3 represents a pore with a coordination number of 2. The coordination number in a pore network represents the average number of pore bodies connected to a specific pore (Sahimi, 2012). A coordination number of 2 is representative of some types of carbonate reservoirs (Hashemi et al., 2021) but not of sandstone reservoirs.

## 3.2. Slip condition and mesh

A two-dimensional (2D) axisymmetric model is used (see Fig. 1). A partial slip condition is imposed at the fluid-wall interface, allowing a non-zero tangential velocity at the wall while maintaining zero normal velocity. The tangential slip velocity can be expressed either by using a slip coefficient multiplied by the shear rate of the fluid at the wall or through the concept of a slip length, as proposed by Wörner (2012). The latter approach, used in this work, is given by the following equation:

$$u_0 = l_s \frac{\partial u_z}{\partial r} \qquad (13)$$

where $\partial u_z/\partial r$, $u_0$ and $l_s$ represent the radial velocity gradient, the velocity at the wall, and the slip length, respectively. The slip length $l_s$ is defined as the distance between the solid wall and the point where the velocity vanishes, by extrapolating the hydrodynamic flow inside the solid. The slip velocity $u_0$ can be expressed in terms of the velocity of the adjacent cell $u_1$ as:

$$u_0 = l_s \frac{u_1 - u_0}{h/2} \qquad (14)$$

where $h$ is the distance of the adjacent cell. By expressing $u_0$ in terms of $u_1$, the Eq. (14) becomes:

$$u_0 = \frac{2l_s}{h + 2l_s} u_1 \qquad (15)$$

The relationship between velocities is defined in terms of the partial slip coefficient $f$, defined as:

$$u_0 = (1 - f)u_1 \qquad (16)$$

Thus, by analogy, the partial slip coefficient $f$ is defined as:

$$f = \frac{h}{h + 2l_s} \qquad (17)$$

A value of $l_s = 0$ leads to $f = 1$, and consequently $u_0 = 0$ (no-slip condition). Conversely, the value of $l_s = +\infty$ leads to $f = 0$, resulting in $u_0 = u_1$ (full-slip condition). For a partial slip condition, the value of $f$ lies between 0 and 1, with $f$ related to the slip length $l_s$ and the size of the adjacent cell $h$. Since $h$ is a characteristic of the mesh, the choice of $l_s$ is crucial. In most applications, the physical slip length is much smaller than the realistic mesh size for numerical simulations. Renardy et al. (2001) introduced an artificial slip length, with a value on the order of the mesh size, and this approach is adopted in the present study. Similar assumptions have been employed in subsequent studies as a practical way to regularize the stress singularity at the moving contact line (Malgarinos et al., 2014; Qian et al., 2005). Studies on slip length are numerous and diverse. For example, Lauga et al. (2007) conducted a seminal review of slip lengths, discussing their physical values and conditions for deviation from the no-slip boundary. They summarized experimental methods and results, presented insights from molecular dynamics simulations, and interpreted findings through both molecular and continuum models. According to these authors the factors affecting slip length include wetting conditions, shear rate, pressure, surface charge, surface roughness and dissolved gas. The authors demonstrated that apparent slip lengths, reported across numerous experimental studies encompassing a wide range of surface wettabilities, vary over several orders of magnitude, from molecular scales up to hundreds of nanometers. Panton et al. (2007) stated that the slip length is on the order of 100 nm, which consistent with the slip length of 277 nm (cell size) taken in our simulations. While the actual slip length has to be determined experimentally, many numerical studies have adopted $l_s \approx h$ in the absence of precise experimental data (Renardy et al., 2001; Qian et al., 2006; Malgarinos et al., 2014) with $h$ representing the approximate mesh size at the first point near the wall (in our case, $h = 277$ nm).

At the wall, there is also a wettability condition defined by the contact angle $\theta$. The fundamental law of wetting (Young's law), expresses the static contact angle $\theta_s$ that the liquid edge makes when it rests on a solid surface. This angle is given by the surface tension involved in the problem: solid-liquid tension ($\sigma_{sl}$), solid-vapor tension ($\sigma_{sv}$), and liquid-vapor tension ($\sigma_{lv}$), according to the Young's law (Laplace, 1805), which assumes equilibrium at the triple point. To describe the behavior of the moving contact line, several studies (De Gennes et al., 1990; Shi et al., 2020; Malgarinos et al., 2014; Hashemi et al., 2021) have shown that when the triple point moves at a certain velocity, the value of $\theta$ deviates from $\theta_s$. The dependence of the dynamic contact angle, named $\theta_d$, on $\theta_s$ has been described by the well-known Cox-Voinov's theory for low contact angles. A line force then exists that moves the triple line, and it is given by the following equation:

$$F = \sigma_{lv}(\cos\theta_d - \cos\theta_s) \qquad (18)$$

where $\theta_d$ is related to the capillary number $C_a$ by the Cox-Voinov model (Voinov, 1977):

$$\theta_d^3 = \theta_s^3 + 144C_a \qquad (19)$$

The capillary number ($C_a$) is defined in terms of the slip velocity:

$$C_a = \frac{\mu_2 u_0}{\sigma} \qquad (20)$$

where $\mu_2$ is the dynamic viscosity of the liquid.
The parameters in the dynamic contact angle model (Eq. 19) were chosen based on literature values, focusing on (Gennes et al., 1990; Malgarinos et al., 2014; Voinov, 1977; Wörner, 2012), and experimental trends summarized by Malgarinos et al. (2014), ensuring a reasonable representation of dynamic wetting behavior.

The computational domain with the mesh used in this study is shown in Fig. 2 (the example corresponds to geometry 2). The domain is modeled as a wedge and discretized with a uniform mesh consisting of 23,850 cells, with a cell size of $h = 0.277$ µm. Key boundaries are defined, including the inlet, the outlet, the walls and the axis of symmetry.

### 3.3. Fluid properties, boundary and initial conditions

The boundary conditions applied for velocity, pressure, and temperature are detailed in Table 2 for both incompressible and compressible cases. For the compressible case, a mass flow rate is imposed at the inlet to ensure controlled hydrogen injection. Imposing pressure instead would have required exceeding the capillary pressure to displace brine; otherwise, due to wettability, the flow could reverse. The mass flow rate condition thus offers better control for this scenario. To reflect realistic conditions, the injection flow rate is selected based on the INERIS Group report (INERIS, 2021), which outlines operating parameters during hydrogen storage. The report provides the total injection flow rate for the devised application. Since the current simulations are at the pore scale, a calculation is performed to estimate the flow rate through a single pore, based on pore size and medium porosity. To ensure a consistent comparison between the compressible and incompressible regimes, the mass flow rate used in the compressible case is converted into an equivalent volumetric flow rate for the incompressible case considering a constant fluid density. The choice of atmospheric pressure (1 bar) as the outlet pressure ($P_{outlet}$) is deliberate, as it reflects typical low-pressure laboratory conditions—such as micromodel (Lysyy et al., 2023, 2022) or core flooding experiments using atmospheric backpressure regulator—and near-wellbore pressures during fluid withdrawal. It also represents the lower bound of pressure ranges reported in previous studies (Ali et al., 2021; Iglauer et al., 2021). This choice enables isolation of the effects studied here from those induced by the higher pressures characteristic of subsurface storage system. It is also noted that previous experimental research showed that hydrogen saturation in the pores after drainage was independent of the pressure within the 1–40 bar range (Lysyy et al., 2023). The fluid properties of hydrogen and brine are shown in Table 3.

The brine-H₂ interface tension was taken as $\sigma = 0.051$ N/m, and the equilibrium contact angle as $\theta_s = 60^{\circ}$ (Hashemi et al., 2021).

<table><caption>**Table 2**
Boundary conditions.</caption><thead><tr><th></th><th>Flow rate</th><th>Pressure</th><th>Temperature</th><th>α</th></tr></thead><tbody><tr><td>Inlet</td><td>Compressible<br>7.92 × 10<sup>-11</sup>kg/m<sup>3</sup><br>Incompressible<br>8.80 × 10<sup>-10</sup>m<sup>3</sup>/s</td><td>ZeroGradient</td><td>Compressible<br>300 K<br>Incompressible<br>N/A</td><td>1</td></tr><tr><td>Outlet</td><td>ZeroGradient</td><td>Compressible<br>$P_{outlet}$ = 10<sup>5</sup> Pa<br>Incompressible<br>$P_{outlet}$ = 0</td><td>Compressible<br>ZeroGradient<br>Incompressible<br>N/A</td><td>ZeroGradient</td></tr><tr><td>wall</td><td>partialSlip</td><td>ZeroGradient</td><td>ZeroGradient</td><td>$\theta_d$</td></tr></tbody></table>

### 3.4. Test case: Air compression

A test case (Caltagirone et al., 2011) with an analytical solution is used to validate the numerical results obtained in the present work. The simplified 2D geometry illustrated in Fig. 3a is used in the test case. The compression chamber has a height of 500 mm and a wall-to-wall distance of 30 mm. The mass of air $m_0$ is kept constant, while its volume gradually decreases due to the injection of water at a constant velocity $v_0 = 0.05$ m/s. The total simulated time is 5 s. Simulations are conducted using different mesh resolutions, with a focus on tracking the evolution of air density, temperature, and pressure over time. Simulations are performed using three different meshes: 60 × 10, 100 × 30, and 300 × 70, with four-time steps $\Delta t$: 0.05 s, 0.01 s, 0.008 s, and 0.005 s. The effects of gravity and surface tension are neglected. The aim is to verify mesh convergence and to evaluate discrepancies between numerical and analytical results. In our model, air is treated as an ideal gas undergoing adiabatic compression, while water is assumed to be incompressible. The analytical expressions for air density, pressure, and temperature as a function of time are given by the following equations:

$$ \rho(t) = \frac{m_0}{V(t)} = \frac{\rho_0}{1 - \frac{v_0 t}{L_0}} \qquad (21) $$

$$ P(t) = P_0 \left( \frac{1}{1 - \frac{v_0 t}{L_0}} \right)^{\gamma} \qquad (22) $$

$$ T(t) = T_0 \left( \frac{1}{1 - \frac{v_0 t}{L_0}} \right)^{1-\gamma} \qquad (23) $$

where $P_0 = 10^5 Pa$ is the initial pressure, $\rho_0 = 1.2 Kg/m^3$ is the initial

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Geometry 2: Mesh of the computational domain for 2D axisymmetric. The close-up views show the mesh at the start (a) middle (b) and end (c) with boundary conditions applied to the domain.</figcaption></figure>
<table><caption>**Table 3**
Fluid properties.</caption><thead><tr><th>Fluids</th><th>Compressibility</th><th>Dynamic Viscosity (Pa.s)</th><th>Density (Kg/m³)</th><th>C<sub>V</sub>-(J/(Kg.K))</th><th>Molar Mass (g/mol)</th></tr></thead><tbody><tr><td>Hydrogen</td><td>Compressible<br>or Incompressible</td><td>8.76 × 10<sup>-6</sup></td><td>0.09(Incompressible case)</td><td>5060</td><td>2.06</td></tr><tr><td>Brine</td><td>Incompressible</td><td>1.20 × 10<sup>-3</sup></td><td>1000.5</td><td>3850</td><td>58.44</td></tr><tr><td>Air</td><td>Compressible</td><td>1.84 × 10<sup>-5</sup></td><td>1.2</td><td>718</td><td>28.96</td></tr><tr><td>Water</td><td>Incompressible</td><td>10<sup>-3</sup></td><td>1000</td><td>4195</td><td>18</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Panel (a) shows a schematic of the considered geometry and boundary conditions. Panels (b) to (d) display the relative errors in density, temperature, and pressure, respectively, plotted as a function of the inverse time step $1/\Delta t$, for three mesh resolutions: 60 × 10, 100 × 30, and 300 × 70.</figcaption></figure>

density of air, $T_0 = 293K$ is the initial temperature and $\gamma$ is the ratio of specific heats.Fsigs. 3b, 3c and 3d, respectively present the relative error $||E||_\phi$ where $\phi$ represents density, temperature, and pressure as a function of the inverse time step $1/\Delta t$, for the different mesh resolutions. The error $||E||_\phi$ is defined as:

$$ ||E||_{\phi}(t) = 100\max \left| \frac{\phi_{N}(t) - \phi_{A}(t)}{\phi_{A}} \right| \qquad (24) $$

where $\phi_N$ is the numerical solution and $\phi_A$ is the analytical solution. This analysis allows assessing the impact of spatial and temporal refinement on numerical accuracy. The convergence analysis (Fig. 3) shows a consistent decrease in the $||E||_\phi$ error as the mesh transitions from coarse (60 × 10) to fine (300 × 70), and as the time step is decreased (i.e., increasing $1/\Delta t$). For the finest mesh, temporal discretization errors become negligible, indicating that spatial resolution becomes the dominant factor influencing accuracy. The results demonstrate that CompressibleInterFoam exhibits good numerical convergence and accurately captures air compression, confirming the importance of mesh refinement and small time steps for reliable compressible flow simulations.

## 4. Results

### 4.1. Immiscible displacement of brine by hydrogen for the compressible case

Fig. 4 shows the simulated phase distribution over time and the corresponding evolution of inlet pressure $P_i$ for the V2 geometry defined

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Evolution of inlet pressure $P_i$ and phase distribution over time for the V2 geometry (The red color represents compressible hydrogen, while the blue color represents brine).</figcaption></figure>
in subsection 3.1 and the considered boundary conditions (subsection 3.3). Due to the compressibility of hydrogen, the interface between hydrogen and brine remains nearly stationary during the initial stages of injection. As a result, hydrogen accumulates within an almost constant volume, leading to a continuous increase in its density and a corresponding rise in pressure (see Fig. 4, for $t < 7.45 \times 10^{-6}$ s). For $7.45 \times 10^{-6}$ s < $t < 3 \times 10^{-5}$ s, the interface moves as if in a conduit with no cross-sectional variation. Beyond this point, the constriction affects the flow, accelerating the interface due to the reduced cross-section. At $t = 4 \times 10^{-5}$ s, the interface surpasses the midpoint of the conduit, reaching its maximum velocity. The interface reaches the outlet of the conduit, marking the breakthrough time at $t = 5.55 \times 10^{-5}$ s.

Shortly after breakthrough, the first gas bubble forms at $t = 5.95 \times 10^{-5}$ s. This bubble formation is caused by the increase in flow velocity within the constricted section leading to a pressure drop, which in turn promotes hydrogen bubble formation. Previous studies by Cai et al. (2017) highlighted the significant impact of the pore-to-constriction diameter ratio on bubble formation in incompressible flows. Their results showed that higher diameter ratios facilitate bubble formation, while ratios close to one tend to suppress it, especially at low injection velocities. However, their study did not address the behavior of compressible fluids.

At $t = 9 \times 10^{-5}$ s, the hydrogen phase once again reaches the conduit outlet, followed by the formation of another hydrogen bubble at $t = 9.33 \times 10^{-5}$ s. This cycle repeats periodically, with a frequency of one bubble every $3.46 \times 10^{-5}$ s. The effect of bubble formation on inlet pressure evolution can also be observed in Fig. 4, where a marked drop in pressure occurs when a continuous hydrogen film forms between the inlet and the outlet of the pore, followed by a sharp increase after bubble formation (e.g., between $5.86$ and $6 \times 10^{-5}$s).

For the considered boundary conditions (see subsection 3.3) and geometrical variant V2, the curve representing inlet pressure $P_i$ as a function of time can be divided into three distinct phases, as shown in Fig. 4:

* A rising phase from 0 to $7.45 \times 10^{-6}$s,

* A declining phase from $7.45 \times 10^{-6}$ to $5.86 \times 10^{-5}$ s,

* And a fluctuation phase for $t > 5.86 \times 10^{-5}$ s.

During the rising phase, $P_i$ gradually increases, as the movement of the fluid-fluid interface is slow and hydrogen accumulates at the inlet. This also leads to an increase in hydrogen density. Once $P_i$ reaches a threshold value sufficient to overcome both the viscous resistance to brine flow and the capillary pressure, the interface begins to move significantly. This is explained by the time required for the brine to accelerate, during which more and more hydrogen is injected, causing the pressure to increase.

In the declining phase, $P_i$ decreases, and the interface advances more rapidly. This phase starts when the displacement flow rate of the brine matches the imposed hydrogen flow rate. The drop in $P_i$ occurs because the brine is displaced by a much less viscous fluid. Within this phase, two distinct behaviors are observed: (i) a steady decrease with a constant slope before the hydrogen phase reaches the constriction at $t = 2.35 \times 10^{-5}$ s, and (ii) a more gradual decrease once the hydrogen reaches the constriction, where the slope becomes significantly lower due to the geometric narrowing.

The fluctuation phase marks the onset of hydrogen bubble formation. The formation of these bubbles increases the value of $P_i$, as it disrupts the continuity of the hydrogen phase. This phenomenon repeats continuously over time. This suggests that it is possible to predict the formation

<figure><img src="image_8.png" /><figcaption>Fig. 5. (a) Velocity magnitude, (b) pressure, and (c) density fields at $t = 4.7 \times 10^{-5}$ s (d-f) the corresponding profiles along the central horizontal axis of the conduit as a function of the axial position. These results correspond to the compressible case with hydrogen injection in the V2 geometry.</figcaption></figure>
of hydrogen bubbles by measuring the evolution of $P_i$ (cf. Fig. 4). These observations are absent in the incompressible case, as well as in studies that do not account for compressibility effects, where the hydrogen inlet pressure progressively decreases over time (see section 4.2).

Fig. 5-a-c shows an example of the velocity magnitude, pressure, and density maps computed with the current numerical simulations. At the time considered in this figure ($t = 4.7 \times 10^{-5}$ s), the hydrogen phase has already passed through the constriction. The maximum velocity is reached at the constriction, as expected. Density variations (Fig. 5c) further illustrate the evolution of hydrogen density through the constriction, where a noticeable decrease is observed. This behavior confirms that hydrogen undergoes compression and expansion as it flows through the geometry. Fig. 5-d and 5-e represent the velocity magnitude and pressure along the central horizontal axis of the conduit as a function of the axial position at the same time $t = 4.7 \times 10^{-5}$ s. A minimum velocity is observed at the hydrogen/brine interface ($x \approx 81 \mu m$), indicating an interaction between the two phases. It can be noted that, while the brine velocity remains roughly constant outside the constriction, the hydrogen velocity varies, emphasizing the behavioral differences between an incompressible and a compressible fluid. Moreover, the pressure decreases linearly along the pore (Fig. 5-e), with a steeper gradient in the brine phase, which can be attributed to its higher viscosity. A sharp pressure drop is observed at the interface position due to capillary pressure. Additionally, the hydrogen phase exhibits noticeable density variation along the flow direction, as shown in Fig. 5-f, highlighting the effect of gas compressibility within the constricted geometry.

## 4.2. Comparison between compressible and incompressible cases

### 4.2.1. Hydrogen phase advancement in a drainage process

Using the inlet conditions from Table 2, Fig. 6 shows a marked difference in immiscible displacement process between compressible and

incompressible cases. Initially, hydrogen occupies the same position in both cases. However, at $t = 10^{-5}$s, the interface advances further in the incompressible case (Fig. 6a'-g') than in the compressible one (Fig. 6a-g). Under constant inlet mass flow, hydrogen accumulates in the compressible case due to an inlet pressure $P_i$ lower than that required to displace brine at the imposed flow rate. As a result, the pressure field must settle before the fluid flows. In contrast, the incompressible case shows immediate interface movement, leading to a faster displacement. At $t = 3.03 \times 10^{-5}$ s, a bubble forms in the incompressible case before the hydrogen phase reaches the outlet (see Fig. 6f'). This bubble later merges with the hydrogen phase at $t = 3.03 \times 10^{-5}$ s, as shown in Fig. 6g'. At $t = 4.02 \times 10^{-5}$ s, the hydrogen phase reaches the outlet in the incompressible case (cf. Fig. 6h'), whereas in the compressible case, it has only traveled ¾ of the conduit (see Fig. 6h). The breakthrough time was $4.02 \times 10^{-5}$ s for the incompressible case compared to $5.55 \times 10^{-5}$ s for the compressible case, corresponding to a relative underestimation of 27.6 % when compressibility is neglected. This demonstrates that neglecting hydrogen compressibility can obscure key physical phenomena such as the advancement of hydrogen and the formation of hydrogen bubbles in the conduit. In the compressible case, pressure variations lead to the formation of multiple bubbles (as discussed in section 4.1), whereas only a single bubble appears in the incompressible case under identical conditions (cf. Fig. 6f'). Indeed, generating multiple bubbles in the incompressible case would require a significantly lower inlet flow rate, which may not be representative of realistic operating conditions. This comparison shows that the compressible hydrogen model better reflects the real behavior, capturing pressure build-up, delayed interface movement, and bubble formation which correspond to key aspects in applications such as drainage and reservoir storage. These single-pore observations are consistent with the findings of Bagheri et al. (2023), who reported bubble formation exclusively in the compressible case within the sandstone-like 2D pore network used in their simulations. Additionally, in the incompressible case, the

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Evolution of volume fraction over time: on the left, compressible hydrogen; on the right, incompressible hydrogen. Red represents hydrogen and blue represents brine.</figcaption></figure>
movement of the triple line along the conduit walls is slower as compared to the compressible case. This is because incompressible fluids do not allow local volume reduction, which would otherwise create space for interface movement. In contrast, the compressible case adapts better to pressure changes, enabling smoother interface displacement.

### 4.2.2. Comparison of pore occupancy by hydrogen in compressible and incompressible cases

By comparing the results obtained for V2, Fig. 7 shows that incompressible hydrogen reaches the outlet of the conduit more rapidly, with a breakthrough time of $4.02 \times 10^{-5}$ s, compared to $5.55 \times 10^{-5}$ s for the compressible hydrogen. This breakthrough time corresponds to the moment when a continuous hydrogen stream from the inlet reaches the outlet. As explained in the previous subsection, this behavior is explained by the fact that, in the compressible case, the gradual buildup of pressure compresses the hydrogen before it can effectively displace the brine, resulting in a delayed breakthrough time. In contrast, in the incompressible case, hydrogen displaces the brine immediately without undergoing compression, leading to a faster advancement. Another key difference is seen in the fraction of the pore volume occupied by hydrogen at breakthrough time, denoted as $S_{H2, BT}$.

$$S_{H2, BT} = \frac{V_{H2, BT}}{V_p} \qquad (25)$$

where $V_{H2, BT}$ is the volume occupied by hydrogen at breakthrough time and $V_p$ is the total pore volume.

In the compressible case, brine displacement is more efficient, obtaining $S_{H2, BT} = 0.903$ compared to $S_{H2, BT} = 0.742$ in the incompressible case. Therefore, omitting hydrogen compressibility results in a relative underestimation of 17.8 % of the pore storage capacity.

Table 4 lists the position of the advancing front, the values of hydrogen saturation, and $P_i$ for the compressible and incompressible cases at all times for the example shown in Fig. 6. The relative error obtained when neglecting compressibility for each parameter ($P_i$, hydrogen saturation, and the position of the advancing front) at each time step are also included.

### 4.3. Influence of injection conditions

The INERIS group report (INERIS, 2021), highlights injection flow rate as a key factor influencing hydrogen transport in porous media. For this reason, a set of simulations using different injection flow rates $7.92 \times 10^{-11}$, $6.33 \times 10^{-11}$, $4.312 \times 10^{-11}$, and $1.86 \times 10^{-11}$ kg/s is conducted to better understand their effect on phase distribution at breakthrough time in the compressible case. Fig. 8 illustrates the relationship between the injected flow rate and $S_{H2,BT}$ for three different geometrical configurations (V2, V3, and V4). The results reveal a clear trend: as the injection flow rate increases, the proportion of stored hydrogen decreases across all configurations. This occurs because, at lower flow rates, hydrogen conforms more effectively to the corners of the geometry. Once a brine pocket becomes surrounded by hydrogen, it remains trapped. Similar to viscous fingering, once hydrogen advances along a given path, it tends to keep penetrating it due to its lower viscosity compared to brine. This effect becomes more pronounced at higher flow rates.

Among the tested configurations, V2 consistently demonstrates the

highest $S_{H2, BT}$ value, followed by V3 and then V4. This can be explained by the geometry of the constrictions: V2 has the smallest grain size, resulting in a less pronounced constriction, which reduces brine trapping. In contrast, V4 exhibits the lowest storage efficiency due to its larger grain size, creating narrower pore throats that promote faster breakthrough time and reduce the volume of hydrogen retained. These findings highlight the critical influence of pore geometry and injection rate on hydrogen retention, emphasizing the need to optimize these parameters to enhance the performance of underground hydrogen storage systems.

A low injection flow rate reduces the maximum hydrogen pressure, thereby minimizing the risks of overpressure and damage to porous structures. Moderate-pressure storage also reduces hydrogen flammability risks, a particularly critical aspect for industrial and energy applications. Moreover, the current results indicate that for a low injection flow rate (see Fig. 9), $S_{H2, BT}$ reaches 0.960, compared to lower values observed at higher flow rates. These single-pore observations align with the 2D pore network results of Bagheri et al. (2023) and Yan et al. (2025), which also show higher hydrogen saturations during drainage. From a technical standpoint, these results demonstrate that adjusting the injection flow rate not only enhances storage safety but also affects hydrogen pore occupancy, thereby influencing storage efficiency.

### 4.4. Influence of the geometry

The pressure field is significantly influenced by the geometry, particularly by the size of the constriction under the conditions defined in Table 2. The narrower the constriction, the higher the inlet pressure (see Fig. 10). In the case of a geometry without constriction (geometry V1), the pressure behavior indicates the absence of hydrogen bubble formation and the pressure decreases monotonically along the conduit.

For the V4 geometry, the inlet pressure over time is significantly higher compared to the other configurations with a single constriction. This indicates that upstream of the constriction, the pressure builds up substantially, while downstream the pressure remains relatively low. In such a scenario, the advancement of the brine could potentially be impeded if a pressure gradient were applied instead of a mass flow rate. Indeed, for the gas phase to overcome the constriction and displace the brine, the applied pressure must exceed the capillary pressure and viscous stresses threshold imposed by the narrower pore throat. Pressure is not only influenced by the size of the constrictions but also by their number in the conduit. As shown in Fig. 11 for V5, an increase in the number of constrictions increases the pressure losses, which in turn requires a higher pressure to overcome the flow resistance, as also demonstrated by Althoey et al. (2018), in order to maintain a constant flow rate. Fewer constrictions mean lower pressure, while a higher number of constrictions leads to higher pressure.

The efficiency of storage, assessed through the proportion of occupied space and breakthrough time, is strongly influenced by the geometry of the porous medium and the size of the constrictions, as shown in Fig. 11. An inverse correlation is observed between these parameters: the narrower the constrictions, the lower the proportion of stored hydrogen. Additionally, the breakthrough time, which corresponds to the time required to reach saturation, is also affected (see Fig. 11). Porous media with larger constrictions saturate more slowly than those with narrower passages due to increased flow resistance, which slows fluid progression. In an industrial context, these insights are key to

<figure><img src="image_10.png" /><figcaption>Fig. 7. Hydrogen-filled pore volume percentage at breakthrough time: a) Compressible case; b) incompressible case for geometry V2. Red represents hydrogen and blue represents brine.</figcaption></figure>
<caption>**Table 4**</caption>

Position of the advancing front, the values of hydrogen saturation, and $P_i$ for compressible and incompressible cases at all times for the example in Fig. 6. The relative error associated with neglecting compressibility is also reported for each parameter.

<table><thead><tr><th rowspan="2">Time (s)</th><th colspan="3">Maximum x-coordinate position of the fluid-fluid interface (m)</th><th colspan="3">Hydrogen saturation</th><th colspan="3">P<sub>i</sub>(Pa)</th></tr><tr><th>Comp.</th><th>Incomp.</th><th>Error (%)</th><th>Comp.</th><th>Incomp.</th><th>Error (%)</th><th>Comp.</th><th>Incomp.</th><th>Error (%)</th></tr></thead><tbody><tr><td>10<sup>-5</sup></td><td>1.7 × 10<sup>-5</sup></td><td>3.2 × 10<sup>-5</sup></td><td>88.2</td><td>0.150</td><td>0.248</td><td>65.3</td><td>167,925</td><td>278,101</td><td>65.6</td></tr><tr><td>1.5 × 10<sup>-5</sup></td><td>2.5 × 10<sup>-5</sup></td><td>4.3 × 10<sup>-5</sup></td><td>72.0</td><td>0.223</td><td>0.325</td><td>45.7</td><td>147,059</td><td>255,561</td><td>73.8</td></tr><tr><td>2 × 10<sup>-5</sup></td><td>3.5 × 10<sup>-5</sup></td><td>6.5 × 10<sup>-5</sup></td><td>85.7</td><td>0.302</td><td>0.461</td><td>52.6</td><td>133,347</td><td>79,498</td><td>40.4</td></tr><tr><td>3.35 × 10<sup>-5</sup></td><td>6.3 × 10<sup>-5</sup></td><td>9.6 × 10<sup>-5</sup></td><td>52.3</td><td>0.508</td><td>0.640</td><td>26.0</td><td>124,626</td><td>65,438</td><td>47.4</td></tr><tr><td>4.02 × 10<sup>-5</sup></td><td>8.2 × 10<sup>-5</sup></td><td>11 × 10<sup>-5</sup></td><td>34.1</td><td>0.652</td><td>0.742</td><td>13.8</td><td>120,556</td><td>19,100</td><td>84.1</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Proportion of stored hydrogen $S_{H2, BT}$ as a function of injected flow rate for different geometries.</figcaption></figure>

optimizing hydrogen storage, where the permeability and porosity are critical. Optimizing the operating parameters by accounting for the constrictions.

### 4.5. Comparison between static and dynamic contact angles

The variation of inlet pressure over time is analyzed for the compressible case by comparing static (60°) and dynamic contact angle. Fig. 12 shows that the overall system evolution is similar in both configurations; however, $P_i$ is slightly higher in the case of a static contact angle as compared to a dynamic one. This difference can be explained by the fact that the dynamic angle allows for better adaptation of the fluid-fluid interface to local variations in capillary forces, facilitating the movement of the interface and reducing constraints on the fluid.

It is also observed that, for a dynamic contact angle, the hydrogen phase progresses slightly faster than in the static angle case (see subset in Fig. 12). This higher velocity can be attributed to a decrease in capillary resistance, making the flow less constrained and enabling smoother hydrogen advancement through the conduit. Furthermore, the behavior of the triple point (the interface between hydrogen, brine, and the wall) differs depending on the nature of the contact angle. With a static contact angle, the transition between the hydrogen and brine phases is more constrained, often leading to local pressure buildup and a slower advancement of the hydrogen front. In contrast, the use of a dynamic contact angle facilitates a more adaptive gas-liquid interface, promoting better spreading of the hydrogen phase and attenuating peak pressures. Additionally, the breakthrough time defined as the moment the hydrogen reaches the outlet is slightly shorter in the dynamic case (5.55

$\times 10^{-5}$s) compared to the static case ($5.60 \times 10^{-5}$s), reflecting improved hydrogen mobility. Notably, the proportion of pore volume occupied by hydrogen at breakthrough time is marginally higher in the dynamic case ($S_{H2, BT} = 0.903$) than in the static case ($S_{H2, BT} = 0.899$), suggesting enhanced displacement efficiency when dynamic wetting conditions are considered.

Table 5 presents the temporal variation of the advancing dynamic contact angle for the case analyzed in Fig. 12. The values listed in Table 5 are comparable to the experimental dynamic advancing contact angle taken as $\theta = 58°-85°$ by Hashemi et al. (2021) for a brine-H$_2$ system.

### 5. Conclusions

This study provides a detailed numerical investigation of hydrogen flow in brine-saturated model pores using OpenFOAM, with the goal of understanding the relationship between operating conditions for underground hydrogen storage and process efficiency. A key contribution lies in the explicit consideration of hydrogen compressibility an often-overlooked factor in pore-scale simulations which enables more accurate predictions of pressure and density variations. Neglecting compressibility can lead to underestimating pressure buildup and misrepresenting displacement dynamics, potentially compromising the safety and reliability of storage system design. This study also highlights the importance of incorporating the dynamic contact angle to accurately capture the interface motion between hydrogen and brine. Unlike a constant angle, the dynamic contact angle reflects wettability changes with flow conditions, affecting capillary pressure, pressure values, interface stability, and the rate of front propagation and the formation of bubbles. Including this factor enhances the physical realism of

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Inlet pressures $P_i$ over time for different geometries.</figcaption></figure>

$$ (a) \quad t = 1.5 \times 10^{-5} s, \quad S_{H2, BT} = 0.903 $$

$$ (b) \quad t = 2.8 \times 10^{-5} s, \quad S_{H2, BT} = 0.960 $$

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

<figure><img src="image_14.png" /><figcaption>**Fig. 9.** Storage ratio and breakthrough time for two injection conditions: (a) $7.92 \times 10^{-11}$ kg/s and (b) $1.89 \times 10^{-11}$ kg/s.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 11.** Breakthrough time and residual brine saturation for the compressible case in geometries V3, V4, and V5. Red represents H₂, and blue represents brine.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 12.** Comparison of Pᵢ between static and dynamic contact angles. An inset shows details of the inlet pressure field.</figcaption></figure>

<table><caption>**Table 5**
Temporal variation of the advancing dynamic contact angle for the case analyzed in Fig. 12.</caption><thead><tr><td>Time (s)</td><td>Dynamic contact angle (°)</td></tr></thead><tbody><tr><td>10<sup>-5</sup></td><td>70.7</td></tr><tr><td>1.5×10<sup>-5</sup></td><td>70.1</td></tr><tr><td>2×10<sup>-5</sup></td><td>69.5</td></tr><tr><td>3.35×10<sup>-5</sup></td><td>69.1</td></tr><tr><td>4.02×10<sup>-5</sup></td><td>70.9</td></tr></tbody></table>

simulations and better represents complex pore-scale interactions during gas injection and trapping. The results also illustrate how the geometrical factors, including the presence and size of constrictions, significantly influence pressure drops and phase distributions, thereby affecting hydrogen pore occupancy and breakthrough times. Although this study focuses on brine-saturated siliciclastic media, other strategies and lithologies have been investigated in the literature. Adsorptive storage using PAN-based nanofibers has shown promising hydrogen uptake at ambient conditions (Yu et al., 2023), while imbibition dynamics in carbonate and basaltic rocks exhibit distinct storage characteristics compared to sandstones (Hosseini et al., 2024). Incorporating such lithology-dependent effects into pore-scale simulations represents a valuable direction for future work.

Additionally, lower injection flow rates enhance hydrogen storage

efficiency and reduce the risk of overpressurization. Building on the current results obtained with idealized pore geometry, future work should aim to extend these simulations by incorporating additional physical parameters such as temperature effects, chemical reactions, and bacterial activity, as well as more complex pore geometries and topologies. This work assumes a rigid porous matrix, thus neglecting possible solid deformation effects. Recent studies (e.g., Yu et al., 2025) suggest that such coupling may alter flow pathways and storage efficiency, and its integration is left for future work. The present work provides application-oriented insights that support and validate the use of OpenFOAM pore-scale modeling for hydrogen storage. While the current study relies on a 2D axisymmetric pore geometry to provide fundamental insights in a controlled setting, we acknowledge that three-dimensional effects such as tortuosity and pore connectivity and the shape of the geometry could alter displacement dynamics, bubble formation, and trapping efficiency (Marafini et al. 2020). Extending this framework to 3D geometries is a key objective of future work, with the aim of providing more representative predictions for complex porous media and enabling upscaling toward pore-network and Darcy-scale models. Further studies should also investigate the effects of coordination number on flow and transport behavior.

CRediT authorship contribution statement

Bourama Diarra: Writing – original draft, Visualization, Validation,
Software, Methodology, Investigation, Formal analysis, Conceptualization. **Antonio Rodríguez de Castro:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Project administration, Methodology, Formal analysis, Conceptualization. **Sakir Amiroudine:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Project administration, Methodology, Formal analysis, Conceptualization.

Declaration of competing interest

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

Data availability

Data will be made available on request.

References

Abdelsamie, A., Fru, G., Oster, T., Dietzsch, F., Janiga, G., Thévenin, D., 2016. Towards direct numerical simulations of low-Mach number turbulent reacting and two-phase flows using immersed boundaries. Comput. Fluids 131, 123–141. https://doi.org/10.1016/j.comfluid.2016.03.017

Ali, M., Jha, N.K., Al-Yaseri, A., Zhang, Y., Iglauer, S., Sarmadivaleh, M., 2021. Hydrogen wettability of quartz substrates exposed to organic acids; Implications for hydrogen geo-storage in sandstone reservoirs. J. Pet. Sci. Eng. 207, 109081. https://doi.org/10.1016/j.petrol.2021.109081

Althoey, F., Wisner, B., Kontsos, A., Farnam, Y., 2018. Cementitious materials exposed to high concentration of sodium chloride solution: Formation of a deleterious chemical phase change. Constr. Build. Mater. 167, 543–552. https://doi.org/10.1016/j.conbuildmat.2018.02.066

Bagheri, M., Mahani, H., Ayatollahi, S., Zivar, D., 2023. Direct pore-scale simulation of the effect of capillary number and gas compressibility on cyclic underground hydrogen storage & production in heterogeneous aquifers. Adv. Water Resour. 181, 104547. https://doi.org/10.1016/j.advwatres.2023.104547

Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modeling surface tension. J. Comput. Phys. 100, 335–354. https://doi.org/10.1016/0021-9991(92)90240-Y

Cai, X., Wörner, M., Marschall, H., Deutschmann, O., 2017. CFD simulation of liquid Back Suction and Gas Bubble Formation in a Circular Tube with Sudden or Gradual expansion. Emiss. Control Sci. Technol. 3, 289–301. https://doi.org/10.1007/s40825-017-0073-3

Caiden, R., Fedkiw, R.P., Anderson, C., 2001. A Numerical Method for Two-phase Flow consisting of Separate Compressible and Incompressible Regions. J. Comput. Phys. 166, 1–27. https://doi.org/10.1006/jcph.2000.6624

Caltagirone, J.-P., Vincent, S., Caruyer, C., 2011. A multiphase compressible model for the simulation of multiphase flows. Comput. Fluids 50, 24–34. https://doi.org/10.1016/j.comfluid.2011.06.011

Cerne, G., Petelin, S., Tiselj, I., 2002. Numerical errors of the volume-of-fluid interface tracking algorithm. Int. J. Numer. Meth. Fluids 38, 329–350. https://doi.org/10.1002/fld.228

Chen, Q., Xu, J., Sun, D., Cao, Z., Xie, J., Xing, F., 2013. Numerical simulation of the modulated flow pattern for vertical upflows by the phase separation concept. Int. J. Multiph. Flow 56, 105–118. https://doi.org/10.1016/j.jmultiphaseflow.2013.05.014

Courant, R., Friedrichs, K., Lewy, H., 1928. Über die partiellen Differenzengleichungen der mathematischen Physik. Math. Ann. 100, 32–74. https://doi.org/10.1007/BF01448839

Da Riva, E., Del Col, D., 2009. Numerical simulation of churn flow in a vertical pipe. Chem. Eng. Sci. 64, 3753–3765. https://doi.org/10.1016/j.ces.2009.04.049

Dai, Y., Zhou, Z., Lin, J., Han, J., 2017. Modeling of Two-phase Flow in Rough-Walled Fracture using Level Set Method. Geofluids 2017, 1–11. https://doi.org/10.1155/2017/2429796

Dai, J., Shi, Y., Lei, X., Zhou, X., Pan, Z., Shen, X., Pi, S., 2025. Study on applicability of ball-and-stick model in reservoir pore-throat network simulation. Processes 13, 433. https://doi.org/10.3390/pr13020433

Dellacherie, S., 2005. On a diphasic low mach number system. ESAIM. Math. Mod. Numer. Anal. 39, 487–514. https://doi.org/10.1051/m2an:2005020

Derijani, H., James, L.A., Haynes, R.D., 2023. Evaluation of the interFoam solver in the prediction of immiscible two-phase flow in imbibition and drainage on the pore-doublet system. E3S Web Conf. 366, 01017. https://doi.org/10.1051/e3sconf/202336601017

Gennes, P.G.D., Hua, X., Levinson, P., 1990. Dynamics of wetting: local contact angles. J. Fluid Mech. 212, 55. https://doi.org/10.1017/S0022112090001859

Gostick, J.T., Khan, Z.A., Tranter, T.G., Kok, M.D.R., Agnaou, M., Sadeghi, M., Jervis, R., 2019. Numerical modeling of porous media images using PoreSpy: a Python toolkit for quantitative analysis. J. Open Source Softw. 4, 1296. https://doi.org/10.21105/j.oss.01296

Hao, Y., Prosperetti, A., 2004. A numerical method for three-dimensional gas-liquid flow computations. J. Comput. Phys. 196, 126–144. https://doi.org/10.1016/j.jcp.2003.10.032

Hashemi, L., Blunt, M., Hajibeygi, H., 2021. Pore-scale modelling and sensitivity analyses of hydrogen-brine multiphase flow in geological porous media. Sci. Rep. 11, 8348. https://doi.org/10.1038/s41598-021-87490-7

Hoseini, M., Ali, M., Fahimpour, J., Keshavarz, A., Iglauer, S., 2024. Energy storage in carbonate and basalt reservoirs: investigating secondary imbibition in H2 and CO2 systems. Adv. Geo-Energy Res. 11, 132–140. https://doi.org/10.46690/ager.2024.02.05.025

Huang, T., Moridis, G.J., Blasingame, T.A., Abdulkader, A.M., Yan, B., 2023. Compositional reservoir simulation of underground hydrogen storage in depleted gas reservoirs. Int. J. Hydrog. Energy 48, 36035–36050. https://doi.org/10.1016/j.ijhydene.2023.05.355

Iglauer, S., Ali, M., Keshavarz, A., 2021. Hydrogen Wettability of Sandstone Reservoirs: Implications for Hydrogen Geo-Storage. Geophys. Res. Lett. 48. https://doi.org/10.1029/2020gl090814

Ineris, 2021. Etat Des Connaissances Sur Le Stockage De L'hydrogène En Cavité Saline et Apport Du Projet ROSTOCK-H No. Ineris-159826-2702809-v1.0

InterFoam, 2024.

Lake, L., Johns, R.T., Rossen, W.R., Pope, G.A., 2014. Fundamentals of Enhanced Oil Recovery. Society of Petroleum Engineers. https://doi.org/10.2118/9781613993286

Lauga, E., Brenner, M. P., & Stone, H. A. *Microfluidics: The no-slip boundary condition.* In Handbook of Experimental Fluid Dynamics, pp. 1219–1240 (2007).

Li, B., Benson, S.M., 2015. Influence of small-scale heterogeneity on upward CO2 plume migration in storage aquifers. Adv. Water Resour. 83, 389–404. https://doi.org/10.1016/j.advwatres.2015.07.010

Lysyy, M., Ersland, G., Fernø, M., 2022. Pore-scale dynamics for underground porous media hydrogen storage. Adv. Water Resour. 163, 104167. https://doi.org/10.1016/j.advwatres.2022.104167

Lysyy, M., Liu, N., Solstad, C.M., Fernø, M.A., Ersland, G., 2023. Microfluidic hydrogen storage capacity and residual trapping during cyclic injections: Implications for underground storage. Int. J. Hydrog. Energy 48, 31294–31304. https://doi.org/10.1016/j.ijhydene.2023.04.253

Malgarinos, I., Nikolopoulos, N., Marengo, M., Antonini, C., Gavaises, M., 2014. VOF simulations of the contact angle dynamics during the drop spreading: Standard models and a new wetting force model. Adv. Colloid Interface Sci. 212, 1–20. https://doi.org/10.1016/j.cis.2014.07.004

Marafini, E., La Rocca, M., Fiori, A., Battiato, I., Prestininzi, P., 2020. Suitability of 2D modelling to evaluate flow properties in 3D porous media. Transp. Porous Media 134, 315–329. https://doi.org/10.1007/s11242-020-01447-4

Mowla, A., Ioanidisid, M.A., 2020. Effect of particle wettability on the hydrodynamics of three-phase fluidized beds subject to foaming. Powder Technol. 374, 58–69. https://doi.org/10.1016/j.powtec.2020.06.077

Nazari, M., Mahani, H., Ayatollahi, S., 2025. Pore-Scale Insights on Mixed-Wettability and its Impact on Underground Hydrogen Storage in Aquifers. Adv. Water Resour. 105044. https://doi.org/10.1016/j.advwatres.2025.105044

Panton, R. L., Kheirandish, S., et Wagner, M. H. (2007) *Experiment as a Boundary-Value Problem*, in Tropea, C., Yarin, A. L., et Foss, J. F. (éd.), *Handbook of Experimental Fluid Mechanics*, Springer, pp. 3–30.

Qian, T., Wang, X.-P., Sheng, P., 2005. Molecular hydrodynamics of the moving contact line in two-phase immiscible flows.

Rabbani, H.S., Joekar-Niasar, V., Pak, T., Shokri, N., 2017. New insights on the complex dynamics of two-phase flow in porous media under intermediate-wet conditions. Sci. Rep. 7, 4584. https://doi.org/10.1038/s41598-017-04545-4

Renardy, M., Renardy, Y., Li, J., 2001. Numerical simulation of moving Contact Line Problems using a Volume-of-Fluid Method. J. Comput. Phys. 171, 243–263. https://doi.org/10.1006/jcph.2001.6785

Rocco, S., Woods, A.W., Harrington, J., Norris, S., 2017. An experimental model of episodic gas release through fracture of fluid confined within a pressurized elastic reservoir. Geophys. Res. Lett. 44, 751–759. https://doi.org/10.1002/2016GL071546

Rodríguez de Castro, A., Goyeau, B., 2021. A pore network modelling approach to investigate the interplay between local and Darcy viscosities during the flow of shear-thinning fluids in porous media. J. Colloid Interface Sci. 590, 446–457. https://doi.org/10.1016/j.cjis.2021.01.081

Rodríguez de Castro, A., Goyeau, B., 2022. Numerical analysis of the pore-scale mechanisms controlling the efficiency of immiscible displacement of a pollutant phase by a shear-thinning fluid. Chem. Eng. Sci. 251, 117462. https://doi.org/10.1016/j.ces.2022.117462

Rodríguez de Castro, A., Chabanon, M., Goyeau, B., 2023a. Numerical analysis of the fluid-solid interactions during steady and oscillatory flows of non-Newtonian fluids through deformable porous media. Chem. Eng. Res. Des. 193, 38–53.

Rodríguez de Castro, A., Agnaou, M., Gostick, J., 2023b. Predicting shear-thinning fluid flows in porous media using pore network modeling: Simulations and experimental validation. Transp. Porous Media 149, 453–478. https://doi.org/10.1007/s11242-023-01968-8

Sahimi, M., 2012. Flow and transport in porous media and fractured rock: from classical methods to modern approaches. John Wiley & Sons.

Shi, Z., Jessen, K., Tsotsis, T.T., 2020. Impacts of the subsurface storage of natural gas and hydrogen mixtures. Int. J. Hydrog. Energy 45, 8757–8773. https://doi.org/10.1016/j.ijhydene.2020.01.044

Laplace, T., 1805. An essay on the cohesion of fluids philosophical transactions of Royal society of London. Typ Crapelet.
Tarkowski, R., 2019. Underground hydrogen storage: Characteristics and prospects. *Renew. Sustain. Energy Rev.* 105, 86–94. https://doi.org/10.1016/j.rser.2019.01.051.

Voinov, O.V., 1977. Hydrodynamics of wetting. *Fluid Dyn.* 11, 714–721. https://doi.org/10.1007/BF01012963.

Wang, J., Yang, Y., Cai, S., Yao, J., Xie, Q., 2023. Pore-scale modelling on hydrogen transport in porous media: Implications for hydrogen storage in saline aquifers. *Int. J. Hydrog. Energy* 48, 13922–13933. https://doi.org/10.1016/j.ijhydene.2022.11.299.

Wardle, K.E., Weller, H.G., 2013. Hybrid Multiphase CFD Solver for coupled Dispersed/Segregated Flows in Liquid-Liquid Extraction. *Int. J. Chem. Eng.* 2013, 1–13. https://doi.org/10.1155/2013/128936.

Weller, H.G., 2008. A new approach to VOF-based interface capturing methods for incompressible and compressible flow. *OpenCFD Ltd Rep. TRHGW* 4, 35.

Wörner, M., 2012. Numerical modeling of multiphase flows in microfluidics and micro process engineering: a review of methods and applications. *Microfluid. Nanofluidics* 12, 841–886. https://doi.org/10.1007/s10404-012-0940-8.

Yan, C., Yan, Z., Huo, T., Zhu, L., Wang, S., Li, X., 2025. MD-CFD investigations of hydrogen storage performance in aquifers considering media deformation. *Chem. Eng. Sci.* 304, 121069. https://doi.org/10.1016/j.ces.2024.121069.

Yin, Z., Khurana, M., Tan, H.K., Linga, P., 2018. A review of gas hydrate growth kinetic models. *Chem. Eng. J.* 342, 9–29. https://doi.org/10.1016/j.cej.2018.01.120.

Yu, J., Liang, W., Zhao, J., 2025. Enhancing dynamic modeling of porous media with compressible fluid: a THM material point method with improved fractional step formulation. *Comput. Methods Appl. Mech. Eng.* 444, 118100. https://doi.org/10.1016/j.cma.2025.118100.

Yu, J., Lin, T., Li, J., Zhang, W., Bao, W., Zhu, B., 2023. Construction of PAN-based activated carbon nanofibers for hydrogen storage under ambient pressure. *Capillarity* 6, 49–56. https://doi.org/10.46690/capi.2023.03.02.

Zhao, X., Li, Y., Zhang, H., Wang, L., Liu, J., 2023. Underground hydrogen storage: a recovery prediction using pore network modeling and machine learning. *Fuel* 357, 130051. https://doi.org/10.1016/j.fuel.2023.130051.