# The solubility of H2 in NaCl brine at high pressures and high temperatures: Molecular simulation study and thermodynamic modeling

Halla Kerkache, Hai Hoang, Pierre Cézac, Guillaume Galliéro, Salaheddine Chabab

► To cite this version:

Halla Kerkache, Hai Hoang, Pierre Cézac, Guillaume Galliéro, Salaheddine Chabab. The solubility of H2 in NaCl brine at high pressures and high temperatures: Molecular simulation study and thermodynamic modeling. Journal of Molecular Liquids, 2024, 400, pp.124497. 10.1016/j.molliq.2024.124497. hal-04623890

HAL Id: hal-04623890

https://univ-pau.hal.science/hal-04623890

Submitted on 25 Jun 2024

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
# The solubility of H2 in NaCl brine at high pressures and high temperatures: Molecular simulation study and thermodynamic modeling

Halla Kerkache, Hai Hoang, Pierre Cézac, Guillaume Galliéro, Salaheddine Chabab

► To cite this version:

Halla Kerkache, Hai Hoang, Pierre Cézac, Guillaume Galliéro, Salaheddine Chabab. The solubility of H2 in NaCl brine at high pressures and high temperatures: Molecular simulation study and thermodynamic modeling. Journal of Molecular Liquids, 2024, 400, pp.124497. 10.1016/j.molliq.2024.124497. hal-04623890

HAL Id: hal-04623890

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

Submitted on 25 Jun 2024

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>

# The solubility of H₂ in NaCl brine at high pressures and high temperatures:
Molecular simulation study and thermodynamic modeling

Halla Kerkache<sup>a</sup>, Hai Hoang<sup>c,d</sup>, Pierre Cézac<sup>b</sup>, Guillaume Galliéro<sup>a,*</sup>, Salaheddine Chabab<sup>a,b,*</sup>

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

## ARTICLE INFO

### Keywords:

H₂ solubility in brine

Underground hydrogen storage

H₂ natural emissions

Monte Carlo simulation

Molecular dynamics

Thermodynamic modeling

## ABSTRACT

To address the energy demands of a growing population and mitigate carbon emissions, it is imperative to transition from fossil fuels to renewable energy sources. However, the intermittency of renewable energies poses a significant challenge. To address this issue, deep saline aquifers have emerged as a viable option for large-scale energy storage, particularly through hydrogen (H₂) storage post Power-to-Gas process. Moreover, natural H₂ emissions have been documented worldwide, and the potential for underground accumulations presents promising zero-carbon energy sources. However, in these different contexts, the interaction between gas, brine, and rock can lead to physico-chemical and biochemical phenomena that can directly impact the mobility and stability of H₂. Therefore, understanding the thermophysical behavior of the involved fluids is essential for the development Underground Hydrogen Storage in porous media and for exploring natural H₂ reserves. However, despite recent advancements, there is still a lack of experimental data on thermophysical properties of hydrogen in contact with brine. This study investigates the equilibrium of the H₂/brine system using Continuous Fractional Component Monte Carlo molecular simulation through two methods: the Gibbs ensemble method and the isobaric-isothermal simulation based on Henry's law. Different force fields for H₂, water and salt (NaCl) ions were evaluated. Through a comparative analysis, two model combinations, Marx-TIP4P/EP-KBF and Marx-TIP4P/2005-Madrid, were identified as providing the most accurate results, albeit requiring a constant binary interaction for enhanced H₂ solubility quantification in brine. After adjustment to some limited experimental data from literature, the simulations were extended to higher temperatures (up to 453 K), pressures (up to 100 MPa), and NaCl salinities (up to 6 mol/kgw). Finally, the newly generated data facilitated the refinement of a thermodynamic model using the Krichevsky and Kasarnovsky approach, improving estimations of H₂ dissolution losses, caprock sealing capacity, and insights into natural H₂ production and accumulation underground.

## 1. Introduction

Reducing CO₂ emissions into the atmosphere is considered one of the most important environmental policies to combat global warming. The use of renewable energies is the primary solution for the decarbonization of energy systems, which are the major source of anthropogenic CO₂ emissions. However, the intermittency, seasonal fluctuations, and geographical constraints associated with these energies lead to instability between energy supply and demand [1]. To manage this instability, excess electricity can be stored in different forms, depending on capacity and storage time, such as flywheels, batteries, compressed air energy storage (CAES), and pumped hydropower stations [2]. However,

these electricity storage technologies are not suitable for storing large amounts of electricity (~MWh) over long periods. Therefore, the combination of power-to-gas (PtG) and underground hydrogen storage (UHS) appears to be a promising option for large-scale (~GWh) energy storage [3].

For more than 30 years, high-purity hydrogen storage in salt caverns has been successfully practiced by the chemical sector. Currently, there are four storage sites in operation, including a site in the UK (Teesside) and three sites on the US Gulf Coast, at the Chevron Phillips Clemens Terminal in Texas (Clemens Dome, Moss Bluff, Spindletop) [4]. On the other hand, for energy purposes, several projects have recently emerged (HyPSTER, HyUnder, HyStock, Rostock H, STOPIL H2, etc.) [5]. However, to reach GWh energy storage, other types of underground storage

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

<footnote><sup>b</sup> Université de Pau et des Pays de l'Adour, Laboratoire de Thermique, Energetique et Procédés (LaTEP), Rue Jules Ferry BP 7511, 64075 Pau Cedex, 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>* Corresponding authors at: Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LFCR, UMR5150, Pau, France (S. Chahab). E-mail addresses: guillaume.galliero@univ-pau.fr (G. Galliéro), salaheddine.chabab@univ-pau.fr (S. Chahab).</footnote>
<figure><img src="image_2.png" /></figure>

# The solubility of H₂ in NaCl brine at high pressures and high temperatures:
Molecular simulation study and thermodynamic modeling

Halla Kerkache<sup>a</sup>, Hai Hoang<sup>c,d</sup>, Pierre Cézac<sup>b</sup>, Guillaume Galliéro<sup>a,*</sup>, Salaheddine Chabab<sup>a,b,*</sup>

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

## ARTICLE INFO

### Keywords:

H₂ solubility in brine

Underground hydrogen storage

H₂ natural emissions

Monte Carlo simulation

Molecular dynamics

Thermodynamic modeling

## ABSTRACT

To address the energy demands of a growing population and mitigate carbon emissions, it is imperative to transition from fossil fuels to renewable energy sources. However, the intermittency of renewable energies poses a significant challenge. To address this issue, deep saline aquifers have emerged as a viable option for large-scale energy storage, particularly through hydrogen (H₂) storage post Power-to-Gas process. Moreover, natural H₂ emissions have been documented worldwide, and the potential for underground accumulations presents promising zero-carbon energy sources. However, in these different contexts, the interaction between gas, brine, and rock can lead to physico-chemical and biochemical phenomena that can directly impact the mobility and stability of H₂. Therefore, understanding the thermophysical behavior of the involved fluids is essential for the development Underground Hydrogen Storage in porous media and for exploring natural H₂ reserves. However, despite recent advancements, there is still a lack of experimental data on thermophysical properties of hydrogen in contact with brine. This study investigates the equilibrium of the H₂/brine system using Continuous Fractional Component Monte Carlo molecular simulation through two methods: the Gibbs ensemble method and the isobaric-isothermal simulation based on Henry's law. Different force fields for H₂, water and salt (NaCl) ions were evaluated. Through a comparative analysis, two model combinations, Marx-TIP4P/EP-KBF and Marx-TIP4P/2005-Madrid, were identified as providing the most accurate results, albeit requiring a constant binary interaction for enhanced H₂ solubility quantification in brine. After adjustment to some limited experimental data from literature, the simulations were extended to higher temperatures (up to 453 K), pressures (up to 100 MPa), and NaCl salinities (up to 6 mol/kgw). Finally, the newly generated data facilitated the refinement of a thermodynamic model using the Krichevsky and Kasarnovsky approach, improving estimations of H₂ dissolution losses, caprock sealing capacity, and insights into natural H₂ production and accumulation underground.

## 1. Introduction

Reducing CO₂ emissions into the atmosphere is considered one of the most important environmental policies to combat global warming. The use of renewable energies is the primary solution for the decarbonization of energy systems, which are the major source of anthropogenic CO₂ emissions. However, the intermittency, seasonal fluctuations, and geographical constraints associated with these energies lead to instability between energy supply and demand [1]. To manage this instability, excess electricity can be stored in different forms, depending on capacity and storage time, such as flywheels, batteries, compressed air energy storage (CAES), and pumped hydropower stations [2]. However,

these electricity storage technologies are not suitable for storing large amounts of electricity (~MWh) over long periods. Therefore, the combination of power-to-gas (PtG) and underground hydrogen storage (UHS) appears to be a promising option for large-scale (~GWh) energy storage [3].

For more than 30 years, high-purity hydrogen storage in salt caverns has been successfully practiced by the chemical sector. Currently, there are four storage sites in operation, including a site in the UK (Teesside) and three sites on the US Gulf Coast, at the Chevron Phillips Clemens Terminal in Texas (Clemens Dome, Moss Bluff, Spindletop) [4]. On the other hand, for energy purposes, several projects have recently emerged (HyPSTER, HyUnder, HyStock, Rostock H, STOPIL H2, etc.) [5]. However, to reach GWh energy storage, other types of underground storage

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

<footnote><sup>b</sup> Université de Pau et des Pays de l'Adour, Laboratoire de Thermique, Energetique et Procédés (LaTEP), Rue Jules Ferry BP 7511, 64075 Pau Cedex, 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>* Corresponding authors at: Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, LFCR, UMR5150, Pau, France (S. Chahab). E-mail addresses: guillaume.galliero@univ-pau.fr (G. Galliéro), salaheddine.chabab@univ-pau.fr (S. Chahab).</footnote>
Nomenclature

g radial distribution function

H Henry's law constant (MPa)

k binary interaction parameter

K salting-out effect coefficient (kg_w/mol)

$k_B$ Boltzmann constant ($k_B = 1.380649 \times 10^{-23} \text{ m}^2 \cdot \text{kg} \cdot \text{s}^{-2} \cdot \text{K}^{-1}$)

m molality (mol/kg_w)

n number of moles (mol)

N number of molecules

P pressure (MPa)

p Boltzmann probability distribution

R ideal gas constant ($R = 8.314462618 \text{ kg} \cdot \text{m}^2 \cdot \text{s}^{-2} \cdot \text{K}^{-1} \cdot \text{mol}^{-1}$)

r distance

T temperature (K)

v molar volume (m³/mol)

V volume (m³)

x liquid molar fraction

x' salt-free liquid molar fraction

y vapor molar fraction

Z compressibility factor

Greek letters

ρ density (kg/m³) or number density

μ chemical potential

φ fugacity coefficient

Π Poynting factor

λ coupling parameter for CFCMC sampling method

δ Dirac function

γ activity coefficient

ε cross interaction energy

Subscripts

i, j compound or phase

min minimum

max maximum

s salt

w water

Superscripts

sat saturation

∞ infinite dilution

Abbreviations

CFC Continuous Fractional Component

EoS Equation of State

Exp Experimental

FF Force Fields

GE Gibbs Ensemble

HKF Helgeson-Kirkham-Flowers EoS

IUPAC International Union of Pure and Applied Chemistry

KBF Force fields based on Kirkwood-Buff integrals

KK Krichevsky and Kasarnovsky model

LB Lorentz-Berthelot combining rules

MC Monte Carlo

MD Molecular Dynamics

MS Molecular Simulation

NIST National Institute of Standards and Technology

RDF Radial Distribution Function

REFPROP NIST Reference Fluid Thermodynamic and Transport Properties Database

UHS Underground Hydrogen Storage

VLE Vapor-Liquid Equilibria

should be envisaged. Among them, saline aquifers are one of the most attractive types of storage due to: (i) their large storage capacity (100 million standard cubic meters) compared to salt caverns (hundreds of thousands of standard cubic meters); (ii) their long-term storage capacity (months) providing sufficient storage for seasonal energy demands, in contrast to the short-to-medium-term storage capacity of salt caverns [3]; and (iii) their geographic availability (located in all sedimentary basins around the world) compared to salt caverns, and depleted gas reservoirs [6]. Using this type of reservoir would be economically viable and environmentally friendly, as it does not require the injection of fresh water or the release of brine [7]. However, experiences and data on the underground storage of hydrogen in porous geological formations (depleted gas fields and saline aquifers) are limited to the historical storage of town gas – a mixture primarily composed of hydrogen (50–60 %), methane (10–33 %), nitrogen (<30 %), and carbon monoxide and dioxide (12–20 %) – as found in France (Beynes), Czechoslovakia (Lobodice), and Germany (Engelborstel, Bad Lauchstaedt, Kiel) [8].

In another context, many locations around the world have exhibited underground hydrogen emanations [9]. This hydrogen, known as natural hydrogen (or native hydrogen), is produced by the earth and is considered as a new potential source of low-carbon energy. Therefore, exploring the potential geological accumulation of native hydrogen offers an additional pathway for zero-carbon energy sources [10,11]. However, many open questions remain regarding its origin, its availability, its underground behavior and its possible subsurface accumulation [11].

Considering the limited knowledge in these different contexts (UHS and native H₂), it is important to study the different phenomena that can

have an impact on the behavior of H₂ along its migration path in this type of underground environment. The contact between gas, brine and rock can lead to a series of physico-chemical and bio-chemical processes that can directly impact the mobility and stability of hydrogen (diffusion/diffusion, reactivity, etc.). Among others, the knowledge of thermophysical properties such as gas/brine solubility is required for the study and development of such infrastructure [1]. These properties are required to determine the impact of bio-geo-physical-chemical phenomena on the efficiency of UHS, and to understand the processes involved in the generation and possible accumulation of native H₂.

For reasons of complexity (very low solubility of H₂ in brine), cost and hazardousness (due to high pressure, high flammability of H₂, and the presence of salt which can cause risk of corrosion and thus leakage), there are very few studies of phase equilibria of H₂/brine systems [5,12–15]. Furthermore, they sometimes exhibit non-negligible differences in terms of results. Therefore, further experimental measurements are needed to evaluate the existing data, validate and/or adjust the different models used in classical and statistical thermodynamics, as recently partially covered, from an experimental perspective in our group [15]. Moreover, experimental measurements would benefit from being complemented by additional approaches such as molecular simulations, which is the main purpose of this work (Fig. 1), in order to: i) Better understand the different physical mechanisms occurring; ii) Produce vapor-liquid equilibrium data under severe conditions (very high temperature, pressure and salinity conditions); iii) Facilitate the development of consistent and physically based thermodynamic models.

Classical Monte Carlo molecular simulation has become a popular tool for predicting the fluid phase equilibria of systems of scientific and industrial interest [16]. It is based on a stochastic approach, where force
Nomenclature

g radial distribution function

H Henry's law constant (MPa)

k binary interaction parameter

K salting-out effect coefficient (kg_w/mol)

$k_B$ Boltzmann constant ($k_B = 1.380649 \times 10^{-23} \text{ m}^2 \cdot \text{kg} \cdot \text{s}^{-2} \cdot \text{K}^{-1}$)

m molality (mol/kg_w)

n number of moles (mol)

N number of molecules

P pressure (MPa)

p Boltzmann probability distribution

R ideal gas constant ($R = 8.314462618 \text{ kg} \cdot \text{m}^2 \cdot \text{s}^{-2} \cdot \text{K}^{-1} \cdot \text{mol}^{-1}$)

r distance

T temperature (K)

v molar volume (m³/mol)

V volume (m³)

x liquid molar fraction

x' salt-free liquid molar fraction

y vapor molar fraction

Z compressibility factor

Greek letters

ρ density (kg/m³) or number density

μ chemical potential

φ fugacity coefficient

Π Poynting factor

λ coupling parameter for CFCMC sampling method

δ Dirac function

γ activity coefficient

ε cross interaction energy

Subscripts

i, j compound or phase

min minimum

max maximum

s salt

w water

Superscripts

sat saturation

$\infty$ infinite dilution

Abbreviations

CFC Continuous Fractional Component

EoS Equation of State

Exp Experimental

FF Force Fields

GE Gibbs Ensemble

HKF Helgeson-Kirkham-Flowers EoS

IUPAC International Union of Pure and Applied Chemistry

KBF Force fields based on Kirkwood-Buff integrals

KK Krichevsky and Kasarnovsky model

LB Lorentz-Berthelot combining rules

MC Monte Carlo

MD Molecular Dynamics

MS Molecular Simulation

NIST National Institute of Standards and Technology

RDF Radial Distribution Function

REFPROP NIST Reference Fluid Thermodynamic and Transport Properties Database

UHS Underground Hydrogen Storage

VLE Vapor-Liquid Equilibria

should be envisaged. Among them, saline aquifers are one of the most attractive types of storage due to: (i) their large storage capacity (100 million standard cubic meters) compared to salt caverns (hundreds of thousands of standard cubic meters); (ii) their long-term storage capacity (months) providing sufficient storage for seasonal energy demands, in contrast to the short-to-medium-term storage capacity of salt caverns [3]; and (iii) their geographic availability (located in all sedimentary basins around the world) compared to salt caverns, and depleted gas reservoirs [6]. Using this type of reservoir would be economically viable and environmentally friendly, as it does not require the injection of fresh water or the release of brine [7]. However, experiences and data on the underground storage of hydrogen in porous geological formations (depleted gas fields and saline aquifers) are limited to the historical storage of town gas – a mixture primarily composed of hydrogen (50–60 %), methane (10–33 %), nitrogen (<30 %), and carbon monoxide and dioxide (12–20 %) – as found in France (Beynes), Czechoslovakia (Lobodice), and Germany (Engelborstel, Bad Lauchstaedt, Kiel) [8].

In another context, many locations around the world have exhibited underground hydrogen emanations [9]. This hydrogen, known as natural hydrogen (or native hydrogen), is produced by the earth and is considered as a new potential source of low-carbon energy. Therefore, exploring the potential geological accumulation of native hydrogen offers an additional pathway for zero-carbon energy sources [10,11]. However, many open questions remain regarding its origin, its availability, its underground behavior and its possible subsurface accumulation [11].

Considering the limited knowledge in these different contexts (UHS and native H₂), it is important to study the different phenomena that can

have an impact on the behavior of H₂ along its migration path in this type of underground environment. The contact between gas, brine and rock can lead to a series of physico-chemical and bio-chemical processes that can directly impact the mobility and stability of hydrogen (diffusion/diffusion, reactivity, etc.). Among others, the knowledge of thermophysical properties such as gas/brine solubility is required for the study and development of such infrastructure [1]. These properties are required to determine the impact of bio-geo-physical-chemical phenomena on the efficiency of UHS, and to understand the processes involved in the generation and possible accumulation of native H₂.

For reasons of complexity (very low solubility of H₂ in brine), cost and hazardousness (due to high pressure, high flammability of H₂, and the presence of salt which can cause risk of corrosion and thus leakage), there are very few studies of phase equilibria of H₂/brine systems [5,12–15]. Furthermore, they sometimes exhibit non-negligible differences in terms of results. Therefore, further experimental measurements are needed to evaluate the existing data, validate and/or adjust the different models used in classical and statistical thermodynamics, as recently partially covered, from an experimental perspective in our group [15]. Moreover, experimental measurements would benefit from being complemented by additional approaches such as molecular simulations, which is the main purpose of this work (Fig. 1), in order to: i) Better understand the different physical mechanisms occurring; ii) Produce vapor-liquid equilibrium data under severe conditions (very high temperature, pressure and salinity conditions); iii) Facilitate the development of consistent and physically based thermodynamic models.

Classical Monte Carlo molecular simulation has become a popular tool for predicting the fluid phase equilibria of systems of scientific and industrial interest [16]. It is based on a stochastic approach, where force
<figure><img src="image_4.png" /><figcaption>Fig. 1. The overall research approach of the study.</figcaption></figure>

fields are used to describe the molecular model and the molecular in-
teractions. The evolution of the system to reach the equilibrium is
achieved thanks to statistical mechanics theories (e.g. the Metropolis
method). This approach has been used to calculate the solubility of
various gases (CH₄, CO₂, H₂S, SO₂, noble gases, etc.) in brine [17–21].
However, its application to the computation of H₂ solubility in brine is
rather limited and relatively recent [13,14], probably due to the low
solubility value leading to statistical difficulties.

In estimating gas solubility in solvents, two Monte-Carlo molecular
simulation methods have been widely employed in the literature. The
first method calculates the molar fraction of gas in a solvent directly by
establishing equilibrium conditions between gas and liquid phases.
However, this method is computationally intensive and yields results
with high statistical uncertainties at low gas solubilities. To address
these limitations when dealing with H₂/brine systems, Rahbari et al.
[22,23] recently proposed a method based on the computation of
chemical potentials dependencies on the compositions.

The second method relies on estimation of Henry's law constant
through Monte Carlo simulations of the saturated liquid phase, along
with technique for calculating the chemical potential of gas in the sol-
vent under infinitely dilute conditions, such as Widom's Test Particle
Insertion or the Continuous Fractional Component (CFC) method
[24,25]. The gas molar fraction in the solvent is then indirectly deter-
mined from the Henry's law constant in combination with the thermo-
dynamics equations. Lopez-Lazaro et al. 2019 [13] applied this indirect
method to calculate the Henry's constant of H₂ in brine and calibrated
an equation of state (Søreide and Whitson 1992 [26]) to predict solu-
bilities at high pressures. van Rooijen et al. 2023 [14] also used this
method to estimate H₂ solubility for pressures up to 10 MPa (beyond

that pressure Henry's law does not apply without correction). In this
work, both methods were employed to calculate the H₂ solubility in
water and brine.

This article is structured as follows. First, an overview of the avail-
able vapor liquid equilibrium properties data for H₂/NaCl brine and H₂
solubility correlations is presented (section 2). Then, equilibrium data
for H₂/brine systems using Monte Carlo (MC) simulations are provided,
with a detailed analysis of the impact of the chosen force fields on sol-
ubility predictions (section 3). Finally, H₂ solubility data in brine (H₂O
+ NaCl) are used to refine an easy to use model applicable across a wide
range of typical underground thermodynamic conditions (section 4).

## 2. Review of available H₂ solubility data in water and brine

In literature, solubility data are given in different forms such as
Bunsen, Kuenen, Ostwald coefficients, molar composition, molality,
molarity or Henry constant. In this work all collected data were con-
verted to Henry's constant (Ln H) using the IUPAC guidelines [27,28] for
low-pressure data and the Krichevsky and Kasarnovsky approach (see
sections 3 and 4 for technical details) for high-pressure data [29].

In Fig. 2, the experimental data of H₂ Henry constants in pure water
taken from the literature (see Table 1) are illustrated along with some
existing correlations. The chosen correlations are: (1) the formulation
recommended by Harvey (1996) [30] which was further optimized by
Fernández-Prini (2003) [31]; (2) the hydrogen Henry's constant in
oxygenated solvents developed by Trinh et al. (2016) [32] that was
optimized by Lopez Lazaro et al. (2019) [13] for the hydrogen Henry's
constant in water; (3) the Akinfiev and Diamond (2002) [33] approach
which relied on a virial-like equation of state (EoS) for describing the
<figure><img src="image_4.png" /><figcaption>Fig. 1. The overall research approach of the study.</figcaption></figure>

fields are used to describe the molecular model and the molecular in-
teractions. The evolution of the system to reach the equilibrium is
achieved thanks to statistical mechanics theories (e.g. the Metropolis
method). This approach has been used to calculate the solubility of
various gases (CH₄, CO₂, H₂S, SO₂, noble gases, etc.) in brine [17–21].
However, its application to the computation of H₂ solubility in brine is
rather limited and relatively recent [13,14], probably due to the low
solubility value leading to statistical difficulties.

In estimating gas solubility in solvents, two Monte-Carlo molecular
simulation methods have been widely employed in the literature. The
first method calculates the molar fraction of gas in a solvent directly by
establishing equilibrium conditions between gas and liquid phases.
However, this method is computationally intensive and yields results
with high statistical uncertainties at low gas solubilities. To address
these limitations when dealing with H₂/brine systems, Rahbari et al.
[22,23] recently proposed a method based on the computation of
chemical potentials dependencies on the compositions.

The second method relies on estimation of Henry's law constant
through Monte Carlo simulations of the saturated liquid phase, along
with technique for calculating the chemical potential of gas in the sol-
vent under infinitely dilute conditions, such as Widom's Test Particle
Insertion or the Continuous Fractional Component (CFC) method
[24,25]. The gas molar fraction in the solvent is then indirectly deter-
mined from the Henry's law constant in combination with the thermo-
dynamics equations. Lopez-Lazaro et al. 2019 [13] applied this indirect
method to calculate the Henry's constant of H₂ in brine and calibrated
an equation of state (Søreide and Whitson 1992 [26]) to predict solu-
bilities at high pressures. van Rooijen et al. 2023 [14] also used this
method to estimate H₂ solubility for pressures up to 10 MPa (beyond

that pressure Henry's law does not apply without correction). In this
work, both methods were employed to calculate the H₂ solubility in
water and brine.

This article is structured as follows. First, an overview of the avail-
able vapor liquid equilibrium properties data for H₂/NaCl brine and H₂
solubility correlations is presented (section 2). Then, equilibrium data
for H₂/brine systems using Monte Carlo (MC) simulations are provided,
with a detailed analysis of the impact of the chosen force fields on sol-
ubility predictions (section 3). Finally, H₂ solubility data in brine (H₂O
+ NaCl) are used to refine an easy to use model applicable across a wide
range of typical underground thermodynamic conditions (section 4).

## 2. Review of available H₂ solubility data in water and brine

In literature, solubility data are given in different forms such as
Bunsen, Kuenen, Ostwald coefficients, molar composition, molality,
molarity or Henry constant. In this work all collected data were con-
verted to Henry's constant (Ln H) using the IUPAC guidelines [27,28] for
low-pressure data and the Krichevsky and Kasarnovsky approach (see
sections 3 and 4 for technical details) for high-pressure data [29].

In Fig. 2, the experimental data of H₂ Henry constants in pure water
taken from the literature (see Table 1) are illustrated along with some
existing correlations. The chosen correlations are: (1) the formulation
recommended by Harvey (1996) [30] which was further optimized by
Fernández-Prini (2003) [31]; (2) the hydrogen Henry's constant in
oxygenated solvents developed by Trinh et al. (2016) [32] that was
optimized by Lopez Lazaro et al. (2019) [13] for the hydrogen Henry's
constant in water; (3) the Akinfiev and Diamond (2002) [33] approach
which relied on a virial-like equation of state (EoS) for describing the
<figure><img src="image_5.png" /><figcaption>Fig. 2. Experimental data (symbols) and correlations (lines) of H₂ Henry's constant in pure water.</figcaption></figure>

<table><caption>**Table 1**

Literature experimental (Exp) and molecular simulations (MS) data for H₂ solubility in pure water.</caption><thead><tr><th>Reference</th><th>Year</th><th>Data type</th><th>Tmin/ K</th><th>Tmax/ K</th><th>Pmin/ MPa</th><th>Pmax/ MPa</th></tr></thead><tbody><tr><td>Wiebe et al. [36]</td><td>1932</td><td>Exp</td><td>298.15</td><td>298.15</td><td>2.5331</td><td>101.325</td></tr><tr><td>Wiebe & Gaddy [37]</td><td>1934</td><td>Exp</td><td>323.15</td><td>373.15</td><td>2.5331</td><td>101.325</td></tr><tr><td>Morrison & Billett [38]</td><td>1952</td><td>Exp</td><td>285.45</td><td>344.85</td><td>0.10276</td><td>0.13490</td></tr><tr><td>Crozier & Yamamoto [39]</td><td>1974</td><td>Exp</td><td>274.60</td><td>302.47</td><td>0.10130</td><td>0.10130</td></tr><tr><td>Gordon et al. [40]</td><td>1977</td><td>Exp</td><td>273.29</td><td>302.40</td><td>0.10130</td><td>0.10130</td></tr><tr><td>Gillespie & Wilson [41]</td><td>1980</td><td>Exp</td><td>310.93</td><td>588.71</td><td>0.34500</td><td>13.800</td></tr><tr><td>Alvarez et al. [42]</td><td>1988</td><td>Exp</td><td>318.9</td><td>475.6</td><td>0.43600</td><td>3.5740</td></tr><tr><td>Kling & Maurer [43]</td><td>1991</td><td>Exp</td><td>323.15</td><td>423.15</td><td>3.1800</td><td>15.370</td></tr><tr><td>Rahbari et al. [23]</td><td>2019</td><td>MS</td><td>283</td><td>423</td><td>1.0000</td><td>100.000</td></tr><tr><td>Chabab et al. [5]</td><td>2020</td><td>Exp</td><td>323.18</td><td>372.71</td><td>2.9272</td><td>12.1706</td></tr><tr><td>Torín-Ollarves et Trusler [12]</td><td>2021</td><td>Exp</td><td>323.15</td><td>323.15</td><td>15.800</td><td>32.970</td></tr><tr><td>Chabab et al. [15]</td><td>2023</td><td>Exp</td><td>298.15</td><td>373.15</td><td>10.000</td><td>20.000</td></tr><tr><td>This work</td><td>2023</td><td>MS</td><td>298.15</td><td>453.15</td><td>0.1000</td><td>100.000</td></tr></tbody></table>

thermodynamic properties of aqueous nonelectrolytes at infinite dilution (Henry constant, partial molar volume, etc) and proposed parameters for H₂/H₂O system; 4) the polynomial correlations of Li et al. (2018) [34] and Torín-Ollarves and Trusler (2021) [12] which were

obtained by adjusting simultaneously Henry constant and partial molar volume on high pressure H₂ solubility data in water based on Krichevsky and Kasarnovsky approach [29]. Interestingly, as shown in Fig. 2, there is a significant discrepancy between these correlations providing H₂ Henry constant in pure water and the available experimental data at high temperatures (up to 5 % at 100 °C and up to 20 % at 200 °C).

Due to its development from solid conventional thermodynamic relations, and its ability to predict all thermodynamic properties of solute at infinite dilution with property-independent parameters, the model of Akinfiev and Diamond 2002 [33,35] has been selected as a reference in this work. Moreover, it is probably the best model to represent the data at both low (where most data are available, see Table 1) and high temperatures, see Fig. 2.

Actually, in aquifers, brines are a mixture of salts dissolved in water, however, in general, the main species considered are sodium and chloride, and other impurities having an insignificant influence on solubility. This is why, in most cases, and this only concerns thermodynamic properties, underground brines are assimilated to NaCl brine (neglecting the other species or converting them to NaCl equivalents [44]). However, in terms of chemical reactivity, real existing species must be considered and not assimilated to NaCl. Since our focus in this work is the H₂ solubility in brine, therefore, the system H₂O + NaCl is considered as a reservoir brine model.

In NaCl brine, the review of available experimental and molecular simulation data on high-pressure hydrogen solubility (see Table 2) revealed very limited data, with only three experimental studies [5,12,15] and two recent molecular simulation studies [13,14] that exhibit remarkable discrepancies as illustrated in Fig. 3 in terms of salt free molar fraction (x'H₂) at a molality of 1 m (m: mol of NaCl/kg of water), a temperature of 323 K and a pressure up to 20 MPa (since all literature experimental data in NaCl brine are in terms of molar fraction).

All of this reinforces the need for further studies to better constrain the values of H₂ solubility in brine under underground conditions.
<figure><img src="image_5.png" /><figcaption>Fig. 2. Experimental data (symbols) and correlations (lines) of H₂ Henry's constant in pure water.</figcaption></figure>

<table><caption>**Table 1**

Literature experimental (Exp) and molecular simulations (MS) data for H₂ solubility in pure water.</caption><thead><tr><th>Reference</th><th>Year</th><th>Data type</th><th>Tmin/ K</th><th>Tmax/ K</th><th>Pmin/ MPa</th><th>Pmax/ MPa</th></tr></thead><tbody><tr><td>Wiebe et al. [36]</td><td>1932</td><td>Exp</td><td>298.15</td><td>298.15</td><td>2.5331</td><td>101.325</td></tr><tr><td>Wiebe & Gaddy [37]</td><td>1934</td><td>Exp</td><td>323.15</td><td>373.15</td><td>2.5331</td><td>101.325</td></tr><tr><td>Morrison & Billett [38]</td><td>1952</td><td>Exp</td><td>285.45</td><td>344.85</td><td>0.10276</td><td>0.13490</td></tr><tr><td>Crozier & Yamamoto [39]</td><td>1974</td><td>Exp</td><td>274.60</td><td>302.47</td><td>0.10130</td><td>0.10130</td></tr><tr><td>Gordon et al. [40]</td><td>1977</td><td>Exp</td><td>273.29</td><td>302.40</td><td>0.10130</td><td>0.10130</td></tr><tr><td>Gillespie & Wilson [41]</td><td>1980</td><td>Exp</td><td>310.93</td><td>588.71</td><td>0.34500</td><td>13.800</td></tr><tr><td>Alvarez et al. [42]</td><td>1988</td><td>Exp</td><td>318.9</td><td>475.6</td><td>0.43600</td><td>3.5740</td></tr><tr><td>Kling & Maurer [43]</td><td>1991</td><td>Exp</td><td>323.15</td><td>423.15</td><td>3.1800</td><td>15.370</td></tr><tr><td>Rahbari et al. [23]</td><td>2019</td><td>MS</td><td>283</td><td>423</td><td>1.0000</td><td>100.000</td></tr><tr><td>Chabab et al. [5]</td><td>2020</td><td>Exp</td><td>323.18</td><td>372.71</td><td>2.9272</td><td>12.1706</td></tr><tr><td>Torín-Ollarves et Trusler [12]</td><td>2021</td><td>Exp</td><td>323.15</td><td>323.15</td><td>15.800</td><td>32.970</td></tr><tr><td>Chabab et al. [15]</td><td>2023</td><td>Exp</td><td>298.15</td><td>373.15</td><td>10.000</td><td>20.000</td></tr><tr><td>This work</td><td>2023</td><td>MS</td><td>298.15</td><td>453.15</td><td>0.1000</td><td>100.000</td></tr></tbody></table>

thermodynamic properties of aqueous nonelectrolytes at infinite dilution (Henry constant, partial molar volume, etc) and proposed parameters for H₂/H₂O system; 4) the polynomial correlations of Li et al. (2018) [34] and Torín-Ollarves and Trusler (2021) [12] which were

obtained by adjusting simultaneously Henry constant and partial molar volume on high pressure H₂ solubility data in water based on Krichevsky and Kasarnovsky approach [29]. Interestingly, as shown in Fig. 2, there is a significant discrepancy between these correlations providing H₂ Henry constant in pure water and the available experimental data at high temperatures (up to 5 % at 100 °C and up to 20 % at 200 °C).

Due to its development from solid conventional thermodynamic relations, and its ability to predict all thermodynamic properties of solute at infinite dilution with property-independent parameters, the model of Akinfiev and Diamond 2002 [33,35] has been selected as a reference in this work. Moreover, it is probably the best model to represent the data at both low (where most data are available, see Table 1) and high temperatures, see Fig. 2.

Actually, in aquifers, brines are a mixture of salts dissolved in water, however, in general, the main species considered are sodium and chloride, and other impurities having an insignificant influence on solubility. This is why, in most cases, and this only concerns thermodynamic properties, underground brines are assimilated to NaCl brine (neglecting the other species or converting them to NaCl equivalents [44]). However, in terms of chemical reactivity, real existing species must be considered and not assimilated to NaCl. Since our focus in this work is the H₂ solubility in brine, therefore, the system H₂O + NaCl is considered as a reservoir brine model.

In NaCl brine, the review of available experimental and molecular simulation data on high-pressure hydrogen solubility (see Table 2) revealed very limited data, with only three experimental studies [5,12,15] and two recent molecular simulation studies [13,14] that exhibit remarkable discrepancies as illustrated in Fig. 3 in terms of salt free molar fraction (x'H₂) at a molality of 1 m (m: mol of NaCl/kg of water), a temperature of 323 K and a pressure up to 20 MPa (since all literature experimental data in NaCl brine are in terms of molar fraction).

All of this reinforces the need for further studies to better constrain the values of H₂ solubility in brine under underground conditions.
<table><caption>**Table 2**

Literature experimental (Exp) and molecular simulations (MS) data for H₂ solubility in NaCl brine.</caption><thead><tr><th>Reference</th><th>Year</th><th>Data type</th><th>Max Molality/ mol/kg<sub>w</sub></th><th>Tmin/ K</th><th>Tmax/ K</th><th>Pmin/ MPa</th><th>Pmax/ MPa</th></tr></thead><tbody><tr><td>Steiner [45]</td><td>1894</td><td>Exp</td><td>5.3</td><td>286.32</td><td>286.95</td><td>0.1028</td><td>0.10290</td></tr><tr><td>Braun [46]</td><td>1900</td><td>Exp</td><td>1.1</td><td>278.15</td><td>298.15</td><td>0.1022</td><td>0.10450</td></tr><tr><td>Gerecke & Bittrich [47]</td><td>1971</td><td>Exp</td><td>4.3</td><td>288.15</td><td>298.15</td><td>0.1030</td><td>0.10450</td></tr><tr><td>Crozier & Yamamoto [39]</td><td>1974</td><td>Exp</td><td>0.5</td><td>274.03</td><td>301.51</td><td>0.1013</td><td>0.10130</td></tr><tr><td>Lopez-Lazaro et al. [13]</td><td>2019</td><td>MS</td><td>2.0</td><td>270.00</td><td>570.00</td><td>Psat</td><td>Psat</td></tr><tr><td>Chabab et al. [5]</td><td>2020</td><td>Exp</td><td>5.0</td><td>323.00</td><td>372.00</td><td>2.000</td><td>20.000</td></tr><tr><td>Torín-Ollarves et Trusler [12]</td><td>2021</td><td>Exp</td><td>2.5</td><td>323.00</td><td>323.00</td><td>12.000</td><td>39.400</td></tr><tr><td>van Rooijen et al. [14]</td><td>2023</td><td>MS</td><td>5.93</td><td>298.00</td><td>363.00</td><td>0.1000</td><td>100.00</td></tr><tr><td>Chabab et al. [15]</td><td>2023</td><td>Exp</td><td>4</td><td>298.15</td><td>373.15</td><td>10.000</td><td>20.000</td></tr><tr><td>This work</td><td>—</td><td>MS</td><td>6</td><td>298.15</td><td>453.15</td><td>0.1000</td><td>100.00</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** High-pressure hydrogen solubility in brine at 1 m and at 323 K: experimental and molecular simulation data.</figcaption></figure>

# 3. Molecular simulation

## 3.1. Simulation approaches

### 3.1.1. Monte Carlo

#### (A) Gibbs ensemble method

One of the best-known methods to achieve the liquid-vapor equilibria is to carry out Monte Carlo simulations in the Gibbs ensemble introduced by Panagiotopoulos [48]. In this method, the simulation systems contain two separated boxes to model vapor and liquid phases in equilibrium. Among the advantages of the Gibbs ensemble in the study of phase equilibria is the ability to exchange molecules between the two phases to reach equilibrium without including an interface. When two phases are separated by an interface, a large number of molecules is required to reduce the influence of the interface on the properties of the system and thus a large CPU time [49]. The MC simulations in Gibbs ensemble consist of MC moves in each box to reach thermal equilibrium, volume change moves to achieve pressure equilibrium and molecular exchange moves to attain chemical equilibrium.

This method provides accurate results (density, composition, ...) for simple systems (small molecules) far from the critical point [49,50]. Furthermore, it is important to emphasize that this method is difficult to apply directly for very dense phases.

#### (B) Henry's law method

The solubility of a gas in a solvent, at low pressure and low dilution of the solute, is well described by the Henry's law:

$$x_i = \frac{P_i}{H_i} = \frac{y_i P}{H_i} \qquad (1)$$

where $x_i$ is the mole fraction of the solute i in the solvent, $P_i$ is the partial pressure of the solute in the gas phase, which is equal to its vapor-phase mole fraction $y_i$ multiplied by the total pressure $P$, and $H_i$ is the Henry's law constant.

The Henry's law constant can be calculated from MC simulations using the following relation [51,52]:

$$H_i = \rho k_B T \exp(\mu_i^{\text{excess}} / k_B T) \qquad (2)$$

where $\rho$ is the density of the saturated liquid phase, $k_B$ is the Boltzmann constant, $T$ is the temperature, and $\mu_i^{\text{excess}}$ is the excess chemical potential of the solute i at infinite dilution. To perform MC simulation of the saturated liquid phase, there are two methods: the first relies on the Gibbs ensemble, while the other one, used in this work, relies on the isobaric-isothermal (NPT) ensemble for which the saturation pressure is given a priori.

For high pressure applications, as in underground conditions, Henry's law should be modified to introduce (1) the fugacity $\varphi_i$ and the Poynting factor $\Pi$ which are used to correct the effect of pressure in the gas and liquid phases respectively, and (2) the activity coefficient $\gamma_i$ for correcting the non-ideality of the liquid phase [53]. Using Krichevsky Kasarnovsky like method, and assuming that the solute i is not compressible in the liquid phase at infinite dilution, $x_i$ can be calculated from [53]:
<table><caption>**Table 2**

Literature experimental (Exp) and molecular simulations (MS) data for H₂ solubility in NaCl brine.</caption><thead><tr><th>Reference</th><th>Year</th><th>Data type</th><th>Max Molality/ mol/kg<sub>w</sub></th><th>Tmin/ K</th><th>Tmax/ K</th><th>Pmin/ MPa</th><th>Pmax/ MPa</th></tr></thead><tbody><tr><td>Steiner [45]</td><td>1894</td><td>Exp</td><td>5.3</td><td>286.32</td><td>286.95</td><td>0.1028</td><td>0.10290</td></tr><tr><td>Braun [46]</td><td>1900</td><td>Exp</td><td>1.1</td><td>278.15</td><td>298.15</td><td>0.1022</td><td>0.10450</td></tr><tr><td>Gerecke & Bittrich [47]</td><td>1971</td><td>Exp</td><td>4.3</td><td>288.15</td><td>298.15</td><td>0.1030</td><td>0.10450</td></tr><tr><td>Crozier & Yamamoto [39]</td><td>1974</td><td>Exp</td><td>0.5</td><td>274.03</td><td>301.51</td><td>0.1013</td><td>0.10130</td></tr><tr><td>Lopez-Lazaro et al. [13]</td><td>2019</td><td>MS</td><td>2.0</td><td>270.00</td><td>570.00</td><td>Psat</td><td>Psat</td></tr><tr><td>Chabab et al. [5]</td><td>2020</td><td>Exp</td><td>5.0</td><td>323.00</td><td>372.00</td><td>2.000</td><td>20.000</td></tr><tr><td>Torín-Ollarves et Trusler [12]</td><td>2021</td><td>Exp</td><td>2.5</td><td>323.00</td><td>323.00</td><td>12.000</td><td>39.400</td></tr><tr><td>van Rooijen et al. [14]</td><td>2023</td><td>MS</td><td>5.93</td><td>298.00</td><td>363.00</td><td>0.1000</td><td>100.00</td></tr><tr><td>Chabab et al. [15]</td><td>2023</td><td>Exp</td><td>4</td><td>298.15</td><td>373.15</td><td>10.000</td><td>20.000</td></tr><tr><td>This work</td><td>-</td><td>MS</td><td>6</td><td>298.15</td><td>453.15</td><td>0.1000</td><td>100.00</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** High-pressure hydrogen solubility in brine at 1 m and at 323 K: experimental and molecular simulation data.</figcaption></figure>

# 3. Molecular simulation

## 3.1. Simulation approaches

### 3.1.1. Monte Carlo

#### (A) Gibbs ensemble method

One of the best-known methods to achieve the liquid-vapor equilibria is to carry out Monte Carlo simulations in the Gibbs ensemble introduced by Panagiotopoulos [48]. In this method, the simulation systems contain two separated boxes to model vapor and liquid phases in equilibrium. Among the advantages of the Gibbs ensemble in the study of phase equilibria is the ability to exchange molecules between the two phases to reach equilibrium without including an interface. When two phases are separated by an interface, a large number of molecules is required to reduce the influence of the interface on the properties of the system and thus a large CPU time [49]. The MC simulations in Gibbs ensemble consist of MC moves in each box to reach thermal equilibrium, volume change moves to achieve pressure equilibrium and molecular exchange moves to attain chemical equilibrium.

This method provides accurate results (density, composition, ...) for simple systems (small molecules) far from the critical point [49,50]. Furthermore, it is important to emphasize that this method is difficult to apply directly for very dense phases.

#### (B) Henry's law method

The solubility of a gas in a solvent, at low pressure and low dilution of the solute, is well described by the Henry's law:

$$x_i = \frac{P_i}{H_i} = \frac{y_i P}{H_i} \qquad (1)$$

where $x_i$ is the mole fraction of the solute i in the solvent, $P_i$ is the partial pressure of the solute in the gas phase, which is equal to its vapor-phase mole fraction $y_i$ multiplied by the total pressure $P$, and $H_i$ is the Henry's law constant.

The Henry's law constant can be calculated from MC simulations using the following relation [51,52]:

$$H_i = \rho k_B T \exp(\mu_i^{\text{excess}} / k_B T) \qquad (2)$$

where $\rho$ is the density of the saturated liquid phase, $k_B$ is the Boltzmann constant, $T$ is the temperature, and $\mu_i^{\text{excess}}$ is the excess chemical potential of the solute i at infinite dilution. To perform MC simulation of the saturated liquid phase, there are two methods: the first relies on the Gibbs ensemble, while the other one, used in this work, relies on the isobaric-isothermal (NPT) ensemble for which the saturation pressure is given a priori.

For high pressure applications, as in underground conditions, Henry's law should be modified to introduce (1) the fugacity $\varphi_i$ and the Poynting factor $\Pi$ which are used to correct the effect of pressure in the gas and liquid phases respectively, and (2) the activity coefficient $\gamma_i$ for correcting the non-ideality of the liquid phase [53]. Using Krichevsky Kasarnovsky like method, and assuming that the solute i is not compressible in the liquid phase at infinite dilution, $x_i$ can be calculated from [53]:
$$x_i = \frac{y_i P \rho_i}{H_i \Pi} \qquad (3)$$

with

$$\Pi = \exp\left(\frac{\theta_i^\infty P - P_{\text{solvent}}^{\text{sat}}}{RT}\right) \qquad (4)$$

where $P_{\text{solvent}}^{\text{sat}}$ is the saturation pressure of the solvent and $\theta_i^\infty$ is the partial molar volume of the solute i at infinite dilution.

### (C) Sampling methods

In the case of dense solvents, as in this work, the choice of an adequate sampling method in the Monte Carlo simulation is an important factor to improve the sampling efficiency. Thus, in this work, the efficient Continuous Fractional Component Monte Carlo (CFCMC) sampling method introduced by Shi and Maginn [24,25] has been used.

In the CFCMC method, molecules are gradually inserted/deleted into the system. During the insertion/deletion process, the interactions of these molecules with their surroundings are modulated. This modulation is controlled by a lambda ($\lambda$) coupling parameter, with $\lambda \in [0, 1]$. A $\lambda_i$ value of 0 indicates that the fractional molecule of type i behaves as an ideal gas molecule, exhibiting no interaction with surrounding molecules. Conversely, a $\lambda_i$ value of 1 signifies that the fractional molecule of type i interacts fully, akin to a whole molecule [54]. In the CFCMC Gibbs ensemble, while fractional molecule may reside in one of the phases, its position and simulation box remain constant during the scaling of the $\lambda$ parameter. It is important to note that the extent of $\lambda$ variation might differ between phases. For example, the maximum change in $\lambda$ could be larger in a gas phase than that in a liquid phase [54].

Thus, in addition to the classical MC trial moves such as translation, rotation, volume change, the CFCMC simulations include moves associated with the fractional parameter $\lambda$ as follows: (1) lambda trial move, where the value of the fractional parameter $\lambda$ is randomly changed; (2) swap trial moves, where the fractional molecule is removed from one simulation box and reinserted at a random position with a random orientation in the same (NVT/NPT ensemble) or in another (Gibbs ensemble) simulation box leaving $\lambda$ unchanged and (3) identity change trial moves, where the fractional molecule is transformed into an entire molecule while at the same time a randomly chosen entire molecule of the same component is transformed into a fractional molecule with the same fractional parameter $\lambda$. The combination of the swap and identity change trial moves (named hybrid CFC moves) allows for an improved sampling.

An advantage of the CFCMC approach is its efficiency to provide rather accurate partial molar properties. In particular, the excess chemical potential that is used to compute the Henry's law constant, equation (5), can be determined by:

$$\mu_i^{ex} = -k_B T \ln \left( \frac{p(\lambda_{ij} = 1)}{p(\lambda_{ij} = 0)} \right) \qquad (5)$$

where $\lambda_{ij}$ is the fractional parameter of the fractional molecule i in the phase j and $p(\lambda_{ij})$ is the Boltzmann probability distribution of $\lambda_{ij}$. The partial molar volume used to calculate the Poynting factor, equation (6), is given by:

$$v_i^\infty = \langle V(\lambda_{ij} = 1) \rangle - \langle \frac{1}{V} (\lambda_{ij} = 0) \rangle^{-1} \qquad (6)$$

where $V(\lambda_{ij})$ is the volume of solvent corresponding to the fractional parameter of $\lambda_{ij}$.

### 3.1.2. Radial distribution function

In this work, to provide insight into the microscopic structure of the H₂-brine system at the molecular level, radial distribution functions (RDF) were calculated using the following relation [55]:

$$g_{AB}(r) = \frac{\rho_{AB}(r)}{\rho_{0B}} = \left\langle \frac{V}{N_A} \times \frac{1}{N_B} \sum_{i \in A}^{N_A} \sum_{j \in B}^{N_B} \frac{\delta(r_{ij} - r)}{4\pi r^2} \right\rangle \qquad (7)$$

where, $g_{AB}(r)$ is the radial distribution function of component B around component A at a distance $r$, $\rho_{AB}(r)$ is the number density of component B at a distance $r$ from a reference particle of component A, $\rho_{0B}$ is the overall number density of component B in the system, $\delta$ is the Dirac function, $r_{ij}$ is the distance between molecules i of component A and particle j of component B, $N_A$ is the number of molecules of component A.

It should be noticed, that due to the low solubility of H₂ in brine, obtaining RDF involving H₂ with good statistical uncertainties using Monte Carlo simulations would have required impractically long simulations. Therefore, to improve the accuracy on the computed RDF, molecular dynamics (MD) simulations, which are efficiently parallelized, were performed to reduce computational expense.

## 3.2. Force fields

### Hydrogen

Partially the description of hydrogen molecules using classical semi-empirical force fields may be challenging because of its quantum effects [56], explaining why several force fields for the hydrogen molecules exist in the literature. These force fields include different representations: single site [57–60], two sites [61,62], and multiple sites with a quadrupole moment [63–65]. According to recent studies carried out in the context of UHS [13,23], the multiple-site models with quadrupole moment are probably the best options to accurately predict the bulk thermodynamic properties of hydrogen as well as the Vapor-Liquid Equilibria (VLE) of H₂-H₂O and H₂-H₂O-NaCl systems.

In this work, three models have been selected, one two sites model and two multiple-sites models. The two sites model is that of Yang & Zhong [62] which describe H₂ using two Lennard-Jones spheres without charge. The multiple sites models are the ones of Marx [64] and of Alavi/2005 [63], which represent the hydrogen molecule by one Lennard-Jones sphere and three fixed charges. The corresponding force fields parameters of the three H₂ models are provided in the supplementary materials.

### Water

Although real water molecules are flexible and polarizable, they have been widely modeled as rigid non-polarizable molecule because of the simplicity and computational efficiency of such models. According to Vega et al (2011) [66], among all non-polarizable rigid force fields available at the time of this study, the TIP4P/2005 model was the best model to consistently predicts various thermophysical properties (density, surface tension, self-diffusion coefficient, etc.) of pure water, followed by the SPC/E model [67] which is widely used in the geoscience community because of its simplicity. In electrolyte solutions containing gas, as studied in this work, it is a priori crucial to describe well the density and dielectric constant of the solvent. Since the SPC/E model represents more accurately the dielectric constant of water than the TIP4P/2005 model [68], both models were selected for this study. In addition, the TIP4P/EP model developed by Fuentes-Azcatl and Alejandre (2014) [69] was also selected for its good representation of both density and dielectric constant. The corresponding force fields parameters of the three selected water models are provided in the supplementary materials.

### Salt

For aqueous NaCl solution, Creton et al. [17] have found that the combination of the TIP4P/2005 model for pure water and the OPLS model [70] for Na⁺ and Cl⁻ allows to well predict densities and osmotic pressures over a wide range of temperature and salinity. This set of models has been subsequently used for the calculation of the solubility of H₂ in NaCl brine at different molalities [13]. However, it was found that the Henry's constants of H₂ in NaCl brine were systematically
$$x_i = \frac{y_i P \rho_i}{H_i \Pi} \qquad (3)$$

with

$$\Pi = \exp\left(\frac{\theta_i^\infty P - P_{\text{solvent}}^{\text{sat}}}{RT}\right) \qquad (4)$$

where $P_{\text{solvent}}^{\text{sat}}$ is the saturation pressure of the solvent and $\theta_i^\infty$ is the partial molar volume of the solute $i$ at infinite dilution.

### (C) Sampling methods

In the case of dense solvents, as in this work, the choice of an adequate sampling method in the Monte Carlo simulation is an important factor to improve the sampling efficiency. Thus, in this work, the efficient Continuous Fractional Component Monte Carlo (CFCMC) sampling method introduced by Shi and Maginn [24,25] has been used.

In the CFCMC method, molecules are gradually inserted/deleted into the system. During the insertion/deletion process, the interactions of these molecules with their surroundings are modulated. This modulation is controlled by a lambda ($\lambda$) coupling parameter, with $\lambda \in [0, 1]$. A $\lambda_i$ value of 0 indicates that the fractional molecule of type $i$ behaves as an ideal gas molecule, exhibiting no interaction with surrounding molecules. Conversely, a $\lambda_i$ value of 1 signifies that the fractional molecule of type $i$ interacts fully, akin to a whole molecule [54]. In the CFCMC Gibbs ensemble, while fractional molecule may reside in one of the phases, its position and simulation box remain constant during the scaling of the $\lambda$ parameter. It is important to note that the extent of $\lambda$ variation might differ between phases. For example, the maximum change in $\lambda$ could be larger in a gas phase than that in a liquid phase [54].

Thus, in addition to the classical MC trial moves such as translation, rotation, volume change, the CFCMC simulations include moves associated with the fractional parameter $\lambda$ as follows: (1) lambda trial move, where the value of the fractional parameter $\lambda$ is randomly changed; (2) swap trial moves, where the fractional molecule is removed from one simulation box and reinserted at a random position with a random orientation in the same (NVT/NPT ensemble) or in another (Gibbs ensemble) simulation box leaving $\lambda$ unchanged and (3) identity change trial moves, where the fractional molecule is transformed into an entire molecule while at the same time a randomly chosen entire molecule of the same component is transformed into a fractional molecule with the same fractional parameter $\lambda$. The combination of the swap and identity change trial moves (named hybrid CFC moves) allows for an improved sampling.

An advantage of the CFCMC approach is its efficiency to provide rather accurate partial molar properties. In particular, the excess chemical potential that is used to compute the Henry's law constant, equation (5), can be determined by:

$$\mu_i^{\text{ex}} = -k_B T \ln \left( \frac{p(\lambda_{ij} = 1)}{p(\lambda_{ij} = 0)} \right) \qquad (5)$$

where $\lambda_{ij}$ is the fractional parameter of the fractional molecule $i$ in the phase $j$ and $p(\lambda_{ij})$ is the Boltzmann probability distribution of $\lambda_{ij}$. The partial molar volume used to calculate the Poynting factor, equation (6), is given by:

$$v_i^\infty = \langle V(\lambda_{ij} = 1) \rangle - \langle \frac{1}{V} (\lambda_{ij} = 0) \rangle^{-1} \qquad (6)$$

where $V(\lambda_{ij})$ is the volume of solvent corresponding to the fractional parameter of $\lambda_{ij}$.

### 3.1.2. Radial distribution function

In this work, to provide insight into the microscopic structure of the H₂-brine system at the molecular level, radial distribution functions (RDF) were calculated using the following relation [55]:

$$g_{AB}(\mathbf{r}) = \frac{\rho_{AB}(\mathbf{r})}{\rho_{0B}} = \left\langle \frac{V}{N_A} \times \frac{1}{N_B} \sum_{i \in A}^{N_A} \sum_{j \in B}^{N_B} \frac{\delta(r_{ij} - \mathbf{r})}{4\pi r^2} \right\rangle \qquad (7)$$

where, $g_{AB}(\mathbf{r})$ is the radial distribution function of component B around component A at a distance $r$, $\rho_{AB}(\mathbf{r})$ is the number density of component B at a distance $r$ from a reference particle of component A, $\rho_{0B}$ is the overall number density of component B in the system, $\delta$ is the Dirac function, $r_{ij}$ is the distance between molecules $i$ of component A and particle $j$ of component B, $N_A$ is the number of molecules of component A.

It should be noticed, that due to the low solubility of H₂ in brine, obtaining RDF involving H₂ with good statistical uncertainties using Monte Carlo simulations would have required impractically long simulations. Therefore, to improve the accuracy on the computed RDF, molecular dynamics (MD) simulations, which are efficiently parallelized, were performed to reduce computational expense.

## 3.2. Force fields

### Hydrogen

Partially the description of hydrogen molecules using classical semi-empirical force fields may be challenging because of its quantum effects [56], explaining why several force fields for the hydrogen molecules exist in the literature. These force fields include different representations: single site [57–60], two sites [61,62], and multiple sites with a quadrupole moment [63–65]. According to recent studies carried out in the context of UHS [13,23], the multiple-site models with quadrupole moment are probably the best options to accurately predict the bulk thermodynamic properties of hydrogen as well as the Vapor-Liquid Equilibria (VLE) of H₂-H₂O and H₂-H₂O-NaCl systems.

In this work, three models have been selected, one two sites model and two multiple-sites models. The two sites model is that of Yang & Zhong [62] which describe H₂ using two Lennard-Jones spheres without charge. The multiple sites models are the ones of Marx [64] and of Alavi/2005 [63], which represent the hydrogen molecule by one Lennard-Jones sphere and three fixed charges. The corresponding force fields parameters of the three H₂ models are provided in the supplementary materials.

### Water

Although real water molecules are flexible and polarizable, they have been widely modeled as rigid non-polarizable molecule because of the simplicity and computational efficiency of such models. According to Vega et al (2011) [66], among all non-polarizable rigid force fields available at the time of this study, the TIP4P/2005 model was the best model to consistently predicts various thermophysical properties (density, surface tension, self-diffusion coefficient, etc.) of pure water, followed by the SPC/E model [67] which is widely used in the geoscience community because of its simplicity. In electrolyte solutions containing gas, as studied in this work, it is a priori crucial to describe well the density and dielectric constant of the solvent. Since the SPC/E model represents more accurately the dielectric constant of water than the TIP4P/2005 model [68], both models were selected for this study. In addition, the TIP4P/EP model developed by Fuentes-Azcatl and Alejandre (2014) [69] was also selected for its good representation of both density and dielectric constant. The corresponding force fields parameters of the three selected water models are provided in the supplementary materials.

### Salt

For aqueous NaCl solution, Creton et al. [17] have found that the combination of the TIP4P/2005 model for pure water and the OPLS model [70] for Na⁺ and Cl⁻ allows to well predict densities and osmotic pressures over a wide range of temperature and salinity. This set of models has been subsequently used for the calculation of the solubility of H₂ in NaCl brine at different molalities [13]. However, it was found that the Henry's constants of H₂ in NaCl brine were systematically
overestimated (the higher the salinity, the larger the overestimation). The authors suggested that this was due to the bad modeling of the cross interactions between hydrogen and ionic species using the Lorentz-Berthelot (LB) combining rules [13], adjusting them to improve the computations. More recently, van Rooijen et al. [14] proposed another option yielding good result on H₂ solubility in brine based on the TIP4P/2005 model for water combined with the Madrid force field for Na⁺ and Cl⁻ ions salts [71,72]. A last option, studied in this work, is the one proposed by Loche et al [73]. They consider that, before including polarization (e.g., as indirectly considered by the two previous approaches), the Lennard-Jones parameters of the existing non-polarizable ionic force fields should be optimized first and provided a new set of parameters adjusted on the experimental solvation free energy, activity coefficient at finite salt concentrations and the Kirkwood-Buff integrals. This force field is called KBF in the following. Thus, in this work, OPLS, Madrid and KBF force fields for for Na⁺ and Cl⁻ have been evaluated. Their molecular parameters are given in the supplementary materials.

### 3.3. Simulation details

#### 3.3.1. Monte Carlo simulations

Monte Carlo simulations were performed using the open-source Brick-CFCMC package. The Lennard Jones interaction potentials were computed by truncating at a cutoff radius of 12 Å, and long rang corrections were included. The Coulomb interaction potentials were computed using Ewald summation. To determine the Lennard-Jones interaction parameters of unlike atoms, the classical Lorentz-Berthelot combining rules have been used [74]. Periodic boundary conditions were imposed in all directions. To estimate the error bars, six independent MC simulations were carried out for a given thermodynamic condition. In figures where there is no error bar, note that the statistical uncertainties are comparable to or smaller than the symbols. Therefore, the error bars have not been displayed for clarity.

#### Gibbs ensemble Monte Carlo simulations

The simulation systems for the CFC-GE method, used to compute H₂ solubility in pure water, were initially set up so that the liquid phase consists of 830 water molecules, and the gas phase contains of 600 hydrogen molecules. For H₂ solubility in brine, the liquid phase contains 839 water molecules and a varying number of Na⁺ and Cl⁻ ions, ranging from 0 to 90 depending on the molality. The gas phase contains 600 hydrogen molecules. For all simulations, 5 × 10⁵ equilibration cycles were carried out, followed by 2 × 10⁶ production cycles. The number of trial moves per cycle was equal to the number of molecules. In the CFC-GE simulations, the trial moves were selected with the following probabilities: 1 % for volume changes, 33 % for translations, 33 % for rotations, 17 % for λ changes, 8.5 % for reinsertions of fractional molecules at randomly selected positions, and 8.5 % for identity changes of fractional molecules. For the hybrid CFC moves (reinsertions and identity changes), with λ < λ<sub>swap</sub> (fixed at 0.3) for implementing a swap trial move, λ > λ<sub>change</sub> (fixed at 0.7) for implementing an identity change trial move. Otherwise, no trial move is performed.

#### Henry's law approach

The Henry's coefficients for H₂ in water and NaCl brine solutions were computed at temperatures ranging from 298.15 K to 453.15 K. These calculations were performed using MC simulations in the CFC-NPT ensemble, where the vapor pressure of pure water was obtained by using the same force field. The simulation systems contain 500 water molecules and a varying number of Na⁺ and Cl⁻ ion pairs depending on the molality (with 9 particles of Na⁺ and Cl⁻ corresponding to a molality of 1 mol/kgH₂O). For all simulations, the systems were equilibrated during 15 × 10⁴ MC cycles, followed by 5 × 10⁵ production cycles. The number of trial moves per cycle is equal to the number of molecules. The trial moves in the CFC-NPT simulations were selected with the following probabilities: 1 % for volume changes, 35 % for translations, 30 % for rotations, 17 % for λ changes and 17 % for reinsertions of fractional

molecules at randomly selected positions. Throughout all simulations, a 50 % acceptance rate was achieved by adjusting the magnitude of translations, rigid rotations, volume changes and changes in λ.

#### 3.3.2. Molecular Dynamics simulations

The Molecular Dynamics simulations in NPT ensemble were performed using LAMMPS package. All simulations consist of two steps. First, the systems were equilibrated during 5 × 10⁶ time steps. Then, a run of 3 × 10⁶ was carried out for the samplings. The simulation systems consist of 1000 water molecules, 1 hydrogen molecule, and a from 0 to 108Na⁺ 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, while the bond length and the bending angle of H₂O were kept constant using the RATTLE algorithm [75]. The bond length of H₂ is fixed using the rigid body function. The Lennard-Jones interactions were computed by truncating at a cutoff of 12 Å, and long-range corrections were included. The Coulomb interactions were computed by using Particle-Particle-Particle-Mesh (PPPM) method [76].

Uncertainties of ensemble averages were computed by performing 20 independent simulations and recording standard deviations.

### 3.4. Simulation results

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

This section initially explores the capability of the different force fields to predict the thermodynamic properties of pure hydrogen and water, followed by an examination of H₂ solubility in pure water. Given that hydrogen exists in a supercritical state under underground conditions, the primary thermodynamic quantity of interest is its compressibility factor. In contrast, for water, which exists in a subcritical state in not too deep underground conditions, the focus is on phase equilibrium calculations.

##### (a) Pure hydrogen gas

To quantify the deviations from the ideal gas law due to underground conditions, the compressibility factor is an interesting quantity. It is defined as:

$$Z = \frac{P\theta}{RT} \qquad (8)$$

where θ is the molar volume of the studied gas and R is the ideal gas constant.

In Fig. 4, the outcomes of the three studied force fields for the compressibility factor of hydrogen at different temperatures and pressures are displayed and compared with the results obtained from the reference equation of state for normal hydrogen, developed by Leachman et al. [77], as implemented in the NIST REFPROP software [78].

At 323 K and 453 K, which corresponds to typical conditions of hydrogen storage or natural hydrogen accumulations, the Marx and Yang & Zhong models effectively provide the compressibility factor for pressure up to 100 MPa (with average absolute deviations (AAD) below 0.5 %, relatively to the reference equation of state). In contrast, the Alavi/2005 model provides a less accurate prediction of the compressibility factor (with an AAD > 1 %) at these temperatures. This suggest that including a quadrupole moment in hydrogen does not always significantly improve the prediction of the bulk properties of hydrogen in the gas phase but can improve the prediction of phase coexistence [23,79].

##### (b) Pure liquid water

The water saturation liquid densities using three rigid models SPC/E, TIP4P/2005, and TIP4P/EP have been computed and are shown in Fig. 5. They are compared against regressed experimental data [80], as implemented in the NIST REFPROP Software [78]. As well-known from the literature, it appears that the TIP4P/2005 and the TIP4P/EP models yield better results than the SPC/E model, with AAD respectively equal
overestimated (the higher the salinity, the larger the overestimation). The authors suggested that this was due to the bad modeling of the cross interactions between hydrogen and ionic species using the Lorentz-Berthelot (LB) combining rules [13], adjusting them to improve the computations. More recently, van Rooijen et al. [14] proposed another option yielding good result on H₂ solubility in brine based on the TIP4P/2005 model for water combined with the Madrid force field for Na⁺ and Cl⁻ ions salts [71,72]. A last option, studied in this work, is the one proposed by Loche et al [73]. They consider that, before including polarization (e.g., as indirectly considered by the two previous approaches), the Lennard-Jones parameters of the existing non-polarizable ionic force fields should be optimized first and provided a new set of parameters adjusted on the experimental solvation free energy, activity coefficient at finite salt concentrations and the Kirkwood-Buff integrals. This force field is called KBF in the following. Thus, in this work, OPLS, Madrid and KBF force fields for for Na⁺ and Cl⁻ have been evaluated. Their molecular parameters are given in the supplementary materials.

### 3.3. Simulation details

#### 3.3.1. Monte Carlo simulations

Monte Carlo simulations were performed using the open-source Brick-CFCMC package. The Lennard Jones interaction potentials were computed by truncating at a cutoff radius of 12 Å, and long rang corrections were included. The Coulomb interaction potentials were computed using Ewald summation. To determine the Lennard-Jones interaction parameters of unlike atoms, the classical Lorentz-Berthelot combining rules have been used [74]. Periodic boundary conditions were imposed in all directions. To estimate the error bars, six independent MC simulations were carried out for a given thermodynamic condition. In figures where there is no error bar, note that the statistical uncertainties are comparable to or smaller than the symbols. Therefore, the error bars have not been displayed for clarity.

#### Gibbs ensemble Monte Carlo simulations

The simulation systems for the CFC-GE method, used to compute H₂ solubility in pure water, were initially set up so that the liquid phase consists of 830 water molecules, and the gas phase contains of 600 hydrogen molecules. For H₂ solubility in brine, the liquid phase contains 839 water molecules and a varying number of Na⁺ and Cl⁻ ions, ranging from 0 to 90 depending on the molality. The gas phase contains 600 hydrogen molecules. For all simulations, 5 × 10⁵ equilibration cycles were carried out, followed by 2 × 10⁶ production cycles. The number of trial moves per cycle was equal to the number of molecules. In the CFC-GE simulations, the trial moves were selected with the following probabilities: 1 % for volume changes, 33 % for translations, 33 % for rotations, 17 % for λ changes, 8.5 % for reinsertions of fractional molecules at randomly selected positions, and 8.5 % for identity changes of fractional molecules. For the hybrid CFC moves (reinsertions and identity changes), with λ < λ<sub>swap</sub> (fixed at 0.3) for implementing a swap trial move, λ > λ<sub>change</sub> (fixed at 0.7) for implementing an identity change trial move. Otherwise, no trial move is performed.

#### Henry's law approach

The Henry's coefficients for H₂ in water and NaCl brine solutions were computed at temperatures ranging from 298.15 K to 453.15 K. These calculations were performed using MC simulations in the CFC-NPT ensemble, where the vapor pressure of pure water was obtained by using the same force field. The simulation systems contain 500 water molecules and a varying number of Na⁺ and Cl⁻ ion pairs depending on the molality (with 9 particles of Na⁺ and Cl⁻ corresponding to a molality of 1 mol/kgH₂O). For all simulations, the systems were equilibrated during 15 × 10⁴ MC cycles, followed by 5 × 10⁵ production cycles. The number of trial moves per cycle is equal to the number of molecules. The trial moves in the CFC-NPT simulations were selected with the following probabilities: 1 % for volume changes, 35 % for translations, 30 % for rotations, 17 % for λ changes and 17 % for reinsertions of fractional

molecules at randomly selected positions. Throughout all simulations, a 50 % acceptance rate was achieved by adjusting the magnitude of translations, rigid rotations, volume changes and changes in λ.

#### 3.3.2. Molecular Dynamics simulations

The Molecular Dynamics simulations in NPT ensemble were performed using LAMMPS package. All simulations consist of two steps. First, the systems were equilibrated during 5 × 10⁶ time steps. Then, a run of 3 × 10⁶ was carried out for the samplings. The simulation systems consist of 1000 water molecules, 1 hydrogen molecule, and a from 0 to 108Na⁺ 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, while the bond length and the bending angle of H₂O were kept constant using the RATTLE algorithm [75]. The bond length of H₂ is fixed using the rigid body function. The Lennard-Jones interactions were computed by truncating at a cutoff of 12 Å, and long-range corrections were included. The Coulomb interactions were computed by using Particle-Particle-Particle-Mesh (PPPM) method [76].

Uncertainties of ensemble averages were computed by performing 20 independent simulations and recording standard deviations.

### 3.4. Simulation results

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

This section initially explores the capability of the different force fields to predict the thermodynamic properties of pure hydrogen and water, followed by an examination of H₂ solubility in pure water. Given that hydrogen exists in a supercritical state under underground conditions, the primary thermodynamic quantity of interest is its compressibility factor. In contrast, for water, which exists in a subcritical state in not too deep underground conditions, the focus is on phase equilibrium calculations.

##### (a) Pure hydrogen gas

To quantify the deviations from the ideal gas law due to underground conditions, the compressibility factor is an interesting quantity. It is defined as:

$$Z = \frac{P\theta}{RT} \qquad (8)$$

where θ is the molar volume of the studied gas and R is the ideal gas constant.

In Fig. 4, the outcomes of the three studied force fields for the compressibility factor of hydrogen at different temperatures and pressures are displayed and compared with the results obtained from the reference equation of state for normal hydrogen, developed by Leachman et al. [77], as implemented in the NIST REFPROP software [78].

At 323 K and 453 K, which corresponds to typical conditions of hydrogen storage or natural hydrogen accumulations, the Marx and Yang & Zhong models effectively provide the compressibility factor for pressure up to 100 MPa (with average absolute deviations (AAD) below 0.5 %, relatively to the reference equation of state). In contrast, the Alavi/2005 model provides a less accurate prediction of the compressibility factor (with an AAD > 1 %) at these temperatures. This suggest that including a quadrupole moment in hydrogen does not always significantly improve the prediction of the bulk properties of hydrogen in the gas phase but can improve the prediction of phase coexistence [23,79].

##### (b) Pure liquid water

The water saturation liquid densities using three rigid models SPC/E, TIP4P/2005, and TIP4P/EP have been computed and are shown in Fig. 5. They are compared against regressed experimental data [80], as implemented in the NIST REFPROP Software [78]. As well-known from the literature, it appears that the TIP4P/2005 and the TIP4P/EP models yield better results than the SPC/E model, with AAD respectively equal
<figure><img src="image_7.png" /><figcaption>Fig. 4. Compressibility factor of pure hydrogen in the gas phase at different temperatures and pressures.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Water saturation liquid densities in liquid at different temperatures.</figcaption></figure>

to 0.48 %, 0.40 % and 1.40 %. However, at high temperatures, all models begin to underestimate water densities consistently with the fact that they all underestimate the critical temperature of water [69,81].

$$ (c) H_2 + H_2O $$

First, before comparing the studied force fields, the results of the Henry's law method were compared with those from the Gibbs ensemble method. Some results are presented in Fig. 6 for the Marx-TIP4P/2005 models at a temperature of 323 K and for pressure up to 20 MPa. To do so, the Krichevsky-Kasarnovsky approach is applied, see section 4.1, using a partial molar volume of hydrogen in water at infinite dilution set at 19 cm³/mol (see discussion in section 4.1. H₂ partial molar volumes) and a fugacity coefficient calculated from the reference equation of state for normal hydrogen developed by Leachman et al. [77] (using REFPROP software).

As shown in Fig. 6, similar results are obtained with both methods with an AAD less than 3 %. Hence, in the following of this section the

presented results are coming from the use of the Henry's law method which is by far more efficient than the Gibbs ensemble method when focusing on the H₂ solubility, as it is the case here.

To find the most appropriate combination of force fields for determining the Henry's law constant of H₂ in pure water, different combinations have been evaluated: Alavi/2005-SPC/E, Alavi/2005-TIP4P/2005, Marx-TIP4P/2005, Yang & Zhong-TIP4P/2005, Marx-TIP4P/EP. All results are summarized in Fig. 7 and compared with the Akinfiev and Diamond 2002 correlation [33].

From the results it appears obviously that the Alavi/2005-SPC/E combination is the least effective both qualitatively and quantitatively. Qualitatively, the minimum solubility temperature (i.e. the maximum Henry's law constant) is equal to ~ 310 K instead of ~ 330 K from the reference correlation. According to equation (2), the Henry's constant depends mainly on the density of the solvent. Thus, as SPC/E yields increasingly underestimates of water density with temperature,
<figure><img src="image_7.png" /><figcaption>Fig. 4. Compressibility factor of pure hydrogen in the gas phase at different temperatures and pressures.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Water saturation liquid densities in liquid at different temperatures.</figcaption></figure>

to 0.48 %, 0.40 % and 1.40 %. However, at high temperatures, all models begin to underestimate water densities consistently with the fact that they all underestimate the critical temperature of water [69,81].

$$ (c) H_2 + H_2O $$

First, before comparing the studied force fields, the results of the Henry's law method were compared with those from the Gibbs ensemble method. Some results are presented in Fig. 6 for the Marx-TIP4P/2005 models at a temperature of 323 K and for pressure up to 20 MPa. To do so, the Krichevsky-Kasarnovsky approach is applied, see section 4.1, using a partial molar volume of hydrogen in water at infinite dilution set at 19 cm³/mol (see discussion in section 4.1. H₂ partial molar volumes) and a fugacity coefficient calculated from the reference equation of state for normal hydrogen developed by Leachman et al. [77] (using REFPROP software).

As shown in Fig. 6, similar results are obtained with both methods with an AAD less than 3 %. Hence, in the following of this section the

presented results are coming from the use of the Henry's law method which is by far more efficient than the Gibbs ensemble method when focusing on the H₂ solubility, as it is the case here.

To find the most appropriate combination of force fields for determining the Henry's law constant of H₂ in pure water, different combinations have been evaluated: Alavi/2005-SPC/E, Alavi/2005-TIP4P/2005, Marx-TIP4P/2005, Yang & Zhong-TIP4P/2005, Marx-TIP4P/EP. All results are summarized in Fig. 7 and compared with the Akinfiev and Diamond 2002 correlation [33].

From the results it appears obviously that the Alavi/2005-SPC/E combination is the least effective both qualitatively and quantitatively. Qualitatively, the minimum solubility temperature (i.e. the maximum Henry's law constant) is equal to ~ 310 K instead of ~ 330 K from the reference correlation. According to equation (2), the Henry's constant depends mainly on the density of the solvent. Thus, as SPC/E yields increasingly underestimates of water density with temperature,
<figure><img src="image_9.png" /><figcaption>Fig. 6. Henry method vs Gibbs ensemble method: solubility data of H₂ in pure water and brine at 1 m and at 323 K and up to 30 MPa (experimental vs. molecular simulation data).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Reference (line) and molecular simulation (symbols) data of Henry constant of H₂ in pure water versus temperature.</figcaption></figure>

see Fig. 5, one could attribute the bad qualitative evaluation of Henry's law constant to this SPC/E weakness. This interpretation is confirmed by performing simulations at fixed densities equal to the experimental densities [78], as shown in Fig. 7 by hollow triangles. This point is further confirmed by the results obtained with the Alavi/2005-TIP4P/2005 combination which provides a good estimate of the minimum solubility temperature, in line with the better water density predictions obtained with the TIP4P/2005 force field.

From a quantitative point of view, therefore, the systematic overestimation of Henry's constant is strongly connected to the choice of the hydrogen force fields, as seen when comparing the results of the Alavi/2005, Marx and Yang & Zhong models combined with the TIP4P/2005

force field. It seems that choosing a hydrogen model which describe well the compressibility factor in the gas phase (Marx and Yang & Zhong models) and contains only a quadrupole moment (Marx force field), to correctly account for cross-interactions between water molecule and hydrogen and so the correct excess chemical potential [79], is, by far, the best option, see Fig. 7.

As a last test in a predictive manner, the Marx-TIP4P/EP combination was used to evaluate the influence of having a good dielectric constant of the water model. The results, shown in Fig. 7, are very similar to those obtained using the Marx-TIP4P/2005 combination, even if generally slightly worse. This indicate that, in pure water, the improvement on the dielectric constant provided by the TIP4P/EP
<figure><img src="image_9.png" /><figcaption>Fig. 6. Henry method vs Gibbs ensemble method: solubility data of H₂ in pure water and brine at 1 m and at 323 K and up to 30 MPa (experimental vs. molecular simulation data).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Reference (line) and molecular simulation (symbols) data of Henry constant of H₂ in pure water versus temperature.</figcaption></figure>

see Fig. 5, one could attribute the bad qualitative evaluation of Henry's law constant to this SPC/E weakness. This interpretation is confirmed by performing simulations at fixed densities equal to the experimental densities [78], as shown in Fig. 7 by hollow triangles. This point is further confirmed by the results obtained with the Alavi/2005-TIP4P/2005 combination which provides a good estimate of the minimum solubility temperature, in line with the better water density predictions obtained with the TIP4P/2005 force field.

From a quantitative point of view, therefore, the systematic overestimation of Henry's constant is strongly connected to the choice of the hydrogen force fields, as seen when comparing the results of the Alavi/2005, Marx and Yang & Zhong models combined with the TIP4P/2005

force field. It seems that choosing a hydrogen model which describe well the compressibility factor in the gas phase (Marx and Yang & Zhong models) and contains only a quadrupole moment (Marx force field), to correctly account for cross-interactions between water molecule and hydrogen and so the correct excess chemical potential [79], is, by far, the best option, see Fig. 7.

As a last test in a predictive manner, the Marx-TIP4P/EP combination was used to evaluate the influence of having a good dielectric constant of the water model. The results, shown in Fig. 7, are very similar to those obtained using the Marx-TIP4P/2005 combination, even if generally slightly worse. This indicate that, in pure water, the improvement on the dielectric constant provided by the TIP4P/EP
model over the TIP4P/2005 model is not sufficient on its own for predicting the H₂ solubility.

Thus, even if very satisfying, the Marx-TIP4P/2005 and Marx-TIP4P/EP combinations yield systematic overestimations of the Henry constant of H₂ in pure water, with AADs of 6 % and 9 %, respectively. This may come from the bad evaluation of the cross-interaction parameters provide by the Lorentz-Berthelot rules, as reported by several authors for similar gas/water systems [82–84]. These studies converge on the fact that LB combining-rule underestimates gas/water polarizability. Hence introducing a positive deviation to the cross-energy Lennard-Jones interactions correct indirectly the polarizability, and leads to a better prediction of the excess chemical potential of gas in water, and consequently its solubility. As a test, cross interaction energies $ε_{ij}$ between H₂ and water, were multiplied by a constant $k_{ij}$ factor that were found to be optimal when equal to 1.05 for both Marx-TIP4P/2005 and Marx-TIP4P/EP. This optimized set of parameters yields AAD below 2 % and 4 % respectively and allows to reproduce perfectly the smoothed experimental data of the Henry constant (Akinfiev and Diамond's correlation) of H₂ in pure water over a very wide range of temperature.

### 3.4.2. H₂ + H₂O + NaCl system

(a) H₂O + NaCl system

As discussed in the previous section on the pure water, the solvent density has a significant impact on the hydrogen solubility. Therefore, to assess the Na⁺ and Cl⁻ ion pairs force fields, their impact on solvent density is crucial. For this purpose, the densities at T = 298 K and P = 0.1 MPa across various NaCl molalities (ranging from 1 to 5 mol/kg<sub>w</sub>) were computed using different molecular model combinations: KBF-TIP4P/EP, OPLS-TIP4P/2005, OPLS-TIP4P/EP and Madrid-TIP4P/2005. All results are presented in Fig. 8 and compared with reference data [85].

Initially, the OPLS ion pairs model, in its original parameterization [70], was selected along with two water models: TIP4P/2005 and TIP4P/EP. This choice was made because this ions pair model has not been fitted to any of these two water models, unlike the KBF and MADRID ion pairs models, which were fitted to the TIP4P/EP [73] and TIP4P/2005 [71,72] water models, respectively. As shown in Fig. 8, the

OPLS-TIP4P/2005 model accurately represent density over the entire molality range, whereas the OPLS-TIP4P/EP model systematically underestimates it, with AAD being equal to 0.3 % and 2.5 %, respectively.

Unsurprisingly, both the KBF-TIP4P/EP and, to an even greater extent, the Madrid-TIP4P/2005 combinations predict brine densities very accurately, with respective AAD of 0.7 % and 0.05 %. These findings reinforce the notion that due of the complexity of the studied systems and the simplicity of the water molecule models used, ions in solution are challenging to model purely predictively. Ions in solution can form specific complexes or ionic associations, which are not fully captured by the dielectric constant alone. This is evident in the Madrid-TIP4P/2005 model [71], which successfully reproduces the thermophysical properties of NaCl brine solution by incorporating specific ingredients: scaled charge for Na⁺ and Cl⁻ to overcome the failures of TIP4P/2005 to represent the dielectric constant, optimized ion-ion and ion-solvent interactions to correct the Lorentz-Berthelot combining rules [86,87].

Therefore, for NaCl aqueous solutions, three model combinations were selected in the following: KBF-TIP4P/EP, OPLS-TIP4P/2005 and, Madrid-TIP4P/2005.

(b) H₂ + H₂O + NaCl system

Similarly to what was done in section 3.4.1.C, the results of the Henry constant method were compared with those of the Gibbs ensemble method for the Marx-TIP4P/2005-OPLS model at a temperature of 323 K, pressure up to 20 MPa and a salinity of 1 mol of NaCl per 1 kg of water, see Fig. 6. To compare the results, the partial molar volume of hydrogen in water at infinite dilution was fixed at 19 cm³/mol and was assumed to be independent of temperature and pressure [88], the fugacity coefficient was calculated from the reference equation of state for normal hydrogen developed by Leachman et al. [77] (using REFPROP software) (see discussion in section 4.1. H₂ partial molar volumes).

As found in pure water, very similar results were obtained using both methods, see Fig. 6, with an AAD = 4 %. This further confirms that the two approaches yield consistent results. Hence, for the same reason than in pure water, the following results will be expressed in terms of Henry's constant only, using the corresponding approach.

As a purely predictive test, Henry constant of hydrogen in NaCl brine

<figure><img src="image_11.png" /><figcaption>Fig. 8. Density of aqueous NaCl solution as a function of molality m (mol NaCl/kg H₂O) at a pressure of 0.1 MPa and temperatures of 298 K.</figcaption></figure>
model over the TIP4P/2005 model is not sufficient on its own for predicting the H₂ solubility.

Thus, even if very satisfying, the Marx-TIP4P/2005 and Marx-TIP4P/EP combinations yield systematic overestimations of the Henry constant of H₂ in pure water, with AADs of 6 % and 9 %, respectively. This may come from the bad evaluation of the cross-interaction parameters provide by the Lorentz-Berthelot rules, as reported by several authors for similar gas/water systems [82–84]. These studies converge on the fact that LB combining-rule underestimates gas/water polarizability. Hence introducing a positive deviation to the cross-energy Lennard-Jones interactions correct indirectly the polarizability, and leads to a better prediction of the excess chemical potential of gas in water, and consequently its solubility. As a test, cross interaction energies $ε_{ij}$ between H₂ and water, were multiplied by a constant $k_{ij}$ factor that were found to be optimal when equal to 1.05 for both Marx-TIP4P/2005 and Marx-TIP4P/EP. This optimized set of parameters yields AAD below 2 % and 4 % respectively and allows to reproduce perfectly the smoothed experimental data of the Henry constant (Akinfiev and Diамond's correlation) of H₂ in pure water over a very wide range of temperature.

### 3.4.2. H₂ + H₂O + NaCl system

(a) H₂O + NaCl system

As discussed in the previous section on the pure water, the solvent density has a significant impact on the hydrogen solubility. Therefore, to assess the Na⁺ and Cl⁻ ion pairs force fields, their impact on solvent density is crucial. For this purpose, the densities at T = 298 K and P = 0.1 MPa across various NaCl molalities (ranging from 1 to 5 mol/kg_w) were computed using different molecular model combinations: KBF-TIP4P/EP, OPLS-TIP4P/2005, OPLS-TIP4P/EP and Madrid-TIP4P/2005. All results are presented in Fig. 8 and compared with reference data [85].

Initially, the OPLS ion pairs model, in its original parameterization [70], was selected along with two water models: TIP4P/2005 and TIP4P/EP. This choice was made because this ions pair model has not been fitted to any of these two water models, unlike the KBF and MADRID ion pairs models, which were fitted to the TIP4P/EP [73] and TIP4P/2005 [71,72] water models, respectively. As shown in Fig. 8, the

OPLS-TIP4P/2005 model accurately represent density over the entire molality range, whereas the OPLS-TIP4P/EP model systematically underestimates it, with AAD being equal to 0.3 % and 2.5 %, respectively.

Unsurprisingly, both the KBF-TIP4P/EP and, to an even greater extent, the Madrid-TIP4P/2005 combinations predict brine densities very accurately, with respective AAD of 0.7 % and 0.05 %. These findings reinforce the notion that due of the complexity of the studied systems and the simplicity of the water molecule models used, ions in solution are challenging to model purely predictively. Ions in solution can form specific complexes or ionic associations, which are not fully captured by the dielectric constant alone. This is evident in the Madrid-TIP4P/2005 model [71], which successfully reproduces the thermophysical properties of NaCl brine solution by incorporating specific ingredients: scaled charge for Na⁺ and Cl⁻ to overcome the failures of TIP4P/2005 to represent the dielectric constant, optimized ion-ion and ion-solvent interactions to correct the Lorentz-Berthelot combining rules [86,87].

Therefore, for NaCl aqueous solutions, three model combinations were selected in the following: KBF-TIP4P/EP, OPLS-TIP4P/2005 and, Madrid-TIP4P/2005.

(b) H₂ + H₂O + NaCl system

Similarly to what was done in section 3.4.1.C, the results of the Henry constant method were compared with those of the Gibbs ensemble method for the Marx-TIP4P/2005-OPLS model at a temperature of 323 K, pressure up to 20 MPa and a salinity of 1 mol of NaCl per 1 kg of water, see Fig. 6. To compare the results, the partial molar volume of hydrogen in water at infinite dilution was fixed at 19 cm³/mol and was assumed to be independent of temperature and pressure [88], the fugacity coefficient was calculated from the reference equation of state for normal hydrogen developed by Leachman et al. [77] (using REFPROP software) (see discussion in section 4.1. H₂ partial molar volumes).

As found in pure water, very similar results were obtained using both methods, see Fig. 6, with an AAD = 4 %. This further confirms that the two approaches yield consistent results. Hence, for the same reason than in pure water, the following results will be expressed in terms of Henry's constant only, using the corresponding approach.

As a purely predictive test, Henry constant of hydrogen in NaCl brine

<figure><img src="image_11.png" /><figcaption>Fig. 8. Density of aqueous NaCl solution as a function of molality m (mol NaCl/kg H₂O) at a pressure of 0.1 MPa and temperatures of 298 K.</figcaption></figure>
were computed for various temperatures and NaCl molalities for the three combinations, Marx-TIP4P/2005-OPLS, Marx-TIP4P/EP-KBF and Marx-TIP4P/2005-Madrid, without any adjustments and compared with the most recent dataset [15]. As shown in Fig. 9A for two molalities, all models capture qualitatively well the salting out effect, i.e. the reduction of H₂ solubility (i.e. the increase of Henry constant) with salinity. It is interesting to point out that this salting out effect is both due to density and excess chemical potential contributions, see equation (2), the latter accounting for 80 % of the increase of the Henry's constant. Quantitatively, the Marx-TIP4P/2005-Madrid combination yields very good results while both Marx-TIP4P/2005-OPLS and Marx-TIP4P/EP-KBF tends to increasingly overestimate the H₂ Henry constant with salinity, yielding AAD relatively to experimental results equal to 4 %, 40 % and 30 %, respectively.

To further improve the results, as shown in Fig. 9B, one could use the $k_{ij}$ optimized on the H₂-pure water system, see section 3.4.1.C. By doing so, all results are improved compared to experimental data. In particular, the Marx-TIP4P/2005-Madrid combination yields excellent results, correcting its slight inadequacy at low salinity. Hence, in the next section, this model will be used to produce pseudo-experimental data of the Henry law constant of H₂ in brine for a large set of temperature and salinities.

Interestingly, even with an optimized *kij* between H₂ and water model both Marx-TIP4P/2005-OPLS and Marx-TIP4P/EP-KBF combinations still overestimate the Henry constant, see Fig. 9B. This indicate that the change in excess chemical potential of H₂ due to the salt in equation (2) is not well captured by these models as the density is well reproduced. To further elucidate this, the radial distribution function between hydrogen and ions have been computed and are shown in Fig. 10. Among others, the RDF clearly indicates that hydrogen has a much stronger affinity towards Cl⁻ ion than towards Na⁺ ion. Furthermore, the Marx-TIP4P/2005-Madrid combination exhibits a higher first peak value compared to the two other combinations. Hence, one option

to improve the results of Marx-TIP4P/2005-OPLS and Marx-TIP4P/EP-KBF would be to increase the Lennard-Jones cross interaction between hydrogen and Cl⁻, see Fig. 9B. It should be emphasized that a similar adjustment was proposed in the work of Lopez-Lazaro et al. [13], but applying a binary interaction parameter both with Cl⁻ and Na⁺. As shown in Fig. 9B and 9C, the Marx-TIP4P/EP-KBF combination used with a LJ cross-diameter parameter between hydrogen and Cl⁻ multiplied by 0.82 is able to provide results as very close to those obtained using the Marx-TIP4P/2005-Madrid combination. The salting-out effect due to the presence of NaCl (from pure water to saturated NaCl brine) and the effect of temperature (from 298 K to 453 K) are well captured by the both of models (Fig. 9C).

# 4. Thermodynamic modeling

In this section, the aim is to correlate solubility data from the literature and those obtained in this work, providing a model that can be used directly to calculate the solubility of H₂ in brine over a wide range of temperature, pressure and salinity.

## 4.1. Modeling approach

The proposed model uses an asymmetric approach (Krichevsky-Kasarnovsky type) to calculate the liquid-vapor equilibrium (gamma-phi), using a Henry's law constant (equation (1)) and an activity coefficient (γ) for the liquid phase, and a fugacity coefficient (phi) for the vapor phase. This approach was chosen in order to exploit directly the Henry constants predicted by molecular simulations from this work, and also because of its simplicity to implement (non-iterative for single-gas systems).

Hence, knowing these properties and Henry's constant $H_i$, the composition of solute i in the liquid phase $x_i$ can be determined. However, when dealing with high-pressure, the simple form of Henry's law

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** H₂ Henry coefficient at different salinities (A, B and C) and its solubility in terms of true molar fraction (D) in aqueous NaCl solutions (experimental vs. molecular simulation data).</figcaption></figure>

<figure><img src="image_13.png" /></figure>
were computed for various temperatures and NaCl molalities for the three combinations, Marx-TIP4P/2005-OPLS, Marx-TIP4P/EP-KBF and Marx-TIP4P/2005-Madrid, without any adjustments and compared with the most recent dataset [15]. As shown in Fig. 9A for two molalities, all models capture qualitatively well the salting out effect, i.e. the reduction of H₂ solubility (i.e. the increase of Henry constant) with salinity. It is interesting to point out that this salting out effect is both due to density and excess chemical potential contributions, see equation (2), the latter accounting for 80 % of the increase of the Henry's constant. Quantitatively, the Marx-TIP4P/2005-Madrid combination yields very good results while both Marx-TIP4P/2005-OPLS and Marx-TIP4P/EP-KBF tends to increasingly overestimate the H₂ Henry constant with salinity, yielding AAD relatively to experimental results equal to 4 %, 40 % and 30 %, respectively.

To further improve the results, as shown in Fig. 9B, one could use the $k_{ij}$ optimized on the H₂-pure water system, see section 3.4.1.C. By doing so, all results are improved compared to experimental data. In particular, the Marx-TIP4P/2005-Madrid combination yields excellent results, correcting its slight inadequacy at low salinity. Hence, in the next section, this model will be used to produce pseudo-experimental data of the Henry law constant of H₂ in brine for a large set of temperature and salinities.

Interestingly, even with an optimized *kij* between H₂ and water model both Marx-TIP4P/2005-OPLS and Marx-TIP4P/EP-KBF combinations still overestimate the Henry constant, see Fig. 9B. This indicate that the change in excess chemical potential of H₂ due to the salt in equation (2) is not well captured by these models as the density is well reproduced. To further elucidate this, the radial distribution function between hydrogen and ions have been computed and are shown in Fig. 10. Among others, the RDF clearly indicates that hydrogen has a much stronger affinity towards Cl⁻ ion than towards Na⁺ ion. Furthermore, the Marx-TIP4P/2005-Madrid combination exhibits a higher first peak value compared to the two other combinations. Hence, one option

to improve the results of Marx-TIP4P/2005-OPLS and Marx-TIP4P/EP-KBF would be to increase the Lennard-Jones cross interaction between hydrogen and Cl⁻, see Fig. 9B. It should be emphasized that a similar adjustment was proposed in the work of Lopez-Lazaro et al. [13], but applying a binary interaction parameter both with Cl⁻ and Na⁺. As shown in Fig. 9B and 9C, the Marx-TIP4P/EP-KBF combination used with a LJ cross-diameter parameter between hydrogen and Cl⁻ multiplied by 0.82 is able to provide results as very close to those obtained using the Marx-TIP4P/2005-Madrid combination. The salting-out effect due to the presence of NaCl (from pure water to saturated NaCl brine) and the effect of temperature (from 298 K to 453 K) are well captured by the both of models (Fig. 9C).

# 4. Thermodynamic modeling

In this section, the aim is to correlate solubility data from the literature and those obtained in this work, providing a model that can be used directly to calculate the solubility of H₂ in brine over a wide range of temperature, pressure and salinity.

## 4.1. Modeling approach

The proposed model uses an asymmetric approach (Krichevsky-Kasarnovsky type) to calculate the liquid-vapor equilibrium (gamma-phi), using a Henry's law constant (equation (1)) and an activity coefficient (γ) for the liquid phase, and a fugacity coefficient (phi) for the vapor phase. This approach was chosen in order to exploit directly the Henry constants predicted by molecular simulations from this work, and also because of its simplicity to implement (non-iterative for single-gas systems).

Hence, knowing these properties and Henry's constant $H_i$, the composition of solute i in the liquid phase $x_i$ can be determined. However, when dealing with high-pressure, the simple form of Henry's law

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** H₂ Henry coefficient at different salinities (A, B and C) and its solubility in terms of true molar fraction (D) in aqueous NaCl solutions (experimental vs. molecular simulation data).</figcaption></figure>

<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 10. Radial distribution function of H₂O-NaCl (A) and H₂-H₂O-NaCl (B) systems at 323 K and 20 MPa.</figcaption></figure>

(equation (1) is no longer applicable. In 1935, Krichevsky and Kasarnovsky (KK) solved this problem by introducing the Poynting factor (equations (3) and (4), which corrects the effect of pressure variation from the solvent vapor pressure at $x_i = 0$ to a pressure P at $x_i > 0$. Hence, the KK approach is used herein to convert the Henry's constant to molar composition (see section 3.1.1.B) and also to establish a reliable correlation model for the solubility of H₂ in water and brine (this section). Moreover, to consider the effect of the electrolyte (NaCl) on gas solubility, an activity coefficient $\gamma_{H_2}$ is also considered. The resulting model for H₂ as solute is as follows:

$$x_{H_2}(\mathbf{T}, \mathbf{P}, m_s) = \frac{y_{H_2} P \varphi_{H_2}}{H_{H_2}(T, P_w^{sat}, m_s) \Pi_{H_2}} = \frac{y_{H_2} P \varphi_{H_2}}{H_{H_2}^0(T, P_w^{sat}) \gamma_{H_2}(T, m_s) \Pi_{H_2}} \quad (9)$$

$$\Pi_{H_2} = \exp\left(\frac{\vartheta_{H_2, w}^{\infty}}{RT}(P - P_w^{sat})\right) \quad (10)$$

where T, P and $m_s$ are the temperature (K), total pressure (bars) and salt molality $m_s$ (molNaCl/kgwater), respectively; $x_{H_2}$ and $y_{H_2}$ are the molar composition of H₂ in the liquid and vapor phase, respectively. $\vartheta_{H_2}$ is its fugacity coefficient in the vapor phase, $H_{H_2}^0$ is the Henry's constant of H₂

in pure water, $\gamma_{H_2}$ is its activity coefficient in the liquid phase, $\Pi_{H_2}$ is the Poynting factor (pressure correction), $P_w^{sat}$ is the vapor pressure of the solvent (water) and $\vartheta_{H_2, w}^{\infty}$ is the partial molar volume of H₂ in pure water at infinite dilution.

### H₂ fugacity coefficient

Since the water content in the vapor phase is very low, the fugacity coefficient of hydrogen in the vapor phase is replaced by that of pure H₂. It's calculated from the reference equation of state for normal hydrogen developed by Leachman et al. [77] (using REFPROP software), but other simpler models can be used instead (e.g. Peng-Robinson EoS [89], or virial-type EoS [90]). The H₂ composition in the gas phase is calculated using $y_{H_2} = 1 - y_{H_2O}$, and the water content in the gas phase ($y_{H_2O}$) is estimated using Raoult's law ($y_{H_2O} = \frac{p_{sat}}{P}$) where saturated vapor pressure of water ($P_w^{sat}$) is calculated using the IAPWS formulation [80].

### H₂ activity coefficient

In non-electrolytic solutions, as H₂ solubility is always low, its activity coefficient in pure water is very close to 1. On the other hand, the addition of NaCl, considerably changes the structure of the solution through the solvation of Na⁺ and Cl⁻ ions by water molecules, and consequently hydrogen's activity increases and its solubility decreases.
<figure><img src="image_14.png" /><figcaption>Fig. 10. Radial distribution function of H₂O-NaCl (A) and H₂-H₂O-NaCl (B) systems at 323 K and 20 MPa.</figcaption></figure>

(equation (1) is no longer applicable. In 1935, Krichevsky and Kasarnovsky (KK) solved this problem by introducing the Poynting factor (equations (3) and (4), which corrects the effect of pressure variation from the solvent vapor pressure at $x_i = 0$ to a pressure P at $x_i > 0$. Hence, the KK approach is used herein to convert the Henry's constant to molar composition (see section 3.1.1.B) and also to establish a reliable correlation model for the solubility of H₂ in water and brine (this section). Moreover, to consider the effect of the electrolyte (NaCl) on gas solubility, an activity coefficient $γ_{H₂}$ is also considered. The resulting model for H₂ as solute is as follows:

$$x_{H_2}(\mathbf{T}, \mathbf{P}, m_s) = \frac{y_{H_2} P \varphi_{H_2}}{H_{H_2}(T, P_w^{sat}, m_s) \Pi_{H_2}} = \frac{y_{H_2} P \varphi_{H_2}}{H_{H_2}^0(T, P_w^{sat}) \gamma_{H_2}(T, m_s) \Pi_{H_2}} \quad (9)$$

$$\Pi_{H_2} = \exp\left(\frac{\vartheta_{H_2, w}^{\infty}}{RT}(P - P_w^{sat})\right) \quad (10)$$

where T, P and $m_s$ are the temperature (K), total pressure (bars) and salt molality $m_s$ (molNaCl/kgwater), respectively; $x_{H_2}$ and $γ_{H_2}$ are the molar composition of H₂ in the liquid and vapor phase, respectively. $γ_{H_2}$ is its fugacity coefficient in the vapor phase, $H_{H_2}^0$ is the Henry's constant of H₂

in pure water, $γ_{H_2}$ is its activity coefficient in the liquid phase, $\Pi_{H_2}$ is the Poynting factor (pressure correction), $P_w^{sat}$ is the vapor pressure of the solvent (water) and $\vartheta_{H_2, w}^{\infty}$ is the partial molar volume of H₂ in pure water at infinite dilution.

### H₂ fugacity coefficient

Since the water content in the vapor phase is very low, the fugacity coefficient of hydrogen in the vapor phase is replaced by that of pure H₂. It's calculated from the reference equation of state for normal hydrogen developed by Leachman et al. [77] (using REFPROP software), but other simpler models can be used instead (e.g. Peng-Robinson EoS [89], or virial-type EoS [90]). The H₂ composition in the gas phase is calculated using $y_{H_2} = 1 - y_{H_2O}$, and the water content in the gas phase ($y_{H_2O}$) is estimated using Raoult's law ($y_{H_2O} = \frac{p_{sat}}{P}$) where saturated vapor pressure of water ($P_w^{sat}$) is calculated using the IAPWS formulation [80].

### H₂ activity coefficient

In non-electrolytic solutions, as H₂ solubility is always low, its activity coefficient in pure water is very close to 1. On the other hand, the addition of NaCl, considerably changes the structure of the solution through the solvation of Na⁺ and Cl⁻ ions by water molecules, and consequently hydrogen's activity increases and its solubility decreases.
In this work, a simple linear dependence is considered between the activity coefficient and the salt molality, as follows:

$$ \gamma_{H_2} = \exp(K_s m_s) \qquad (11) $$

$$ K_s (\text{kg.mol}^{-1}) = A_0 T + \frac{A_1}{T} + A_2 \frac{T}{m_s} + A_3 + A_4 m_s T \qquad (12) $$

This Sechenov-type formulation for expressing salting-out effect on gas solubility, which is also similar to the first-order Pitzer model, has been used for the same system (H₂-H₂O-NaCl) by Li et al. 2018 [34] and Torín-Ollarves and Trusler 2021 [12]. The adjustable coefficient A0-A4 are listed in Table 3.

### H₂ Henry constant model

The Henry's constant $H_{H_2}$ is calculated using the Akinfiev and Diamond model [35] (equation (13)):

$$ \ln H_{iAkinfiev} (T, P_w^{sat}) = (1 - \varepsilon) \ln f_w^0 + \varepsilon \ln \left( \frac{RT}{M_w \cdot \rho_w^0} \right) + 2 \Delta B_{\rho} \cdot \rho_w^0 \\ 2 \Delta B_{\rho} = a + b \left( \frac{10^3}{T} \right)^{0.5} \qquad (13) $$

where: $f_w^0 = P_w \cdot \rho_w^0$ and $\rho_w^0$ are the fugacity (bars) and density (g/cm³) of the pure solvent (H₂O) at given temperature and pressure conditions. $\rho_w^0$ is calculated using the NIST REFPROP software. $\Delta B = \Delta B_{\rho}$. M = $B_{\text{solute/water}} - B_{\text{water/water}}$ (cm³/mol) represents the difference in interaction between different molecules and between identical solvent molecules. R is the ideal gas constant R = 83.1441 (cm³·bar/mol·K). $\varepsilon$, a (cm³/g) and b (cm³·K⁰.⁵/g) are the empirical parameters of the equation of state and are given respectively: 0.3090, -8.4596 and 10.8301 for hydrogen. A simplified polynomial form given by the equation (14) was developed and is used to calculate Henry's constant (in MPa) instead of using equation (13). The parameters $B_0 - B_4$ are listed in Table 3.

$$ \ln H_{H_2Akinfiev} (T) = B_0 T^4 + B_1 T^3 + B_2 T^2 + B_3 T + B_4 \qquad (14) $$

From equation (9), the H₂ Henry's constant in brine can be calculated by the following equation:

$$ H_{H_2} (T, P_w^{sat}, m_s) = H_{H_2Akinfiev}^0 (T, P_w^{sat}) \gamma_{H_2} = H_{H_2Akinfiev}^0 (T, P_w^{sat}) \exp(K_s m_s) \qquad (15) $$

### H₂ partial molar volumes $\delta_{H_2,w}^{\infty}$

There are various methods to measure directly or indirectly the partial molar volumes $\delta_{H_2,w}^{\infty}$ of gases and liquids, including determining the pressure dependence of the Henry constants (Poynting's factor), measuring the density, and using dilatometric methods. In the literature, there are only a few experimental data on the partial molar volume of H₂ in pure water under ambient conditions of pressure and temperature with some discrepancies. As mentioned earlier in section 3.1.1.C, the partial molar volume at infinite dilution can also be obtained by molecular simulation, but this type of calculation is time-consuming. In this work, some data have been predicted just to get the order of magnitude and the temperature dependence, while optimization of these results is outside the scope of this work. The results obtained are illustrated in Fig. 11 and compared with experimental data from the literature and some models and correlations.

In 1988, Tanger and Helgeson developed the Helgeson-Kirkham-Flowers (HKF) equation of state, which provides excellent predictions of the thermodynamic properties of aqueous ions, including the partial molar volume at infinite dilution over a wide range of conditions. However, when extended to neutral aqueous species (Shock et al., 1989), the HKF model fails to correctly predict the behavior of non-electrolyte solutes in the near-critical and supercritical regions of water. To address this problem, Akinfiev and Diamond 2002 proposed a new virial-type equation of state (EoS) to describe the thermodynamic properties of aqueous non-electrolytes at infinite dilution (Henry's law constant, partial molar volume, entropy, apparent heat capacity, etc). The partial molar volume using this model is given by the following equation:

$$ \delta_i^{\infty} = \frac{\partial \rho_{water}^0}{\partial P} |_T = \delta_{water}^{\infty} (1 - \varepsilon) + \varepsilon R T \frac{1}{\rho_w^0} \frac{\partial \rho_w^0}{\partial P} + R \left( a T \frac{\partial \rho_w^0}{\partial P} + b \cdot 10^{1.5} \cdot T^{0.5} \frac{\partial \rho_w^0}{\partial P} \right) \qquad (16) $$

Note that the derivative of the saturated liquid density with respect to pressure is calculated using the reference Helmholtz equation of state for water of Wagner and Pruss (2002) available in the REFPROP software. As shown in Fig. 11, the results predicted from molecular simulations, the values obtained from the EOS of Akinfiev and Diamond and experimental literature data are in good agreement. As observed in Henry's constant, the partial molar volume also has an optimum (minimum) at around 330 K. This trend may be explained by the fact that, in liquid water at low temperature, all water molecules are effectively bonded to neighboring water molecules with no free water molecules present. The hydrogen bonds in liquid water form chains that connect the water molecules separated by large distances, although the tetrahedral grouping is only locally present and decreases with increasing temperature [88].

On the other hand, Li et al. 2018 [34], developed a temperature- and pressure-dependent correlation for H₂ partial molar volume applicable for temperatures up to 373 K. Torín-Ollarves and Trusler 2021 [12] developed a temperature-dependent correlation applicable up to 423 K and without pressure dependence. As shown in Fig. 11, these two correlations give results of the order of magnitude of the limited experimental data and molecular simulation predictions. However, from a fundamental point of view, the assumption made in these works as well as that the activity coefficient of H₂ in pure water is equal to one ($\gamma_{H_2-H_2O} = 1$), is not entirely valid especially at very high pressure where the solubility of hydrogen becomes relatively non-negligible, i. e. at 20 MPa / 323 K, the H₂ activity coefficient is estimated by the KK approach (considering experimental solubility in pure water) to be 1.012 and 0.960 for H₂ partial molar volumes of 19 and 26 cm³/mol, respectively. Hence, this assumption can generate errors on the partial molar volume, since it is obtained by adjustment on data of solubilities at very high pressure (the error due to the assumption that $\gamma_{H_2-H_2O} = 1$ is compensated in the partial molar volume obtained). From a practical point of view, given that very few data exist on the partial molar volume of H₂ in water, and no experimental data exist above 298 K, we are "constrained" to fix a value of $\delta_{H_2,w}^{\infty}$, or to consider a dependence on temperature and possibly pressure, by adjusting it on solubility data, if the aim of the model is rather to correlate solubility data, which is the case in their works [12,34] as well as in this work. Based on this discussion, an adjusted constant value of the partial molar volume at infinite dilution $\delta_{H_2,w}^{\infty}$ of 19 cm³/mol was obtained and considered in the

<table><caption>**Table 3**

Optimized coefficients of the developed model (Equations 9–12, 14–15) for the calculation of H₂ solubility in NaCl brine.</caption><thead><tr><th>Parameters</th><th>A<sub>0</sub></th><th>A<sub>1</sub></th><th>A<sub>2</sub></th><th>A<sub>3</sub></th><th>A<sub>4</sub></th><th>B<sub>0</sub></th><th>B<sub>1</sub></th><th>B<sub>2</sub></th><th>B<sub>3</sub></th><th>B<sub>4</sub></th><th>δ<sup>∞</sup><sub>H<sub>2,w</sub></sub>(cm<sup>3</sup>/mol)</th></tr></thead><tbody><tr><td>Value</td><td>1.59054<br>E-03</td><td>2.63936<br>+ 02</td><td>3.79192<br>E-05</td><td>-1.14604<br>+ 00</td><td>-1.18419<br>E-05</td><td>-4.40347<br>E-10</td><td>8.33536<br>E-07</td><td>-6.05506<br>E-04</td><td>1.90837<br>E-01</td><td>-1.27991<br>+ 01</td><td>19</td></tr></tbody></table>
In this work, a simple linear dependence is considered between the activity coefficient and the salt molality, as follows:

$$ \gamma_{H_2} = \exp(K_s m_s) \quad (11) $$

$$ K_s (\text{kg.mol}^{-1}) = A_0 T + \frac{A_1}{T} + A_2 \frac{T}{m_s} + A_3 + A_4 m_s T \quad (12) $$

This Sechenov-type formulation for expressing salting-out effect on gas solubility, which is also similar to the first-order Pitzer model, has been used for the same system (H₂-H₂O-NaCl) by Li et al. 2018 [34] and Torín-Ollarves and Trusler 2021 [12]. The adjustable coefficient A0-A4 are listed in Table 3.

### H₂ Henry constant model

The Henry's constant $H_{H_2}$ is calculated using the Akinfiev and Diamond model [35] (equation (13)):

$$ \ln H_{iAkinfiev} (T, P_w^{sat}) = (1 - \varepsilon) \ln f_w^0 + \varepsilon \ln \left( \frac{RT}{M_w \cdot \rho_w^0} \right) + 2 \Delta B_{\rho} \cdot \rho_w^0 \\ 2 \Delta B_{\rho} = a + b \left( \frac{10^3}{T} \right)^{0.5} \quad (13) $$

where: $f_w^0 = P_w \cdot \rho_w^0$ and $\rho_w^0$ are the fugacity (bars) and density (g/cm³) of the pure solvent (H₂O) at given temperature and pressure conditions. $\rho_w^0$ is calculated using the NIST REFPROP software. $\Delta B = \Delta B_{\rho}$. M = $B_{\text{solute/water}} - B_{\text{water/water}}$ (cm³/mol) represents the difference in interaction between different molecules and between identical solvent molecules. R is the ideal gas constant R = 83.1441 (cm³·bar/mol·K), $\varepsilon$, a (cm³/g) and b (cm³·K⁰.⁵/g) are the empirical parameters of the equation of state and are given respectively: 0.3090, -8.4596 and 10.8301 for hydrogen. A simplified polynomial form given by the equation (14) was developed and is used to calculate Henry's constant (in MPa) instead of using equation (13). The parameters $B_0 - B_4$ are listed in Table 3.

$$ \ln H_{H_2Akinfiev} (T) = B_0 T^4 + B_1 T^3 + B_2 T^2 + B_3 T + B_4 \quad (14) $$

From equation (9), the H₂ Henry's constant in brine can be calculated by the following equation:

$$ H_{H_2} (T, P_w^{sat}, m_s) = H_{H_2Akinfiev}^0 (T, P_w^{sat}) \gamma_{H_2} = H_{H_2Akinfiev}^0 (T, P_w^{sat}) \exp(K_s m_s) \quad (15) $$

### H₂ partial molar volumes $\delta_{H_2,w}^{\infty}$

There are various methods to measure directly or indirectly the partial molar volumes $\delta_{H_2,w}^{\infty}$ of gases and liquids, including determining the pressure dependence of the Henry constants (Poynting's factor), measuring the density, and using dilatometric methods. In the literature, there are only a few experimental data on the partial molar volume of H₂ in pure water under ambient conditions of pressure and temperature with some discrepancies. As mentioned earlier in section 3.1.1.C, the partial molar volume at infinite dilution can also be obtained by molecular simulation, but this type of calculation is time-consuming. In this work, some data have been predicted just to get the order of magnitude and the temperature dependence, while optimization of these results is outside the scope of this work. The results obtained are illustrated in Fig. 11 and compared with experimental data from the literature and some models and correlations.

In 1988, Tanger and Helgeson developed the Helgeson-Kirkham-Flowers (HKF) equation of state, which provides excellent predictions of the thermodynamic properties of aqueous ions, including the partial molar volume at infinite dilution over a wide range of conditions. However, when extended to neutral aqueous species (Shock et al., 1989), the HKF model fails to correctly predict the behavior of non-electrolyte solutes in the near-critical and supercritical regions of water. To address this problem, Akinfiev and Diamond 2002 proposed a new virial-type equation of state (EoS) to describe the thermodynamic properties of aqueous non-electrolytes at infinite dilution (Henry's law constant, partial molar volume, entropy, apparent heat capacity, etc). The partial molar volume using this model is given by the following equation:

$$ \delta_i^\infty = \frac{\partial \rho_{water}^0}{\partial \mathbf{P}}|_T = \delta_{water}^\infty (1 - \varepsilon) + \varepsilon R T \frac{1}{\rho_w^0} \frac{\partial \rho_w^0}{\partial \mathbf{P}} + R \left( a \mathbf{T} \frac{\partial \rho_w^0}{\partial \mathbf{P}} + b \cdot 10^{1.5} \mathbf{T}^{0.5} \frac{\partial \rho_w^0}{\partial \mathbf{P}} \right) \quad (16) $$

Note that the derivative of the saturated liquid density with respect to pressure is calculated using the reference Helmholtz equation of state for water of Wagner and Pruss (2002) available in the REFPROP software. As shown in Fig. 11, the results predicted from molecular simulations, the values obtained from the EOS of Akinfiev and Diamond and experimental literature data are in good agreement. As observed in Henry's constant, the partial molar volume also has an optimum (minimum) at around 330 K. This trend may be explained by the fact that, in liquid water at low temperature, all water molecules are effectively bonded to neighboring water molecules with no free water molecules present. The hydrogen bonds in liquid water form chains that connect the water molecules separated by large distances, although the tetrahedral grouping is only locally present and decreases with increasing temperature [88].

On the other hand, Li et al. 2018 [34], developed a temperature- and pressure-dependent correlation for H₂ partial molar volume applicable for temperatures up to 373 K. Torín-Ollarves and Trusler 2021 [12] developed a temperature-dependent correlation applicable up to 423 K and without pressure dependence. As shown in Fig. 11, these two correlations give results of the order of magnitude of the limited experimental data and molecular simulation predictions. However, from a fundamental point of view, the assumption made in these works as well as that the activity coefficient of H₂ in pure water is equal to one ($\gamma_{H_2-H_2O} = 1$), is not entirely valid especially at very high pressure where the solubility of hydrogen becomes relatively non-negligible, i. e. at 20 MPa / 323 K, the H₂ activity coefficient is estimated by the KK approach (considering experimental solubility in pure water) to be 1.012 and 0.960 for H₂ partial molar volumes of 19 and 26 cm³/mol, respectively. Hence, this assumption can generate errors on the partial molar volume, since it is obtained by adjustment on data of solubilities at very high pressure (the error due to the assumption that $\gamma_{H_2-H_2O} = 1$ is compensated in the partial molar volume obtained). From a practical point of view, given that very few data exist on the partial molar volume of H₂ in water, and no experimental data exist above 298 K, we are "constrained" to fix a value of $\delta_{H_2,w}^{\infty}$, or to consider a dependence on temperature and possibly pressure, by adjusting it on solubility data, if the aim of the model is rather to correlate solubility data, which is the case in their works [12,34] as well as in this work. Based on this discussion, an adjusted constant value of the partial molar volume at infinite dilution $\delta_{H_2,w}^{\infty}$ of 19 cm³/mol was obtained and considered in the

<table><caption>**Table 3**

Optimized coefficients of the developed model (Equations 9–12, 14–15) for the calculation of H₂ solubility in NaCl brine.</caption><thead><tr><th>Parameters</th><th>A<sub>0</sub></th><th>A<sub>1</sub></th><th>A<sub>2</sub></th><th>A<sub>3</sub></th><th>A<sub>4</sub></th><th>B<sub>0</sub></th><th>B<sub>1</sub></th><th>B<sub>2</sub></th><th>B<sub>3</sub></th><th>B<sub>4</sub></th><th>δ<sup>∞</sup><sub>H<sub>2</sub>w</sub>(cm<sup>3</sup>/mol)</th></tr></thead><tbody><tr><td>Value</td><td>1.59054<br>E-03</td><td>2.63936<br>+ 02</td><td>3.79192<br>E-05</td><td>-1.14604<br>+ 00</td><td>-1.18419<br>E-05</td><td>-4.40347<br>E-10</td><td>8.33536<br>E-07</td><td>-6.05506<br>E-04</td><td>1.90837<br>E-01</td><td>-1.27991<br>+ 01</td><td>19</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>Fig. 11. Partial molar volume of H₂ in water at infinite dilution versus temperature at 0.1 MPa.</figcaption></figure>

proposed model.

## 4.2. Modeling results

In equation (11) and (12), the Sechenov salting-out coefficient, $K_s$, depends on temperature and salinity, and is obtained by fitting on H₂ solubility data in brine from literature and molecular simulation obtained in this work. By considering a dependence only on temperature, the results are less accurate, which is the reason why a dependence on salinity has also been added, consistently with some observations made for other gases (CO₂ [91] and CH₄ [92]) where salting-out was found to be also dependent on salinity, especially at high salinities. The resulting model is valid for temperatures ranging from 298.15 K to 453.15 K, pressures up to 100 MPa, and NaCl salinities up to 6 mol/kgw. The optimized model parameters are listed in Table 3.

The comparison of computed results from the model with experimental data from the literature and molecular simulation data from this work, is illustrated in Fig. 12. The effects of temperature and salinity on H₂ Henry's constant are very well captured by the model presented here, as shown in Fig. 12A. Henry's constant increases with temperature until the minimum solubility temperature ($T_{X/H_2 \text{min}}$) is reached, then decreases with increasing temperature. This effect is probably due to the hydrogen bonds (as discussed above for the partial molar volume) formed between water molecules, which are more important at lower temperatures. These bonds decrease as the temperature increases, thus the Henry's constant decreases (and solubility increases) at temperatures higher than $T_{X/H_2 \text{min}}$.

Introducing salt into the water decrease the H₂ solubility (salting-out effect) in proportion to the salt concentration, up to saturation of the solution. In Fig. 12B, 12C and 12D, the solubility in terms of mole fraction is plotted as a function of temperature at fixed pressure and salinity, as a function of pressure at fixed temperature and salinity, and as a function of salinity at fixed temperature and pressure, respectively. The model reproduces the solubility of H₂ in pure water and up to saturated NaCl brine (6 mol/kgw) with an AAD of less than 5 %.

According to the model and the data, one can observe that the minimum solubility temperature moves towards lower temperatures with increasing salinity, until there is no longer any minimum at high salinities.

In addition, the model was assessed on data that were not included in the parameterization, especially the data on H₂ solubility in pure water at high temperature measured by Torín-Ollarves and Trusler 2021 [12]. As shown in Fig. 13, at 423 K, the model agrees very well with the experimental data, especially at high pressure, where the experimental uncertainty is less significant (the higher the solubility, the easier the measurement). Finally, as the calculated solubilities are in very good agreement with experimental data from the literature and molecular simulation data (in pure water and brine), smoothed H₂ solubility values in water and brine (in terms of Henry constant and molar fraction) were generated at different temperatures (298–453 K), pressures (up to 100 MPa) and NaCl molalities (up to 6 mol NaCl/kg water) and are listed in the supplementary materials.

## 5. Summary and conclusions

In this investigation, new equilibrium data for the H₂/brine system under high-pressure conditions (up to 1000 bar), temperatures from 298 to 453 K, and salinities up to 6 mol/kgw of NaCl were computed using Monte Carlo simulation. The study was aimed at understanding the influence of bulk and excess properties representation on H₂ solubility predictions. Key findings include:

The investigation began with an assessment of the dielectric constant's impact by employing various salt models alongside different water models to compute brine density. Despite the critical role of the dielectric constant in electrolytes, it alone was insufficient to model the thermodynamic properties of electrolytic solutions accurately, possibly due to the Lorentz-Berthelot combining rules for ion-ion and ion-solvent interactions.

A comparative analysis of different combinations of force fields for H₂ and NaCl brine was undertaken. Notably, three combinations,
<figure><img src="image_15.png" /><figcaption>Fig. 11. Partial molar volume of H₂ in water at infinite dilution versus temperature at 0.1 MPa.</figcaption></figure>

proposed model.

## 4.2. Modeling results

In equation (11) and (12), the Sechenov salting-out coefficient, $K_s$, depends on temperature and salinity, and is obtained by fitting on H₂ solubility data in brine from literature and molecular simulation obtained in this work. By considering a dependence only on temperature, the results are less accurate, which is the reason why a dependence on salinity has also been added, consistently with some observations made for other gases (CO₂ [91] and CH₄ [92]) where salting-out was found to be also dependent on salinity, especially at high salinities. The resulting model is valid for temperatures ranging from 298.15 K to 453.15 K, pressures up to 100 MPa, and NaCl salinities up to 6 mol/kgw. The optimized model parameters are listed in Table 3.

The comparison of computed results from the model with experimental data from the literature and molecular simulation data from this work, is illustrated in Fig. 12. The effects of temperature and salinity on H₂ Henry's constant are very well captured by the model presented here, as shown in Fig. 12A. Henry's constant increases with temperature until the minimum solubility temperature ($T_{X/H_2_{min}}$) is reached, then decreases with increasing temperature. This effect is probably due to the hydrogen bonds (as discussed above for the partial molar volume) formed between water molecules, which are more important at lower temperatures. These bonds decrease as the temperature increases, thus the Henry's constant decreases (and solubility increases) at temperatures higher than $T_{X/H_2_{min}}$.

Introducing salt into the water decrease the H₂ solubility (salting-out effect) in proportion to the salt concentration, up to saturation of the solution. In Fig. 12B, 12C and 12D, the solubility in terms of mole fraction is plotted as a function of temperature at fixed pressure and salinity, as a function of pressure at fixed temperature and salinity, and as a function of salinity at fixed temperature and pressure, respectively. The model reproduces the solubility of H₂ in pure water and up to saturated NaCl brine (6 mol/kgw) with an AAD of less than 5 %.

According to the model and the data, one can observe that the minimum solubility temperature moves towards lower temperatures with increasing salinity, until there is no longer any minimum at high salinities.

In addition, the model was assessed on data that were not included in the parameterization, especially the data on H₂ solubility in pure water at high temperature measured by Torín-Ollarves and Trusler 2021 [12]. As shown in Fig. 13, at 423 K, the model agrees very well with the experimental data, especially at high pressure, where the experimental uncertainty is less significant (the higher the solubility, the easier the measurement). Finally, as the calculated solubilities are in very good agreement with experimental data from the literature and molecular simulation data (in pure water and brine), smoothed H₂ solubility values in water and brine (in terms of Henry constant and molar fraction) were generated at different temperatures (298–453 K), pressures (up to 100 MPa) and NaCl molalities (up to 6 mol NaCl/kg water) and are listed in the supplementary materials.

## 5. Summary and conclusions

In this investigation, new equilibrium data for the H₂/brine system under high-pressure conditions (up to 1000 bar), temperatures from 298 to 453 K, and salinities up to 6 mol/kgw of NaCl were computed using Monte Carlo simulation. The study was aimed at understanding the influence of bulk and excess properties representation on H₂ solubility predictions. Key findings include:

The investigation began with an assessment of the dielectric constant's impact by employing various salt models alongside different water models to compute brine density. Despite the critical role of the dielectric constant in electrolytes, it alone was insufficient to model the thermodynamic properties of electrolytic solutions accurately, possibly due to the Lorentz-Berthelot combining rules for ion-ion and ion-solvent interactions.

A comparative analysis of different combinations of force fields for H₂ and NaCl brine was undertaken. Notably, three combinations,
<figure><img src="image_16.png" /><figcaption>**Fig. 12.** H₂ solubility in water and NaCl brine (Experimental and molecular simulation data vs. developed model): (A) H₂ Henry constant at different temperatures and salinity; (B) Solubility versus temperature; (C) Solubility versus pressure and (D) Solubility versus salinity.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 13.** Prediction of H₂ solubility in water and NaCl brine at 323 K (A) and 423 K (B) versus pressure: Experimental and molecular simulation data vs. developed model (Equations 9–12, 14–15).</figcaption></figure>

namely Marx-TIP4P/2005-OPLS, Marx-TIP4P/EP-KBF, and Marx-TIP4P/2005-Madrid, were able to reproduce the solubility of H₂ in water and brine. The representation of H₂ solubility in water (Marx-

TIP4P/2005 and Marx-TIP4P/EP) and brine (Marx-TIP4P/2005-OPLS and Marx-TIP4P/EP-KBF) was enhanced through the utilization of the radial distribution function, assessing the solvation structure of the H₂/
<figure><img src="image_16.png" /><figcaption>**Fig. 12.** H₂ solubility in water and NaCl brine (Experimental and molecular simulation data vs. developed model): (A) H₂ Henry constant at different temperatures and salinity; (B) Solubility versus temperature; (C) Solubility versus pressure and (D) Solubility versus salinity.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 13.** Prediction of H₂ solubility in water and NaCl brine at 323 K (A) and 423 K (B) versus pressure: Experimental and molecular simulation data vs. developed model (Equations 9–12, 14–15).</figcaption></figure>

namely Marx-TIP4P/2005-OPLS, Marx-TIP4P/EP-KBF, and Marx-TIP4P/2005-Madrid, were able to reproduce the solubility of H₂ in water and brine. The representation of H₂ solubility in water (Marx-

TIP4P/2005 and Marx-TIP4P/EP) and brine (Marx-TIP4P/2005-OPLS and Marx-TIP4P/EP-KBF) was enhanced through the utilization of the radial distribution function, assessing the solvation structure of the H₂/
brine mixtures. This structure is predominantly influenced by H₂-H₂O and H₂-Cl⁻ interactions. Consequently, a constant binary interaction parameter on energy for H₂-H₂O in the Marx-TIP4P/2005 and Marx-TIP4P/EP models, and on volume for H₂-Cl⁻ in the Marx-TIP4P/EP-KBF, was implemented. Following this adjustment, both model combinations provided excellent representations of H₂ solubility in NaCl brine.

Despite both models' accurate estimations of hydrogen solubility in brine, the predictive model Marx-TIP4P/2005-Madrid was chosen to generate data, as it was not directly adjusted to the ternary system. This model underwent validation against experimental literature data for H₂ solubility in brine (for pressures up to 200 bar, temperatures between 298 K and 373 K, and salinities up to 4 mol/kg NaCl), extending the data for H₂/brine to higher temperatures (up to 453 K), pressures (up to 1000 bar), and salinities (up to 6 mol/kg NaCl).

Both experimental and molecular simulations data were leveraged to refine parameters of a thermodynamic model. The asymmetric (gamma-phi) approach of the Krichevsky-Kasarnovsky type, chosen for its non-iterative simplicity in single-gas systems, facilitated the calculation of the liquid-vapor equilibrium. The model showcased very good results for H₂ solubility in water and brine across varying temperatures (298–453 K), pressures (up to 1000 bar), and NaCl molalities (up to 6 mol NaCl/kg water), proving invaluable for evaluating H₂ storage capacity, dissolution losses, biotic and abiotic reactivity, and caprock integrity.

Future work may delve into the interactions of H₂ with other sub-surface gases (CO₂, CH₄, etc.) or the impact of clay on its solubility, building on the current findings.

### 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, Writing – original draft, Validation, Software, Methodology, Investigation, Formal analysis. **Pierre Cézac:** Writing – review & editing, Supervision, Resources, Project administration, Funding acquisition. **Guillaume Galliéro:** Writing – review & editing, Writing – original draft, Supervision, Resources, Project administration, Methodology, Funding acquisition, Formal analysis, Conceptualization. **Salaheddine Chabab:** Writing – review & editing, Writing – original draft, Supervision, Project administration, Methodology, Investigation, 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.

### 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 grant 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 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.mollicq.2024.124497.

### References

[1] N. Heinemann, J. Alcalde, J.M. Miocic, S.J. Hangx, J. Kallmeyer, C. Ostertag-Henning, A. Hassanpouryouzband, E.M. Thaysen, G.J. Strobel, C. Schmidt-Hattenberger, Energ. Environ. Sci. 14 (2021) 853.

[2] P. Schulze, J. Holstein, A. van den Noort, J. Knijp, DNV GL, 2017.

[3] S.M.J. Raad, Y. Leonenko, H. Hassanzadeh, Renew. Sustain. Energy Rev. 168 (2022) 112846.

[4] M. Panfilov, Compendium of hydrogen energy, Elsevier (2016) 91–115.

[5] S. Chabab, P. Theveneau, C. Coquelet, J. Corvisier, P. Paricaud, Int. J. Hydrogen Energy 45 (2020) 32206.

[6] N.S. Muhammed, B. Haq, D. Al Shehri, A. Al-Ahmed, M.M. Rahman, E. Zaman, Energy Rep. 8 (2022) 461.

[7] V.K. Singh, (2010).

[8] GaffneyCline, Underground Hydrogen Storage (UHS). https://www.gaffneycline.com/ (2022).

[9] V. Zgonnik, Earth Sci. Rev. 203 (2020) 103140.

[10] E.C. Gaucher, Elements: an International Magazine of Mineralogy, Geochemistry, and Petrology 16 (2020) 8.

[11] E. Gaucher, I. Moretti, N. Péllissier, G. Burridge, N. Gonthier, Eur. Geol. (2023) 5.

[12] G.A. Torín-Ollarves, J.M. Trusler, Fluid Phase Equilib. 539 (2021) 113025.

[13] C. Lopez-Lazaro, P. Bachaud, I. Moretti, N. Ferrando, BSGF-Earth Sciences Bulletin 190 (2019) 7.

[14] W. van Rooijen, P. Habibi, K. Xu, P. Dey, T. Vlugt, H. Hajibeygi, O. Moultos, J. Chem. Eng. Data (2023).

[15] S. Chabab, H. Kerkache, I. Bouchkira, M. Poulain, O. Baudouin, É. Moine, M. Ducousso, H. Hoang, G. Galliéro, P. Cézac, Int. J. Hydrogen Energy 50 (2024) 648.

[16] P. Ungerer, B. Tavitian, A. Boutin, Applications of molecular simulation in the oil and gas industry: Monte Carlo methods, Editions Technip, 2005.

[17] B. Creton, C. Nieto-Draghi, T. De Bruin, V. Lachet, E. El Ahmar, A. Valtz, C. Coquelet, S. Lasala, R. Privat, J.-N. Jaubert, Fluid Phase Equilib. 461 (2018) 84.

[18] H. Jiang, I.G. Economou, A.Z. Panagiotopoulos, Acc. Chem. Res. 50 (2017) 751.

[19] E.S. Tsai, H. Jiang, A.Z. Panagiotopoulos, Fluid Phase Equilib. 407 (2016) 262.

[20] J. Vorholz, V. Harismiadiis, A.Z. Panagiotopoulos, B. Rumpf, G. Maurer, Fluid Phase Equilib. 226 (2004) 237.

[21] Y. Liu, T. Lafitte, A.Z. Panagiotopoulos, P.G. Debenedetti, AIChE J 59 (2013) 3514.

[22] A. Rahbari, R. Hens, I. Nikolaidis, A. Poursaeidesfahani, M. Ramdin, I. Economou, O. Moultos, D. Dubbeldam, T. Vlugt, Mol. Phys. 116 (2018) 3331.

[23] A. Rahbari, J. Brenkman, R. Hens, M. Ramdin, L.J. Van Den Broeke, R. Schoon, R. Henkes, O.A. Moultos, T.J. Vlugt, J. Chem. Eng. Data 64 (2019) 4103.

[24] W. Shi, E.J. Maginn, J. Chem. Theory Comput. 3 (2007) 1451.

[25] W. Shi, E.J. Maginn, J. Comput. Chem. 29 (2008) 2520.

[26] I. Søreide, C.H. Whitson, Fluid Phase Equilib. 77 (1992) 217.

[27] R. Sander, W.E. Acree Jr, A. De Visscher, S.E. Schwartz, T.J. Wallington, Pure Appl. Chem. 94 (2022) 71.

[28] C. Yong, Pergamon, 1981.

[29] I. Krichevsky, J. Kasarnovsky, J. Am. Chem. Soc. 57 (1935) 2168.

[30] A.H. Harvey, AIChE J 42 (1996) 1491.

[31] R. Fernández-Prini, J.L. Alvarez, A.H. Harvey, J. Phys. Chem. Ref. Data 32 (2003) 903.

[32] T.-K.-H. Trinh, J.-C. de Hemptinne, R. Lugo, N. Ferrando, J.-P. Passarello, J. Chem. Eng. Data 61 (2016) 19.

[33] N.N. Akinfiev, L.W. Diamond, Geochim. Cosmochim. Acta 67 (2003) 613.

[34] D. Li, C. Beyer, S. Bauer, Int. J. Hydrogen Energy 43 (2018) 512.

[35] E.F. Bazarkina, I.-M. Chou, A.F. Goncharov, N.N. Akinfiev, Elements: an International Magazine of Mineralogy, Geochemistry, and Petrology 16 (2020) 33.

[36] R. Wiebe, Ind. Eng. Chem. 24 (1932) 823.

[37] R. Wiebe, V. Gaddy, J. Am. Chem. Soc. 56 (1934) 76.

[38] T. Morrison, F. Billett, Journal of the Chemical Society (resumed) (1952) 3819.

[39] T.E. Crozier, S. Yamamoto, J. Chem. Eng. Data 19 (1974) 242.

[40] L.I. Gordon, Y. Cohen, D.R. Standley, Deep-Sea Res. 24 (1977) 937.

[41] P. Gillespie, G. Wilson, Tulsa, OK (1980).

[42] J. Alvarez, R. Crovetto, R. Fernández-Prini, Ber. Bunsen. Phys. Chem 92 (1988) 935.

[43] G. Kling, G. Maurer, J. Chem. Thermodyn. 23 (1991) 531.

[44] S. Chabab, J.L. Cruz, M. Poulain, M. Ducousso, F. Contamine, J.P. Serin, P. Cézac, Energies 14 (2021) 5239.

[45] P. Steiner, Ann. Phys. 288 (1894) 275.

[46] L. Braun, Z. Phys. Chem. 33 (1900) 721.

[47] J. Gerecke, H. Bittrich, Wiss Z Tech Hochsch Chem Carl Shorlemmer Leuna Merseburg 13 (1971) 115.

[48] A.Z. Panagiotopoulos, Mol. Phys. 61 (1987) 813.

[49] B. Smit, D. Frenkel, Mol. Phys. 68 (1989) 951.

[50] A. Rahbari, A. Poursaeidesfahani, A. Torres-Knoop, D. Dubbeldam, T.J. Vlugt, Mol. Simul. 44 (2018) 405.

[51] S. Murad, S. Gupta, Chem. Phys. Lett. 319 (2000) 60.

[52] T. Schnabel, J. Vrabec, H. Hasse, Fluid Phase Equilib. 233 (2005) 134.

[53] K. Lekvam, P.R. Bishnoi, Fluid Phase Equilib. 131 (1997) 297.
brine mixtures. This structure is predominantly influenced by H₂-H₂O and H₂-Cl⁻ interactions. Consequently, a constant binary interaction parameter on energy for H₂-H₂O in the Marx-TIP4P/2005 and Marx-TIP4P/EP models, and on volume for H₂-Cl⁻ in the Marx-TIP4P/EP-KBF, was implemented. Following this adjustment, both model combinations provided excellent representations of H₂ solubility in NaCl brine.

Despite both models' accurate estimations of hydrogen solubility in brine, the predictive model Marx-TIP4P/2005-Madrid was chosen to generate data, as it was not directly adjusted to the ternary system. This model underwent validation against experimental literature data for H₂ solubility in brine (for pressures up to 200 bar, temperatures between 298 K and 373 K, and salinities up to 4 mol/kg NaCl), extending the data for H₂/brine to higher temperatures (up to 453 K), pressures (up to 1000 bar), and salinities (up to 6 mol/kg NaCl).

Both experimental and molecular simulations data were leveraged to refine parameters of a thermodynamic model. The asymmetric (gamma-phi) approach of the Krichevsky-Kasarnovsky type, chosen for its non-iterative simplicity in single-gas systems, facilitated the calculation of the liquid-vapor equilibrium. The model showcased very good results for H₂ solubility in water and brine across varying temperatures (298–453 K), pressures (up to 1000 bar), and NaCl molalities (up to 6 mol NaCl/kg water), proving invaluable for evaluating H₂ storage capacity, dissolution losses, biotic and abiotic reactivity, and caprock integrity.

Future work may delve into the interactions of H₂ with other sub-surface gases (CO₂, CH₄, etc.) or the impact of clay on its solubility, building on the current findings.

### 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, Writing – original draft, Validation, Software, Methodology, Investigation, Formal analysis. **Pierre Cézac:** Writing – review & editing, Supervision, Resources, Project administration, Funding acquisition. **Guillaume Galliéro:** Writing – review & editing, Writing – original draft, Supervision, Resources, Project administration, Methodology, Funding acquisition, Formal analysis, Conceptualization. **Salaheddine Chabab:** Writing – review & editing, Writing – original draft, Supervision, Project administration, Methodology, Investigation, 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.

### 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 grant 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 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.mollicq.2024.124497.

### References

[1] N. Heinemann, J. Alcalde, J.M. Miocic, S.J. Hangx, J. Kallmeyer, C. Ostertag-Henning, A. Hassanpouryouzband, E.M. Thaysen, G.J. Strobel, C. Schmidt-Hattenberger, Energ. Environ. Sci. 14 (2021) 853.

[2] P. Schulze, J. Holstein, A. van den Noort, J. Knijp, DNV GL, 2017.

[3] S.M.J. Raad, Y. Leonenko, H. Hassanzadeh, Renew. Sustain. Energy Rev. 168 (2022) 112846.

[4] M. Panfilov, Compendium of hydrogen energy, Elsevier (2016) 91–115.

[5] S. Chabab, P. Theveneau, C. Coquelet, J. Corvisier, P. Paricaud, Int. J. Hydrogen Energy 45 (2020) 32206.

[6] N.S. Muhammed, B. Haq, D. Al Shehri, A. Al-Ahmed, M.M. Rahman, E. Zaman, Energy Rep. 8 (2022) 461.

[7] V.K. Singh, (2010).

[8] GaffneyCline, Underground Hydrogen Storage (UHS). https://www.gaffneycline.com/ (2022).

[9] V. Zgonnik, Earth Sci. Rev. 203 (2020) 103140.

[10] E.C. Gaucher, Elements: an International Magazine of Mineralogy, Geochemistry, and Petrology 16 (2020) 8.

[11] E. Gaucher, I. Moretti, N. Pélissier, G. Burridge, N. Gonthier, Eur. Geol. (2023) 5.

[12] G.A. Torín-Ollarves, J.M. Trusler, Fluid Phase Equilib. 539 (2021) 113025.

[13] C. Lopez-Lazaro, P. Bachaud, I. Moretti, N. Ferrando, BSGF-Earth Sciences Bulletin 190 (2019) 7.

[14] W. van Rooijen, P. Habibi, K. Xu, P. Dey, T. Vlugt, H. Hajibeygi, O. Moultos, J. Chem. Eng. Data (2023).

[15] S. Chabab, H. Kerkache, I. Bouchkira, M. Poulain, O. Baudouin, É. Moine, M. Ducousso, H. Hoang, G. Galliéro, P. Cézac, Int. J. Hydrogen Energy 50 (2024) 648.

[16] P. Ungerer, B. Tavitian, A. Boutin, Applications of molecular simulation in the oil and gas industry: Monte Carlo methods, Editions Technip, 2005.

[17] B. Creton, C. Nieto-Draghi, T. De Bruin, V. Lachet, E. El Ahmar, A. Valtz, C. Coquelet, S. Lasala, R. Privat, J.-N. Jaubert, Fluid Phase Equilib. 461 (2018) 84.

[18] H. Jiang, I.G. Economou, A.Z. Panagiotopoulos, Acc. Chem. Res. 50 (2017) 751.

[19] E.S. Tsai, H. Jiang, A.Z. Panagiotopoulos, Fluid Phase Equilib. 407 (2016) 262.

[20] J. Vorholz, V. Harismiadiis, A.Z. Panagiotopoulos, B. Rumpf, G. Maurer, Fluid Phase Equilib. 226 (2004) 237.

[21] Y. Liu, T. Lafitte, A.Z. Panagiotopoulos, P.G. Debenedetti, AIChE J 59 (2013) 3514.

[22] A. Rahbari, R. Hens, I. Nikolaidis, A. Poursaeidesfahani, M. Ramdin, I. Economou, O. Moultos, D. Dubbeldam, T. Vlugt, Mol. Phys. 116 (2018) 3331.

[23] A. Rahbari, J. Brenkman, R. Hens, M. Ramdin, L.J. Van Den Broeke, R. Schoon, R. Henkes, O.A. Moultos, T.J. Vlugt, J. Chem. Eng. Data 64 (2019) 4103.

[24] W. Shi, E.J. Maginn, J. Chem. Theory Comput. 3 (2007) 1451.

[25] W. Shi, E.J. Maginn, J. Comput. Chem. 29 (2008) 2520.

[26] I. Søreide, C.H. Whitson, Fluid Phase Equilib. 77 (1992) 217.

[27] R. Sander, W.E. Acree Jr, A. De Visscher, S.E. Schwartz, T.J. Wallington, Pure Appl. Chem. 94 (2022) 71.

[28] C. Yong, Pergamon, 1981.

[29] I. Krichevsky, J. Kasarnovsky, J. Am. Chem. Soc. 57 (1935) 2168.

[30] A.H. Harvey, AIChE J 42 (1996) 1491.

[31] R. Fernández-Prini, J.L. Alvarez, A.H. Harvey, J. Phys. Chem. Ref. Data 32 (2003) 903.

[32] T.-K.-H. Trinh, J.-C. de Hemptinne, R. Lugo, N. Ferrando, J.-P. Passarello, J. Chem. Eng. Data 61 (2016) 19.

[33] N.N. Akinfiev, L.W. Diamond, Geochim. Cosmochim. Acta 67 (2003) 613.

[34] D. Li, C. Beyer, S. Bauer, Int. J. Hydrogen Energy 43 (2018) 512.

[35] E.F. Bazarkina, I.-M. Chou, A.F. Goncharov, N.N. Akinfiev, Elements: an International Magazine of Mineralogy, Geochemistry, and Petrology 16 (2020) 33.

[36] R. Wiebe, Ind. Eng. Chem. 24 (1932) 823.

[37] R. Wiebe, V. Gaddy, J. Am. Chem. Soc. 56 (1934) 76.

[38] T. Morrison, F. Billett, Journal of the Chemical Society (resumed) (1952) 3819.

[39] T.E. Crozier, S. Yamamoto, J. Chem. Eng. Data 19 (1974) 242.

[40] L.I. Gordon, Y. Cohen, D.R. Standley, Deep-Sea Res. 24 (1977) 937.

[41] P. Gillespie, G. Wilson, Tulsa, OK (1980).

[42] J. Alvarez, R. Crovetto, R. Fernández-Prini, Ber. Bunsen. Phys. Chem 92 (1988) 935.

[43] G. Kling, G. Maurer, J. Chem. Thermodyn. 23 (1991) 531.

[44] S. Chabab, J.L. Cruz, M. Poulain, M. Ducousso, F. Contamine, J.P. Serin, P. Cézac, Energies 14 (2021) 5239.

[45] P. Steiner, Ann. Phys. 288 (1894) 275.

[46] L. Braun, Z. Phys. Chem. 33 (1900) 721.

[47] J. Gerecke, H. Bittrich, Wiss Z Tech Hochsch Chem Carl Shorlemmer Leuna Merseburg 13 (1971) 115.

[48] A.Z. Panagiotopoulos, Mol. Phys. 61 (1987) 813.

[49] B. Smit, D. Frenkel, Mol. Phys. 68 (1989) 951.

[50] A. Rahbari, A. Poursaeidesfahani, A. Torres-Knoop, D. Dubbeldam, T.J. Vlugt, Mol. Simul. 44 (2018) 405.

[51] S. Murad, S. Gupta, Chem. Phys. Lett. 319 (2000) 60.

[52] T. Schnabel, J. Vrabec, H. Hasse, Fluid Phase Equilib. 233 (2005) 134.

[53] K. Lekvam, P.R. Bishnoi, Fluid Phase Equilib. 131 (1997) 297.
[54] R. Hens, A. Rahbari, S.n. Caro-Ortíz, N. Dawass, M.t. Erdős, A. Poursaeidesfahani, H.S. Salehi, A.T. Celebi, M. Ramdin, O.A. Moultos, J. Chem. Inform. Model. 60 (2020) 2678.

[55] M.E. Tuckerman, Statistical mechanics: theory and molecular simulation, Oxford University Press, 2023.

[56] H. Nagashima, S. Tsuda, N. Tsuboi, M. Koshi, K. Hayashi, T. Tokumasu, J. Chem. Phys. 140 (2014).

[57] V. Buch, J. Chem. Phys. 100 (1994) 7610.

[58] J.O. Hirschfelder, C.F. Curtiss, R.B. Bird, The molecular theory of gases and liquids, John Wiley & Sons, 1964.

[59] A. Köster, M. Thol, J. Vrabec, J. Chem. Eng. Data 63 (2018) 305.

[60] S. Mondal, S. Ghosh, P. Chattaraj, J. Mol. Model. 19 (2013) 2785.

[61] R.F. Cracknell, PCCP 3 (2001) 2091.

[62] Q. Yang, C. Zhong, J. Phys. Chem. B 109 (2005) 11862.

[63] S. Alavi, J. Ripmeester, D. Klug, J. Chem. Phys. 123 (2005).

[64] D. Marx, P. Nielaba, Phys. Rev. A 45 (1992) 8968.

[65] J.L. Belof, A.C. Stern, B. Space, J. Chem. Theory Comput. 4 (2008) 1332.

[66] C. Vega, J.L. Abascal, M. Conde, J. Aragones, Faraday Discuss. 141 (2009) 251.

[67] H.J. Berendsen, J.R. Grigera, T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269.

[68] S.P. Kadaoluwa Pathirannahalage, N. Meftahi, A. Elbourne, A.C. Weiss, C. F. McConville, A. Padua, D.A. Winkler, M. Costa Gomes, T.L. Greaves, T.C. le, J. Chem. Inf. Model. 61 (2021) 4521.

[69] R.l. Fuentes-Azcatl, J. Alejandre, J. Phys. Chem. B 118 (2014) 1263.

[70] J. Chandrasekhar, D.C. Spellmeyer, W.L. Jorgensen, J. Am. Chem. Soc. 106 (1984) 903.

[71] A. Benavides, M. Portillo, V. Chamorro, J. Espinosa, J. Abascal, C. Vega, J. Chem. Phys. 147 (2017).

[72] I. Zeron, J. Abascal, C. Vega, J. Chem. Phys. 151 (2019).

[73] P. Loche, P. Steinbrunner, S. Friedowitz, R.R. Netz, D.J. Bonthuis, J. Phys. Chem. B 125 (2021) 8581.

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

[75] H.C. Andersen, J. Comput. Phys. 52 (1983) 24.

[76] B.A. Luty, W.F. van Gunsteren, J. Phys. Chem. 100 (1996) 2581.

[77] J.W. Leachman, R.T. Jacobsen, S. Penoncello, E.W. Lemmon, J. Phys. Chem. Ref. Data 38 (2009) 721.

[78] E.W. Lemmon, M.L. Huber, M.O. McLinden, Reference fluid thermodynamic and transport properties (REFPROP), Version 9 (2010).

[79] R. Sun, S. Lai, J. Dubessy, Fluid Phase Equilib. 390 (2015) 23.

[80] W. Wagner, A. Pruß, J. Phys. Chem. Ref. Data 31 (2002) 387.

[81] V. Vins, D. Celný, B. Planková, T. Němec, M. Duska, J. Hrubý, EPJ Web of Conferences, EDP Sciences (2016) 02136.

[82] J. Delhommelle, P. Millié, Mol. Phys. 99 (2001) 619.

[83] H. Docherty, A. Galindo, C. Vega, E. Sanz, J. Chem. Phys. 125 (2006).

[84] P.J. Dyer, H. Docherty, P.T. Cummings, J. Chem. Phys. 129 (2008).

[85] S. Al Ghafri, G.C. Maitland, J.M. Trusler, J. Chem. Eng. Data 57 (2012) 1288.

[86] J.L. Aragones, M. Rovere, C. Vega, P. Gallo, J. Phys. Chem. B 118 (2014) 7680.

[87] J. Aragones, E. Sanz, C. Vega, J. Chem. Phys. 136 (2012).

[88] W.A. Gerth, J. Solution Chem. 12 (1983) 655.

[89] D.-Y. Peng, D.B. Robinson, Ind. Eng. Chem. Fundam. 15 (1976) 59.

[90] N. Sakoda, K. Shindo, K. Motomura, K. Shinzato, M. Kohno, Y. Takata, M. Fujii, Int. J. Thermophys. 33 (2012) 381.

[91] D. Koschel, J.-Y. Coxam, L. Rodier, V. Majer, Fluid Phase Equilib. 247 (2006) 107.

[92] T.D. O'Sullivan, N.O. Smith, J. Phys. Chem. 74 (1970) 1460.
[54] R. Hens, A. Rahbari, S.n. Caro-Ortíz, N. Dawass, M.t. Erdős, A. Poursaeidesfahani, H.S. Salehi, A.T. Celebi, M. Ramdin, O.A. Moultos, J. Chem. Inform. Model. 60 (2020) 2678.

[55] M.E. Tuckerman, Statistical mechanics: theory and molecular simulation, Oxford University Press, 2023.

[56] H. Nagashima, S. Tsuda, N. Tsuboi, M. Koshi, K. Hayashi, T. Tokumasu, J. Chem. Phys. 140 (2014).

[57] V. Buch, J. Chem. Phys. 100 (1994) 7610.

[58] J.O. Hirschfelder, C.F. Curtiss, R.B. Bird, The molecular theory of gases and liquids, John Wiley & Sons, 1964.

[59] A. Köster, M. Thol, J. Vrabec, J. Chem. Eng. Data 63 (2018) 305.

[60] S. Mondal, S. Ghosh, P. Chattaraj, J. Mol. Model. 19 (2013) 2785.

[61] R.F. Cracknell, PCCP 3 (2001) 2091.

[62] Q. Yang, C. Zhong, J. Phys. Chem. B 109 (2005) 11862.

[63] S. Alavi, J. Ripmeester, D. Klug, J. Chem. Phys. 123 (2005).

[64] D. Marx, P. Nielaba, Phys. Rev. A 45 (1992) 8968.

[65] J.L. Belof, A.C. Stern, B. Space, J. Chem. Theory Comput. 4 (2008) 1332.

[66] C. Vega, J.L. Abascal, M. Conde, J. Aragones, Faraday Discuss. 141 (2009) 251.

[67] H.J. Berendsen, J.R. Grigera, T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269.

[68] S.P. Kadaoluwa Pathirannahalage, N. Meftahi, A. Elbourne, A.C. Weiss, C. F. McConville, A. Padua, D.A. Winkler, M. Costa Gomes, T.L. Greaves, T.C. le, J. Chem. Inf. Model. 61 (2021) 4521.

[69] R.l. Fuentes-Azcatl, J. Alejandre, J. Phys. Chem. B 118 (2014) 1263.

[70] J. Chandrasekhar, D.C. Spellmeyer, W.L. Jorgensen, J. Am. Chem. Soc. 106 (1984) 903.

[71] A. Benavides, M. Portillo, V. Chamorro, J. Espinosa, J. Abascal, C. Vega, J. Chem. Phys. 147 (2017).

[72] I. Zeron, J. Abascal, C. Vega, J. Chem. Phys. 151 (2019).

[73] P. Loche, P. Steinbrunner, S. Friedowitz, R.R. Netz, D.J. Bonthuis, J. Phys. Chem. B 125 (2021) 8581.

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

[75] H.C. Andersen, J. Comput. Phys. 52 (1983) 24.

[76] B.A. Luty, W.F. van Gunsteren, J. Phys. Chem. 100 (1996) 2581.

[77] J.W. Leachman, R.T. Jacobsen, S. Penoncello, E.W. Lemmon, J. Phys. Chem. Ref. Data 38 (2009) 721.

[78] E.W. Lemmon, M.L. Huber, M.O. McLinden, Reference fluid thermodynamic and transport properties (REFPROP), Version 9 (2010).

[79] R. Sun, S. Lai, J. Dubessy, Fluid Phase Equilib. 390 (2015) 23.

[80] W. Wagner, A. Pruß, J. Phys. Chem. Ref. Data 31 (2002) 387.

[81] V. Vins, D. Celný, B. Planková, T. Němec, M. Duska, J. Hrubý, EPJ Web of Conferences, EDP Sciences (2016) 02136.

[82] J. Delhommelle, P. Millié, Mol. Phys. 99 (2001) 619.

[83] H. Docherty, A. Galindo, C. Vega, E. Sanz, J. Chem. Phys. 125 (2006).

[84] P.J. Dyer, H. Docherty, P.T. Cummings, J. Chem. Phys. 129 (2008).

[85] S. Al Ghafri, G.C. Maitland, J.M. Trusler, J. Chem. Eng. Data 57 (2012) 1288.

[86] J.L. Aragones, M. Rovere, C. Vega, P. Gallo, J. Phys. Chem. B 118 (2014) 7680.

[87] J. Aragones, E. Sanz, C. Vega, J. Chem. Phys. 136 (2012).

[88] W.A. Gerth, J. Solution Chem. 12 (1983) 655.

[89] D.-Y. Peng, D.B. Robinson, Ind. Eng. Chem. Fundam. 15 (1976) 59.

[90] N. Sakoda, K. Shindo, K. Motomura, K. Shinzato, M. Kohno, Y. Takata, M. Fujii, Int. J. Thermophys. 33 (2012) 381.

[91] D. Koschel, J.-Y. Coxam, L. Rodier, V. Majer, Fluid Phase Equilib. 247 (2006) 107.

[92] T.D. O'Sullivan, N.O. Smith, J. Phys. Chem. 74 (1970) 1460.