Predicting the phase behavior of hydrogen in NaCl brines by molecular simulation for geological applications

Hydrogen (H₂) is currently taking an increasing role in the energy mix. If hydrogen remains today a main raw material for the industry, its capacity for long-term storage of decarbonized energy is highlighted by numerous authors and the H₂ council. First pilot projects are taking place and use electrolyzers to convert solar- and wind-farm electricity to H₂ which is then stored (Darras et al., 2012; Perez et al., 2016). It can then be converted again to electricity (such as the Methycentre project [ADEME, 2018]), injected in the natural gas network (such as the Grhyd project in Dunkerque, France [ENGIE, 2018]) or used as fuel for vehicles. But the development of a large-scale use of hydrogen implies increasing storage needs, and injection in geological formations, used extensively in the gas and compressed air energy industries, is an attractive solution to reach sizable storage volumes (Lord et al., 2014). Various options that already proved their efficiency and economic viability for other gases, such as injection in aquifers, depleted oil and gas fields, and salt caverns, are currently under study (Le Duigou et al., 2017).

In addition, numerous native H₂ emissions have been detected worldwide, in the mid oceanic ridge, in active fault zones, in ophiolitic contexts or in intracratonic basins (Charlou et al., 2002; Charlou et al., 2010; Deville and Prinzhofer 2016; Larin et al., 2015; Moretti et al., 2018; Prinzhofer and Deville, 2015; Satake et al., 1984; Sato et al., 1986; Vacquand et al., 2018). Accumulations have even been found in Kansas (Guélard et al., 2017) or Mali (Prinzhofer et al., 2018), where natural hydrogen is used to produce electricity since 2015. The sources of H₂ and the mechanisms responsible for its formation remain to be clearly identified, even if some chemical reactions are already known to be efficient H₂ producers, such as serpentinization or oxidation of iron-rich rocks (Bachaud et al., 2017; Deville and Prinzhofer 2016; Vacquand et al., 2018). The behavior of this molecule along its migration path is another aspect that needs to be better constrained. The high mobility and reactivity of hydrogen have been highlighted in numerous studies (Carden and Paterson, 1979; Paterson, 1983; Reitenbach et al., 2015), both of which being impacted by hydrogen solubility in formation waters. Understanding the phase behavior of systems involving hydrogen and brines is thus of primary importance for exploration of native H₂ accumulations and to correctly design hydrogen storage in geological formations.

Various equations of state have been already developed to accurately predict solubility of hydrogen in pure water at high temperature and pressure for geological applications (e.g.Akinfiev and Diamond, 2003; Plyasunov and Bazarkina, 2018; Plyasunov and Shock, 2001). Concerning electrolytic aqueous solutions, light gas solubility has been successfully modeled in the past using either activity coefficient models mainly based on the Pitzer theory (e.g.Duan and Sun, 2003; Duan et al., 1992; Li et al., 2018; Rumpf et al., 1998; Xia et al., 2000; Xia et al., 1999) or equations of state such as Virial-based equations (e.g.Duan et al., 1996, 2003), cubic equations (e.g. (Li et al., 2015; Soreide and Whitson, 1992; Sorensen et al., 2002), or more advanced equations of state including association forces as e-CPA (e.g.Carvalho et al., 2015; Courtial et al., 2014; Hemptinne et al., 2006; Kontogeorgis and Folas, 2010; Maribo-Mogensen et al., 2015) or e-SAFT (e.g.Ahmed et al., 2018; Held et al., 2014; Ji et al., 2005; Llovell et al., 2012; Patel et al., 2003; Rozmus et al., 2013; Sun and Dubessy, 2012). Nevertheless, all of these models involve empirical parameters to be fitted on available experimental data. Their application range is thus restricted to the temperature, pressure and salinity ranges of the experimental data used for their parameterization.

Concerning hydrogen solubility in brines, an important issue is the lack of experimental data in thermodynamic conditions relevant for geological applications. As described in the first section of this article, almost all the available solubility data concern temperature less than 303 K, and salinities less than 1 mol/kg_H2O. The acquisition of new experimental data is a challenging and costly task for safety reasons (high flammability of hydrogen, corrosive aspect of brines) and technical aspects (low H₂ solubility). In this context, a promising alternative is molecular simulation. Thanks to robust simulations, Monte Carlo simulation technique has become a very powerful method to predict thermodynamic data and phase equilibrium of systems of industrial interest (Theodorou, 2010; Ungerer et al., 2005). When using an appropriate force field, molecular simulation can be seen as a reliable tool to generate data in complement of classical experiments.

Previous studies have demonstrated the ability of Monte Carlo simulation to predict gas solubility in brines, such as CO₂ (Creton et al., 2018; Jiang et al., 2017; Liu et al., 2013; Tsai et al., 2016; Vorholz et al., 2004), H₂S (Fauve et al., 2017), SO₂ and other diatomic light gases (Creton et al., 2018). However, to the best of our knowledge, no molecular simulation study dealing with the solubility of hydrogen in electrolyte solutions has been published. Hence, the goal of this work is first to validate a force field able to accurately predict hydrogen solubility in aqueous NaCl solutions over a wide range of temperatures and salt concentrations, and second to use the simulation data in complement of the existing experimental data to recalibrate an equation of state.

This paper is organized as follows: in a first section, a literature review is carried out to build a consistent database for H₂ + H₂O and H₂ + H₂O + NaCl systems. In a second section, a force field able to reproduce experimental data for both salt-free and salted systems is proposed. Then, Monte Carlo simulations are carried out at higher temperatures and salinities to generate new simulated data. In a third section, experimental and simulated solubility data are used to recalibrate the Soreide and Whitson equation of state (Soreide and Whitson, 1992) applicable over wide ranges of temperature, pressure and salinity. Finally, some geological applications for this new model are provided in the context of hydrogen storage and naturally-produced subsurface hydrogen. Today, underground gas storage as well as natural gas production refer mainly to hydrocarbon. As a comparison, we conclude on the difference between the methane and the H₂ solubility and so transport and accumulation modes in subsurface.

Mole fractions are then converted in Henry constant using equation (2). The fugacity coefficient of hydrogen in vapor phase is calculated with the Soave-Redlich-Kwong equation of state (Soave, 1972) after having checked that this model is able to accurately reproduce pure hydrogen reference compressibility factors (Younglove, 1982). Solubility data (“TPxy” data) are also converted in Henry constants with equation (2) by taking for a given isotherm the smallest available experimental partial pressure to stay as close as possible of the infinite dilution domain. Finally, Table 1 details the type and the references of the various experimental data selected in this work for the binary system H₂ + H₂O. From these experimental data, the evolution of the Henry constant of H₂ in pure water is shown on Figure 1. The repeatability of the experiments coming from various authors allows to evaluate an average experimental uncertainty of ±10%.

Experimental data exhibit a maximum in the Henry constant curve for a temperature close to 330 K (i.e. about 57 °C), indicating that hydrogen solubility is minimal at this temperature. This is the typical behavior for the solubility of a small molecule in pure water, the minimum-solubility temperature depending on the nature of the solute (Ahmed et al., 2016).

Some experimental data for hydrogen solubility in aqueous NaCl solutions are available in literature, but they concern very limited temperature and salinity ranges. All the data available are provided under Bunsen coefficient form, and have been converted into Henry constant as previously described. Note that in this conversion, the density of the solvent is required, and to this end the model of Rowe and Chou (Rowe, 1970) was used to calculate density of NaCl aqueous solution at the desired temperature and salt concentration. The most significant contribution is the work of Crozier and Yamamoto (1974) who acquired a large amount of hydrogen solubility data in seawater and NaCl aqueous solutions. Among the other data available, those of Braun (1900) have been rejected because they were found inconsistent with all other available data in similar conditions. Finally, a total of 296 consistent data are considered, and references are detailed in Table 1.

It is worth noticing that 95% of these data concern salt concentrations less than 1 mol/kg_H2O (that is, 58 g/L, twice the usual salinity of seawater) and temperatures less than 303 K. As these operating conditions are very far from temperatures and even salinities that can be encountered in geological formations, it fully justifies the need to generate new data at higher temperature and salinity. Figure 2 shows a selection of Henry constants as a function of temperature for various salinity values. As data are mainly restricted to low temperatures, they cannot be correlated as previously done for pure water. Nevertheless some typical behaviors can be emphasized. First, the greater the molality, the greater the value of Henry constant. This is a consequence of the “salting-out” effect which corresponds to the solubility decrease when salt concentration in water increases. Second, as observed for pure water, data seem to exhibit a maximum in the curve comprised between 320 and 350 K.

The molecular model used for hydrogen is the model of Darkrim et al. (1999). This choice is motivated by a previous study led on hydrogen solubility in organic solvents (Ferrando and Ungerer, 2007) which showed that this molecular model is able to accurately reproduce pure hydrogen compressibility factor over a large temperature range and is particularly suitable for phase equilibrium calculation. The hydrogen molecule is represented by a single Lennard-Jones spheres with three point charges aim to mimic the quadrupole moment of the molecule. Details of this model are given in Table 3.

For aqueous NaCl solutions, the combination of the TIP4P/2005 model for pure water (Abascal and Vega, 2005) and the OPLS model for Na⁺ and Cl⁻ (Chandrasekhar et al., 1984) is selected. The TIP4P/2005 model is able to accurately reproduce pure water density and vapor pressures, and it has been shown in a recent study (Creton et al., 2018) that its combination with the OPLS model for ions allows to reproduce correctly densities and osmotic pressures of aqueous NaCl solutions over a large temperature and salinity range. Details of these molecular models are given in Table 3.

Aqueous NaCl solutions form a well-structured liquid phase in which the insertion move can be often rejected. Thus, Widom insertion method could appear in this context less efficient than more advanced techniques to evaluate chemical potential, such as umbrella sampling (Torrie and Valleau, 1977), slow-growth method (Nezbeda and Kolafa, 1991) or thermodynamic integration (Frenkel and Smit, 1996). Nevertheless, a prior study on the solubility of H₂S in NaCl aqueous solution (Fauve et al., 2017) shown that, with a sufficient number of Monte Carlo steps, the Widom insertion method lead to similar results than thermodynamic integration technique within statistical uncertainties.

The Widom insertion method is applied during a Monte Carlo simulation in the NPT ensemble. This ensemble consists in one single box representing the liquid phase, where the total number of molecule N, the temperature T and the pressure P are imposed. For a given temperature, simulations are carried out at a pressure equal to the vapor pressure of the pure solvent (computed with the current force field). A total number of 500 molecules of water is used, and the number of Na⁺ and Cl⁻ ions is chosen according to the desired molality (for example 9 particles of Na⁺ and Cl⁻ for a molality of 1 mol/kg_H2O). A simulation consists typically in an equilibrium run of about 200 to 500 million Monte Carlo steps to achieve equilibrium state, followed by a production run lasting for additional 400 million Monte Carlo steps to compute average properties. Each step generates a new configuration, which is accepted or rejected following the classical Metropolis criterion (Frenkel and Smit, 1996). The excess chemical potential value required in equation (8) is computed each 50 000 steps. The different Monte Carlo moves and their corresponding attempt probabilities used during a simulation are molecular translation (20%), molecular rotation (20%), insertion move (59.5%) with a pre-insertion statistical bias (Mackie et al., 1997) using 10 trial positions, and volume change (0.5%). The amplitude of translations, rigid rotations and volume changes was adjusted during the simulation to achieve an acceptance ratio of 40% for these moves. The simulations are carried out using the GIBBS software jointly developed by IFP Energies nouvelles and the Laboratoire de Chimie-Physique (CNRS-Université Paris-Sud) (Ungerer et al., 2005).

Monte Carlo simulation of the system H₂ + H₂O have been carried out with the three combining rules of Table 4 for the water – hydrogen cross interactions. From the results plotted on Figure 3, it can be observed that the choice of the combining rule has a negligible effect at high temperatures (above 400 K). At lower temperatures, the effect is more pronounced, although results remain close considering statistical uncertainties. The best compromise is obtained when using the Lorentz-Berthelot combining rule. At high temperatures (typically from 380 K, 107 °C), a very good accordance between experimental and simulation results is achieved. At lower temperatures, simulated Henry constants are slightly overestimated, but offer a correct agreement within statistical and experimental uncertainties. Finally, an average absolute deviation of 10% is obtained on the whole range of temperatures using these combining rules. It is interesting to notice that the maximum in the Henry constant curve is predicted at the correct temperature (330 K) without introducing specific binary parameter adjusted on this value. Consequently, the Lorentz-Berthelot combining rules are definitively adopted for this study. We emphasis that this choice is motivated only on hydrogen solubility considerations, which is the focus of this work.

Henry constants of hydrogen in aqueous NaCl solutions have been first computed from a fully predictive manner for NaCl molalities of 1 and 2 mol/kg_H2O and temperatures up to 333 K, and compared to the available experimental data in same conditions, as illustrated on Figure 4.

The effect of the ξ_ij parameter is first investigated by setting the l_ij parameter to 0 and by tuning the ξ_ij parameter to match the experimental data at 1 and 2 mol/kg_H2O NaCl in the same temperature range. However, this optimization was not successful since it was observed that this binary parameter has only a negligible effect on the simulation results. It suggests that the attractive Van der Waals forces do not play an important role for the interactions between hydrogen and ionic species. Thus, rather, it was decided to set the ξ_ij value to 0 and to optimize the l_ij value, which influences the repulsive part of the Lennard-Jones potential. In the case of the H₂O + H₂ + NaCl system, simulation results are very sensitive to this parameter. The optimal l_ij value, which led to the best agreement between experimental data and simulations, is provided in Table 5. This observation suggests that hydrogen solubility is rather driven by repulsion forces and consequently by volume effects. This assumption is corroborated by the recent work of Trinh et al. (2016b) which exhibits a relationship between hydrogen solubility and free volume in various organic oxygenated solvents, showing thus that hydrogen solubility is more dependent on the available free space in the bulk than on the nature of the solvent. According to Trinh et al., since hydrogen has no kernel electron, its electronic cloud becomes highly deformable. Considering that repulsion is driven by electronic cloud overlap, this justifies the use of a non-additive combining rule for cross-diameter.

Finally, Figure 2 compares the experimental Henry constants with the simulation results for NaCl molalities of 1 and 2 mol/kg_H2O, but also for 0.5 mol/kg_H2O (not considered in the l_ij regression).

The average deviation between experimental and simulated data is about 6%. Note that six experimental Henry constants are available at higher NaCl molality (4 mol/kg_H2O) (Gerecke and Bittrich, 1971). A preliminary optimization test including these data concluded that a specific l_ij value should be used at this high molality, different of the unique l_ij value previously tuned, making the approach less predictive. Consequently, additional specific forces should be probably taken into account to model data at high salinities, such as polarization forces. Thus, the force field presented in this work is currently restricted to NaCl molalities up to 2 mol/kg_H2O.

From Figure 2, it can be stated that the temperature for which H₂ Henry constant is maximum (and thus its solubility is minimum) is only slightly modified by the NaCl concentration, and remains close to the 330 K value evidenced for pure water. It can also be seen that the salting-out effect is less pronounced at elevated temperature. In this temperature range (typically above 460 K, 187 °C), hydrogen solubility is not significantly modified by salt concentration. This observation is actually not surprising, since the experimental data of Cramer et al. (Cramer, 1982) exhibit the same behavior for methane in salted water.

Modeling electrolytic solutions with an equation of state is a challenging task. As described in introduction, many approaches exist involving various equation of state families. In this work, we focus on the Soreide and Whitson (SW) equation of state (Soreide and Whitson, 1992). The advantage of this model is how convenient it is for simple geological applications: it is a cubic equation of state in which ions or salts are not explicitly considered as specific species, but the salinity is rather implicitly included in the solvent properties (salted water instead of pure water). But in return, this is a highly correlative approach strictly restricted to the validation range of the empirical parameters of the model. The SW model is specifically suited to model gas solubility in brines: in their original work, Soreide and Whitson modeled accurately the solubility of nitrogen, carbon dioxide, hydrogen sulfide and hydrocarbons up to C4 in a large range of temperature (up to 200 °C), pressure (up to 700 bar) and salinity (up to 5 mol/kg_H2O). It was later successfully used for modeling phase equilibrium of light gas + brine systems for reservoir and geochemical studies (e.g.Nichita et al., 2008; Plennevaux et al., 2013; Qiao et al., 2016; Wei and Zhang 2013; Yan et al., 2011).

It is important to stress that for a given i,j pair, different interaction binary parameters k_ij can be used for aqueous phase and vapor phase (see Eqs. (11) and (12)), making this approach heterogeneous and thus not suitable to calculate phase equilibrium in the vicinity of critical points. In this section, an aqueous and a vapor interaction binary parameter between hydrogen and salted water is adjusted.

Note that this binary interaction parameter is between hydrogen and salted water, since this equation of state does not consider ionic species explicitly. This approach is thus fundamentally different from the molecular simulation one described in Section 3, where the optimized binary parameters (ξ_ij = 0 and l_ij = −0.24) represent explicitly the interactions between hydrogen and the ions Na⁺ and Cl⁻.

These SW empirical parameters are adjusted in order to reproduce experimental Henry constants of hydrogen in pure water up to 573 K (smoothed experimental data by correlation 3 are used), experimental Henry constants of hydrogen in salted water from 273 to 323 K and salinities from 0.5 to 2 mol/kg_H2O, and simulated Henry constant obtained from molecular simulation from 273 to 573 K and salinities from 0.5 to 2 mol/kg_H2O. The optimized parameters are given in Table 6. The resulting equation of state matches the experimental data with an average deviation of 0.3% for pure water, and 3% for salted water, as shown in Figure 5.

After the calibration of the SW interaction binary parameters for the aqueous phase, the k_ij in the vapor phase were adjusted to match experimental gaseous composition in binary H₂ + H₂O system (Gillespie and Wilson, 1980) at given temperatures and pressures, for temperatures ranging from 310 to 440 K and pressures from 3 to 140 bar. The optimized values are given in Table 6 and lead to an average deviation between calculated and experimental data of about 6.5% (Fig. 6).

In order to illustrate the possible applications of the thermodynamic work presented in the previous sections, the calibrated SW equation of state was used in several simple calculations relevant for geological situations involving hydrogen.

The first example concerns hydrogen geological storage. The calibrated Soreide and Whitson equation of state can be used to perform flash calculations to determine both the aqueous and vapor phase compositions. Water content in the produced vapor phase, necessary to design surface facilities, can thus be predicted. Taking as examples hypothetical pure hydrogen storages in the formations of Lobodice, in Czech Republic, and Diadema, in Argentina, and considering the thermodynamic conditions given by Panfilov (2015) (respectively 34 °C/90 bar and 50 °C/10 bar) and seawater salinity for connate waters (0.6 mol/kg_H2O NaCl), the developed model predicts water molar fractions in the vapor phase of, respectively, 5.6 × 10⁻⁴ and 1.171 × 10⁻² (Tab. 7). In other words, considering that the two gaseous components and their mixing have similar molar volume, it means that, for each 100 000 m³ of withdrawn gas, a storage in Lobodice would only produce 56 m³ of water vapor versus 1 171 m³ in Diadema. The H₂ concentration in the water at equilibrium with these gas phases can also be calculated and would be 0.056 mol/kg_H2O in the case of Lobodice and 0.006 mol/kg_H2O in the Diadema Formation. As could be expected, pressure is playing a key role on the phase behavior of this system, and can modify by several orders of magnitude the water content in the produced gas.

The geochemistry of the vent fluids around the middle oceanic ridge, especially the Atlantic one that has been studied by various oceanographic campaigns, have all shown high hydrogen concentrations (e.g.Charlou et al., 2002 for the Acores Sea, 36°N). As a second illustrative application case, the phase behavior of native hydrogen in thermodynamic conditions typical of these Mid-Atlantic ridge hydrothermal systems was explored. Hydrothermal fluids temperatures come from (Charlou et al., 2010) and pressures were deduced from the same data-set considering hydrostatic equilibrium. Pure-H₂ solubility calculations were performed for temperatures varying between 300 and 350 °C and pressures from 270 to 320 bar with a salinity close to 0.65 mol/kg_H2O NaCl. In addition, Lost City (about 30°N on the Atlantic mid oceanic ridge) conditions were also investigated (94 °C and 78 bar), as well as the ones of 3 km-deep seawater (2 °C and 300 bar). The calculated solubilities are indicated in Table 8, and can be compared to the H₂ concentrations measured in the corresponding hydrothermal fluid (Charlou et al., 2010). If we consider a negligible effect of the other dissolved gases (measured in lesser concentrations than H₂), these hydrothermal systems should be purely monophasic, due to hydrogen content up-to 70 times below the solubility limit. According to these results, the presence of gas bubbles, described by (Charlou et al., 2010) throughout the smokers of Ashadze (mid oceanic ridge, 13°N), is, as mentioned by the authors of this study, “unexpected”. The phase behavior of the H₂O-H₂-NaCl system alone cannot explain this observation, which should be investigated with a more comprehensive thermodynamic model considering additional species, such as CO₂ or H₂S.

Natural continental hydrogen emissions have been detected on several geographic locations (Larin et al., 2015; Moretti et al., 2018; Prinzhofer and Deville, 2015; Prinzhofer et al., 2019). In order to investigate the phase behavior of such continental natural emissions, hydrogen solubility was plotted versus depth (Fig. 7). Thermodynamic conditions were estimated considering a surface temperature of 20 °C, a geothermal gradient of 30 °C per km and hydrostatic pressure. Calculations were performed for pure water and brines with salinity of 0.6, 1 and 2 mol/kg_H2O NaCl.

Results show that the solubility of hydrogen tends to increase with depth, to rapidly become significant. As an illustration, H₂ solubility calculated at 1000 m in pure water (about 0.07 mol/kg_H2O) can be compared with the maximum concentration of about 0.01 mol/kg_H2O calculated by fluid-rock interactions modeling in Oman ophiolite (Bachaud et al., 2017). This indicates that, in such systems, transport of hydrogen would occur mostly in the aqueous phase and that formation of gaseous H₂ would be a quite shallow process.

Due to the salting-out effect, the solubility decreases with increasing salinity, which is observed in Figure 7. The impact of NaCl concentration does not change significantly between 0 and 3000 m: the H₂ solubility is about 35% lower in 1 mol/kg_H2O NaCl brines than in pure water at each investigated depth, and about 30% lower in 2 mol/kg_H2O NaCl solutions than in 1 mol/kg_H2O NaCl ones.

Solubility of hydrogen was also compared with the one of methane. Solubility of both components was calculated for depths from 0 to 1500 m considering a salinity of 0.1 mol/kg_H2O NaCl. The model parameters for the H₂O + CH₄ + NaCl system come from the original work of Soreide and Whitson (Soreide and Whitson, 1992). Results are plotted in Figure 8.

The solubility of methane is superior to the hydrogen one down to 1300 m and, below that point, H₂ is actually more soluble than CH₄. This result can be put in perspective to some observations presented by Prinzhofer et al. (2018) on the hydrogen and methane accumulations discovered in the Bourabougou field (Mali). Indeed, it seems that highest methane contents in the gas phase are found below 1500 m. If a solubility effect alone cannot explain the gas reservoirs organization, the results presented on Figure 8 indicate that a given-composition gas mixture of CH₄ and H₂ equilibrated with formation water would tend to contain more H₂ than CH₄ above 1300 m, and more CH₄ than H₂ below.

In this work, Monte Carlo simulations have been carried out to generate new simulated data of hydrogen solubility in aqueous NaCl solutions in temperature and salinity ranges of interest for geological applications, and for which no experimental data are currently available. In order to validate the force field used, a review of existing data for the binary system H₂ + H₂O and the ternary system H₂ + H₂O + NaCl is first proposed, all the data being converted into Henry constants. For the simulations, the molecular models of Darkrim for hydrogen (Darkrim et al., 1999), TIP4P/2005 for water (Abascal and Vega, 2005) and OPLS for Na⁺ and Cl⁻ (Chandrasekhar et al., 1984) have been selected since they are able to accurately reproduce phase property of pure components and brine densities. To model solvent-solutes and solutes-solutes interactions, it was shown that the Lorentz-Berthelot combining rule is the most appropriate to reproduce the experimental hydrogen Henry constants. Nevertheless, it was necessary to correct the cross-diameter of the Lorentz-Berthelot combining rule between H₂ and Na⁺ and Cl⁻ with a constant binary interaction parameter. With this force field, simulation results reproduce the hydrogen Henry constant in salted water with an average deviation of 6%. Simulation results accurately predict the expected behavior such as the salting-out effect and the maximum of hydrogen solubility close to 330 K. This force field was then used in extrapolation to generate Henry constant for temperatures up to 573 K and salinities up to 2 mol/kg_H2O. It is worth noticing that the extrapolation at higher salinities is questionable, and adding additional molecular interactions in the force field, such as polarization, could be useful to keep the predictivity of the approach. Using the experimental data and these new simulated data generated by molecular simulation at higher temperatures and salinities, a binary interaction parameter of the Soreide and Whitson equation of state has been fitted in order to extend the validity range of this thermodynamic model for systems containing hydrogen and salted water. This new parameterization of the equation of state allows a fast and reliable calculations of phase equilibrium for the H₂ + H₂O + NaCl system in a range of temperatures and salinities corresponding to subsurface conditions. Consequently, this model allows to predict hydrogen phase behavior in various geological contexts. To illustrate the possible use of the model, simple situations involving hydrogen in subsurface were investigated. In particular, the results showed that:

• at low depths (typically below 200 m) hydrogen solubility is low, suggesting that, close to the surface, transportation through the aqueous phase should only lead to minor hydrogen migrations;

• in underground storages (for which temperature is usually lower than 60 °C), the key factor governing the phase equilibrium is the pressure rather than salinity;

• when temperature and pressure increase, such as in the hydrothermal vents or at important depths, a larger amount of hydrogen can be dissolved in water (typically around 100 times more at a depth of 1000 meter than in surface conditions), indicating that deep aquifers or hydrothermal fluids could be a major transport mode for hydrogen.

Authors gratefully acknowledge Dr. J.C. de Hemptinne and Dr. E. Brosse (IFPEN) for fruitful discussions. This work is extracted from the Master Thesis of Cristina Lopez Lazaro which has been funded by ENGIE.