Effect of Water Content on Thermodynamic Properties of Compressed Hydrogen.

Force field-based molecular simulations were used to calculate thermal expansivities, heat capacities, and Joule-Thomson coefficients of binary (standard) hydrogen-water mixtures for temperatures between 366.15 and 423.15 K and pressures between 50 and 1000 bar. The mole fraction of water in saturated hydrogen-water mixtures in the gas phase ranges from 0.004 to 0.138. The same properties were calculated for pure hydrogen at 323.15 K and pressures between 100 and 1000 bar. Simulations were performed using the TIP3P and a modified TIP4P force field for water and the Marx, Vrabec, Cracknell, Buch, and Hirschfelder force fields for hydrogen. The vapor-liquid equilibria of hydrogen-water mixtures were calculated along the melting line of ice Ih, corresponding to temperatures between 264.21 and 272.4 K, using the TIP3P force field for water and the Marx force field for hydrogen. In this temperature range, the solubilities and the chemical potentials of hydrogen and water were obtained. Based on the computed solubility data of hydrogen in water, the freezing-point depression of water was computed ranging from 264.21 to 272.4 K. The modified TIP4P and Marx force fields were used to improve the solubility calculations of hydrogen-water mixtures reported in our previous study [Rahbari A.;J. Chem. Eng. Data2019, 64, 4103-4115] for temperatures between 323 and 423 K and pressures ranging from 100 to 1000 bar. The chemical potentials of ice Ih were calculated as a function of pressure between 100 and 1000 bar, along the melting line for temperatures between 264.21 and 272.4 K, using the IAPWS equation of state for ice Ih. We show that at low pressures, the presence of water has a large effect on the thermodynamic properties of compressed hydrogen. Our conclusions may have consequences for the energetics of a hydrogen refueling station using electrochemical hydrogen compressors.

Introduction

To supply the energy demand for a growing population worldwide and to reduce carbon emissions, efforts are being made to switch from fossil-based energy production to sustainable energy production using renewable energy sources.[1−5] Current estimations show that in 2050 ca. 75% of the global energy production will be based on fossil fuels.[1] In recent years, hydrogen is considered as one of the most promising alternatives in the energy transition for combating global warming and reducing the carbon footprint from anthropogenic activities.[4,6−11] Hydrogen is one of the most abundant elements found on earth mainly in water and organic compounds.[3] To produce hydrogen as a free gas, energy is required.[2,12] Hydrogen can be a carbon-free renewable energy source depending on the production pathway,[12]e.g., green hydrogen can be produced from water electrolysis in which the required energy is supplied from renewable sources. The high energy density of hydrogen per unit of mass (higher heating value of 141.8 MJ/kg at 298 K)[12] makes it a flexible off-peak energy carrier, especially for storing intermittent renewable energy at peak times, or as a fuel for transportation.[6,13−15] To date, most of the required hydrogen in refineries is produced from steam methane reforming (SMR) with CO2 as a byproduct.[2,4,8,12,16] About 50% of the produced hydrogen worldwide is used in ammonia synthesis plants. Hydrogen is also used as a hydrocarbon upgrader or feedstock in refineries, e.g., reducing nitrogen oxides, sulfur oxides, and other particles causing smog.[6,17]

The heating value of hydrogen (per mass) is higher than that of most common fossil fuels.[2,12] However, due to the low volumetric density of hydrogen at standard temperature and pressure, 273.15 K and 1 bar (10.79 MJ/m3),[9] and to save space, gaseous hydrogen should be compressed, liquefied, or absorbed by metal hydrides or nanoporous materials such as metal–organic frameworks (MOFs).[2,3,7,17−19] Compressing hydrogen increases the heating value (per unit of volume) compared to the heating value of conventional liquid fossil fuels.[12] Compressed hydrogen can be used in fuel-cell electric vehicles or hydrogen-fueled combustion engines.[1,17,20] To facilitate and increase the use of hydrogen-driven vehicles, hydrogen refueling stations (HRSs) should provide compressed hydrogen to fuel-cell electric vehicles or internal combustion hydrogen engines.[21] Current vehicle technologies allow direct on-board hydrogen storage with pressures up to 700 bar inside the tank.[9,18,20,22,23] This allows a driving range of at least 500 km.[1,9,24−26] To meet the target refueling time (3–4 min),[22,23] a large pressure difference is required when filling the hydrogen tank, e.g., the pressure of compressed hydrogen before entering the vehicle should be as high as 875 bar.[21] Therefore, energy- and cost-efficient hydrogen compression technologies are required for operating refueling stations at these high pressures.

Hydrogen compression technologies can be divided into mechanical and nonmechanical categories.[9] There are well-known limitations associated with mechanical compression such as high mechanical losses, high maintenance, high levels of noise, etc.[20] In addition, the compression work can be as high as 30% of the energy of the stored hydrogen,[27] and the efficiency is low when operating at small scales.[9,28] Alternatives to mechanical compression are (1) cryogenic compression, (2) metal hydride or thermally powered compressors,[29] and (3) electrochemical hydrogen compressors (EHCs), also known as electrochemical hydrogen pumps.[9,20,30] A detailed overview of conventional and state-of-the-art technologies for hydrogen compression technologies is provided in ref (9).

Recently, EHCs have become cost-competitive alternatives to their mechanical counterparts.[16,30,31] Due to the advantages over mechanical compression, EHCs are on the way to be commercialized on a large scale. EHCs can generate a high output pressure (up to 1000 bar), which is suitable for refueling hydrogen vehicles.[20,30] EHCs have a higher efficiency compared to mechanical compression, especially for low volumetric flows.[9,32] Since there are no moving parts in an EHC, there are no maintenance costs associated with rotary parts, and the compressor operates silently. The working principle of an EHC is based on the proton-exchange membrane fuel-cell technology.[4,9,30] A schematic representation of the working principle of an EHC is shown in Figure . On the anode side, electricity is used to split the hydrogen molecule into protons and electrons, and the protons pass through the membrane. Hydrated protons (H2O+ in which n is the electro-osmosis coefficient) are transferred to the cathode side.[4,9,13] At the cathode side, the protons are combined with electrons to form hydrogen molecules.[20] As the number of hydrogen molecules increases on the cathode side, the discharge pressure at the outlet of the EHC increases. In EHCs, the driving force for compression is the electrical potential difference between the anode and the cathode. This means that the discharge pressure can be as large as allowed by material limitations. It is important to note that generating a large differential pressure over a single compression cell lowers the efficiency of the EHC and may cause back-diffusion of hydrogen from the cathode to the anode side.[9,20] To circumvent this issue, one can raise the discharge pressure using a cascaded arrangement of single compression cells.[30,33] The highest differential pressure reported in the literature for a single compression cell is 168 bar.[20] Two types of proton conducting membranes are typically used for electrochemical compression:[16][16] (1) proton electrolyte membranes such as Nafion-type membranes,[16,28,34] which are suitable for low-temperature applications; and (2) protonic-ceramic electrolyte membranes, which are suitable for high-temperature applications (e.g., 600 °C).[16] These membranes are highly selective to permeation of protons as impurities cannot pass the membrane.[33] The proton-exchange membrane requires hydration for its functionality. Therefore, it is important that the membrane inside the compressor remains hydrated with water when operating. While compressing the hydrogen stream, the EHC can simultaneously purify the hydrogen stream such that the compressed hydrogen at the outlet is saturated with water, and no other impurities are present. Recently, HyET BV[30,31] has developed state-of-the-art EHCs operating based on the proton-exchange membrane technology. Currently, each compression stack can compress 50 kg H2/day up to 1000 bar. The EHC capacity can be scaled up using several compression stacks.[31,33] The EHC by HyET BV can achieve the same pressure difference as in mechanical compression with fewer steps and silently. The compressed hydrogen is saturated with water and free from other contaminants.

Figure 1

Schematic representation of the working principle of electrochemical compression of hydrogen.[30] Hydrogen molecules are split at the low-pressure side of the membrane (i.e., using platinum-alloy catalysts). Protons are forced through the membrane and form hydrogen molecules at the high-pressure side.

Here, we consider the application of an EHC to generate compressed hydrogen for refueling a hydrogen vehicle. Although the solubility of water in compressed hydrogen at the outlet of the EHC is very low, this affects the thermodynamic behavior of the gas phase (depending on the water content). As the compressed hydrogen gas enters the on-board storage tank, it follows directly from the first law of thermodynamics of an open system that the temperature of hydrogen inside the tank rises. A similar effect happens for the isenthalpic expansion of hydrogen (e.g., using a throttle) due to the negative Joule–Thomson coefficient of hydrogen at high pressures.[23,35] To protect the mechanical integrity of the tank, the maximum temperature limit of the hydrogen in the process of filling should not exceed 85 °C.[26,35,36] To avoid high temperatures inside the hydrogen tank, precooling is required before refueling to bring the temperature of compressed hydrogen down between −33 and −40 °C.[23,26,35] A complete hydrogen fueling protocol is compiled by the Society of Automotive Engineers (SAE) in “Fueling Protocols for Light Duty Gaseous Hydrogen Surface Vehicles (SAE J2601)”.[22,37]

Due to the presence of water content at the output of an EHC, ice formation may take place during precooling, which may clog the hydrogen delivery system. Experimental data on the solubility and other thermodynamic properties of water–hydrogen systems in the published literature are limited, especially at high pressures.[38,39] Published experimental data for solubilities of water–hydrogen mixtures in the gas or/and liquid phase at high pressures can be found in refs (40−49). A summary of the experimental data in the tabulated form can be found in ref (38). However, little is known about the influence of water on the thermodynamic properties of hydrogen at high pressures. Performing experiments for obtaining thermodynamic data at high pressures is challenging due to safety protocols, material limitations, high costs, etc.[50−54] A robust alternative for high-pressure experiments is force-field-based molecular simulation.[55,56] In this work, we use force-field-based molecular simulation to calculate properties of compressed hydrogen by varying the water content at different temperatures and pressures. Molecular simulation is an approach for calculating macroscopic as well microscopic properties of materials based on the knowledge of their constituent molecules and atoms as well as their interactions. Molecular simulation allows the computation of material properties at high pressures, which would be difficult or expensive to access using experiments. By modeling the VLE of water–hydrogen systems, it is possible to compute the amount of water present in compressed hydrogen at temperatures and pressures corresponding to the freezing point of water. In this way, one can predict water/ice formation when precooling the hydrogen stream. The VLE calculations can also be used to predict the freezing-point depression of water in the presence of hydrogen. In this work, the VLE calculations of hydrogen–water are performed at temperatures and pressures corresponding to the water–ice Ih equilibrium in the pressure range between 100 and 1000 bar. Using molecular simulation, we also calculate thermodynamic properties such as thermal expansivity, heat capacity, and the Joule–Thomson coefficient of the water–hydrogen mixtures.

This paper is organized as follows. Thermodynamic properties of compressed hydrogen–water mixtures are investigated using force-field-based molecular simulation. An overview of the molecular models used in this paper is provided in Table . Thermodynamic properties of mixtures are related to the derivatives of extensive properties of the mixtures.[57] One can calculate the appropriate thermodynamic derivatives (i.e., for calculating heat capacities, thermal expansivities, and Joule–Thomson coefficients) from ensemble fluctuations.[55,56,58] In Section , we provide the appropriate expressions for calculating thermodynamic derivatives to obtain the required thermodynamic properties. In Section , we explain how vapor–liquid equilibria of hydrogen–water mixtures along the melting line of ice can be computed. This allows one to estimate how much water would exit the saturated hydrogen stream during deep cooling at low temperatures. To the best of our knowledge, high-pressure experimental data of water–hydrogen systems at temperatures between 264.2 and 272.4 K are not available in the literature. Volumetric data from the literature are used to compute the chemical potential of ice for a large pressure range. Simulation details are provided in Section . All simulation results are presented and discussed in Section . Our conclusions are summarized in Section . The results show that small concentrations of water dissolved in hydrogen in the gas phase can significantly change the properties of the gas phase (i.e., Joule–Thomson coefficient, heat capacity, and thermal expansivity) compared to pure hydrogen at the same conditions. Our results will therefore have consequences for the energetics of a hydrogen refueling station using EHCs including drying or precooling of the compressed hydrogen stream.

Table 1

Chemical Compounds Used for Molecular Simulation in This Worka

chemical name	chemical formula	CAS number	force field
water	H2O	7732-18-5	TIP3P[100]
water	H2O	7732-18-5	modified TIP4P (this work)
hydrogen	H2	1333-74-0	Marx[99]
hydrogen	H2	1333-74-0	Vrabec[102]
hydrogen	H2	1333-74-0	Cracknell[97]
hydrogen	H2	1333-74-0	Hirschfelder[107]
hydrogen	H2	1333-74-0	Buch[106]

For direct computation of chemical potentials, improving the efficiency, and calculating the solubilites, fractional molecules are used. Interaction parameters of the modified TIP4P force field are provided in the Supporting Information (SI). Mixture compositions of water−hydrogen systems are explicitly provided for every temperature and pressure in Table .

Table 2

Mixture Compositions of Saturated Hydrogen–Water Mixtures Obtained from VLE Simulations Using the Marx[99] and TIP3P[100] Force Fields

T (K)	P (bar)	yH2	yH2O
423.15	50.00	0.862	0.138
423.15	80.00	0.914	0.086
423.15	100.00	0.928	0.072
423.15	300.00	0.975	0.025
423.15	500.00	0.984	0.016
423.15	800.00	0.990	0.010
423.15	1000.00	0.992	0.008
366.48	10.00	0.892	0.108
366.48	50.00	0.978	0.022
366.48	80.00	0.986	0.014
366.48	100.00	0.989	0.011
366.48	300.00	0.996	0.004

Thermodynamic Properties of Mixtures Obtained from Ensemble Fluctuations

To compute thermodynamic properties of compressed hydrogen with and without traces of water using molecular simulations, we use derivatives of volume, internal energy, and enthalpy with respect to temperature and pressure. These derivatives are required to calculate properties such as thermal expansivity, heat capacity, and the Joule–Thomson coefficient.[57,59,60] These thermodynamic derivatives are directly obtained from ensemble fluctuations at constant composition.[56,58,61,62] Lagache et al.[58] showed that the derivative of an extensive property X with respect to β = 1/(kBT) (in which kB is the Boltzmann constant and T the absolute temperature) in the NPT ensemble can be obtained from the ensemble fluctuations as followswhere Ĥ = U + PV is the configurational enthalpy of the system,[58]P is the imposed pressure, and U is the potential energy of the system consisting of an intermolecular contribution Uext and an intramolecular contribution Uint. The mathematical proof for this is provided in the Supporting Information. In a similar manner, one can obtain an expression for the derivative of the ensemble average of an extensive property X with respect to pressure[58]The last term on the right-hand side of eq is nonzero if the extensive property X includes a pressure-dependent contribution, for instance, for the configurational enthalpy Ĥ. The term was not explicitly provided in the corresponding equation in ref (58). The derivation of eq is also provided in the Supporting Information. The extensive property X may be replaced by an intensive property x. Using eqs and 2, one can obtain various thermodynamic properties of mixtures from ensemble fluctuations. For instance, combining the definition of thermal expansivity α with eq (58) leads toTo compute the Joule–Thomson coefficient, one needs to compute the heat capacity of the mixture as well. The extensive heat capacity follows from combining the definition of the heat capacity with eq where the enthalpy H also includes the kinetic energy term (K): H = Uint + Uext + K + PV = Ĥ + K. Jorgensen et al.[63] and Lagache et al.[58] emphasize that the expression for the heat capacity should be split into an ideal contribution and a residual contribution.Including the ideal gas contribution of heat capacity from the force field usually leads to overestimating the heat capacity.[64,65] To calculate Cid(T) and Cres(T, P), Lagache et al.[58] proposed the following split for the enthalpy of the systemin which N is the number of molecules in the systems. The ideal gas heat capacity for a pure component can be obtained from tabulated thermodynamic reference tables.[59,66] Alternatively, Cid can be calculated by taking into account the contribution of translational, rotational, vibrational, and electronic energy levels of an isolated molecule, using quantum mechanical calculations or experimental data.[67,68] The residual heat capacity Cres is obtained by combining eqs and 6The Joule–Thomson coefficient μJT is obtained using

Vapor–Liquid Equilibrium of Hydrogen–Water Mixtures along the Melting Line of Ice

Computation of Solubilities in the Liquid and Gas Phases

In molecular simulation, VLEs are conveniently computed by Monte Carlo simulations[55] in the Gibbs ensemble.[69,70] In the Gibbs ensemble, the liquid and gas phases are simulated in different simulation boxes. During the simulation, standard Monte Carlo trial moves are used to generate configurations in each phase and to exchange molecules between the phases. In simulations of multicomponent systems in the Gibbs ensemble, at each temperature and pressure, chemical equilibrium is obtained when the chemical potentials of both phases are equal. From the few published experimental VLE data of hydrogen–water mixtures at high pressures,[38,42,44,47−49] it is known that the solubility of hydrogen in the liquid phase and the solubility of water in the gas phase are very low. The experimental data in refs (40−49) are tabulated in the Supporting Information of ref (38). Simulating the VLE of hydrogen–water mixtures at such dilute concentrations in the Gibbs ensemble[69,70] would be computationally expensive (requiring large system sizes), and the simulation time would be prohibitively long. It is possible to simulate the VLE of hydrogen–water mixtures using a smaller system size; however, one needs to compensate by performing much longer simulations to obtain reliable statistics.[38] Despite its reliability and small finite-size effects (except in the vicinity of the critical point),[71,72] the Gibbs ensemble is not the preferred method here. The liquid phase is almost pure water at high pressures, and the gas phase is almost pure hydrogen.[38] Therefore, one can simulate the liquid phase and the gas phase independently in the isobaric–isothermal ensemble and obtain the chemical potentials of water and hydrogen at phase equilibrium. Since the concentration of the solute is very low in each phase, the chemical potential of the solvent will not be different from that of the pure solvent.

In Monte Carlo simulations, one of the methods to efficiently compute the chemical potentials, especially in dense systems such as liquid water, is the continuous fractional component Monte Carlo (CFCMC) method.[73−75] To compute the chemical potentials of water and hydrogen in the liquid and the gas phase, we use the CFCMC in the NPT ensemble.[73,76] In the CFCMC method, a so-called fractional molecule (one per component type) is added to the system. The interaction potential of the fractional molecule of type i is scaled using a coupling parameter λ ∈ [0, 1]. By constantly performing a random walk between λ = 0 and 1 during the simulation, one can mimic gradual insertion or deletion of a molecule and compute the chemical potential μ of a certain molecule type[77]where μ0(T) is the reference chemical potential, which contains the intramolecular contributions (rotation, vibration, and translation) and the bond dissociation energy of component i.[67] ρ is the number density of i (i.e., the ratio between the number of molecules of i and the volume of the system). ρ0 = 1 Å–3 is the reference number density to make the argument inside the logarithm dimensionless. The reference chemical potential is important when simulating chemical reactions, taking into account free energy differences of bond breaking/formation. In this work, the choice of μ0(T) is not important since chemical reactions are not considered. The first and the second term on the right-hand side of eq are the ideal gas parts of the chemical potential (μid), and the last term on the right-hand side of eq is the excess chemical potential (μex). The excess chemical potential is related to the interactions between a molecule and the surrounding molecules and is independent of the reference chemical potential μ0. For an ideal gas, βμex = 0. p(λ) is the Boltzmann probability distribution of λ. The excess chemical potential can be directly computed using the values of p(λ) at λ = 0 and 1 or by performing extrapolation or interpolation; see eq .[73,78] Further details regarding the CFCMC method and recent advances of this method are provided in ref (73) and references in there.

By including a fractional molecule of water in the liquid phase, it is straightforward to calculate the chemical potential of water in the liquid phase (μH).[73,76] To calculate the excess chemical potential of water in the gas phase, a fractional molecule of water is added to the gas phase (containing hydrogen molecules). The solubility of water in the gas phase is obtained by imposing equal chemical potentials of water in the gas and liquid phases (μH = μH) and solving the ideal gas contribution of the chemical potential of water in the gas phaseIn a similar way, one can compute the solubility of hydrogen in the liquid phase. Since the solubility of hydrogen is very low in the liquid phase, the excess chemical potential of hydrogen can be computed by adding a single fractional molecule of hydrogen in the liquid phase. The ideal gas chemical potential of hydrogen (related to the number density) in the liquid phase is initially unknown. By imposing equal chemical potentials of hydrogen in the gas phase (almost pure hydrogen) and the liquid phase (μH = μH), the ideal gas chemical potential of hydrogen and the number density of hydrogen in the liquid phase at equilibrium are solved usingEquations and 11 are derived in the Supporting Information.

Computation of the Chemical Potential of the Most Stable Phase

Precooling compressed hydrogen containing water may lead to ice formation during the refueling process for a hydrogen-driven vehicle. To predict ice formation, one can use an empirical equation of state for ice and water or compare the chemical potentials of water and ice at different temperatures and pressures. At equilibrium, the chemical potentials of ice and water are equal, which corresponds to the minimum Gibbs energy of the system. By departing from the equilibrium pressure, the most stable phase (liquid or ice) will be the phase that lowers the Gibbs energy of the system the most.[59,79]

Simulating the molecular structure of ice,[80−84] solid–liquid equilibria, and computing the chemical potential of ice require special simulation techniques such as anisotropic NPT simulations, the Einstein crystal method,[85] the Einstein molecule approach,[86] and/or the Parrinello–Rahman-type sampling.[83] For details on free energy calculations of ice, the reader is referred to the works of Vega et al.[81,86−88] An alternative and a simpler pathway for computing the chemical potential of ice is to compute the chemical potential of liquid water along the experimental solid–liquid equilibrium. It is straightforward to compute the chemical potential of liquid water from simulations in the CFCNPT ensemble.[73,76] Alternatively, thermodynamic tables, or the empirical equation of state for ice and water,[89,90] can be used to compute the chemical potentials of ice/water along the melting line. In this work, we use both molecular simulation of liquid water and the IAPWS equation of state for water,[91] which is implemented in REFPROP.[92] The equations of state used in REFPROP for pure hydrogen, liquid water, and ice are fitted to experimental data. Since the compressibility of both water and ice is low, it is possible to calculate the chemical potentials of liquid water and ice at different melting temperatures by varying the pressure, usingwhere Tm is the melting temperature of ice. The notation for melting pressure Pm(T) is used to emphasize that the melting pressure is a function of the melting temperature. From eqs and 13, it is clear that knowledge of the molar volume of ice and water is required to calculate the chemical potential of the most stable phase at different pressures for each melting temperature. Marion et al.[93] reported experimental values for the molar volume υice of pure ice at P0 = 1 atm, without cracks or gas bubbles, based on crystal lattice parameters for ice Ih[94]where T is in units of kelvin. To calculate the molar volume of ice at P ≠ P0, the isothermal expansivity of ice is required. Marion et al.[93] estimated the isothermal compressibility β of ice usingwhere T is in units of kelvin, and the molar volume of ice at pressure P follows fromThe data for molar volume of liquid water are available in REFPROP.[91,92] By combining eqs and 14–16, one can estimate the chemical potentials of ice and water at different pressures as a function of melting temperature.

Calculating the Freezing-Point Depression of Water with Dissolved Hydrogen

The solubility of hydrogen in liquid water affects the freezing temperature of ice. In the Supporting Information, we show that the freezing-point depression is obtained from[57]where Δh̅fus =h̅ice – h̅water is the enthalpy of fusion of ice, R is the gas constant, and xH is the mole fraction of hydrogen dissolved in liquid water. In eq , it is assumed that the solution is ideal and that the solvent is nonvolatile. The difference between the enthalpy of ice and water is a function of temperature and does not vary significantly with pressure.[95] Osborne and Dickinson[95] performed an experiment to obtain the enthalpy of fusion of ice, which is Δh̅fus = 6010.44 J/mol. Osborne and Dickinson[95] also tabulated the difference in enthalpy between water and ice at different temperatures. As a service to the reader, in Table S1 of the Supporting Information, this enthalpy table is provided.

Simulation Details

In this work, two types of simulations are performed: simulations with fractional molecules and simulations without fractional molecules. To calculate the chemical potentials of water or hydrogen from simulations, fractional molecules are used. For computation of thermodynamic derivatives from ensemble fluctuations, no fractional molecules are used. To calculate the thermodynamic derivatives of water–hydrogen mixtures from ensemble fluctuations, Monte Carlo simulations are performed in the NPT ensemble, as explained in Section . To compute the chemical potentials of water and hydrogen in the gas and liquid phases, independent simulations are performed by expanding the NPT ensemble with fractional molecules of water and hydrogen. All simulations are performed using open-source Brick-CFCMC software.[96] All molecules are considered rigid objects. Intermolecular interactions consist of pairwise Lennard-Jones (LJ) and Coulomb interactions. Force-field parameters for hydrogen and water models are provided in the Supporting Information of ref (38). The interaction parameters of the modified TIP4P force field are provided in the Supporting Information. The Ewald summation with a precision of 10–6 is used to handle electrostatic interactions. Lennard-Jones interactions are truncated at 12 Å, and analytic tail corrections are applied. Periodic boundary conditions are used. The Lorentz–Berthelot mixing rules are applied for cross-interactions between different interaction sites.[55,56] For the temperature range considered in this study, no quantum corrections are needed for hydrogen.[97,98]

For all simulations of pure water and pure hydrogen, 730 molecules of water and 600 molecules of hydrogen are used. The mixture compositions of the water–hydrogen mixtures are based on the VLE calculations of ref (38), using the Marx[99] force field for hydrogen and the TIP3P[100] force field for water. All mixture compositions are provided in Table . In the VLE calculations of water–hydrogen systems, simulations using the TIP3P force field predict the composition of the gas phase in good agreement with available experimental data.[38] Also, compared to other commonly used force fields, the chemical potential of the TIP3P water is in good agreement with the chemical potential of water obtained from the IAPWS equation of state.[91] It is important to note that while the performance of the TIP3P force field is not the best overall among the different force fields of water,[101] it performs reasonably well for calculating the chemical potential of water. As shown in ref (38), this property is crucial for VLE calculations of water–hydrogen systems. An in-depth discussion on the role of the chemical potential of water can also be found in Section S5 of the Supporting Information. The hydrogen force field by Vrabec et al.[102] was also used to investigate the thermodynamic properties of water–hydrogen mixtures.[102]

For simulations including fractional molecules, equilibration runs were performed between 2 × 105 and 5 × 105 cycles until the optimum weight function is obtained. In each cycle, k trial moves are performed to generate new configurations (to be accepted or rejected based on the Metropolis acceptance rules[103]). In Brick-CFCMC software, k is set as the maximum between the number of molecules in the system and the number 20.[96] In the CFCMC algorithm, the scaling of the interactions of the fractional molecules is performed as follows. Starting from the value of λ = 1, first the electrostatic interactions are switched off. This is followed by switching off the van der Waals interactions. Details about scaling the interactions of the fractional molecule can be found in the manual of Brick-CFCMC.[96] A weight function W(λ) was used in the sampling scheme to ensure that the observed probability distribution pobs(λ) is approximately flat. The Boltzmann probability distribution p(λ) is recovered by p(λ) ∼ exp[−W(λ)] × pobs(λ). The weight function is computed using the Wang–Landau algorithm[104] followed by an iterative scheme. Applying a weight function results in a simulation where all different λ states are visited with the same probability. This improves the accuracy of the computed chemical potentials.[105] The probabilities p(λ = 1) and p(λ = 0) are directly sampled using the scheme in ref (105) without using interpolation or extrapolation. For each temperature and pressure, 107 production runs are performed, and the uncertainties of ensemble averages are calculated using block averaging (five blocks). Simulations in the CFCNPT ensemble are performed using trial moves with the following probabilities: 1% volume changes, 35% translations, 30% rotations, 17% λ changes, 8.5% random reinsertions, and 8.5% identity changes. The random reinsertions and identity changes are combined in a hybrid trial move.[76,96]

For simulations of hydrogen–water mixtures without fractional molecules, 2 × 105 equilibration cycles are performed followed by 8 × 106 production cycles. Simulations in the NPT ensemble are performed using trial moves with the following probabilities: 1% volume changes, 49.5% translations, and 49.5% rotations. The uncertainties of ensemble averages are calculated using block averaging (five blocks). The confidence level of all reported uncertainties in this work is 95%.

Results and Discussion

Thermophysical Properties of Pure Hydrogen at High Pressures

Thermodynamic properties of pure hydrogen were computed using different force fields. The Cracknell,[97] Buch,[106] Hirschfelder,[107] Vrabec,[102] and Marx[99] force fields for hydrogen were used to compute the thermal expansivity α, specific heat capacity, c, and the Joule–Thomson coefficient μJT. For pure hydrogen, these thermodynamic properties were also obtained from an empirical equation of state for hydrogen in REFPROP.[92] The results are shown in Figure , and the raw data are provided in Table . The Vrabec force field predicts α and c, in good agreement with the empirical equation of state obtained from REFPROP.[92] This quantitative agreement is consistent with the fact that the hydrogen force field by Vrabec et al. is optimized for thermal properties of hydrogen and the speed of sound.[102] None of the force fields could very accurately predict the Joule–Thompson coefficient from experimental data; however, the force fields can capture the behavior of the Joule–Thomson coefficient qualitatively, i.e., the Joule–Thomson coefficient decreases with the pressure as it is the case for the Joule–Thomson coefficient obtained from experimental values. The results from the Marx and Vrabec force fields show the least deviation. The Joule–Thomson coefficient of TIP3P water is also calculated, and the results are compared to those obtained from REFPROP.[92] For pure water, the simulations were performed along the VLE line at pressures between 10 and 100 bar using experimental densities. The results are shown in Figure S1 of the Supporting Information. Excellent agreement is found between the experiments and the simulations. The Marx and Vrabec force fields, together with the TIP3P force field, are used in this work to further investigate thermal properties of hydrogen–water mixtures.

Figure 2

Comparison of different force fields to predict (a) thermal expansivity, (b) constant pressure heat capacity, and (c) the JT coefficient of pure hydrogen in the gas phase at T = 323 K and pressures ranging between P = 10 and 1000 bar. Experimental data from REFPROP[92] (lines), molecular force fields: Cracknell[97] (circles), Buch[106] (upward-pointing triangles), Hirschfelder[107] (squares), Vrabec[102] (downward-pointing triangles), and Marx[99] (diamonds). Raw simulation data are provided in Table . Error bars are smaller than symbols.

Table 3

Thermodynamic Properties of Pure Hydrogen Obtained from the Empirical Equation of State for Hydrogen in REFPROP,[92] and Different Hydrogen Force Fields: Hirschfelder,[107] Vrabec,[102] Buch,[106] Cracknell[97] and Marx[99],a

(volumetric) thermal expansivity αP (bar–1)
T (K)	P (bar)	REFPROP[92,109]	Cracknell[97]	Buch[106]	Hirschfelder[107]	Vrabec[102]	Marx[99]
323.15	100	0.00297	0.003002(3)	0.002984(4)	0.003012(4)	0.002938(2)	0.003001(3)
323.15	200	0.00283	0.002902(5)	0.002868(5)	0.002905(6)	0.002790(3)	0.002890(7)
323.15	300	0.00270	0.002806(3)	0.002747(5)	0.002795(4)	0.002650(3)	0.002776(4)
323.15	400	0.00256	0.002703(5)	0.002630(2)	0.002681(2)	0.002523(3)	0.002656(7)
323.15	500	0.00244	0.002602(3)	0.002517(2)	0.002579(3)	0.002408(3)	0.002544(2)
323.15	600	0.00233	0.002513(5)	0.002411(3)	0.002476(4)	0.002300(4)	0.002444(7)
323.15	800	0.00213	0.002344(5)	0.002222(3)	0.002290(5)	0.002116(4)	0.002254(9)
323.15	1000	0.001958	0.002193(5)	0.002070(4)	0.002131(3)	0.001963(3)	0.002093(4)

Numbers in brackets indicate uncertainties in the last digit (95% confidence interval).

Thermodynamic Properties of Hydrogen–Water Mixtures at Low Water Concentrations in the Gas Phase

In Figure , the thermal expansivities α are shown for force-field combinations TIP3P–Marx and TIP3P–Vrabec at temperatures T = 366 and 423 K. The composition of the mixture at each temperature and pressure, provided in Table , is based on the solubility calculations of ref (38). The mixture compositions, and raw data of Figure , obtained from molecular simulations are provided in Tables and 5. Thermodynamic properties of hydrogen–water mixtures in Table computed using force-field combinations TIP3P–Marx and TIP3P–Vrabec were used for thermodynamic properties provided in Table . Thermodynamic properties of pure hydrogen in the gas phase obtained from REFPROP are provided in Table S2 of the Supporting Information. To better see the influence of the water content on the properties of compressed hydrogen, the corresponding properties of pure hydrogen from simulations and empirical data are also provided in this figure. The simulation results show that the values of thermal expansivity of the mixture are very similar to the behavior of pure hydrogen, with deviations of ca. 1% or less for pressures higher than 50 bar, at T = 366 K, and for pressures higher than 800 bar, at T = 423 K. In both cases, the mole fraction of water is ca. yH = 0.01. At pressures below 50 bar, the influence of water content on the properties of the mixture becomes significant.

Figure 3

Thermal expansivities of pure hydrogen in the gas phase and hydrogen in the gas phase that is saturated with water. Triangles, pure hydrogen; circles, hydrogen–water mixtures; lines, empirical data for pure hydrogen from REFPROP.[92] Thermal expansivities of mixtures computed using the TIP3P[100] and Marx[99] force fields are shown in subfigures (a) and (b) for T = 366 and 423 K, respectively. Thermal expansivities of mixtures computed using the TIP3P[100] and Vrabec[102] force fields are shown in subfigures (c) and (d) for T = 366 and 423 K, respectively. The composition in the gas phase is obtained from VLE simulations of the water–hydrogen mixture at T = 423 and 366 K. Raw data obtained from molecular simulations (including the composition in the gas phase) are provided in Tables and 5. Thermodynamic properties of pure hydrogen in the gas phase, obtained from REFPROP, are provided in Table S2 of the Supporting Information. Error bars are smaller than symbols.

Table 4

Thermodynamic properties of hydrogen-water mixtures and pure hydrogen in the gas phase obtained from CFCMC simulations in the NPT ensemblea

T (K)	P (bar)	αP (K–1)	cP (J/(mol K))	μJT (K/bar)
Hydrogen–Water
423.15	50	0.00262(1)	37.1(3)	0.208 (8)
423.15	80	0.002475(7)	34.2(2)	0.062 (4)
423.15	100	0.002435(8)	33.36(6)	0.033 (4)
423.15	300	0.002172(3)	31.01(4)	–0.035(1)
423.15	500	0.002013(4)	31.05(3)	–0.041(1)
423.15	800	0.001814(4)	30.99(4)	–0.045(1)
423.15	1000	0.001701(3)	31.01(3)	–0.046(1)
366.48	10	0.00283(1)	32.5(2)	0.36(3)
366.48	50	0.002713(4)	30.09(8)	–0.012(3)
366.48	80	0.002675(5)	29.9(1)	–0.026(2)
366.48	100	0.002656(4)	29.82(1)	–0.029(1)
366.48	300	0.002459(6)	30.23(2)	–0.038(1)
Pure Hydrogen
423.15	50	0.002326(2)	29.37(1)	–0.038(2)
423.15	80	0.002302(4)	29.45(1)	–0.040(2)
423.15	100	0.002287(3)	29.50(1)	–0.040(2)
423.15	300	0.002134(2)	29.94(1)	–0.043(1)
423.15	500	0.001987(4)	30.23(2)	–0.046(1)
423.15	800	0.001802(5)	30.56(3)	–0.047(1)
423.15	1000	0.001694(4)	30.71(2)	–0.047(1)
366.48	10	0.002719(4)	29.18(1)	–0.03(1)
366.48	50	0.002686(2)	29.35(1)	–0.033(1)
366.48	80	0.002661(2)	29.47(1)	–0.033(1)
366.48	100	0.002641(3)	29.53(1)	–0.035(1)
366.48	300	0.002450(4)	30.10(2)	–0.040(1)

Water was simulated using the TIP3P[100] force field, and hydrogen was simulated using the Marx[99] force field. The mixture compositions of the water hydrogen mixtures are obtained from VLE simulations using the Marx[99] and TIP3P[100] force fields,[38] and are provided in Table . The numbers in brackets indicate uncertainties in the last digit (95% confidence interval).

Table 5

Thermodynamic Properties of Hydrogen–Water Mixtures and Pure Hydrogen in the Gas Phase Obtained from CFCMC Simulations in the NPT Ensemblea

T (K)	P (bar)	αP (K–1)	cP (J/(mol K))	μJT (K/bar)
Hydrogen–Water
423.15	50.00	0.00260(1)	37.0(2)	0.196(8)
423.15	80.00	0.00245(1)	34.0(1)	0.049(6)
423.15	100.00	0.00241(1)	33.3(2)	0.021(6)
423.15	300.00	0.002107(3)	30.70(3)	–0.048(1)
423.15	500.00	0.001939(3)	30.73(4)	–0.052(1)
423.15	800.00	0.001729(3)	30.51(3)	–0.054(1)
423.15	1000.00	0.001619(4)	30.49(2)	–0.054(1)
366.48	10.00	0.00284(1)	32.5(1)	0.37(5)
366.48	50.00	0.002689(6)	29.90(3)	–0.030(5)
366.48	80.00	0.002640(3)	29.76(6)	–0.044(2)
366.48	100.00	0.002610(3)	29.69(1)	–0.047(1)
366.48	300.00	0.002368(3)	29.88(1)	–0.053(1)
Pure Hydrogen
423.15	50	0.002309(4)	29.31(1)	–0.057(4)
423.15	80	0.002280(4)	29.37(2)	–0.055(3)
423.15	100	0.002258(2)	29.40(1)	–0.056(1)
423.15	300	0.002070(3)	29.66(1)	–0.057(1)
423.15	500	0.001913(3)	29.87(2)	–0.056(1)
423.15	800	0.001717(3)	30.07(2)	–0.056(1)
423.15	1000	0.001614(2)	30.20(1)	–0.055(1)
366.48	10.00	0.002716(4)	29.17(1)	–0.05(1)
366.48	50.00	0.002663(4)	29.28(1)	–0.051(3)
366.48	80.00	0.002623(4)	29.35(1)	–0.052(2)
366.48	100.00	0.002601(2)	29.40(1)	–0.052(1)
366.48	300.00	0.002364(2)	29.76(1)	–0.054(1)

Water was simulated using the TIP3P[100] force field, and hydrogen was simulated using the Vrabec[102] force field. The mixture compositions of the water hydrogen mixtures are obtained from VLE simulations using the Marx[99] and TIP3P[100] force fields,[38] and are provided in Table . The numbers in brackets indicate uncertainties in the last digit (95% confidence interval).

A similar behavior is observed for the heat capacity of hydrogen–water mixtures at T = 366 and 423 K; see Figure . Compared to pure hydrogen, the heat capacity of hydrogen–water mixtures deviates around 1% or less for pressures higher than 100 bar, at T = 366 K, and for pressures around 800 bar, at T = 423 K. As can be seen in Figure , the presence of water at pressures lower than 50 bar significantly changes the heat capacity of the mixture compared to pure hydrogen. For instance, based on the raw data in Table , one can observe that the presence of around 13% water (molar) in the gas phase results in a 26% deviation of the heat capacity at P = 50 bar and T = 423 K. A similar observation is made based on the raw data in Table . Based on eq , it is expected that the Joule–Thomson coefficient shows a similar behavior as the thermal expansivity and heat capacity. In Figure , it can be seen that at around 9% water content in the gas phase (P = 80 bar), the deviation of the Joule–Thomson coefficient is around 255%, at T = 423 K. Based on the raw data in Table , the deviation of the Joule–Thomson coefficient is around 189%. Due to the presence of water, the Joule–Thomson coefficient has a positive sign for pressures below 50 bar. The deviation from the Joule–Thomson coefficient of pure hydrogen is reduced to 2% at P = 1000 bar and T = 423 K, where the mole fraction of water is around yH = 0.008. At T = 366 K and P = 300 bar, the mole fraction of water is around 0.004 and the Joule–Thomson coefficient of the mixture deviates around 3.6% from the Joule–Thomson coefficient of pure hydrogen. It can be observed that at relatively low pressures, the effect of water on the properties of the gas phase (i.e., the Joule–Thomson coefficient, heat capacity, and thermal expansivity) is significant. However, at pressures higher than 700 bar, the effect of water content on the properties of hydrogen gas is almost negligible. Since the Joule–Thomson coefficient is also a function of molar volume, additional plots of the Joule–Thomson coefficients as a function of molar volume of the mixture are shown in Figure S2 of the Supporting Information.

Figure 4

Heat capacities of hydrogen in the gas phase and hydrogen in the gas phase that is saturated with water. Triangles, pure hydrogen; circles, hydrogen–water mixtures; lines, empirical data for pure hydrogen from REFPROP.[92] Heat capacities of mixtures computed using the TIP3P[100] and Marx[99] force fields are shown in subfigures (a) and (b) for T = 366 and 423 K, respectively. Heat capacities of mixtures computed using the TIP3P[100] and Vrabec[102] force fields are shown in subfigures (c) and (d) for T = 366 and 423 K, respectively. The composition in the gas phase is obtained from VLE simulations of the water–hydrogen mixture at T = 423 and 366 K. Raw data obtained from molecular simulations (including the composition in the gas phase) are provided in Tables and 5. Thermodynamic properties of pure hydrogen in the gas phase, obtained from REFPROP, are provided in Table S2 of the Supporting Information. Error bars are smaller than symbols.

Figure 5

Joule–Thomson coefficients of hydrogen in the gas phase and hydrogen that is saturated with water. Triangles, pure hydrogen; circles, hydrogen–water mixtures; lines, empirical data for pure hydrogen obtained from REFPROP.[92] Joule–Thomson coefficients of mixtures computed using the TIP3P[100] and Marx[99] force fields are shown in subfigures (a) and (b) for T = 366 and 423 K, respectively. Joule–Thomson coefficients of mixtures computed using the TIP3P[100] and Vrabec[102] force fields are shown in subfigures (c) and (d) for T = 366 and 423 K, respectively. The composition in the gas phase is obtained from VLE simulations of the water–hydrogen mixture at T = 423 and 366 K. Raw data obtained from molecular simulations (including the composition in the gas phase) are provided in Tables and 5. Thermodynamic properties of pure hydrogen in the gas phase, obtained from REFPROP, are provided in Table S2 of the Supporting Information. Error bars are smaller than symbols.

Chemical Potential of Ice, Water, and Hydrogen at the Melting Temperature

In Figure , the phase boundary diagram of water and ice Ih is shown. The corresponding data for this figure are obtained from experimental data.[91,92,108] The melting temperatures and pressures in Figure are used for calculating chemical potentials of ice and water, from molecular simulation and empirical data. The IAPWS equation of state in REFPROP[91] is used to calculate the chemical potential of ice and water along the melting line (Tm, Pm(Tm)). The results are shown in Figure , and the raw data are provided in Table . For hydrogen, excellent agreement is observed between the chemical potentials of the Marx force field and the chemical potentials obtained from REFPROP. The chemical potential computed with the TIP3P force field differs by a maximum of ca. 1 kJ/mol compared to the chemical potential obtained from the empirical equation of state of water.[91] This difference does not cause significant deviation between the computed solubilities and the experimental solubilities of water in compressed hydrogen, as was shown in ref (38).

Figure 6

Phase boundary diagram of water and ice Ih taken from experimental literature data.[91,92,108] The open circle indicates the solid–liquid–vapor triple point: P = 611.657 Pa and T = 273.16 K.[92]

Figure 7

Chemical potentials of pure hydrogen and water along the melting line of ice (Tm, Pm) computed from CFCMC simulations using the TIP3P[100] and Marx[99] force fields. Lines indicate the chemical potentials calculated from the empirical equation of states using REFPROP.[92] Raw data are provided in Table . Error bars are smaller than symbols.

Table 6

Densities and Chemical Potentials of Pure Water and Pure Hydrogen along the Melting Line of Watera

pure water (TIP3P[100])		CFCNPT		REFPROP[92]
Tm (K)	Pm (bar)	ρ (kg/m3)	(μ – μ0) (kJ/mol)	ρ (kg/m3)	(μ – μ0) (kJ/mol)
272.4	100	1010.6(4)	–34.4(1)	1004.8	–35.3491
270.79	300	1021.4(5)	–34.15(5)	1014.6	–35.0356
269.06	500	1032.1(4)	–33.80(3)	1024.1	–34.7292
266.73	750	1045.2(3)	–33.44(3)	1035.7	–34.3561
266.24	800	1047.1(5)	–33.40(6)	1038.0	–34.2829
264.21	1000	1056.9(4)	–33.12(5)	1046.9	–33.9943

The numbers in brackets indicate uncertainties in the last digit (95% confidence interval).

By changing the pressure, the equilibrium between ice and water shifts toward the most stable phase depending on the pressure. At different melting temperatures of ice, one can compute the chemical potentials of the most stable phase (liquid water or ice) as a function of pressure using eqs and 14–16. This allows one to calculate the chemical potential of ice, based on empirical data without explicitly simulating the crystal structure of ice. The chemical potential of the most stable phase at each temperature can be calculated as a function of pressure using eqs and 13. An example is provided in Figure where the melting point of ice is at T = 266.73 K and P = 750 bar. It can be seen in Figure that by increasing the pressure, the chemical potential of water would be lower than that of ice, which is the most stable phase for P > 750 bar. Decreasing the pressure would lead to having ice as the most stable phase.

Figure 8

Chemical potentials of ice and water at Tm = 266.73 K as a function of pressure. Lines, crystalline ice; dashed lines, liquid phase. At this temperature, ice and water are in equilibrium at Pm = 750 bar. At P > Pm, the chemical potential of water is smaller than the chemical potential of ice, which indicates the most stable phase for that pressure range. For P < Pm, ice is the most stable phase. Raw data are provided in Table S3 of the Supporting Information, including different melting temperatures 264.21 K < Tm < 272.4 K. Equations and 14–16 and the IAPWS equation of state for water[91] were used to generate the data.

VLE of the Water–Hydrogen Mixtures at Freezing Temperature of Water

The VLE calculations of water–hydrogen mixtures are performed along the melting line of ice; see Figure . The corresponding melting temperatures and pressures for ice Ih are directly obtained from REFPROP.[91,92] One can obtain the same melting curve from ref (108). The melting temperatures selected in this work are between T = 272.4 and 264.21 K, which correspond to pressures between P = 100 and 1000 bar. It is crucial to consider the phase diagram of Figure in precooling of the compressed hydrogen mixed with water (from the output of an EHC) and to calculate the amount of water dissolved in hydrogen at low temperatures. By performing independent simulations of the liquid phase (water) and the gas phase (hydrogen), the chemical potentials of water and hydrogen are computed along the melting line of ice. In Table , we provide an example to emphasize that the chemical potential of the solvent (water) is not affected by the chemical potential of the solute (hydrogen) at low concentrations. Due to the application of the EHC, the pressure range is limited to the range between P = 100 and 1000 bar. Using eqs and 11, the VLEs of water–hydrogen mixtures were computed at the melting temperature of ice as a function of pressure. At each melting temperature, different pressures were considered for the VLE calculations (while remaining in the liquid water region). The results are shown in Figure . The raw data corresponding to the liquid phase are provided in Table , and the raw data corresponding to the gas phase are provided in Table . In Figure , only simulation results are provided, as to the best of our knowledge, no experimental VLE data of water and hydrogen, along the melting line of ice, are reported in the published literature. From simulations, the solubility of water is the highest at ca. 90 ppm (molar) at 272.4 K (and 100 bar) and the lowest at ca. 7 ppm (molar) at 264.21 K (and 1000 bar). Such low solubilities justify our choice of not using the Gibbs ensemble. The solubilities of hydrogen in the liquid phase vary between 0.01 and 0.001, corresponding to 272.4 and 264.21 K, respectively. It is well-known that to obtain solubilities in good agreement with experimental data, the accuracy of the chemical potentials of water and hydrogen obtained from the force fields is crucial. In the Supporting Information (Section S5), we provide an extended discussion on the importance of the chemical potentials (focusing on water) regarding the accuracy of solubility calculations for water–hydrogen systems. Additional VLE calculations of water–hydrogen mixtures were performed using the TIP3P and TIP4P/2005 force fields. For the purpose of this discussion, a modified force field based on the TIP4P/2005 force field is introduced in the Supporting Information, which improves the accuracy of solubility calculations simultaneously in the liquid and the gas phase.

Table 7

Computed Chemical Potentials of Water in Water–Hydrogen Mixtures Where the Mole Fraction of Hydrogen (xH) Is Close to Saturationa

T (K)	P (bar)	NH2	NH2O	xH2	(μ – μ0) (kJ/mol)
268	875	0	730	0	–33.2(2)
268	875	1	730	0.001	–33.3(2)
268	875	3	730	0.004	–33.3(2)
268	875	5	730	0.007	–33.3(2)
268	875	10	730	0.014	–33.1(2)
273	300	0	730	0	–34.1(2)
273	300	1	730	0.001	–34.1(2)
273	300	3	730	0.004	–34.1(2)
273	300	5	730	0.007	–34.1(2)
273	300	10	730	0.014	–34.1(2)
273	500	0	730	0	–33.8(2)
273	500	1	730	0.001	–33.7(2)
273	500	3	730	0.004	–33.7(2)
273	500	5	730	0.007	–33.7(2)
273	500	10	730	0.014	–33.7(3)
273	875	0	730	0	–33.1(2)
273	875	1	730	0.001	–33.1(2)
273	875	3	730	0.004	–33.1(3)
273	875	5	730	0.007	–33.2(2)
273	875	10	730	0.014	–33.1(3)
273	1000	0	730	0	–32.9(2)
273	1000	1	730	0.001	–32.9(2)
273	1000	3	730	0.004	–32.8(3)
273	1000	5	730	0.007	–32.9(2)
273	1000	10	730	0.014	–32.8(3)

The mixture is defined using the TIP3P[100] and Marx[99] force fields.[38]N is the number of molecules of component i. Numbers in brackets indicate uncertainties in the last digit (95% confidence interval).

Figure 9

VLE of water–hydrogen at low temperatures as a function of pressure, obtained from eqs and 11. (a) Mole fraction of water (yH) in the gas phase. (b) Mole fraction of hydrogen (xH) in the liquid phase. For each temperature (Tm), the points on the dashed line correspond to the composition calculated at the melting pressure (Pm) of ice; see Table . For each value of Tm, solubility calculations are performed for pressures in the range between Pm and 1000 bar. The melting temperature and pressure were selected based on the values from the IAPWS equation of state.[91,92] Raw data, including the melting temperatures and pressures, are provided in Tables and 9. Error bars are smaller than symbols.

Table 8

Densities of Water (ρH), Chemical Potentials of Water (μH – μH0), Excess Chemical Potentials of Hydrogen at Infinite Dilution in Liquid Water (μHex), and the Solubilities of Hydrogen in Liquid Water, Obtained from Simulations in the CFCNPT Ensemble[73,76],a

T (K)	P (bar)	ρH2O (kg/m3)	(μH2O – μH2O0) (kJ/mol)	μH2ex (kJ/mol)	xH2/10–3
272.4*	100	1009.8(5)	–34.5(1)	9.18(3)	1.44(2)
272.4	300	1020.3(6)	–34.10(5)	9.58(5)	4.0(1)
272.4	500	1029.8(5)	–33.8(1)	10.01(4)	6.3(1)
272.4	750	1040.7(4)	–33.28(5)	10.46(3)	8.9(1)
272.4	1000	1051.5(5)	–32.89(4)	10.97(5)	11.2(3)
270.79*	300	1021.2(5)	–34.1(1)	9.56(4)	4.01(6)
270.79	500	1030.2(4)	–33.84(8)	9.93(4)	6.4(1)
270.79	750	1041.8(3)	–33.37(5)	10.45(4)	8.8(2)
270.79	1000	1052.5(3)	–32.85(7)	10.96(3)	11.0(2)
269.06*	500	1031.9(5)	–33.84(7)	9.90(4)	6.3(1)
269.06	750	1042.9(4)	–33.4(1)	10.38(4)	8.9(1)
269.06	1000	1053.6(1)	–32.96(6)	10.90(3)	11.0(1)
266.73*	750	1045.0(4)	–33.47(8)	10.35(6)	8.7(3)
266.73	1000	1055.1(6)	–33.1(1)	10.80(4)	11.2(2)
264.21*	1000	1056.9(3)	–33.15(6)	10.76(4)	11.0(2)

Water was simulated using the TIP3P force field, and a fractional molecule of hydrogen using the Marx force field[99] was added to the system. Stars indicate melting temperatures and the corresponding melting pressures for ice. xH is the mole fraction of hydrogen in liquid water. Numbers in brackets indicate uncertainties in the last digit (95% confidence interval).

Table 9

Densities of Hydrogen (ρH), Chemical Potentials of Hydrogen (μH – μH0), Excess Chemical Potentials of Water at Infinite Dilution in Hydrogen (μHex), and the Solubilities of Water in Compressed Hydrogen, Obtained from Simulations in the CFCNPT Ensemble[73,76],a

T (K)	P (bar)	ρH2 (kg/m3)	(μH2 – μH20) (kJ/mol)	μH2Oex (kJ/mol)	yH2O/10–6
272.4*	100	8.3776(6)	–13.31(1)	0.112(5)	87(5)
272.4	300	22.396(3)	–10.54(1)	0.394(9)	31(1)
272.4	500	33.475(2)	–9.10(1)	0.748(8)	20.8(5)
272.4	750	44.391(4)	–7.81(1)	1.24(1)	14.6(7)
272.4	1000	53.080(8)	–6.77(1)	1.746(7)	12.1(2)
270.79*	300	22.517(1)	–10.46(1)	0.400(7)	27.4(9)
270.79	500	33.635(3)	–9.02(1)	0.741(8)	18.7(8)
270.79	750	44.584(4)	–7.74(1)	1.228(9)	12.9(4)
270.79	1000	53.283(7)	–6.71(1)	1.73(1)	10.7(2)
269.06*	500	33.812(2)	–8.95(1)	0.73(1)	16.0(3)
269.06	750	44.788(2)	–7.67(1)	1.22(1)	11.7(4)
269.06	1000	53.505(3)	–6.65(1)	1.714(9)	9.5(5)
266.73*	750	45.072(5)	–7.58(1)	1.195(5)	10.1(3)
266.73	1000	53.808(3)	–6.56(1)	1.70(2)	8.1(3)
264.21*	1000	54.141(6)	–6.47(1)	1.671(5)	6.8(2)

Hydrogen was defined using the Marx force field,[99] and a fractional molecule of water using the TIP3P force field[100] was added to the system. Stars indicate melting temperatures and the corresponding melting pressures for ice. yH is the mole fraction of water in the gas phase. Numbers in brackets indicate uncertainties in the last digit (95% confidence interval).

Based on the solubility of hydrogen in liquid water along the melting line, one can calculate the freezing-point depression of water using eq . The results are shown in Table . The largest change in the freezing point of ice is around ΔTF = 1 K, corresponding to 1.1% dissolved hydrogen in water at P = 1000 bar. At lower pressures, the freezing-point depression is less than 1 K.

Table 10

Freezing-Point Depression (ΔTF) of Water as a Function of the Mole Fraction of Dissolved Hydrogen (xH)[57],a

Tm (K)	Pm (bar)	xH2	ΔTF (K)
272.4	100	0.00144(2)	0.148(2)
270.79	300	0.00401(6)	0.407(6)
269.06	500	0.0063(1)	0.63(1)
266.73	750	0.0087(2)	0.86(2)
264.21	1000	0.011(2)	1.06(2)

Tm and Pm denote the freezing point of pure water. ΔTF is obtained from eq . The mole fractions xH are taken from Table . Numbers in brackets indicate uncertainty in the last digit.

Conclusions

In this work, we investigated the thermodynamic properties of hydrogen–water mixtures using molecular simulation. Using molecular-based modeling, we obtained important thermochemical properties of hydrogen–water mixtures, which are not readily available from experiments. Ensemble fluctuations were used to calculate thermodynamic derivatives from molecular simulation. The chemical potentials of water and ice were calculated along the melting line of ice. Different force fields for hydrogen (Cracknell,[97] Buch,[106] Hirschfelder,[107] Vrabec,[102] and Marx[99]) were used to calculate thermal expansivities, heat capacities, and the Joule–Thomson coefficients at high pressures. The heat capacities and thermal expansivities calculated using the Vrabec force field were in good agreement with experimental data. The results from other force fields agreed qualitatively well with the empirical results (e.g., following similar trends). No force field could very accurately calculate the Joule–Thomson coefficient of hydrogen; however, the qualitative behavior of the thermodynamic properties of hydrogen is correctly captured. The results obtained from the Vrabec and Marx force fields are closest to the experimental data. The TIP3P force field for water was used for simulating hydrogen–water mixtures and performing the VLE calculations, at temperatures and pressures corresponding to the melting line of ice Ih. From the commonly used rigid water models with point charges, the TIP3P force field is one of the best force fields for calculating the solubility calculations without using any binary interaction parameters. To examine the effect of the chemical potential on the VLE of water–hydrogen mixtures, a modified force field based on the TIP4P/2005 force field was introduced. Compared to the TIP4P/2005, the accuracy of the computed chemical potentials was significantly improved for the modified force field, which consequently led to some loss of accuracy for density predictions in the liquid phase. The computed solubilities of water in the gas phase using the modified TIP4P force field were similar to those obtained using the TIP3P force field. However, the modified TIP4P force field outperformed the TIP3P force field for the solubility of hydrogen in liquid water. Simulation results showed that thermodynamic properties of compressed hydrogen significantly change as the mole fraction of water in the gas phase is increased. For example, in sharp contrast to pure hydrogen, the Joule–Thomson coefficient of a saturated hydrogen–water mixture, containing 0.7% water (at T = 423.15 K and P = 100 bar), has a positive value, clearly showing that the effect of water cannot be neglected. At pressures above 700 bar, the solubility of water in hydrogen is low so that the behavior of the hydrogen–water mixture becomes very similar to the behavior of pure hydrogen. Experimental validation of the effect of water on thermodynamic properties of hydrogen would be very helpful and complementary to the simulations results. Unfortunately, to the best of our knowledge, such experimental data of the exact system at similar conditions are not available in the literature. The chemical potentials of water and ice along the melting of ice were obtained from molecular simulation and the IAPWS equation of state. For each melting temperature, the chemical potentials of water and ice were calculated as a function of pressure using experimental volumetric data. From the VLE of water–hydrogen, the compositions in the gas and liquid phase were obtained at low temperatures. The solubility calculations were performed in the pressure range between 100 and 1000 bar. The solubility of water in the gas phase is the highest at 100 bar (272.4 K), which is around 90 ppm (molar), and the solubility of hydrogen in the liquid phase is the highest at 1000 bar (264.21 K), which is around 1%. Based on the solubility data in the liquid phase, the freezing-point depression of water was calculated. The largest freezing-point depression is around 1.1 K corresponding to a pressure of 1000 bar and temperature of 264.21 K. For lower pressures, the freezing-point depression is very small. Our simulation results may have consequences for the energetics of a hydrogen refueling station using EHCs. This includes drying or precooling of the compressed hydrogen stream.