Evaluation of cubic, PC-SAFT, and GERG2008 equations of state for accurate calculations of thermophysical properties of hydrogen-blend mixtures

Hydrogen (H 2 ) is a clean fuel and key enabler of energy transition into green renewable sources and a method of achieving net-zero emissions by 2050. Underground H 2 storage (UHS) is a prominent method offering a permanent solution for a low-carbon economy to meet the global energy demand. However

Hydrogen (H 2 ) is a clean fuel and key enabler of energy transition into green renewable sources and a method of achieving net-zero emissions by 2050. Underground H 2 storage (UHS) is a prominent method offering a permanent solution for a low-carbon economy to meet the global energy demand. However, UHS is a complex procedure where containment security, pore-scale scattering, and largescale storage capacity can be influenced by H 2 contamination due to mixing with cushion gases and reservoir fluids. The literature lacks comprehensive investigations of existing thermodynamic models in calculating the accurate transport properties of H 2 -blend mixtures essential to the efficient design of various H 2 storage processes. This work benchmarks cubic equations of state (EoSs), namely Peng-Robinson (PR) and Soave Redlich-Kwong (SRK) and their modifications by Boston-Mathias (PR-BM) and Schwartzentruber-Renon (SR-RK), for their reliability in predicting the thermophysical properties of binary and ternary H 2 -blend mixtures, including CH 4 , C 2 H 6 , C 3 H 8 , H 2 S, H 2 O, CO 2 , CO, and N 2 , in addition to Helmholtz-energy-based EoSs (i.e., PC-SAFT and GERG2008). The benchmarked models are regressed against the experimental data for vapor-liquid equilibrium (VLE) that covers a wide range of pressures (0.01 to 101 MPa), temperatures (92 K to 367 K), and mole fractions (0.001 to 0.90) of H 2 . The novelty of this work is in benchmarking and optimizing the parameters of the mentioned EoSs to study VLE envelopes, densities, and other critical transport properties, such as heat capacity and the Joule-Thomson coefficient of H 2 mixtures in a wide range of associated conditions. The results highlight the significant effect of the temperature-dependent binary interaction parameters on the calculations of thermophysical properties. The SR-RK EoS demonstrated the highest agreement with VLE data among the cubic EoSs with a low root mean square error and absolute average deviation. The PC-SAFT VLE models demonstrated results comparable to the SR-RK. The sensitivity analysis highlighted the high influence of impurity on changing the thermophysical behavior of H 2 -blend streams during the H 2 storage process.

Introduction
The worldwide main energy supply is achieved by fossil fuels and high-carbon emitters, such as oil and natural gas (International Energy Agency (IEA), 2018b,a, 2019; Mohanty et al., 2021). The dramatic growth of the worldwide population causes an additional challenge for energy demand, which is expected to rise by about 40% by 2040, calling for an urgent need for effective clean energy solutions (Ali, 2018;Ali et al., 2022;Rahbari et al., 2019). Renewable sources, such as wind power and solar energy, have emerged as promising clean alternatives (Gallo et al., 2016). However, their intermittency and seasonal nature depend on availability, geographic location, and atmospheric events, which results in a discrepancy between the energy supply and demand, leading to uncontrollable energy output (Gallo et al., 2016;Moustakas et al., 2020;Rahbari et al., 2019;Saidi and Omri, 2020). To remedy that, hydrogen (H 2 ) is considered a clean energy fuel that can provide a sustainable energy market and overcome intermittent production issues (Mahlia et al., 2014).
Currently, H 2 is produced from various sources, depending on economics and feasibility (Dawood et al., 2020). For instance, green H 2 is produced through an electrolysis process or biomass gasification. Similarly, blue H 2 (the main source of H 2 production) is produced from natural gas through steam methane reforming and coal gasification (Dawood et al., 2020). In addition, H 2 will play a critical role in energy transition and mitigate heavy carbon emissions from petrochemical industries (DNV, 2020;Edelenbosch et al., 2017;Hanley et al., 2018;Lazarou et al., 2018;McPherson et al., 2018). A vast expansion of the H 2 economy requires a massive storage capacity of Giga to Terawatt-scale (which can be achieved via underground H 2 storage (UHS), such as deep aquifers or depleted oil and gas reservoirs) compared to the limited storage and discharge capacity of the existing surface facilities (i.e., tanks and pipelines) (Crotogino et al., 2010;Panfilov, 2010;Pfeiffer et al., 2017;Taylor et al., 1986). The main types of H 2 production and expected operating conditions of the H 2 -based industry are described in Fig. 1.
Since the advent of the first underground storage of natural gas in the 1970s, the closest approximation to H 2 storage is deep aquifers, depleted oil and gas reservoirs, and coal seam gas and salt caverns Pan et al., 2021;Tarkowski, 2019;Zivar et al., 2021). However, recovering H 2 from these reservoirs requires the push of a cushion gas or the presence of impurities in the reservoir, which macerates the mixing zone Hosseini et al., 2022;Iglauer et al., 2021). The behavior of the mixing zone can be studied via various parameters, including the differences in density and viscosity, molecular diffusion, mechanical dispersion, and mobility ratios (Feldmann et al., 2016). The gas mixtures have different thermophysical properties that affect the accuracy and prediction of different scenarios (Ma et al., 2019;Oldenburg and Pan, 2013) but can lead to engineering issues related to toxicity, safety, compression, and transportation in the gas streams (Shao et al., 2014). Many typical impurities, such as CH 4 , N 2 , O 2 , Ar, CO 2 , and H 2 S, might be found in transportation pipelines (Heinemann et al., 2021;Li and Yan, 2009). Therefore, accurate calculations to determine uncertainty in reservoir modeling is crucial.
The knowledge of pure H 2 thermodynamics is well established (Michels and Goudeket, 1941;Seward and Franck, 1981), although the available experimental data for H 2 -blend mixtures do not cover the full range of gas mixtures or the range of operational conditions for UHS in the H 2 economy (Hassanpouryouzband et al., 2020). Therefore, semi-empirical equations of state (EoSs) are employed to calculate the properties over an extended range of pressures and temperatures. Cubic EoSs are widely used in reservoir compositional simulators. Several studies have intensively investigated their reliability (Abbott, 1979;Danesh et al., 1991;Yang et al., 2003). The common conclusion of these studies has demonstrated varying reliability for a wide range of properties, component types, and temperature/pressure conditions. Two famous classical EoSs are Peng-Robinson (PR) (Peng and Robinson, 1976) and Soave Redlich-Kwong (SRK) (Soave, 1972). Both EoSs are typically used in engineering applications to predict compositional phase behavior of gas mixtures and model recovery processes (Carroll, 2010;Li, 2008;Wei et al., 2011).
However, challenges arise when classical cubic EoSs are extended to calculate phase equilibrium and mixture density in elevated conditions (Wilhelmsen et al., 2017).
The thermodynamic modeling of vapor-liquid equilibria (VLE) for H 2 -blend mixtures have been of a recent interest to many researchers (Alanazi et al., 2022;Hassanpouryouzband et al., 2020). Some researchers investigated modeling the VLE using classic cubic (i.e., PR) with and without temperature-dependent binary interaction parameters calculated through a group-contribution approach (Koulocheris et al., 2019;Qian et al., 2013). Other work studied specifically H 2 mixed with n-alkane binary mixtures using a group contribution SAFT EoS with a kij group contribution method (Le thi et al., 2006).
In this study, cubic EoSs (RK, SRK, and PR) and their modifications are evaluated against the perturbed-chain statistical associating fluid theory (PC-SAFT) (Gross and Sadowski, 2001) and GERG2008 Kunz andWagner, 2012 in modeling the thermophysical properties of H 2 mixtures with other gases. Binary interaction parameters are regressed to fit the experimental reference data for binary mixtures and provide optimized coefficients for each EoS. This work assesses the capabilities of the thermodynamic models to provide accurate predictions for various H 2 concentrations in the mixture at high pressure and temperature conditions relevant to UHS (i.e., pressure up to 80 MPa and temperature from 313 K to 415 K) (Hassanpouryouzband et al., 2020). Furthermore, the models evaluate the influence of impurities on the H 2 physical transport for storage applications.
The paper is structured as follows. First, it discusses the methodology and thermodynamic modeling using various EoSs. Then, the second section briefly describes the thermodynamic calculation method of the different types of the applied EoSs. The following section lists the experimental reference data collected and used in the present work. Afterward, the paper analyzes and demonstrates the findings of the approach, starting with the results of running the regression against experimental data, comparing and benchmarking the results of different cubic EoSs. The generated thermodynamic models from cubic EoSs are compared to those of PC-SAFT and GERG2008, including assessing density matching. Finally, the best models from cubic EoSs, PC-SAFT, and GERG2008 are used to run prediction scenarios to demonstrate the expandability to conditions not covered by the experimental work, mostly relevant to UHS, as mentioned earlier.

Methodology
A complete set of experimental data extended over the full range of reservoir conditions and compositions is rare. Hence, developing thermodynamic models to accurately simulate the thermodynamics in a reservoir is necessary. Before building thermodynamic models and running regression scenarios, accurate knowledge of the uncertainties in the experimental data is required for accurate results. Hence, only experimental data with reported uncertainty information are included in this work. Datasets with insufficient data points are excluded.
This work includes only datasets with full information on pressure, temperature, and compositions and information on standard deviations for each measurement. Although many publications have provided many binary and ternary VLE datasets, only a limited subset satisfies these constraints (Skiborowski et al., 2018). Most default values for standard deviations found in the literature assume, for example, an error estimate of 1% for the vapor mole fraction (Finlayson, 2012;Skiborowski et al., 2018). Further, sophisticated objective functions were used to regress thermodynamic parameters against experimental data. Finally, the parameters with the least root mean square error (RMSE) were used to predict different properties for several isothermal systems. For example, the RMSE between the experimental mole fraction of component y i,exp and the calculated mole fraction of component i, y i,cal over the total number of components n is expressed in Eq. (1) (Barker, 1953): This work used ASPEN Plus (v. 12.0) (Aspen Technology Inc., 2015), a popular application, to simulate the thermophysical behavior of a fluid mixture. Parameter regression was performed using the available regression methods in ASPEN properties. The most statistically reliable parameter estimates are obtained using a maximum likelihood function (MLF). Assuming that all measurements are independent of each other and that the measurement noise follows a Gaussian distribution with a zero mean, the MLF can be obtained by a weighed least-squares minimization, for which the weights (w n ) are used to weight each date group. The standard deviation (STD) is used to normalize all measurements within the data group. In ASPEN Plus, the MLF considers all compositions (liquid phase mole fraction (x) and vapor-phase mole fraction (y)), the temperature (T ) and pressure (P), as presented in Eq. (2): where NG is the number of the data group, NP denotes the number of points in each data group, and NC represents the total number of components. The errors in the dependent parameters are minimized by fine-tuning the model parameters.

Cubic equations of state
The EoSs are mathematical correlations that aim to calculate the thermodynamic behavior of a fluid by describing the interrelationship between the pressure (P), temperature (T ), volume (V ), and composition (x i ). Most cubic EoSs used for hydrocarbons are pressure-explicit equations. The volume is commonly solved, then the rest of the properties are derived using thermophysical relationships (Diamantonis et al., 2013). Several developments of existing cubic EoSs have followed the introduction of the first cubic EoS, the van der Waals EoS, and have been reviewed extensively in the literature (Economou, 2010;Valderrama, 2003;Mangold et al., 2019).
Modified versions of the cubic EoS have the same cubic structure but with a more generalized temperature-dependent term (a(T )) (Redlich and Kwong, 1949;Balaji et al., 2011;Mangold et al., 2019;Reid et al., 1987). Daridon et al. (1993) proposed a general formalism of the cubic EoS based on the work by Schmidt and Wenzel (1980), as indicated in Eq. (3): where R is the universal gas constant, V m denotes the molar volume, and u and w are parameters of the generalized EoS. In addition, a and b represent constants that depend on the component, where a represents the attraction between the molecules, whereas b defines the volume of a pure component as a function of critical temperature (T c ) and critical pressure (P c ). The assumption is that b is independent of temperature, whereas a is a function of temperature, allowing a more accurate calculation of the vapor-phase pressure of a pure component (Mathias et al., 1991). Several efforts have been exerted to derive mathematical formulations for the temperature dependence term (a(T )) with a high level of consensus, as presented in Eqs. (4) to (6) (Mangold et al., 2019;Mathias and Copeman, 1983;Twu et al., 1991): where Ω a , Ω b , u, and w are unitless constants that vary based on the developed EoS. The purpose of introducing the α-function (α(T )) into the equation is to achieve better matching with experimental data, especially at declined temperature and critical zone regions. Several studies have focused on modifying the α-function to improve the prediction of various fluid systems (El-Twaty and Prausnitz, 1980;Twu et al., 1996;Valderrama, 2003;Wenchuan and Chongli, 1989). The forms of different α-functions in this work are summarized in Table 1. The work introduced by Boston and Mathias (1980) attempts to improve the developed relations Table 1 The α-functions as a function of reduced temperature (T r = T /T c ) and the parameters (u, w) for the general cubic equations of state investigated in this work.  Graboskl and Daubert (1978), Kabadi and Danner (1985), Soave Abbott (1979), Peng and Robinson (1976) Boston-Mathias Boston and Mathias (1980), Mathias (1983), Mathias and Copeman (1983) Schwartzentruber- Schwartzentruber and Renon (1989) T r denotes the reduced temperature; C d , C r , and d are SR-RK parameters.
to cover highly polar substances, such as water, carbon dioxide, and carbon monoxide, by introducing a polar parameter in the α-function (Mathias, 1983). However, at higher densities, the predictions of the pressure (P) -volume (V ) -temperature (T ) state are still not accurate, particularly, for H 2 -blend mixtures, attributed to the quantization of translational motion and the quantum nature of H 2 (Sadus, 1992). Therefore, Schwartzentruber and Renon (1989) introduced three additional polar parameters (i.e., p o , p 1 , p 2 ) to improve the calculations presented by Mathias (1983). The early EoSs used the mixing rules of the single temperatureindependent binary interaction parameter (k a,ij is not a function of temperature) van der Waals proposed, which allows the extrapolation of the calculations over a wide range of temperatures and pressures (Poling et al., 2001). A basic approach for the determination of the mixture composition as a function of the cross parameters a ij and b ij is the use of quadratic mixing rules (Privat et al., 2016;Sandler, 2006;Valderrama, 2003;Vidal, 1978). The simplest quadratic form, including the volume translation parameter (c ij ) of a mixture, can be given in an arithmetic or geometric mean, as shown in Eqs. (7) to (12): where k a,ij , k b,ij , l ij , and X ij denote binary interaction parameters; c i is the volume translation factor; and a i and b i represent the pure component parameters responsible for molecular interaction and co-volume, respectively. The properties of pure components found in H 2 mixtures in relation to the storage process in the investigated EoSs are listed in Table A.1 in the Appendix. These properties are used to predict the thermodynamics behavior of the mixtures using different EoSs.
For H 2 -hydrocarbon mixtures, some proposed mixing methods are based on improving the calculation of k a,ij (Moysan et al., 1986;Valderrama et al., 1990;Yang et al., 2003). Other methods improve the cubic EoS by combining more complex mixing rules or by adding a second interaction parameter (k b,ij ), commonly added in the co-volume mixing rule (Gao et al., 2003;Huang et al., 1994;Ioannidis and Knox, 1999). The mixing rules and formulations of the EoSs in this work are listed in Table 2. Each of the interaction parameters (k a,ij , k b,ij , l ij , and X ij ) can have one or more coefficient parameters determined by fitting to experimental data.

Perturbed-chain statistical associating fluid theory (PC-SAFT) equation of state
The perturbation theory helps describe the effects of molecular attributes, such as the size, shape, and interactions between molecules, with modifications of the expressions for dispersion forces (Nikolaidis et al., 2021). Similar to other higher-order statistical associating fluid theory (SAFT) EoSs (Chapman et al., 1989;Huang and Radosz, 1991), PC-SAFT is based on statistical mechanics. In addition, the PC-SAFT was developed by Gross and Sadowski in 2001 using the perturbation theory (Tumakaka et al., 2005). The theoretical bases of SAFT models are founded on the first-order perturbation thermodynamic theory by Wertheim (1986Wertheim ( , 1984a, leading to the development of EoSs, such as the one introduced by Jackson et al. (1988) and Chapman et al. (1989). The perturbation-based models are often introduced to Table 2 Expressions of various cubic equations of state (EoSs) in this work with mixing rules and associated binary interaction parameters (k a,ij , k b,ij , and l ij ).

EoS
Function form Mixing rule k ij and l ij Reference Redlich and Kwong (1949) PR Abbott (1979), Peng and Robinson (1976) SRK Graboskl and Daubert (1978), Soave (1972) Boston and Mathias (1980) SR-RK Schwartzentruber and Renon (1989) represent a simplified solution for a given molecular model. In PC-SAFT, the underlying molecular model is described as a coarsegrained representation of the molecules and their intermolecular interactions, as illustrated in Fig. 2. The principal idea of using the perturbation solutions is to split the total intermolecular forces into a reference term representing repulsive interactions and a perturbation or correction term that accounts for attractive forces. The attractive forces are further divided into different contributions. Theoretically, the first term is usually known as a function of temperature, density or pressure, and composition. The SAFT and PC-SAFT EoSs are expressed as a summation of the reduced residual Helmholtz free energy (F res ) for each contributor term that represents the type of intermolecular force that occurs in the system. The residual Helmholtz free energy is the same as the Helmholtz free energy at the same temperature and volume but excluding the ideal gas Helmholtz free energy. Thus, the molecular interaction forces for a specific number of molecules (N i ) of each component, volume (V ), and density (ρ) in PC-SAFT are written in Eq. (13), as follows: where k B is the Boltzmann constant. The right-side expression in Eq. (13) represents the hard-chain reference fluid that characterizes the PC-SAFT. The superscripts in the various Helmholtz energy contributions denote the contribution from the chain formation (hc), hard-sphere repulsion (hs), dispersion (disp), association (assoc), and inter-polar (polar) interactions, respectively. In PC-SAFT, three parameters for each pure component are required that account for the nonassociation components: the number of molecular chain segments (M), dispersion energy between segments (ε), and either the diameter (σ ) or volume of the chain segment (v 00 ). For the pure components with association interactions, two more parameters are included: the association volume (κ AB ) and association energy between sites (the molecules, ε AB ).
Following the methodology adopted in this work, the above parameters in PC-SAFT were adjusted to fit the experimental data for the pure component vapor and liquid saturation pressures.
The PC-SAFT is extended to mixtures by modifying σ mix and ε mix using mixing rules (Chapman et al., 2006(Chapman et al., , 1989 derived from the single-fluid theory by van der Waals, as indicated in Eqs. (14) and (15): The association parameters, such as the dispersion interaction, are calculated using the Lorentz-Berthelot combining rules (Chapman et al., 1989;Gross and Sadowski, 2001). Accordingly, the dispersion cross energy between segments (ε ij ) and the diameter of the chain segment (σ ij ) are given below (Eqs. (16) and (17)): The combining rules incorporate the binary interaction parameter (k ij ) allowing a direct comparison with other EoSs used in this work. In addition, the binary interaction parameter can be extended to include complex temperature-dependent terms with multiple coefficients similar to the one in this study, as revealed in Eq. (18): where a ij , b ij , c ij , d ij , and e ij are equation parameters.

GERG-2008 equation of state
The thermodynamics calculations of fluids using the GERG2008 EoS are based on a multifluid approximation expressed explicitly in Helmholtz free energy α as a function of the temperature T, density ρ, molar composition vector x, and independent variables of a mixture (Kunz and Wagner, 2012). The Helmholtz free energy f (ρ, T , x) has two components: one represents the ideal gas phase of the mixture f o as a function of ρ, T, and x and another component accounts for the residual contribution of the gas mixture f res , including the influence of intermolecular forces, as displayed in Eqs. (19) to (23): where f res is a function of the residual of dimensionless Helmholtz free energy of component i(f res oi ), the composition (x), and the departure function (∆f r ), expressed in Eq. (20): where f res oi is a function of the reduced mixture density (δ r ) and the inverse reduced mixture temperature (τ ), which are expressed in Eqs. (21) and (22), respectively: where ρ r (x), and T r (x) denote the reducing functions for the mixture density and mixture temperature, respectively. The departure function (∆f r ) is expressed as follows (Kunz and Wagner, 2012;Rowland et al., 2016): where f res ij is the part of ∆f r that depends only on δ, and τ .
The schematic diagram of the calculations of the dimensionless Helmholtz free energy (f ) using the GERG2008 EoS is in Fig. 3,

Reference data
The collected experimental data on H 2 binary and ternary mixtures cover a wide range of compositions, temperatures, and pressures, as listed in Tables A.2 and A.3 in Appendix, respectively, including VLE, densities, and other data on transport properties. The composition of H 2 in the binary mixtures ranges from a mole fraction of 0.001 to 0.9 tested at temperatures and pressures that range from 63.2 K to 588.7 K and 0.01 to 138.98 MPa, respectively. However, the data on H 2 ternary mixtures include a mole fraction of 0.01 to 0.27 of H 2 composition conducted over experimental conditions at temperatures of 80 K to 323.15 K and under 3 to 28 MPa of pressure. The models are calibrated by fine-tuning the binary interaction (k ij ) coefficients for each EoS, allowing an overall variation of 10% in the RMSE using the presented objective function (Eq. (2)) in the methodology section. The experimental data are first quality checked, and similar isothermal tests are compared against each other. All temperature ranges are plotted for trend analysis. Data on ternary H 2 mixtures are limited to VLE with little information about the solubility and uncertainty level of the reported experimental results.

Phase equilibrium of H 2 -blend mixtures
The experimental VLE data for eight binary H 2 mixtures containing CH 4 , C 2 H 6 , C 3 H 8 , H 2 S, CO 2 , CO, H 2 O, and N 2 were used to validate the results of various EoSs. The work in the section below benchmarks the VLE calculations generated using the cubic EoSs. First, the binary interaction parameters were optimized by regression against experimental data. Thereafter, the parameters were used to predict the missing data region and were plotted along with various mixtures. Finally, the percentage of average absolute deviation (AAD, %) and RMSE were used to provide a qualitative indication of the matches.

Cubic equations of state and their modifications
The predictions of the H 2 /CO 2 mixtures are plotted in Fig. 4, presenting an example of the validation of the VLE models against the experimental data. The regression was performed using the temperature-dependent coefficients of the binary interaction parameters over different forms depending on the type of classical cubic EoSs, as listed in Table 2. The results reveal a comparable

Table 3
Percentage of average absolute deviation (AAD, %) between experimental data and estimates using the investigated cubic equations of state with the root mean square error (RMSE), fitted using regression over various binary interaction coefficients for H 2 /CO 2 mixtures. prediction of H 2 mixture phase diagrams that can be achieved with the proposed iterative regression process. However, PR, SRK, and BM-PR EoSs fail to accurately match the bubble pressure curves. Hence, a quantitative analysis is performed based on the RMSE and AAD to determine which EoS has the best matching, as presented in Table 3. The lowest RMSE was achieved using the SR-RK EoS with very close matching, particularly to the experimental dew-point pressure, whereas other cubic EoSs failed to accurately match the liquid curves. The better matching achieved using SR-RK is attributed to the number of interaction parameters (i.e., k a,ij , k b,ij , and l ij ). Hence, the SR-RK EoS models are used as a representative of the cubic EoSs in the comparison with noncubic EoSs (PC-SAFT and GERG2008) in subsequent sections.

Cubic EoS (SR-RK) versus PC-SAFT
Comparing the results from cubic EoS, namely the SR-RK, with isothermal phase diagrams of the studied binary H 2 mixtures predicted using PC-SAFT establishes an understanding of the strengths and weaknesses of each approach. Similar to the cubic EoS, the binary interaction parameters in PC-SAFT were  Table 4. The optimized coefficients of the binary interaction parameters regressed for SR-RK and PC-SAFT are provided in Tables 5 and 6, respectively.  Table 6 Optimized coefficients of the binary interaction parameters in SR-RK equations of state for eight H 2 mixtures. The results indicate that SR-RK models are slightly better at predicting bubble pressure curves than PC-SAFT for the studied H 2 mixtures, particularly mixtures that include one or more carbon atoms in one of the components, such as CO 2 and C 2 H 6 , which are plotted in Figs. 4 and 6, respectively. The results of H 2 /N 2 mixture (Fig. 7) exhibit a slight shift in the phase envelop for both EoSs, although SR-RK has a slightly better match to the experimental results. For the H 2 /CH 4 mixture (Fig. 5), the qualitative analysis demonstrates that both EoSs provide reasonably accurate results; hence, the quantitative analysis is mandatory to determine which EoS is more accurate. The quantitative analyses in Table 7 support the results achieved in the figures and reveal that the predictions of the SR-RK EoS are slightly better than PC-SAFT, as indicated by the RMSE (%) and AAD (%) values. It is worth mentioning that SR-RK EoS modeling utilizes more adjustable coefficients in the regression of the binary interaction parameters than the PC-SAFT EoS to achieve the best agreement possible with experimental data.  Table 7 are low, data scarcity does not allow a conclusive opinion on the validity of the phase equilibrium.

Density of H 2 -blend mixtures
The density of an H 2 mixture varies with depth as the thermodynamic parameters, such as pressure and temperature change. Hence, it is a crucial thermo-physical parameter to model using the best EoSs accurately. Classical cubic EoSs, such as PR, can model the densities of H 2 mixtures by incorporating a correction factor to account for the volume translation, called volumetranslated Peng-Robinson (VTPR), as demonstrated in Fig. 10. However, Our results showed that PR accuracy is not as good as the other investigated EoSs (i.e., GERG2008, PC-SAFT, and SR-RK) for the same range of pressures and temperatures, as shown in Figs. 11 to 13 for the H 2 /N 2 mixtures over varying H 2 mixing concentrations of 25%, 50%, and 75%. The variation in the accuracy between the models and actual experimental densities can be attributed to the measured percentage of the mixing gases and increased pressures, as indicated by the previous studies (Sadus, 1992).
Over the various H 2 mixtures, it is notable that the predictions slightly lose accuracy at higher pressures (i.e., above 60 MPa) and higher mixing concentrations of H 2 (e.g., 75%). In contrast, other EoSs have higher precision with slightly varying relative uncertainties similar to H 2 /N 2 mixtures ( Table 8).
The models created using the extended cubic SR-RK EoS present very good density predictions comparable to the reliable GERG2008 models (Hassanpouryouzband et al., 2020). The PC-SAFT models also exhibit high accuracy in density predictions and can be trusted to be used in compositional simulation models. Similarly, the density models for the H 2 /CH 4 mixture were validated and plotted in Fig. 14, which presents a similar conclusion, except that the pressure range does not mimic underground reservoir conditions. Nevertheless, PC-SAFT models provide better matches than SR-RK with a close agreement with experimental data and GERG2008 results. The work was extended to validate the densities over ternary H 2 -blend mixtures demonstrating a match similar to experimental data, as depicted in Fig. 15.
To demonstrate the applicability of the EoSs in predicting both liquid and vapor densities of H 2 mixtures, we considered an H 2 /CO 2 mixture and calculated the phase densities at different concentration and pressure conditions along three lines, as shown in Fig. 16a. The calculated liquid and vapor densities by SR-RK, PC SAFT, and GERG2008 are plotted on Fig. 16b, c, and d, respectively.

Calculations of VLE and transport properties for H 2 -blend mixtures
Prediction of thermophysical properties is considered a key performance measure to mark the reliability of EoSs for UHS (Heinemann et al., 2021). Accurate thermodynamic models should accurately predict critical transport properties at relevant temperatures and pressures. The section below uses the generated models to assess the gas impurity influence and variation in the component concentration in the blend mixtures on VLE, density, heat capacity, and Joule-Thomson (JT) coefficient.

Vapor-liquid equilibria
An accurate prediction of the bubble pressure curve is critical for understanding the behavior of stored H 2 in subsurface geological formations. The VLE profiles using SR-RK, GERG2008, and PC-SAFT for a range of H 2 concentrations of H 2 /CH 4 are depicted in Fig. 17. The bubble pressure curves and regions close to the critical conditions are comparable, whereas the predictions vary dramatically between dew-point pressure lines.

Density
The thermodynamic models were extended using three ratios of H 2 concentrations (i.e., 10%, 50%, and 90%) to calculate densities for up to 450 K and 100 MPa. The results of the calculations are displayed in Fig. 18. The plots present reasonable predictions comparable among the EoSs and over the extended range of temperatures and pressures.

Heat capacity
The calculations of heat capacity values were similarly extended using the thermodynamic models. The results are illustrated in Fig. 19. The predictions of the heat capacity are consistent at H 2 concentrations lower than 50%. At a 50% H 2 concentration, the heat capacity values are high, and the models start to reveal discrepancies due to the various calculations by each applied EoS. Nevertheless, the overall trend is consistent, where the heat capacity increases monotonically until it reaches a fixed value above 50 MPa.

Joule-Thomson effect
The JT throttling effect describes the change in the temperature of a gas mixture as it expands through a constraint, such as casing nozzles, pipelines, tubing, or an orifice at constant enthalpy (H) (Tarom et al., 2018). The JT effect is a critical element in designing H 2 injection and production schemes during the storage process, represented by the JT coefficient (µ JT ), which is implicitly a function of the gas properties and flow rate that lead to either positive or negative values, depending on the change in pressure (Steffensen and Smith, 1973;Tarom and Hossain, 2017), as follows in Eq. (24): The JT coefficient varies with temperature and pressure and becomes negative when the change in pressure (∂p) takes an sign opposite the change in temperature (∂T ). Unlike most gases (including CH 4 and CO 2 ), H 2 exhibits a negative JT coefficient (Fig. 20) in the typical temperature and pressure operating ranges.
A negative µ JT implies that H 2 causes an increase in temperature with decreasing pressure (i.e., gas expansion) and a decrease in temperature with increasing pressure (i.e., gas compression). This trend is opposite that of CH 4 and CO 2 . For instance, the JT effect during CO 2 injection may cause cooling, leading to hydrate or ice formation (Hoteit et al., 2019). However, when H 2 is injected into a reservoir, the gas expands at the bottom hole of the well and across perforations, resembling a throttling effect. Depending on the pressure drop range and injection duration, this mechanism may cause a localized increase in temperature near the wellbore and surrounding zone. The opposite thermal effect (i.e., cooling) may occur during production; however, it is less significant because the pressure drop between the reservoir and wellbore occurs over a larger radial distance. The resulting cooling effect is not localized and, therefore, ineffective in altering the reservoir temperature. The unintentional blending of H 2 with impurity caused by gas mixing with cushion gases or reservoir fluids can influence the JT coefficient differently. The H 2 JT coefficient decreases when the contamination increases, as depicted in Fig. 21.
Efforts have been established to evaluate the JT effect for natural gases using cubic EoSs and compositional simulation models (Tarom et al., 2018;Wei et al., 2011). Comparing the studied EoS calculations against the experimental data in this work indicates that better matches are provided by GERG2008, supported by the calculated uncertainties (Fig. 23). However, the calculations from SR-RK and PC-SAFT are comparable to GERG2008 and the experimental data. (See Fig. 22.)

Conclusions and remarks
The present work benchmarks the reliability of various cubic EoSs and their modifications (i.e., PR, SRK, BM-PR, VTPR, and SR-RK) and noncubic EoSs (i.e., PC-SAFT and GERG2008) for predicting accurate phase equilibria and thermophysical properties of binary and ternary H 2 -blend mixtures commonly found during the H 2 geo-storage process. The EoSs are regressed against VLE experimental data covering pressures up to 101 MPa, temperatures up to 367 K, and mole fractions from 0.001 to 0.90 of H 2 . Optimized binary interaction parameters were acquired to study VLE envelopes, densities, and other critical transport properties, such as the heat capacity and JT coefficient for the associated reservoir conditions. The results highlight the significant effect of the temperaturedependent binary interaction parameters on the calculations of thermophysical properties. The SR-RK EoS demonstrated the highest agreement with VLE data among the cubic EoSs. Hence, SR-RK was compared to noncubic EoSs (PC-SAFT). The results revealed that the predictions of PC-SAFT and SR-RK EoSs for VLE are comparable, although SR-RK has a slightly better match with the experimental results. The PC-SAFT models have high accuracy in density predictions for binary and ternary H 2 -blend mixtures and can be trusted to be used in compositional simulation models. The PC-SAFT models agree better with the experimental data and GERG2008 than SR-RK. Therefore, any of these EoSs can be applied adequately in compositional simulators and engineering applications, with GERG2008 having higher popularity in the industry because it was specifically optimized for 21 components over various gases with H 2 -blend mixtures.
The sensitivity analysis highlighted the high influence of impurity on changing the thermophysical behavior of H 2 -blend streams during the H 2 storage process. The sensitivity of the impurity influence on the thermodynamic and transport properties was evaluated by extending the calibrated models over a wider range of temperatures, pressures, and compositions using the Density models of H 2 + CH 4 binary mixtures compared to experimental data at 323 K and 348 K: (a) 26% H 2 + 74% CH 4 (b) 26% H 2 + 74% CH 4 , and (c) 26% H 2 + 74% CH 4 . Fig. 15. Thermodynamic density models of H 2 + CH 4 + N 2 ternary mixtures compared against experimental data at 323 K and 348 K: (a) 28% H 2 , 40% CH 4 , and 32% N 2 and (b) 25% H 2 , 20% CH 4 , 55% N 2 .
three EoSs. The VLE envelopes changed drastically with the H 2 concentrations, where the two-phase region sharply shrinks as the H 2 concentration decreases. The bubble pressure curves in the VLE envelopes and near-critical region are comparable, except for the VLE predictions of the dew-point pressure lines, which vary significantly. The temperature variation near the wellbore or deep into the reservoir might not be significant, though it varies substantially when other gases are mixed. The impurity influence on the density was assessed by varying the H 2 concentrations. The density of the gas increases as impurities increase. Hence, the heat capacity, which varies with density, was studied. The heat capacity increased monotonically as the H 2 concentration in the mixture increased, indicating that H 2 -rich mixtures have higher heat capacity. The calculations of the JT coefficient using GERG2008 had superior matches, but the results from the SR-RK and PC-SAFT were quite comparable.  The presented data are significantly important for large-scale H 2 geo-storage projects, and reservoir modeling software should account for these observations to accurately determine the feasibility of H 2 geo-storage projects. A future work can be dedicated to evaluate the effect of using a volume translation on improving the thermodynamics calculations of classical cubic EoSs of H 2blend mixtures.

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