# Assessment of H2 diffusivity in water and brine for underground storage: A molecular dynamics approach

Halla Kerkache, Hai Hoang, Tu Khai Nam Nguyen, Antoine Geoffroy-Neveux, Carlos Nieto-Draghi, Pierre Cézac, Salaheddine Chabab, Guillaume Galliéro

► To cite this version:

Halla Kerkache, Hai Hoang, Tu Khai Nam Nguyen, Antoine Geoffroy-Neveux, Carlos Nieto-Draghi, et al. Assessment of H2 diffusivity in water and brine for underground storage: A molecular dynamics approach. International Journal of Hydrogen Energy, 2025, 128, pp.279-290. 10.1016/j.ijhydene.2025.04.157. hal-05075832

HAL Id: hal-05075832

https://ifp.hal.science/hal-05075832v1

Submitted on 20 May 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>

# Assessment of H₂ diffusivity in water and brine for underground storage: A molecular dynamics approach

Halla Kerkache<sup>a, b</sup>, Hai Hoang<sup>c, d</sup>, Tu Khai Nam Nguyen<sup>c, d</sup>, Antoine Geoffroy-Neveux<sup>e</sup>, Carlos Nieto-Draghi<sup>c, f</sup>, Pierre Cézac<sup>b</sup>, Salaheddine Chabab<sup>a, b, **</sup>, Guillaume Galliéro<sup>a, *</sup>

## ARTICLE INFO

Handling Editor: Dr F Gallucci

**Keywords:**

Hydrogen

Water

Brine

Diffusion

Non-equilibrium molecular dynamics

Excess entropy

## ABSTRACT

Hydrogen (H₂) is poised to play a key role in the future energy landscape, but its efficient storage remains a challenge due to its small size and high mobility. Underground hydrogen storage in aquifers offers potential, yet risks arise from H₂ diffusivity in brine, which can lead to migration or chemical reactions. In this study, H₂ diffusivity in water and NaCl brine is studied using Equilibrium and Non-Equilibrium Molecular Dynamics simulations. A detailed benchmarking of force fields for H₂, water, and Na⁺/Cl⁻ ions is performed to ensure accuracy and reliability. These results provide diffusivity data for H₂ across a broad range of underground conditions (298 K–473 K, 1–1000 bar, 0–6 mol/kg NaCl), analyzed within an excess entropy framework and modeled using the Stokes-Einstein relationship and an Arrhenius-type equation, offering a framework for predicting H₂ diffusion behavior in water and brine in various conditions.

## 1. Introduction

Hydrogen (H₂) is expected to play a major role in reducing greenhouse gas emissions during the energy transition. For its deployment, in the short term, hydrogen can be stored in small quantities in solid, liquid (liquefied or transformed into H₂ carrier molecules) or gaseous form in metallic reservoirs [1]. In the medium/long term, and for large-scale storage, underground geological reservoirs, historically used for storing town gas and currently for natural gas, are promising candidates. These include salt caverns, depleted natural gas fields, and aquifers. Storage in salt caverns is a well-established technology that offers considerable flexibility (with storage cycle ranging from monthly to even daily) and high safety due to the high impermeability of salt layers [2]. In contrast, porous reservoirs, such as natural gas fields and aquifers, are more suited for seasonal and strategic storage. They offer the advantage of massive storage volumes and better geographical availability (especially aquifers) [3], but their development is not yet mature and requires further research and technological advancement.

While experience from town gas storage (where H₂ is present in

significant quantities) has shown that the composition of stored gas can change significantly [4], the behavior of H₂ in aquifers or depleted gas fields remains poorly understood. When hydrogen is injected into this underground reservoir, it can disturb the chemical equilibrium between the formation water, dissolved gases and rock minerals, potentially affecting the stability and mobility of the hydrogen and the integrity of the storage system. More specifically, in such reservoirs, a water-saturated caprock typically acts as a physical barrier to the migration of hydrogen to the surface. As long as capillary pressure is not exceeded and mechanical fracturing does not occur, hydrogen losses through the caprock will primarily occur via molecular diffusion [5]. Therefore, one of the crucial properties to evaluate is the diffusion coefficient of H₂ in brine under storage conditions (temperature, pressure and salinity). Furthermore, accurate knowledge of the diffusion coefficient of H₂ in reservoir brine is important for estimating hydrogen losses by dissolution and biogeochemical reactivity in the aquifer, as hydrogen is only reactive in the aqueous phase under storage conditions.

Several studies in the literature provide experimental data on H₂ diffusion in pure water at low temperature and pressure [6–14]. These

<footnote><sup>a</sup> Université de Pau et des Pays de L'Adour, CNRS, LFCR, UMR5150, Pau, France</footnote>

<footnote><sup>b</sup> Université de Pau et des Pays de L'Adour, LaTEP, Pau, France</footnote>

<footnote><sup>c</sup> Institute of Fundamental and Applied Sciences, Duy Tan University, Tran Nhat Duat Street, District 1, Ho Chi Minh City, Viet Nam</footnote>

<footnote><sup>d</sup> Faculty of Environmental and Natural Sciences, Duy Tan University, 03 Quang Trung Street, Da Nang, Viet Nam</footnote>

<footnote><sup>e</sup> IFP Energies Nouvelles, 1 et 4 Avenue de Bois Préau, 92852, Rueil-Malmaison, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>** Corresponding author. Université de Pau et des Pays de l'Adour, CNRS, LFCR, UMR5150, Pau, France.</footnote>

<footnote>E-mail addresses: salaheddine.chabab@univ-pau.fr (S. Chabab), guillaume.galliero@univ-pau.fr (G. Galliéro).</footnote>
measurements are typically conducted at atmospheric pressure and within the temperature range of 278.15-333.15 K using methods such as the rate of collapse of small bubbles, permeability of a stagnant liquid layer in the quasi-steady state, and diaphragm cell, with the exception of a recent study Wang et al., 2023 [13], which explores diffusion at higher temperatures (up to 423 K) and pressures (30 MPa) using Taylor Dispersion technique. However, given the complexity of such measurements, these studies exhibit up to 70 % variability in their results, according to Tsimpanogiannis et al. [15]. Regarding brine, to the best of the author's knowledge, there is only one experimental set of data on H₂ diffusion in NaCl brine reported recently by Piszko et al., 2024 [16] (at temperature of 293K and salinity up to 0.34 mol of NaCl/kg of water). This study highlighted the difficulties of carrying out these measurements at higher temperatures, citing analytical difficulties and issues with the thermal homogenization in diffusion cells.

Molecular simulation offers a valuable tool to complement experimental data at subsurface conditions relevant to hydrogen storage and native hydrogen production. While some molecular simulation studies exist on H₂ diffusion in pure water [15,17-25], only two recent studies have focused on H₂ diffusion in brine under subsurface storage conditions [24,25]. However, these studies employed less realistic models for H₂, where H₂ was represented without a quadrupole moment. In addition, they used the equilibrium molecular dynamics approach, which can introduce notable uncertainties systems at infinite dilution, such as H₂-H₂O and H₂-H₂O-NaCl systems [25].

As H₂ is highly diluted in water and brine, its Maxwell-Stefan and Fick diffusion coefficients can be well approximated to by its self-diffusion coefficient [26,27]. Hence, to advance the current state of the art, in this work, the self-diffusion coefficient of H₂ in brine was calculated using a more realistic model for H₂ that accurately reproduces H₂ solubility in water and brine. Additionally, has been employed a different methodological approach compared to those used in recent studies [24,25], ensuring improved accuracy of the reported results, namely the external-field non-equilibrium molecular dynamic (NEMD) [28]. This method significantly reduces the uncertainty compared to those obtained using the Equilibrium Molecular Dynamic (EMD) method that was used previously. The EMD method is also used here to cross-validate the results obtained with the NEMD method [28]. Furthermore, to provide insights into the microscopic structure of the H₂-brine system at the molecular level, radial distribution functions (RDF) were calculated and were used to estimate the two-body excess entropy so as to analyze the results within the entropy scaling framework [29]. Finally, the H₂ diffusivity results were used to develop an easy-to-use model applicable to a wide range of underground thermodynamic conditions.

This article is structured as follows. First, the methods and simulation details are described in section 2. Then, in section 3, the results of simulations using different molecular models and computational approaches are presented, discussed and correlated in section 4 to facilitate the use of the data under various temperature, pressure and salinity conditions.

## 2. Molecular simulation methods

Molecular simulation is a powerful and well-established tool for computing diffusion coefficients in pure fluids and mixtures, serving as a valuable complement to experimental data and theoretical models. With this numerical tool, it is possible to compute both self-diffusion and mutual (also known as transport or collective) diffusion. Self-diffusivity characterizes the mobility of individual molecules in pure fluids or mixture driven by the random walk (Brownian) motion of molecules in the absence of any external driving force for diffusion [30-32]. In contrast, mutual diffusivity is only defined for mixtures and describe the net transport of a group of molecules in response to a driving force. Mutual diffusion can be represented within both Maxwell-Stefan (MS) and Fick frameworks. MS diffusivity quantifies the mass transport due to

a gradient in the chemical potential while Fick diffusivity quantifies the mass transport due to a gradient in the concentration [28,33]. As both frameworks describe the same physical process, they are related by the so-called Thermodynamic Factor, which accounts for the non-ideality of mixtures [34,35].

In molecular dynamic, the thermodynamic factor can be determined by both direct and indirect approaches. Indirect approaches involve fitting simulation data to activity models (obtained from Widom's Test Particle Insertion (WTPI)) or equations of state. Direct approaches include several techniques such as:1) Kirkwood-Buff integrals (i.e. it converges for large systems) [36], 2) Simulations in the grand-canonical ensemble [37,38], 3) The Permuted Widom Test Particle Insertion (PWTPI) method [39], and 4) The Continuous Fractional Component Monte Carlo (CFCMC) with the adopted PWTPI method [40].

For systems where the solute is infinitely diluted in a solvent, as in the case of hydrogen - water/brine system studied in this work, the MS, Fick, and self-diffusivities are identical. Therefore, in this work, the self-diffusion coefficient (SDC) of H₂ in water/brine has been computed.

### 2.1. Force fields

For water molecules, the rigid non-polarizable TIP4P/2005 force field (FF) was selected due to its simplicity, its ability to accurately predict thermophysical properties and computational efficiency [41]. Indeed, the density of water is very well reproduced by TIP4P/2005 FF throughout the range of application of this work (from 1 to 1000 bar and from 273 to 473 K) [15]. In the presence of NaCl salt and in combination with the TIP4P/2005 model, Na⁺ and Cl⁻ ions can be represented with different models that were investigated in our previous study on hydrogen solubility in water and NaCl brine [42]. To ensure a good representation of the dielectric constant of the brine, the Madrid model [43,44], that uses scaled charge of 0.85 for Na⁺ and -0.85 for Cl⁻, was used. This adjustment allows accurate predictions of NaCl brine density and hydrogen solubility in brine. Furthermore, to improve the prediction of transport properties of NaCl solutions, such as viscosities and diffusion coefficients, Blazquez et al. [45] improved the Madrid model (denoted here as "MadridTransport" model) by slightly modifying the Lennard-Jones (LJ) parameters and further scaling the Na⁺ and Cl⁻ ion charges to 0.75 and -0.75, respectively [45]. In this work, both Madrid and MadridTransport force fields were tested to represent NaCl brine.

For the H₂ molecule, which is a light gas, its quantum nature can influence its properties. However, for the temperatures considered in this study, above 298 K, these effects are negligible (Marx and Nielaba [46]; Gunn et al. [47]). Three models were selected to describe the H₂ molecule: a single site non-charged model "Vrabec" [48], a two site non-charged model "Yang & Zhong" [49] and a multiple-sites charged model "Marx" [46]. According to the literature, all these force fields are able to correctly reproduce the thermophysical properties of pure H₂ [15]. The Yang and Zhong's model, which was used in our previous work on the solubility of H₂ in brine [42], did not perform well in predicting solubility due to the absence of charges, an important factor in predicting phase equilibria. However, this force field has been used in several studies and is able to predict well the diffusivity of H₂ in bulk (water/brine) [15,18,19,24,25] and in clay [23]. For H₂ solubility in water and brine, the Marx force [24,42,50,51] gives very good results compared to experimental data, whereas the Vrabec model exhibits significant deviations (around 50 %) according to van Roojen et al. [24]. Regarding H₂ self-diffusivity in pure water, Tsimpanogiannis et al. reported that both the Vrabec and Marx models, when combined with the TIP4P/2005 water model, yielded good result over a wide range of temperatures and pressures [15].

The classical Lorentz-Berthelot (LB) [30,52,53] combining rules have been used to determine the Lennard-Jones interaction parameters of unlike atoms. However, as indicated in our recent study [42], correcting the cross-interaction energies between H₂ and water by multiplying them by an optimal constant factor equal to 1.05 for the
Marx-TIP4P/2005 force fields give a good representation of the experimental data of the solubility of H₂ in pure water and brine over a very wide range of temperatures. Thus, in this work, this correction was considered for Marx-TIP4P/2005 force fields.

Details of the force fields parameters used in this work are provided in Tables 1, 2 and 3 in the Supplementary Materials (SM).

## 2.2. Methodology

All Molecular Dynamics simulations (Equilibrium and Non-Equilibrium) were performed using the LAMMPS package [54] and the initial configurations were generated using the Towhee software [55].

### 2.2.1. Self-diffusion coefficient from EMD simulations

The self-diffusion coefficient of H₂ in water/brine is usually computed from EMD simulations using Einstein's relation [30]:

$$D_{H_2} = \frac{1}{6N_{H_2}} \lim_{t \to \infty} \frac{\left( \sum_{i=1}^{N_{H_2}} \left[ \mathbf{r}_{H_2}^i(t) - \mathbf{r}_{H_2}^i(0) \right]^2 \right)}{t} \quad (1)$$

where $N_{H_2}$ is the number of H₂ molecules, and $\mathbf{r}_{H_2}^i(t)$ is the position of H₂ molecule at the time t. This relation is derived from the fluctuation-dissipation theorem at thermodynamic equilibrium, but it often results in a slow convergence and a low signal-to-noise ratio. This is especially true for systems with very high H₂ dilution, like those studied in this work, where only a single solute molecule is present in the simulation box, leading to significant statistical uncertainties [56,57].

To improve accuracy and reliability, it is common to average the results over multiple initial time points. There are two methods for this averaging: “post-processing”, which occurs after the simulation is completed, and “on-the-fly”, which takes place during the simulation itself. The latter method is generally preferred as it reduces memory usage and computational time. However, implementing the “on-the-fly” method can be challenging. Recently, the OCTP plugin [58] was developed to compute the self-diffusion coefficient using the “on-the-fly” method in LAMMPS. While the OCTP package is very effective, it presents challenges when applied to multi-component systems that require different simulation algorithms for each component. The plugin calculates the self-diffusion coefficient for all groups defined in the LAMMPS input script. To obtain the self-diffusion coefficient of a specific compound, a group must be defined that includes only one atom of that compound. This requirement can be inconvenient, as multiple groups need to be defined to accommodate different simulation algorithms, making the use of the OCTP plugin less practical in these multi-component cases.

### 2.2.2. Self-diffusion coefficient from NEMD simulations

An alternative approach to determine the self-diffusion coefficient of H₂ in water/brine is to employ NEMD simulations by applying external forces to the simulation system [28]. The H₂ self-diffusion coefficient is then calculated using the theory of thermodynamics of irreversible processes (TIP) as follows [59]. The diffusive flux of H₂ in a solvent $J_{H_2}$, as described by TIP, is given as:

$$J_{H_2} = c_{H_2} \times \mathbf{u}_{H_2} = c_{H_2} \frac{L_{H_2 H_2}}{RT} \mathbf{F}_{H_2} + \sum_{i \in \text{Solvent}} c_i \frac{L_{H_2}}{RT} \mathbf{F}_i \quad (2)$$

where c is the concentration, u is the molar diffusive velocity, R is the gas constant, T is the temperature, F is the driving force acting on species (in Kcal/Å/mol), and L are the Onsager's coefficients. The molar diffusive velocity of H₂, $\mathbf{u}_{H_2}$, is defined as:

$$\mathbf{u}_{H_2} = \frac{1}{N_{H_2}} \sum_{i=1}^{N_{H_2}} \mathbf{v}_{H_2}^i - \mathbf{V} \quad (3)$$

where $\mathbf{v}_{H_2}^i$ is the velocity of H₂ molecule i, and V is the drift velocity of the whole system. Onsager's coefficient $L_{H_2,i}$ can be expressed in terms of displacement of the center of mass of species by using Einstein's relation as:

$$L_{H_2,i} = \frac{N_{H_2}}{6} \lim_{t \to \infty} \frac{\left( \mathbf{r}_{H_2}^{\text{com}}(t) - \mathbf{r}_{H_2}^{\text{com}}(0) \right) \times \left( \mathbf{r}_i^{\text{com}}(t) - \mathbf{r}_i^{\text{com}}(0) \right)}{t} \quad (4)$$

where $\mathbf{r}_i^{\text{com}}(t)$ is the center of mass of species i at the time t. This expression leads to $c_i \times L_{H_2,i} = c_{H_2} \times L_{i,H_2}$. Therefore, the molar diffusive velocity of H₂ reduces to:

$$\mathbf{u}_{H_2} = \frac{L_{H_2 H_2}}{RT} \mathbf{F}_{H_2} + \sum_{i}^{L_{i,H_2}} \mathbf{F}_i \quad (5)$$

For solutions at very high H₂ dilution, as here where the simulation systems contain only one H₂ molecule, the H₂ self-diffusion coefficient $D_{H_2}$ is physically identical to the Onsager's coefficient $L_{H_2,H_2}$.

Using Equation (5), the self-diffusion coefficient of infinitely diluted H₂ in water/brine can be determined by performing NEMD simulations as follows. An external force $\mathbf{F}_{H_2}^{\text{Ext}}$ is applied to the H₂ molecule, and another external force $F_{\text{Solv}}^{\text{Ext}} = c_{H_2} \times F_{H_2}^{\text{Ext}} / \sum_i c_i$ is applied to all solvent molecules. For such solutions, where $c_{H_2} \ll \sum_i c_i$, the external force acting on the solvent molecules $F_{\text{Solv}}^{\text{Ext}}$ is negligible compared to the external force acting on the H₂ molecule $F_{H_2}^{\text{Ext}}$. Therefore, the molar diffusive velocity of H₂ derived from equation (5), can be simplified as:

$$\mathbf{u}_{H_2} = \frac{L_{H_2 H_2}}{RT} \mathbf{F}_{H_2}^{\text{Ext}} \quad (6)$$

Then, the H₂ self-diffusion coefficient is simply given by:

$$D_{H_2} = L_{H_2 H_2} = RT \frac{\mathbf{u}_{H_2}}{\mathbf{F}_{H_2}^{\text{Ext}}} = \frac{RT}{\mathbf{F}_{H_2}^{\text{Ext}}} \lim_{t \to \infty} \frac{\mathbf{r}_{H_2}(t)}{t} \quad (7)$$

As the position of the hydrogen molecule, $\mathbf{r}_{H_2}(t)$, is known during the simulation, calculating $D_{H_2}$ with the NEMD approach is rather straightforward and efficient. However, it is crucial to choose appropriate values for the applied external force $F_{H_2}^{\text{Ext}}$. It should be neither too small to ensure good statistical uncertainties, nor too large to keep the simulation system in the linear response regime. Thus, in this work, $F_{H_2}^{\text{Ext}}$ applied have been taken between 0.20 and 0.45 in reduced units.

### 2.2.3. Excess entropy calculation

The self-diffusion coefficient of a solute in a solvent is linked to its excess entropy [29], s, a physical quantity that reflects the fluid microscopic structure [60]. More precisely, the self-diffusion coefficient of a given species scales exponentially with its excess entropy (D increase with s). Therefore, to provide insights into how the force fields, salinity and thermodynamic conditions affect the diffusion coefficient of H₂ in water/brine, the two-body partial excess entropy of H₂ was determined as follows [29]:

$$s_{H_2} = -\frac{1}{2} \sum_{\nu} 4\pi r^2 \rho_{\nu} \int \left\{ g_{H_2-\nu}(r) \ln g_{H_2-\nu}(r) - [g_{H_2-\nu}(r) - 1] \right\} dr \quad (8)$$

where, $\rho_{\nu}$ is the density of component $\nu$ and $g_{H_2-\nu}(r)$ represents the radial distribution function (RDF) of component $\nu$ surrounding H₂ at a distance r. Physically, a more negative value of $s_{H_2}$ corresponds to a more structured solvent around H₂.

In practice, the RDFs $g_{H_2-\nu}(r)$, relatively to the center of mass of H₂ molecule, were computed by performing EMD simulations as follows [61]:

$$g_{H_2-\nu}(r) = \frac{\rho_{H_2-\nu}(r)}{\rho_{\nu}} = \left\langle \frac{1}{N_{H_2}} \times \frac{V}{N_{\nu}} \sum_{i \in H_2}^{N_{H_2}} \sum_{j \in \nu}^{N_{\nu}} \frac{\delta(r_{ij} - r)}{4\pi r^2} \right\rangle \quad (9)$$
where, $\rho_{H_2-\nu}(r)$ denotes the number density of component $\nu$ at a distance $r$ from $H_2$, $\delta$ stands for the Dirac function, $r_{ij}$ is the distance between molecules i of $H_2$ and particle j of component $\nu$, $V$ is the volume, and $N_\nu$ indicates the number of molecules of component $\nu$.

### 2.2.4. Simulation details

Both the EMD and NEMD simulations to compute the $H_2$ self-diffusion coefficient were conducted in the NVT ensemble. To determine the density of the solutions for these simulations, EMD simulations in the NPT ensemble were first performed in two steps. In the equilibration step, the systems were equilibrated during $2 \times 10^6$ time steps. This was followed by the production step, where a run of $5 \times 10^6$ was carried out for the density computation.

The NVT-EMD simulations to calculate $H_2$ self-diffusion coefficient were also performed in two steps. First, the equilibration step consisted in a run of $2 \times 10^6$ time steps, followed by a production run of $2 \times 10^6$ for samplings. The NVT-NEMD simulations involved three steps: First, the simulation system was relaxed for $2 \times 10^6$ time steps. Then, during the transient step, the external forces were applied to the molecules within the simulation systems over $1 \times 10^7$ time steps to reach the steady state. Finally, a production run of $5 \times 10^6$ time steps was conducted to perform samplings, exceptions noted.

All simulation systems consist of 1000 water molecules, 1 hydrogen molecule, and between 0 and 108 Na+ and Cl-ion pairs, depending on the salt molality. Periodic boundary conditions were imposed in all directions. The Newtonian equations of motion were integrated using the Verlet-velocity algorithm with a time step of 1 fs. In the NPT-EMD simulations, the temperature and pressure were kept constant using the Nose-Hoover algorithm, which was applied on all three directions [30,62]. For the NVT-MD simulations without applying the external forces, i.e. NVT-EMD and the first step of NVT-NEMD, the Nose-Hoover thermostat was applied to all three components of atomic velocities. Otherwise, when external forces were applied in the subsequent two steps of NVT-NEMD simulations, the thermostat was applied only to the two velocity components perpendicular to the external forces to remove the effect of streaming velocity along the direction of the applied forces [63].

To maintain the bond length and the bending angle of $H_2O$, as well as the bond length of $H_2$ using Yang and Zhong's model, the RATTLE algorithm was used [64]. For $H_2$, using the Marx model, which is a 3-site linear model where RATTLE cannot be applied, molecular integrity was maintained using the rigid body function in LAMMPS. As a consequence, no thermostat was applied directly to the $H_2$ solute molecule in simulations using the Marx model (its thermalization occurred via collisions with the surrounding solvent molecules). The Lennard-Jones interactions were computed by truncating at a cutoff of $12 \AA^\circ$, and long-range corrections were included. The Coulomb interactions were computed by using Particle-Particle-Particle-Mesh (PPPM) method with a relative error of $10^{-5}$ [65].

All self-diffusion coefficients were corrected for finite-size effects using the Yeh-Hummer equation [66]. This correction, valid for cubic simulation boxes, allows extrapolation of diffusion coefficient for infinite size systems ($D_i^\infty$) from results obtained with finite-size systems ($D_i$), using the following expression:

$$D_i^\infty = D_i + \frac{k_B T_\xi}{6\pi\eta L} \qquad (10)$$

where $k_B$ in the Boltzmann constant, $T$ is temperature, $\eta$ is the solvent viscosity, calculated from Mao and Duan (2009) model [67], $L$ is the simulation box length, and $\xi = 2.837298$. In this work, the magnitude of the finite size correction ranged from 4 to 8 % of the calculated self-diffusion coefficients, depending on the thermodynamic conditions.

To improve the statistical uncertainties of the RDF computed from the EMD simulation results, 20 independent simulations were performed. The error bars represent the standard deviation of these results.

For the NEMD simulation results, the error bars were computed using the block method by dividing the sampling step into six blocks. For the clarity, the error bars have been omitted on some of the figures but are provided in the Supplementary Materials.

Simulations details, including specific thermodynamic conditions used in our study, and results are provided in the SM (Tables 4–6).

## 3. Results and discussions

### 3.1. Existing data and correlation of $H_2$ diffusion in pure water

As previously noted, according to Tsimpanogiannis et al. [15], the available $H_2$ diffusivity data in pure water prior to 2021 were all at atmospheric pressure and temperatures no higher than 333 K. These data exhibit deviations up to 70 %. To address this variability, the authors proposed two distinct Arrhenius-type correlations, shown in Fig. 2, using $D = D_0 \exp(\alpha/T)$; for correlation 1: $D_0 = 4.931 \times 10^{-6} \text{ m}^2/\text{s}$ and $\alpha = 1990 \text{ K}$; and for correlation 2: $D_0 = 1.691 \times 10^{-6} \text{ 1}$ and $\alpha = 1810 \text{ K}$. Based on their EMD simulation results, the authors concluded that correlation 2, with its weaker temperature dependence, is the most consistent and reliable one. More recently, Wang et al. [13] extended the experimental conditions to higher temperatures (up to 423 K) and two pressures (5 and 280 bar). As shown in Fig. 2, the data from Wang et al., 2023 also align well with correlation 2. This agreement supports the idea that correlation 2 is the most consistent correlation for $H_2$ self-diffusion in pure water at low pressures.

### 3.2. Preliminary results: brine density

Before computing $H_2$ diffusion in water/brine, it was essential to validate the correct implementation of the selected force fields and to assess the accuracy of the predicted solvent density, particularly for NaCl brine. Thus, Fig. 2 presents the results of brine density computations using the TIP4P/2005-Madrid and TIP4P/2005-MadridTransport models. The correlation proposed by Al Ghafri et al., 2012 [68], obtained by fitting their experimental density data for the $H_2O + \text{NaCl}$ system, was used as a reference.

As shown in Fig. 2 and in Table 4 of the SM, both force fields provide the brine density very accurately at different salinities. This results is consistent with literature [24,45] and confirms the successful implementation of the selected force fields for water and ions.

### 3.3. EMD vs. NEMD comparison

As discussed in Section 2.2, calculating the diffusion coefficient of $H_2$ in water and brine using the NEMD simulations offers a potentially advantageous alternative to the usual EMD simulations approach. Thus, in this subsection, the efficiency of the simulation using these two methods is evaluated.

For this purpose, a combination of force fields was selected: Marx for $H_2$, TIP4P/2005 for water and MadridTransport for NaCl. To ensure a fair comparison between the two methods, the production step was performed using the same number of time steps for both simulation types. More precisely, the production step was divided into six blocks, with each block consisting of $2 \times 10^6$ timesteps.

Fig. 3a and 3b show the relationship between the mean square displacement (MSD) of $H_2$ ($=\langle|r_{H_2}(t) - r_{H_2}(0)|^2\rangle$) and $t$, as used in Equation (1) to determine self-diffusion coefficient of $H_2$ in water and brine at $T = 323\text{K}$ and $P = 200\text{Bar}$, with $m = 3 \text{ M}$ for brine. These results were obtained from EMD simulations using the “post-processing” method. Each line represents the simulation results of individual blocks. The results show significant deviations between the lines, indicating large statistical uncertainties in the simulation outcomes. Furthermore, as $t$ increases, these uncertainties tend to become more pronounced. This is because the number of samples used to calculate the average value of the MSD of $H_2$ decreases as the time windows used for the computation
increases. The time window directly affects the determination of the slope in the linear relationship between MSD and $t$ in Equation (1).

For instance, the diffusion coefficient of H₂ in water so computed is $D_{H_2} = 5.76 \pm 0.94 (\times 10^{-9}m^2/s)$ for a time window of 100 ps, and $D_{H_2} = 5.38 \pm 1.21 (\times 10^{-9}m^2/s)$ for a time window of 200 ps. Time windows of this magnitude are often considered as sufficient to ensure decorrelation of velocity fluctuations when calculating diffusion coefficients in dense fluids. Similarly, for H₂ in brine, the diffusion coefficient computed with such EMD simulations is $D_{H_2} = 4.50 \pm 0.34 (\times 10^{-9}m^2/s)$ for a time window of 100 ps, and $D_{H_2} = 4.88 \pm 0.52 (\times 10^{-9}m^2/s)$ for a time window of 200 ps.

Fig. 3c and 3d presents the relationship between $r_{H_2}(t)$ and $t$ for H₂ in water and brine, used to compute self-diffusion of H₂, see Equation (7). These results were obtained from NEMD simulations under the same conditions as the EMD simulations. In contrast to the EMD, the lines representing individual blocks in the NEMD simulations nearly overlap. This suggests that the diffusion coefficients calculated from NEMD simulations exhibit smaller statistical uncertainties compared to those from EMD. Specifically, the diffusion coefficient of H₂ in water and brine are $D_{H_2} = 6.18 \pm 0.44 (\times 10^{-9}m^2/s)$ and $D_{H_2} = 4.89 \pm 0.16 (\times 10^{-9}m^2/s)$, respectively. These results exhibit error bars that are approximately two or three times smaller than those obtained from EMD simulations for an equivalent computational time. In addition, it is worth pointing out that the results obtained from both methods are consistent with each other, even if not perfectly aligned.

Based on these findings and other results for different thermodynamic conditions provide in Table 5 of the SM, NEMD simulations were used to calculate the diffusion coefficient of H₂ in water and brine under various thermodynamic conditions in the subsequent section.

### 3.4. Force fields comparison and selection

#### 3.4.1. H₂ + H₂O system

As a benchmark exercise, the force field combinations recently reported in the literature [15,24,25] for calculating the diffusivity of H₂ in pure water by the EMD method (Vrabeč-TIP4P/2005, Yang&Zhong-TIP4P/2005) were evaluated in this work using NEMD. In addition, the best combination of force fields for thermodynamic properties of H₂ and

water (density and Henry's law constant) reported in our previous study [42], i.e., Marx-TIP4P/2005, was investigated as well.

In Figs. 1 and 4, the results of the NEMD calculations of the diffusivity of H₂ in pure water at low (Fig. 1) and high (Fig. 4) pressures are compared with recent EMD results from the literature [15,24,25], as well as with the two correlations previously discussed in section 3.1. Interestingly, these figures confirm that our NEMD simulation provides results consistent with those coming from EMD simulations from the literature. In addition, when compared with reference results (experimental data or correlation 2), it appears that all three force field combinations effectively capture the temperature dependence of H₂ diffusion coefficients. This is likely because all these simulations use the same water force field, TIP4P/2005, which predicts very accurately the viscosity of water and its temperature dependence. As the diffusion coefficient is inversely proportional to viscosity (cf. the Stokes-Einstein relation [27]), the consistency between the force fields predictions and experimental data can be attributed to the accurate representation of

<figure><img src="image_4.png" /><figcaption>Fig. 2. Density of aqueous NaCl solution as a function of molality $m$ (mol NaCl/kg H₂O) at a pressure of 1 bar and a temperature of 298 K. Uncertainties are smaller than 0.2 %.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 1. Variation of H₂ self-diffusion in pure water at different thermodynamic conditions. Comparison of the MD results from this work with experimental data, existing correlations, and MD data from the literature at low pressures (1–10 bar). The dashed and solid lines are Arrhenius-type correlation [15].</figcaption></figure>
<figure><img src="image_6.png" /></figure>

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

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

<figure><img src="image_9.png" /><figcaption>Fig. 3. EMD vs. NEMD results. Mean Square Displacements of the hydrogen molecule in water and brine (a and b, respectively) and evolution of the position of the H₂ molecule *r*<sub>H₂</sub> in the water and brine (c and d, respectively) as a function of time at 323K, 200 bar and 3 M for brine.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 4. Temperature dependence of the H₂ diffusion in pure water. Comparison of the MD results from this work with experimental data and MD data from the literature at high pressure (≈200 bar).</figcaption></figure>

water viscosity. Furthermore, our results confirm that the H₂ diffusion is weakly dependent on pressure, which can also be explained by the relatively small change in water viscosity under the pressure conditions explored in this study.

Although the choice of the water model appears to be the dominant factor, the H₂ force field also plays a role, as shown in Figs. 2 and 4. Among the tested combinations, the Marx-TIP4P/2005 combination

shows the smallest deviation from the experimental data of Zhang et al., 2023 (average absolute relative deviation (AAD) of 5.2 %), followed by the Vrabeč-TIP4P/2005 combination (AAD of 8.8 %) and the Yang&Z-hong-TIP4P/2005 combination (AAD of 17.3 %).

This influence of the H₂ force field on its self-diffusion coefficient in water can be analyzed in terms of its excess entropy, which reflects the local structure around the H₂ molecule. It was found that its two-body
excess entropy, $s_{H_2}$, is approximately the same for the Marx and Vrabec force fields, whereas the Yong & Zhang force field show a less negative value, see Table 1. This indicates that the water structure around $H_2$ is less organized with the Yang & Zhong force field compared to the other two force fields, leading to a slightly higher $H_2$ self-diffusion coefficient for the Yang & Zhong force field.

### 3.4.2. $H_2 + H_2O + NaCl$ system

Fig. 5 presents the predictions of $H_2$ diffusion in NaCl brine using the two NaCl force fields, Madrid (Fig. 5-a) and MadridTransport (Fig. 5-b) combined with the three $H_2$ force fields (Vrabec, Yang&Zhong, and Marx) tested earlier. The NEMD results from this work are comparable with those of Kalati et al. (EMD with Yang&Zhong) and van Rooijen et al. (EMD with Vrabec). Similar trends are observed against salinity, although small deviations (of about $\pm 0.5 \times 10^{-9} \text{ m}^2/\text{s}$), within the range of uncertainties, are noted. In terms of experimental comparison, the trends are consistent, showing a decrease of $H_2$ diffusion with increasing salinity. However, the comparison remains qualitative as, to the best our knowledge, only one set of experimental data of $H_2$ diffusion is available [16]. and that data is limited to low salinity brine and is reported with large uncertainties.

Fig. 5-a and 5-b also show that the Vrabec and Marx force fields produce almost identical $H_2$ diffusion values over the entire salinity range, consistently with their behavior in pure water. The Marx force field, which is the only that accurately captures both the solubility (Henry's law constant), see our previous work [42], and diffusion of $H_2$ in brine, was therefore selected for further data production and correlation, as discussed in the following.

As mentioned earlier, since the diffusion coefficient is inversely

<table><caption>**Table 1**
Excess entropies.</caption><thead><tr><th></th><th>S<sup>Oa</sup><sub>Hm</sub></th><th>S<sup>Na</sup><sub>Hm</sub></th><th>S<sup>Cl</sup><sub>Hm</sub></th><th>S<sub>H2</sub></th><th>D<sub>H2</sub><br>[m<sup>2</sup>/s]</th></tr></thead><tbody><tr><td colspan="6"><strong>Effect of H<sub>2</sub> force fields</strong></td></tr><tr><td>Marx_TIP4P/2005, T = 323K and P = 200 bar</td><td>-1.779</td><td>-</td><td>-</td><td>-1.779</td><td>6.48E-09</td></tr><tr><td>Vrabec_TIP4P/2005, T = 323K and P = 200 bar</td><td>-1.760</td><td>-</td><td>-</td><td>-1.760</td><td>6.75E-09</td></tr><tr><td>Yang & Zhong_TIP4P/2005, T = 323K and P = 200 bar</td><td>-1.675</td><td>-</td><td>-</td><td>-1.675</td><td>7.55E-09</td></tr><tr><td colspan="6"><strong>Effect of ion force fields</strong></td></tr><tr><td>Marx_TIP4P/2005_MadridTransport</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>m = 1, T = 323K and P = 200 bar</td><td>-1.752</td><td>-0.027</td><td>-0.039</td><td>-1.817</td><td>6.09E-09</td></tr><tr><td>m = 3, T = 323K and P = 200 bar</td><td>-1.701</td><td>-0.070</td><td>-0.114</td><td>-1.885</td><td>5.36E-09</td></tr><tr><td>m = 6, T = 323K and P = 200 bar</td><td>-1.617</td><td>-0.130</td><td>-0.215</td><td>-1.962</td><td>4.25E-09</td></tr><tr><td>Marx_TIP4P/2005_Madrid</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>m = 1, T = 323K and P = 200 bar</td><td>-1.759</td><td>-0.033</td><td>-0.040</td><td>-1.831</td><td>6.13E-09</td></tr><tr><td>m = 3, T = 323K and P = 200 bar</td><td>-1.718</td><td>-0.085</td><td>-0.116</td><td>-1.918</td><td>4.90E-09</td></tr><tr><td>m = 6, T = 323K and P = 200 bar</td><td>-1.644</td><td>-0.143</td><td>-0.232</td><td>-2.019</td><td>3.48E-09</td></tr><tr><td colspan="6"><strong>Effect of temperature</strong></td></tr><tr><td>Marx_TIP4P/2005_MadridTransport</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>m = 1, T = 373K and P = 200 bar</td><td>-1.524</td><td>-0.025</td><td>-0.035</td><td>-1.584</td><td>1.19E-08</td></tr><tr><td>m = 1, T = 423K and P = 200 bar</td><td>-1.344</td><td>-0.024</td><td>-0.033</td><td>-1.400</td><td>1.92E-08</td></tr><tr><td>m = 1, T = 473K and P = 200 bar</td><td>-1.186</td><td>-0.024</td><td>-0.030</td><td>-1.240</td><td>2.82E-08</td></tr><tr><td colspan="6"><strong>Effect of pressure</strong></td></tr><tr><td>Marx_TIP4P/2005_MadridTransport</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>m = 1, T = 323K and P = 600 bar</td><td>-1.815</td><td>-0.026</td><td>-0.040</td><td>-1.881</td><td>6.00E-09</td></tr><tr><td>m = 1, T = 323K and P = 1000 bar</td><td>-1.876</td><td>-0.026</td><td>-0.041</td><td>-1.943</td><td>5.74E-09</td></tr></tbody></table>

proportional to viscosity (according to the Stokes-Einstein relation), one can expect that a good description of brine viscosity should imply a good prediction of the diffusion coefficient. Fig. 5-c and 5-d illustrate the representation of the viscosity of NaCl brine (solvent) and the diffusion of $H_2$ (solute) by the two ions force fields, Madrid and MadridTransport, respectively. Up to 1 mol/kg, brine viscosity (and by extension $H_2$ diffusion) is well represented by both force fields. However, at higher salinities, the Madrid model leads to deviation from experimental viscosity data, whereas the MadridTransport model maintains its accuracy. As a result, the Marx-TIP4P/2005-MadridTransport combination emerges as the most consistent for representing $H_2$ diffusion in brine over the full range of temperature, pressure and salinity, making it suitable for various applications such as underground $H_2$ storage, native $H_2$ exploration, deep geothermal, etc. Therefore, the data generated for correlation in Fig. 6, are produced by this combination of force fields.

Several general observations can be made regarding the trends in $H_2$ diffusion with brine characteristics. For example, the decrease in the $H_2$ self-diffusion coefficient with increasing molality can be attributed to the enhanced brine structure (and resulting increase of viscosity) due to ion hydration. This trend is consistent with the variation in $s_{H_2}$, which becomes more negative as molality increases, see Table 1. Furthermore, a deeper analysis of the impact of molality on $H_2$ diffusion reveals that the contributions to $s_{H_2}$ change with molality. Results indicate that $s_{H_2}^v$ becomes less negative for water but more negative for ions when molality increases. This suggests that the decrease in the $H_2$ self-diffusion coefficient with increasing molality is also related to the more structured arrangement of ions around $H_2$.

In terms of temperature, as it increases, $s_{H_2}$ becomes less negative, indicating a weakening of the solvent structure around $H_2$, see Table 1. This leads to an increase in the $H_2$ diffusion with increasing temperature. Furthermore, the temperature dependence of $s_{H_2}$ is primarily driven by the contribution from the water component, $s_{H_2}^{H_2O}$, as shown in Table 1. This suggests that the temperature dependence of the $H_2$ self-diffusion coefficient in brine is also related to the changes in the structure of $H_2O$ around $H_2$ with temperature.

With respect to pressure, the small variation of the $H_2$ diffusion coefficient follows a pattern similar to that of $s_{H_2}$. More precisely, as pressure increases, the structure of the solvent around $H_2$ becomes slightly more ordered, leading to a small decrease in the $H_2$ self-diffusion coefficient, see Table 1. This suggests that the variation of the $H_2$ self-diffusion coefficient with pressure is also related to the structural changes of $H_2O$ around $H_2$.

## 4. $H_2$ diffusion modeling

In this last part, due to the lack of experimental data, extensive $H_2$ self-diffusion data in water and brine were generated by MD. These 150 points of data were generated for conditions ranging from 298 to 473 K, 1–1000 bars and 0–6 mol/kg of NaCl molality. All data are provided in Table 6 of the SM.

These data were then correlated to make this property more accessible for different applications. To this end, both a Stokes-Einstein-type model and an Arrhenius-type correlation were tested.

### 4.1. Stokes-Einstein model

Semi-empirical correlations based on the Stokes-Einstein (SE) theory have been extensively used for predicting the self-diffusivity of gases in liquids. By applying the SE hydrodynamic theory [27,69,70], it is possible to relate the diffusion coefficient of a rigid, spherical Brownian particle in a dilute solution directly to the fluid's shear viscosity $\eta$ as follows:

$$D = k_B T / (n_{SE} \pi \eta a) \qquad (11)$$

where, $k_B$ is the Boltzmann constant, $\eta$ is the solvent viscosity, and $a$ is
<figure><img src="image_11.png" /><figcaption>Fig. 5. Variation of H₂ diffusion with the NaCl molality using different FF combinations and simulation methods (EMD, NEMD) (a, b), variation of the brine viscosity (c) and diffusivity (d) with the NaCl molality using two different FF for NaCl.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 6. Compilation of MD simulation data of the H₂ diffusion in water/brine (a) at different temperatures and salinities at 200 bar, and (b) at different temperatures and pressure at one specific salinity (1 molNaCl/kgw).</figcaption></figure>

the hydrodynamic radius of the solute [71]. $n_{SE}$ is the Stokes–Einstein number, which depends on the boundary conditions at the particle/fluid interface. When assuming no-slip conditions (typically for large spherical molecules in low molecular weight solvents or suspended particles), $n_{SE}$ is equal to 6. For complete slip conditions, where the solute and the solvent are similar, $n_{SE}$ is equal to 4 [69].

However, deviations from the simple SE relationship may occur in

certain systems due to factors such as particle shape, specific interactions between particles, and non-ideal behavior of the solution. To address these deviations, various modified forms of the SE equation have been proposed, including the Wilke-Chang [72], Scheibel [73], Reddy-Doraiswamy [74], Lusis-Ratcliff [75] models. These semi-theoretical models can generally be expressed as:
$$D = A \frac{T}{\eta} \quad \text{or equivalently: } \ln(D) = \ln(A) + \ln\left(\frac{T}{\eta}\right) \qquad (12)$$

where A is the “intercept”, which may be a function of various parameters of the system.

While previous research has often focused on re-evaluating the parameter A in these semi-theoretical approaches based on SE, Magalhaes et al. [76] showed that in the case of real fluids, this assumption is not always valid. In such cases, the fractional Stokes-Einstein (FSE) relation should be applied as:

$$D \propto \left(\frac{T}{\eta}\right)^s \qquad (13)$$

or equivalently:

$$\ln(D) = \ln(A) + s \ln\left(\frac{T}{\eta}\right)$$

where the exponent s is determined from the slope of the corresponding log-log plot. When s = 1, the SE model holds true, while for s ≠ 1, the FSE relation applies. According to Tsimpanogiannis et al., 2023 [77], the SE relation does not adequately describe self-diffusion of hydrogen in pure water, suggesting that the FSE relation is more appropriate in this context.

## 4.2. Arrhenius-type correlation

An empirical Arrhenius-type correlation may be useful for rapid estimates of the H₂ diffusion, as it does not depend on other properties like viscosity, which is required in the case of SE and FSE. In addition, this correlation that takes into account the effect of temperature can be adapted [24] to account for the effect of pressure and salinity as well, using the following equation:

$$\ln D = \ln D_o + \frac{E_D}{RT} + k_s m_s + k_p P \qquad (14)$$

where D is in m²/s, T in K and P in bar. The parameters D₀, E₀, k_s and k_p, obtained by fitting on our data, are given in Table 2.

## 4.3. Modeling results

First, the validity of the SE relation was tested. The results shown in Fig. 7 indicate that the slope differs from 1, in agreement with the findings of Tsimpanogiannis et al., 2023 [77], and is equal to about 0.83. Therefore, the FSE relation was adopted in this work, utilizing only two fixed parameters (one slope and one intercept) for all temperatures, pressures and salinities (see Table 2). The viscosity of the solvent (water or NaCl brine) in the FSE model was calculated using the Mao and Duan (2009) [67] correlation. All data are provided in Table 6 of the SM.

Fig. 8 shows the results of the FSE modeling in pure water and brine at different temperatures, pressures and salinities. Up to 373k, the FSE model accurately captures the values of the H₂ diffusion coefficient as a function of salinity and pressure. Beyond this, the model continues to perform well at low salinities but starts to deviate for high salinities. Overall, the FSE model reproduces the data with an average absolute relative deviation of 5.6 %.

The results of the empirical Arrhenius-type correlation, which includes four adjustable parameters (see Table 2) linked to temperature,

<figure><img src="image_13.png" /><figcaption>Fig. 7. Check of the validity of the SE model.</figcaption></figure>

pressure and salinity dependencies, are also shown in Fig. 8. Interestingly, this Arrhenius-type correlation reproduces the H₂ diffusion data with a relative deviation of 3.9 %, representing a slight improvement over the FSE model, in particular at high temperatures and salinities. These results confirm that both the FSE and Arrhenius-type models can provide reasonable estimates of H₂ diffusion in water and brine across a wide range of temperature and salinities. For extrapolation purposes, the FSE model may be preferable due to its strong physical basis, which relates diffusivity to solvent viscosity and involves fewer parameters.

## 5. Conclusions

In this study, the diffusion coefficient of H₂ in water and NaCl brine was investigated using two simulation approaches: Equilibrium Molecular Dynamics and External-Field Non-Equilibrium Molecular Dynamics. Both methods yielded similar results, validating the protocols used. However, the latter method provided lower uncertainties for the same computational time and was subsequently chosen for the systematic computation of H₂ diffusion coefficients. Three force fields were tested for H₂: a 1-site model (Vrabec), a 2-site model (Yang&Zhong) and a 3-site model with charges (Marx). For NaCl brine, two force field combinations were tested, both using the same water force field (TIP4P/2005). One combination was designed for thermodynamic properties (Madrid), while the other was tailored to transport properties (MadridTransport).

The results show that, regardless of the H₂ model used, the predicted H₂ diffusion coefficients in pure water at different temperatures and pressures were qualitatively consistent. This is most likely due to the use of TIP4P/2005 force field for water, which accurately predicts the viscosity of water and its dependence on temperature. As diffusion coefficient is inversely proportional to viscosity, this is in line with the Stokes-Einstein relation.

Quantitatively, when compared with “reference” literature data, the 3-site H₂ force field with charges (Marx), is the best at predicting H₂ diffusion, closely followed by the 1-site model (Vrabec). The analysis of H₂ excess entropy revealed that the self-diffusion coefficient of H₂ in water is influenced as well by the local structure around the H₂ molecules. Both the Marx and Vrabec force fields yielded similar excess entropy values, while the Yong&Zhang force field produced a less negative

<table><caption>**Table 2**

Optimized model parameters.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="2">FSE model</th><th colspan="4">Arrhenius-type correlation</th></tr><tr><th>s</th><th>ln A</th><th>D<sub>0</sub></th><th>E<sub>D</sub></th><th>k<sub>s</sub></th><th>k<sub>p</sub></th></tr></thead><tbody><tr><td></td><td>0.8296</td><td>-24.105</td><td>7.55020E-07</td><td>-1.26738E+04</td><td>-7.06360E-02</td><td>-1.11910E-04</td></tr><tr><td>AAD (%)</td><td>5.6</td><td></td><td>3.9</td><td></td><td></td><td></td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>Fig. 8. Variation of the H₂ diffusivity at different temperature and salinity conditions: FSE model (solid lines); Arrhenius-type correlation (dashed lines).</figcaption></figure>

excess entropy, indicating a less structured water environment around H₂ and consequently, a higher self-diffusion coefficient.

Regarding the ion model, the MadridTransport model for NaCl proved to be the most consistent choice across the full range of temperatures, pressures, and salinities relevant to H₂ in underground applications. Its superiority lies in its ability to more accurately predict the evolution of brine viscosity with temperature, outperforming the Madrid model.

Finally, the best combination of force fields (Marx, TIP4P/2005 and MadridTransport) was selected to generate H₂ self-diffusion data in water and brine using NEMD simulations across a wide range of thermodynamic conditions (298 K–473 K, 1–1000 bar, and 0–6 mol/kg NaCl molality). These data were subsequently modeled using the fractional Stokes-Einstein relationship and correlated with an Arrhenius-type equation, showing a good ability to describe the H₂ diffusion in water and brine.

These correlations offer a reliable framework for predicting H₂ diffusion behavior under varying conditions, which is crucial for underground hydrogen storage and other related applications. However, experimental validation of some of these predictions would be highly beneficial, despite the significant challenges associated with conducting such experiments.

### CRediT authorship contribution statement

**Halla Kerkache:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation. **Hai Hoang:** Writing – review & editing, Software, Methodology, Investigation, Formal analysis. **Tu Khai Nam Nguyen:** Software, Methodology, Investigation. **Antoine Geoffroy-Neveux:** Writing – original draft, Software, Methodology, Investigation. **Carlos Nieto-Draghi:** Writing – review & editing, Software, Methodology, Investigation, Formal analysis. **Pierre Cézac:** Writing – review & editing, Project administration, Investigation, Funding acquisition. **Salaheddine Chabab:** Writing – review & editing, Writing – original draft, Supervision, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization. **Guillaume Galléro:** Writing – review & editing, Writing – original draft, Supervision, Project administration, Methodology, Investigation, Formal analysis,

Conceptualization.

### Declaration of interests

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.

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Galliero Guillaume reports financial support was provided by TotalEnergies SE. 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.

### Acknowledgments

Financial support from the "HYDR" project partners is gratefully acknowledged. The S.C. "HYDR" Chair and H.K. PhD grant are funded by the ISITE E2S program, supported by ProSim SA, ANR PIA, IFPEN and Région Nouvelle-Aquitaine. We are grateful to TOTALENERGIES S.A. for partial funding of the PhD project of H.K. and for in-depth discussions with Sabine Delahaye and Magali Pujol. The Institut Carnot "ISIFOR" is gratefully acknowledged for financing the H2Stock-2 Project. We also thank the Pau University and the MCIA for providing computational facilities.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.ijhydene.2025.04.157.

### References

[1] Rivard E, Trudeau M, Zaghib K. Hydrogen storage for mobility: a review. Materials 2019;12(12):1973.

[2] Blanco H, Faaij A. A review at the role of storage in energy systems with a focus on Power to Gas and long-term storage. Renew Sustain Energy Rev 2018;81:1049–86.

[3] Heinemann N, et al. Enabling large-scale hydrogen storage in porous media—the scientific challenges. Energy Environ Sci 2021;14(2):853–64.

[4] Cline G. Underground hydrogen storage. https://www.gaffneycline.com/; 2022.
[5] Benazzouk BK, et al. Diffusive transport of gases in saturated nanopores: caprock leakage from a molecular simulation perspective. J Nat Gas Sci Eng 2022;98: 104383.

[6] Wise DL, Houghton G. The diffusion coefficients of ten slightly soluble gases in water at 10-60 C. Chem Eng Sci 1966;21(11):999-1010.

[7] De Blok W, Fortuin J. Method for determining diffusion coefficients of slightly soluble gases in liquids. Chem Eng Sci 1981;36(10):1687-94.

[8] Verhallen P, et al. The diffusion coefficients of helium, hydrogen, oxygen and nitrogen in water determined from the permeability of a stagnant liquid layer in the quasi-s. Chem Eng Sci 1984;39(11):1535-41.

[9] Jähne B, Heinz G, Dietrich W. Measurement of the diffusion coefficients of sparingly soluble gases in water. J Geophys Res: Oceans 1987;92(C10):10767-76.

[10] Gertz K, Loeschcke H. Bestimmung der Diffusions-Koeffizienten von H2, O2, N2, und He in Wasser und Blutserum bei konstant gehaltener Konvektion. Z Naturforsch B Chem Sci 1954;9(1):1-9.

[11] Akgerman A, Gainer JL. Predicting gas-liquid diffusivities. J Chem Eng Data 1972; 17(3):372-7.

[12] Baird M, Davidson J. Annular jets—II: gas absorption. Chem Eng Sci 1962;17(6): 473-80.

[13] Wang S, et al. Diffusion coefficients of N2O and H2 in water at temperatures between 298.15 and 423.15 K with pressures up to 30 MPa. J Chem Eng Data 2023; 68(6):1313-9.

[14] Vivian JE, King CJ. Diffusivities of slightly soluble gases in water. AIChE J 1964;10 (2):220-1.

[15] Tsimpanogiannis IN, et al. Engineering model for predicting the intradiffusion coefficients of hydrogen and oxygen in vapor, liquid, and supercritical water based on molecular dynamics simulations. J Chem Eng Data 2021;66(8):3226-44.

[16] Piszko M, et al. Diffusivities in water or aqueous solutions of sodium chloride with dissolved hydrogen and methane by dynamic light scattering. Int J Hydrogen Energy 2024;68:499-510.

[17] Liu J, et al. Hydrogen diffusion in clay slit: implications for the geological storage. Energy & Fuels 2022;36(14):7651-60.

[18] Zhao X, Jin H. Correlation for self-diffusion coefficients of H2, CH4, CO, O2 and CO2 in supercritical water from molecular dynamics simulation. Appl Therm Eng 2020;171:114941.

[19] Zhao X, et al. Numerical study of H2, CH4, CO, O2 and CO2 diffusion in water near the critical point with molecular dynamics simulation. Comput Math Appl 2021; 81:759-71.

[20] Zhao X, Jin H. Investigation of hydrogen diffusion in supercritical water: a molecular dynamics simulation study. Int J Heat Mass Tran 2019;133:718-28.

[21] Zhao X, Luo T, Jin H. A predictive model for self-, Maxwell-Stefan, and Fick diffusion coefficients of binary supercritical water mixtures. J Mol Liq 2021;324: 114735.

[22] Kallikragas DT, Plugatyr AY, Svishchev IM. High temperature diffusion coefficients for O2, H2, and OH in water, and for pure water. J Chem Eng Data 2014;59(6): 1964-9.

[23] Ghasemi M, et al. Molecular dynamics simulation of hydrogen diffusion in water-saturated clay minerals; implications for Underground Hydrogen Storage (UHS). Int J Hydrogen Energy 2022;47(59):24871-85.

[24] van Rooijen W, et al. Interfacial tensions, solubilities, and transport properties of the H2/H2O/NaCl system: a molecular simulation study. J. Chem Eng Data 2023; 69(2):307-19.

[25] Kalati SS, et al. Molecular dynamics simulation of hydrogen diffusion into brine: implications for underground hydrogen storage. Int J Hydrogen Energy 2024;53: 17-28.

[26] Krishna R, Van Baten J. The darken relation for multicomponent diffusion in liquid mixtures of linear alkanes: an investigation using molecular dynamics (MD) simulations. Ind Eng Chem Res 2005;44(17):6939-47.

[27] Poling BE, Prauszinit JM, O'connell JP. The properties of gases and liquids, vol. 5. New York: Mcgraw-hill; 2001.

[28] Liu X, et al. Diffusion coefficients from molecular dynamics simulations in binary and ternary mixtures. Int J Thermophys 2013;34:1169-96.

[29] Samanta A, Ali SM, Ghosh SK. Universal scaling laws of diffusion in a binary fluid mixture. Phys Rev Lett 2001;87(24):245901.

[30] Frenkel D, Smit B. Understanding molecular simulation: from algorithms to applications. San Diego: Academic Press; 2002.

[31] Bian X, Kim C, Karniadakis GE. 111 years of Brownian motion. Soft Matter 2016;12 (30):6331-46.

[32] Chempath S, Krishna R, Snurr RQ. Nonequilibrium molecular dynamics simulations of diffusion of binary mixtures containing short n-alkanes in faujasite. J Phys Chem B 2004;108(35):13481-91.

[33] Krishna R, Wesselingh JA. The Maxwell-Stefan approach to mass transfer. Chemical engineering science 1997;52(6):861-911.

[34] Krishna R, van Baten JM. Describing diffusion in fluid mixtures at elevated pressures by combining the Maxwell-Stefan formulation with an equation of state. Chem Eng Sci 2016;153:174-87.

[35] Schnell SK, et al. Calculating thermodynamic properties from fluctuations at small scales. J Phys Chem B 2011;115(37):10911-8.

[36] Dawass N, et al. Kirkwood-Buff integrals from molecular simulation. Fluid Phase Equilib 2019;486:21-36.

[37] Schnell SK, et al. How to apply the Kirkwood-Buff theory to individual species in salt solutions. Chem Phys Lett 2013;582:154-7.

[38] Balaji SP, Schnell SK, Vlugt TJ. Calculating thermodynamic factors of ternary and multicomponent mixtures using the Permuted Widom test particle insertion method. Theor Chem Acc 2013;132:1-8.

[39] Balaji SP, et al. A direct method for calculating thermodynamic factors for liquid mixtures using the Permuted Widom test particle insertion method. Mol Phys 2013; 111(2):287-96.

[40] Hulikal Chakrapani T, et al. Calculating thermodynamic factors for diffusion using the continuous fractional component Monte Carlo method. J Chem Theor Comput 2023;20(1):333-47.

[41] Vega C, et al. What ice can teach us about water interactions: a critical comparison of the performance of different water models. Faraday discussions 2009;141: 251-76.

[42] Kerkache H, et al. The solubility of H2 in NaCl brine at high pressures and high temperatures: molecular simulation study and thermodynamic modeling. J Mol Liq 2024;124497.

[43] Benavides A, et al. A potential model for sodium chloride solutions based on the TIP4P/2005 water model. The Journal of chemical physics 2017;147(10).

[44] Zeron I, Abascal J, Vega C. A force field of Li+, Na+, K+, Mg2+, Ca2+, Cl-, and SO42- in aqueous solution based on the TIP4P/2005 water model and scaled charges for the ions. The Journal of chemical physics 2019;151(13).

[45] Blazquez S, Conde M, Vega C. Scaled charges for ions: an improvement but not the final word for modeling electrolytes in water. J Chem Phys 2023;158(5).

[46] Marx D, Nielaba P. Path-integral Monte Carlo techniques for rotational motion in two dimensions: quenched, annealed, and no-spin quantum-statistical averages. Phys Rev 1992;45(12):8968.

[47] Gunn R, Chueh P, Prauszinit J. Prediction of thermodynamic properties of dense gas mixtures containing one or more of the quantum gases. AIChE J 1966;12(5): 937-41.

[48] Koster A, Thol M, Vrabec J. Molecular models for the hydrogen age: hydrogen, nitrogen, oxygen, argon, and water. J Chem Eng Data 2018;63(2):305-20.

[49] Yang Q, Zhong C. Molecular simulation of adsorption and diffusion of hydrogen in metal- organic frameworks. J Phys Chem B 2005;109(24):11862-4.

[50] Rahbari A, et al. Solubility of water in hydrogen at high pressures: a molecular simulation study. J Chem Eng Data 2019;64(9):4103-15.

[51] Lopez-Lazaro C, et al. Predicting the phase behavior of hydrogen in NaCl brines by molecular simulation for geological applications. BSGF-Earth Sciences Bulletin 2019;190(1):7.

[52] Lorentz HA. Ueber die Anwendung des Satzes vom Virial in der kinetischen Theorie der Gase. Ann Phys 1881;248(1):127-36.

[53] Berthelot D. Sur le mélange des gaz. Compt. Rendus 1898;126(3):15.

[54] Plimpton S. Fast parallel algorithms for short-range molecular dynamics. J Comput Phys 1995;117(1):1-19.

[55] Martin MG. MCCCS Towhee: a tool for Monte Carlo molecular simulation. Mol Simul 2013;39(14-15):1212-22.

[56] Hoang H, et al. On elemental and isotopic fractionation of noble gases in geological fluids by molecular diffusion. Geochem Cosmochim Acta 2021;315:172-84.

[57] Hoang H, et al. Elemental and isotopic fractionation of noble gases in gas and oil under reservoir conditions: impact of thermodiffusion*. The European Physical Journal E 2019;42:1-10.

[58] Jamali SH, et al. OCTP: a tool for on-the-fly calculation of transport properties of fluids with the order-n algorithm in LAMMPS. J Chem Inf Model 2019;59(4): 1290-4.

[59] De Groot SR, Mazur P. Non-equilibrium thermodynamics. Dover Publications; 1984.

[60] Dyre JC. Perspective: excess-entropy scaling. The Journal of chemical physics 2018;149(21).

[61] Tuckerman ME. Statistical mechanics: theory and molecular simulation. Oxford university press; 2023.

[62] Nosé S. A molecular dynamics method for simulations in the canonical ensemble. Mol Phys 1984;52(2):255-68.

[63] J Evans D, Morriss GP. Statistical mechanics of nonequilibrium liquids. ANU Press; 2007.

[64] Andersen HC. Rattle: a "velocity" version of the shake algorithm for molecular dynamics calculations. J Comput Phys 1983;52(1):24-34.

[65] Luty BA, van Gunsteren WF. Calculating electrostatic interactions using the particle-particle-mesh method with nonperiodic long-range interactions. The Journal of Physical Chemistry 1996;100(7):2581-7.

[66] Yeh I-C, Hummer G. System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions. J Phys Chem B 2004;108(40):15873-9.

[67] Mao S, Duan Z. The viscosity of aqueous alkali-chloride solutions up to 623 K, 1,000 bar, and high ionic strength. Int J Thermophys 2009;30:1510-23.

[68] Al Ghafri S, Maitland GC, Trusler JM. Densities of aqueous MgCl2 (aq), CaCl2 (aq), KI (aq), NaCl (aq), KCl (aq), AlCl3 (aq), and (0.964 NaCl+ 0.136 KCl)(aq) at temperatures between (283 and 472) K, pressures up to 68.5 MPa, and molalities up to 6 mol·kg⁻¹. J Chem Eng Data 2012;57(4):1288-304.

[69] Bird, R.B., W.E. Stewai, and E.N. Lightfoot, Phenomena second ed.].

[70] Gainer J, Metzner A. Diffusion in liquids-theoretical analysis and experimental verification. In: AIChE-I. Chem. E. Symp. Ser.; 1965.

[71] Tarjus G, Kivelson D. Breakdown of the Stokes-Einstein relation in supercooled liquids. The Journal of chemical physics 1995;103(8):3071-3.

[72] Wilke C, Chang P. Correlation of diffusion coefficients in dilute solutions. AIChE J 1955;1(2):264-70.

[73] Scheibel EG. Correspondence. Liquid diffusivities. Viscosity of gases. Ind Eng Chem 1954;46(9):2007-8.

[74] Reddy K, Doraiswamy L. Estimating liquid diffusivity. Ind Eng Chem Fundam 1967;6(1):77-9.
[75] Lusis M, Ratcliff C. Diffusion in binary liquid mixtures at infinite dilution. Can J Chem Eng 1968;46(5):385-7.

[76] Magalhães AL, et al. Simple and accurate correlations for diffusion coefficients of solutes in liquids and supercritical fluids over wide ranges of temperature and density. J Supercrit Fluids 2013;76:94-114.

[77] Tsimpanogiannis IN, Moultos OA. Is Stokes-Einstein relation valid for the description of intra-diffusivity of hydrogen and oxygen in liquid water? Fluid Phase Equilib 2023;563:113568.