Modeling Subsurface Hydrogen Storage With Transport Properties From Entropy Scaling Using the PC‐SAFT Equation of State

Hydrogen is a promising alternative to carbon based energy carriers and may be stored in large quantities in subsurface storage deposits. This work assesses the impact of static (density and phase equilibria) and dynamic (viscosity and diffusion coefficients) properties on the pressure field during the injection and extraction of hydrogen in the porous subsurface. In a first step, we derive transport properties for water, hydrogen and their mixture using the Perturbed‐Chain Statistical Associating Fluid Theory equation of state in combination with an entropy scaling approach and compare model predictions to alternative models from the literature. Our model compares excellently to experimental transport coefficients and models from literature with a higher number of adjustable parameters, such as GERG2008, and shows a clear improvement over empirical correlations for transport coefficients of hydrogen. In a second step, we determine the effect of further model reduction by comparing our against a much simpler model applying empirical transport coefficients from the literature. For this purpose, hydrogen is periodically injected into and extracted out of a dome‐shaped porous aquifer under a caprock. Our results show that density and viscosity of hydrogen have the highest impact on the pressure field, and that a thermodynamic model like the new model presented here is essential for modeling the storage aquifer, while keeping the number of coefficients at a minimum. In diffusion‐dominated settings such as the diffusion of hydrogen through the caprock, our developed diffusion coefficients show a much improved dependence on temperature and pressure, leading to a more accurate approximation of the diffusive fluxes.


of 19
aquifer remains below the entry pressure of the caprock. Furthermore, the distribution of water and hydrogen saturation during injection and extraction as well as the fluctuating pressure field are critical for the safe operation of the storage system, and determine the resulting electrical power output. Another important aspect is the high chemical activity of hydrogen which determines the interaction between hydrogen and microbiological life (Hagemann, 2018) as well as other present chemicals such as acids in the soil (Ali et al., 2022), which in turn might influence the storage capabilities of hydrogen. However, studying chemical or biochemical interactions is beyond the scope of this paper.
Experience with underground storage of hydrogen in porous formations is very limited. Therefore, mathematical and numerical methods, like CFD (Computational-Fluid-Dynamics) simulation tools, are required to understand and predict the complex processes in the underground, as well as to quantify capacity, optimize performance and assess risks associated with subsurface hydrogen storage. Amid et al. (2016) studied the feasibility of storing hydrogen in a depleted gas reservoir using a simple static model fitted to field data. Hydrogen losses from dissolution and diffusion were found to be minimal, while the modeled storage site only showed less than half of the energy deliverability and capacity of its natural gas counterpart. This is due to hydrogen being less compressible and having a significantly lower volumetric energy density than natural gas. Hagemann et al. (2015) focused on the spatial hydrodynamic behavior of hydrogen versus methane in a porous aquifer, concluding that unstable displacement with lateral fingering occurs earlier in hydrogen storage. Feldmann et al. (2016) performed a numerical case study in a depleted gas reservoir. Here, gravity override and viscous fingering was found to play a minor role in the gas saturated reservoir, with a stable displacement of the native gas by the injected hydrogen and a highly hydrogen-concentrated area in the vicinity of the operation well. Pfeiffer et al. (2017) numerically investigated a large subsurface hydrogen storage in a sandstone formation between 400 and 3,000 m depth in the North German Basin with a continuous power output of 245 MW for 1 week and a peak power output of 363 MW. The working conditions in the subsurface where estimated to temperatures of 25°C and pressures from 30 to 70 bar. Their model neglects diffusion and dispersion of hydrogen, which is a suitable approximation during the storage cycles.
Hydrogen has distinct physical and chemical properties compared to other geologically stored fluids such as CO 2 or methane (Heinemann et al., 2021). Key ingredients required for the analysis of subsurface hydrogen storage are therefore accurate transport properties, that is, viscosities and diffusion coefficients, next to densities and water/ hydrogen phase equilibria. In this study, we derive precise equilibrium and transport properties for hydrogen, water and their mixture from entropy scaling Loetgering-Lin et al., 2018) using the Perturbed-Chain Statistical Associating Fluid Theory Equation of State (PC-SAFT EoS) (Gross & Sadowski, 2001) for the subsurface storage of hydrogen. The static properties (densities and phase equilibria) and dynamic properties (diffusion coefficients and viscosities) are provided by a single thermodynamic model. This approach has clear advantages over empirical correlations from the literature, which are often adjusted to experimental data with poor extrapolation capabilities, and enables the simulation of an entire hydrogen underground storage with its advection-dominated injection and extraction processes and diffusive processes through the caprock. Furthermore, the newly derived model requires less coefficients than comparable models such as GERG2008, making it easier to use and leading to a more reliable fit. To assess the influence of the derived transport properties and further model reduction on the resulting pressure field during the injection and extraction process, we consider a numerical test case of subsurface hydrogen storage using the CFD simulation software DUMUX (Flemisch et al., 2011). We inject and extract hydrogen under the caprock of a homogeneous dome-shaped aquifer and investigate the differences in the resulting pressure field using equilibrium and transport properties from PC-SAFT and entropy scaling versus alternative simplified approaches from literature. We show that density and viscosity of hydrogen have the biggest impact during the advection-driven injection and extraction cycles. Furthermore, we approximate the diffusive flux of hydrogen through the porous caprock and analyze the behavior of the thermodynamic models with changes in pressure and temperature. The comparison shows that sophisticated thermodynamic models such as the one developed in this study to determine properties of hydrogen at storage conditions are critical to accurately predict the hydrodynamic behavior of the subsurface storage system, while keeping the number of coefficients at a minimum.
This work is structured as follows. In Section 2 we develop static and dynamic transport properties of hydrogen, water and their mixture using the PC-SAFT EoS in combination with an entropy scaling approach. Next, in Section 3 we introduce the fluid-mechanical model that describes multiphase flow in porous media and will 3 of 19 be used to simulate hydrogen storage. In Section 4 we compare the developed static and dynamic properties to alternative simplified models from the literature and GERG2008, and analyze the impact of our newly developed model and literature models on simulation results of hydrogen storage in a porous aquifer.

Static and Dynamic Properties of Hydrogen, Water, and Their Mixture
In this section, we introduce the PC-SAFT EoS and entropy scaling to calculate the required static and dynamic properties for the fluid mechanical model. PC-SAFT provides the densities of pure hydrogen and the binary phase equilibrium of the hydrogen/water mixture. The dynamic properties are calculated by an entropy scaling approach which delivers the self-diffusion coefficients and viscosities for both hydrogen and water.

PC-SAFT Equation of State
The SAFT EoS are a class of modern EoS based on statistical mechanical theories, which assume an underlying molecular model with a (coarse-grained) representation of the inter-and intramolecular interactions. This statistical mechanical basis provides significant advantages over empirical models for the molar excess Gibbs energy g E or cubic EoS which are adjusted to experimental data of limited range and usually extrapolate poorly to conditions outside the scope of the adjusted data set.
The SAFT approach is based on Wertheim's first-order thermodynamic perturbation theory (TPT I) (Werthei m, 1984a(Werthei m, , 1984b(Werthei m, , 1986a(Werthei m, , 1986b and was introduced by Chapman, Jackson, Gubbins, and co-workers (W. G. Chapman et al., 1988Chapman et al., , 1989Jackson et al., 1988). Several SAFT models were proposed, such as SAFT-VR for square-well potentials of variable range (Alejandro & Jackson, 1998;Gil-Villegas et al., 1997), soft-SAFT (Blas & Vega, 1998;Felipe, 1997), SAFT-γ-Mie (Lymperiadis et al., 2007;Papaioannou et al., 2014) and PC-SAFT (Gross & Sadowski, 2000, 2001, 2002. The SAFT EoS are formulated in terms of the Helmholtz energy ({ }, , ) in the canonical variables, number of molecules { } of each individual component, volume V and temperature T. In this work, we apply the PC-SAFT EoS (Gross, 2005;Gross & Sadowski, 2001;Gross & Vrabec, 2006). The molecular model of the PC-SAFT EoS described molecules as chains of tangentially bound spherical segments. The chains are fully flexible, thus, neglecting intramolecular angle-constraints between segments of the same chain, except for the ideal gas contribution. The intermolecular interactions, that is, interactions between segments of different molecular chains, are decomposed into various contributions each leading to individual Helmholtz energy contributions, as =̃i g +̃h s +̃h c +̃d isp +̃p olar +̃a ssoc where is the reduced Helmholtz energy ̃ = ∕ ∑ and k B is the Boltzmann constant. These contributions are illustrated in Figure 1. The superscripts denote hard-sphere repulsion (hs), chain formation (hc), dispersive van-der-Waals attraction (disp), associating (hydrogen-bonding) interactions (assoc), and multipolar contributions (polar), respectively. Dispersive attraction, in this context, refers to the non-directional attraction between individual segments of different molecular chains. The ideal gas contributions ig is exactly known from statistical mechanics (Hansen & McDonald, 1990). The hard-sphere repulsion is described by the EoS of Boublik and Mansoori (Boublík, 1970;Mansoori et al., 1971). The aforementioned TPT I fuses individual hardsphere segments to form tangentially bound chains and was formulated by W. G. Chapman et al. (1988) based on the work of Wertheim (1984aWertheim ( , 1984bWertheim ( , 1986aWertheim ( , 1986b. The Helmholtz energy contributions disp containing the 4 of 19 dispersive interactions, developed by Gross and Sadowski (2001), is derived from a second-order Barker-Henderson perturbation theory (Barker & Henderson, 1967) extended to chain fluids (Gross & Sadowski, 2000). The first four Helmholtz energy contributions in Equation 1 are sufficient to describe non-polar non-associating fluids, like n-alkanes, mono-atomic gases and non-polar polymers. The required pure component parameters, namely, the number of segments per molecule m i , the segment size parameter σ i , and the dispersive energy parameter ɛ i , can be directly connected to physical properties of the molecular model.
Substances with directional interactions, like alcohols or water, require additional Helmholtz energy contributions. Hydrogen bonds are captured by the associating contribution (Huang & Radosz, 1991), which is also based on TPT I. This contribution introduces additional pure component parameters, namely the effective association volume and the association energy between two association sites A and B on molecules i and j, respectively. Multipolar interactions arising from dipole or quadrupole moments are separated into dipole-dipole (Gross & Vrabec, 2006), dipole-quadrupole (Vrabec & Gross, 2008) and quadrupole-quadrupole (Gross, 2005) interactions. Each contribution is based on a third-order perturbation theory and requires either the dipole moment μ i or the quadrupole moment Q i as a pure component parameter. Rosenfeld (1977Rosenfeld ( , 1999 introduced the entropy scaling principle which states that transport properties ∈ { self,0 } , made dimensionless in an appropriate way by a reference ψ ref , are solely a function of the reduced molar residual entropy * ({ } , ) . This reduces the complex temperature and density dependence of the transport coefficients to univariate relations of the residual entropy s*, only. In this work, we utilize entropy scaling to calculate self-diffusion coefficients self, 0 and shear viscosities η i for both hydrogen and water. The dimensionless transport coefficients can then be expressed as

Transport Coefficients From Entropy Scaling
where the self-diffusion coefficient is for a pure component, with self,0 = self ( = 1) and molar fraction x i of component i. The reduced entropy is defined according to Here, res ({ } , ) is the residual entropy, = ∑ is the average number of segments per chain. Dividing the entropy by ensures similar ranges of s* for different substances. The residual entropy s res is defined as the molar entropy minus the ideal gas entropy, as We note the ideal gas contribution ig ) (with p as pressure) by a factor of ln ∕ ∑ and consequently the corresponding residual quantities are different. The entropy scaling principle, as of now, has no convincing statistical mechanical basis, yet it is proven to be an effective and robust method to calculate transport coefficients of pure components and mixtures Loetgering-Lin & Gross, 2015;Loetgering-Lin et al., 2018;Zmpitas & Gross, 2021). In this work, we use nonlinear relations from the literature for the transport coefficients of water and the diffusion coefficient of hydrogen. We propose a new entropy scaling relation for the viscosity of hydrogen, as 5 of 19 The substance-specific parameters a i , b i , c i , and d i are adjusted to experimental transport properties using a least squares algorithm. The transport properties are reduced by the reference Chapman-Enskog self-diffusion coefficient and viscosity, respectively Here (M i /N A ) is the mass per molecule, with M i as the molar mass and N A as Avogadro's constant, Ω (1,1) * and Ω (2,2) * are the dimensionless collision integrals (S. Chapman & Cowling, 1990;Hirschefelder et al., 1954). Quantity Ω (1,1) * H 2 of hydrogen is set to 1, because a better correlation with experimental data was achieved compared to a density dependent collision integral . Our entropy scaling model provides predictions for the self-diffusion coefficient self,0 for molecules of type i in pure species i, whereas a Maxwell-Stefan transport diffusion model requires a transport diffusion coefficient  ( ) in the considered mixture. The relation between the self-diffusion coefficients self,0 , the self-diffusion coefficients at infinite dilution self,∞ , the self-diffusion coefficients in a binary mixture self ( ) and the Maxwell-Stefan binary diffusion coefficients  ( ) are depicted in Figure 2.
where = 0.5 ( + ) and d i are the effective molecular diameters  10.1029/2021WR030885 The liquid density pure of pure component i is calculated at the species normal boiling point, that is, at pressure p ⊖ = 1.013 25 bar and boiling temperature sat ( ⊖ ) . The exponent β is an empirical parameter set to β = 2.86. The self-diffusion coefficient self in binary mixtures is obtained from the relation of Liu et al. (2011), as The Maxwell-Stefan binary diffusion coefficients, according to Darken (1948) and Sridhar (2010), is given by with the thermodynamic factor where R = k B N A is the universal gas constant. The exponent α is an empirical parameter often set to α = 0.64.

Fluid-Mechanical Model
Multiphase flow through porous media can be described on different scales, whereby intermolecular interactions on the molecular scale, captured by thermodynamic models such as, for example, the PC-SAFT EoS, determine the macroscopic fluid properties on larger scales.
It is possible to describe flow through a pore network by a continuum consideration on the pore scale. To describe this flow, the complete Navier-Stokes equation must be solved. This can be prohibitively tedious and costly in larger domains, such as the storage aquifer considered for underground hydrogen storage. Moreover, it is unfeasible since the pore-scale geometry is not known in the required detail in the entire domain. For this reason, the REV scale is introduced as a simplification. Instead of exactly describing the distribution of the liquid phase and the gas phase as well as the grain structure, only the various phases are regarded in volume-averaged (i.e., integrated) form. Figure 3 illustrates the structure of the proposed model, from the molecular scale to the representative elementary volume (REV) scale.
A further and common system simplification on the REV scale is the assumption of local thermodynamic equilibrium between phases within a REV, which states the phases within an REV are in local mechanical, thermal and chemical equilibrium to one another, albeit gradients in field variables may be present that lead to fluxes across the REV boundaries.

Property Model
Hydrogen density H 2 Ideal gas law Water density H 2 O IAPWS (Wagner & Kretzschmar, 2008) Hydrogen viscosity H 2 Chung (Chung et al., 1988) Water viscosity H 2 O IAPWS (Cooper & Dooley, 2008) Diffusion coefficient gas phase Fuller (Fuller et al., 1966)  We consider a two-phase, two-component system consisting of water and hydrogen phase, and water and hydrogen components. Furthermore we assume a rigid and inert porous medium. The mass balance of component ∈ {H2, H2O} with the storage term, the advection and diffusion term follows as with the phase index ∈ {gas, liquid} , porosity ϕ, saturation S α , sink/source term q i , time t, molar density ρ mol,α , and mass density ρ mass,α of phase α, molar fraction x iα , and mass fraction w iα of component i in phase α and the diffusion coefficient in the porous medium pm . The Darcy velocity v α can be expressed by the extended Darcy's law for multiphase flow (momentum balance) with the relative permeability k rα which is an empirical function that depends on saturation, for example, Brooks-Corey (Brooks & Corey, 1964), intrinsic permeability K, phase pressure p α , viscosity η α , and gravitation g. The viscosities of the vapor and liquid phase η α are assumed to be the pure component viscosities of hydrogen and water, respectively.
The diffusion coefficient on the REV scale pm in a porous medium can be expressed by with the tortuosity τ α and the self-diffusion coefficient self , ( ) of component i in phase α, given by Equation 13. The tortuosity τ α is a degree of twisting of the flow paths through a porous medium and can be calculated by the approach of Millington and Quirk (1961), as We implement two versions of the fluid mechanical model: one uses empirical approaches for fluid properties of water, hydrogen and their mixture, as in Table 1. The other determines fluid properties from entropy scaling and the PC-SAFT EoS according to the procedure described in the previous section.
One of the phase pressures is eliminated using the capillary pressure p c = p gas − p liquid , which is an empirical function that depends on saturation, for example, Brooks and Corey (1964). One saturation is eliminated by making use of the closure equation S liquid − S gas = 1. The finite-volume method is used as the spatial discretization, where we apply a two-point flux approximation (Helmig, 1997). The discretization in time is handled by an implicit Euler scheme and we apply the Newton method to solve the system of nonlinear partial differential equations in a robust and consistent manner.

Results and Discussion
In the following, we present the properties for hydrogen, water and their mixture as determined by the PC-SAFT EoS and entropy scaling and compare our results to experimental data and empirical correlations. Next, we simulate hydrogen storage in a dome-shaped aquifer with our fluid-mechanical model, using the newly developed properties for hydrogen, water and their mixture, as well as empirical correlations. Finally, we analytically calculate the diffusive flux of hydrogen through the caprock, depending on pressure and temperature changes. We compare the impact of the various static and dynamic properties on our results and identify best-suited thermodynamic models for the simulation of underground hydrogen storage.

Pure Component PC-SAFT Parameters
The pure component parameters of PC-SAFT are commonly adjusted to experimental data of vapor pressures and liquid densities. The critical temperature of hydrogen ( = 33.145 K) , however, is significantly lower than for most other species. At this temperature, hydrogen does not behave like a classical fluid, which is evident when we compare the thermal de Broglie wavelength ΛH 2 ( ) at the critical temperature to the mean distance L c between hydrogen molecules at the critical point where h is Planck's constant, H 2 is the mass of a hydrogen atom, and ρ c is the number density of hydrogen at the critical point. A fluid i behaves classically if Λ i ≪ L i , whereas for hydrogen at critical conditions (and below) the thermal de Broglie wavelength is of the same order of magnitude as the mean intermolecular distance. As a consequence hydrogen does not follow Maxwell-Boltzmann statistics at low temperatures, instead it has to be treated as a quantum fluid. Employing a classical equation of state, like the PC-SAFT equation of state, for the modeling of the vapor-liquid equilibrium of hydrogen would not yield accurate results. This can be remedied by introducing quantum corrections to classical interaction potentials. The theory was first introduced by Wigner and Kirkwood (Kirkwood, 1933;Wigner, 1997) and subsequently re-derived by Feynmann and Hibbs (Feynman  et al., 2010) using the path-integral method. The molecules experiencing quantum effects can then still be described by classical statistical mechanics, as long as they interact through effective, temperature-dependent potentials. This principle was applied by Wilhelmsen and co-workers who extended the SAFT equation of state for variable-range MIE potentials (SAFT-VR Mie) to also include Feynman-Hibbs corrections. This quantum corrected (SAFT-VRQ Mie) EoS is then applied to the calculation of vapor-liquid coexistence lines of hydrogen, helium, and neon (Aasen et al., 2019) and their mixtures (Aasen et al., 2020).
In this work, with subsurface storage of hydrogen in mind, we consider temperatures well above 200 K and we are not troubled with an accurate description of the vapor-liquid coexistence behavior. We chose to adjust the PC-SAFT parameter to supercritical pVT data with pressures between 1,000 and 3,500 bar and temperatures between 280 and 400 K from the NIST database (Lemmon & Friend, 2020), as experimental hydrogen self-diffusion coefficients are also available at this pressure and temperature ranges. Hydrogen is assumed to be a spherically symmetric molecule, therefore, we set the number of segments H 2 to unity. The resulting PC-SAFT parameters for hydrogen are presented in Table 2. Of course, this parameter set is not suited for the calculation of vapor-liquid coexistence properties and significantly overestimates the critical point. Figure 4 shows a comparison for the hydrogen density for varying temperatures and pressure typically encountered in the subsurface storage of hydrogen. We observe excellent agreement between the calculated densities from PC-SAFT, GERG2008, and reference data from the NIST database over the whole pressure and temperature range. This illustrates the strong extrapolation of the PC-SAFT EoS to conditions not included in the adjusted data set. The GERG2008 EoS is a highly parametrized empirical EoS which was specifically developed for the calculation of thermodynamic properties of natural gases, similar gases, and other mixtures (Kunz & Wagner, 2012). The main areas of applications are processing, transportation through pipelines or by shipping, storage and liquefaction of natural gas, and processes to separate gas components.
The PC-SAFT EoS, however, is a general purpose EoS with applications ranging from simple fluids like argon to associating fluids, like alcohols, and chain-like polymers. The statistical mechanical basis of PC-SAFT significantly reduces the number of required parameter compared to empirical EoS while retaining the high degree of accuracy required for industrial applications. Therefore, the PC-SAFT EoS is our model of choice for the description of hydrogen, water and their mixture. The ideal gas law is also shown to be a good approximation for the pVT behavior of hydrogen for pressures below 40 bar.
In contrast to hydrogen, the accurate prediction of thermodynamic properties of water is a challenging subject, even for modern EoS. The parameter set for water used in this work was proposed by Rehner and Gross (2020) and we chose the non-polar water model in combination with a 2B association scheme. They included surface tensions from predictive density gradient theory in addition to vapor pressures and liquid densities in a multi-objective optimization to ensure that the resulting parameter set is suitable for the calculation of interfacial properties.
In Figure 5, we compare liquid densities from the PC-SAFT EoS, the empirical IAPWS model (Wagner & Kretzschmar, 2008), and the NIST database (Lemmon & Friend, 2020). For most liquids, the density monotonically increases for decreasing temperatures. Liquid water, however, reaches its maximum density at 4°C and expands again at lower temperature, to the point where solid water has a lower density than liquid water. The density anomaly of water is notoriously difficult to reproduce for EoS based on perturbation theories, as water forms higher-order tetrahedral structures between multiple water molecules. These higher-order structures significantly alter  the fluid structure and cannot be represented by a simple hydrogen-bonding (association) models as employed in PC-SAFT. This is evident from the observed deviation between liquid densities from the PC-SAFT EoS and NIST data. The IAPWS relation is explicitly derived and highly parametrized to mimic the density behavior of water. It should be noted that the PC-SAFT parameters of water are not only used to calculate densities, but also for phase equilibria and for the residual entropy in the entropy scaling approach.  Table 3 Parameters of the Entropy Scaling Nonlinear Correlation Functions Figure 10. Dimensionless self-diffusion coefficients and viscosities from entropy scaling and experiments.

Binary Hydrogen/Water Phase Equilibrium
The phase behavior of the binary hydrogen/water mixture requires a binary k ij parameter which acts on the strength of the cross-dispersive interactions ɛ ij between molecular segments of different components Water forms hydration shells around individual gas molecules which significantly alter the solubility of the dissolved gas. The underlying hydrogen-bonding network is also responsible for the density anomaly and can not be captured explicitly by the EoS. Instead a k ij parameter has to be adjusted to binary experimental liquid-vapor equilibrium data. The intricate temperature behavior of the molar fraction H 2 in the liquid phase is presented in Figure 6, exhibiting a parabolic shape with a minimum at 330 K. This curvature requires a temperature-dependent binary k ij parameter. We adjust the parameter individually at different temperatures to vapor-liquid equilibrium data from the Dortmund Data Base (Dortmund Data Bank, 2020) and use a third-order polynomial to interpolate between different temperatures shown in Figure 7, resulting in the following relation The temperature-dependent correlation should not be used outside the temperature range 280-500 K, because the extrapolation of the cubic polyno-  mial is not well-behaved We also include the empirical correlation by Fernández-Prini et al. (2003) in Figure 8 which contains more than 10 parameters adjusted to experiments. The excellent agreement highlights the strength of the statistical mechanical basis of the PC-SAFT EoS, as we are able to accurately calculate the water-hydrogen phase equilibrium with a lower number of adjusted parameters (Figures 8 and 9).

Entropy Scaling
In this section, we present the entropy-scaling models for the shear viscosity and for the self-diffusion coefficients of hydrogen and water. Each data point is assigned a dimensionless reduced entropy, according to Equation 4, calculated from the experimental conditions, that is, temperature and pressure/ density. The experimental transport coefficient is then made dimensionless by their Chapman-Enskog counterpart from Equation 10. The parameters of the non-linear correlation functions (Equations 6-9) are adjusted to these data points using a least squares algorithm. The resulting parameters are listed in Table 3 and Figure 10 shows the correlations of the individual transport properties for hydrogen and water. The self-diffusion coefficients of hydrogen in Figure 10a shows scattering of the experimental data, instead of collapsing onto a single line. This could indicate a possible limitation of the Entropy Scaling approach or hint at inaccuracies of the reported experimental data. The newly proposed correlation for the viscosity of hydrogen yields excellent agreement to experiment, as is evident from Figure 10b. We expect to calculate accurate viscosities of hydrogen over wide temperature and pressure windows. The strong non-ideal fluid behavior of water is reflected by the high absolute values of the residual entropy s*. The steep slope of the curves appear for the supercooled fluid region in Figures 10c and 10d.

Binary Self-and Maxwell-Stefan Diffusion Coefficients
We can now calculate the self-diffusion coefficients self H 2 ( ) and self H 2 O ( ) and the Maxwell-Stefan diffusion coefficient H 2 H 2 O from the pure self-diffusion coefficients, as suggested from the correlations by Liu (Equation 13) and Darkin (Equation 14). Figure 11 shows the relation between the diffusion coefficients, similar to Figure 2, at a pressure of 20 bar and different temperatures between 300 and 450 K. The diffusion coefficients are only defined for compositions outside the vapor-liquid coexistence lines (for given pressure and temperatures). All diffusion coefficients collapse onto a single line, without any discernable difference between them, and are strongly dependent on the composition, as even trace amounts of water significantly reduce the diffusion rate.

Comparison of Models for Transport Coefficients
In this section, we compare transport coefficients from entropy scaling to empirical correlations and experimental data in the relevant temperature and pressure range. The self-diffusion coefficient of hydrogen in Figure 12 is systematically underestimated by the Fuller method for the entire range of temperatures and pressures when compared to results from entropy scaling. Figure 13 shows the self-diffusion coefficients of water.
While the diffusion coefficients after the temperature dependent linear approximation of experimental data (Engineering ToolBox, 2021) in DUMUX do have a linear and flatter dependence on the temperature, the entropy scaling approach delivers a more logarithmic and steep curve and correlates more accurately with the experimental data.
Comparing the viscosities of water in Figure 15, the IAPWS relation is able to reproduce the experimental data very precisely, while the entropy scaling  approach differs from the experimental data points at such low temperatures near the density anomaly of water, where the fluid behavior of water is strong non-ideal.
For the viscosity of hydrogen in Figure 14 it is remarkable that the temperature-dependence from the model of Chung et al. (1988) and from the entropy-scaling model are very similar, but the model of Chung et al. results in much lower viscosities over the whole considered temperature range. Furthermore, the Chung-model does not depend on pressure, whereas the entropy scaling approach is able to represent the pressure-dependence of viscosity and fits the experimental data points of hydrogen rather precisely.

Numerical Case Study of a Dome-Shaped Aquifer
We simulate hydrogen storage in a porous aquifer and assess the qualitative influence of static and dynamic properties like density, phase equilibrium, viscosity and diffusion coefficients on the saturation and pressure field. Furthermore, we compare simulation results with fluid properties determined by entropy scaling and the PC-SAFT EoS versus simplified, empirical approaches from literature, and analyze the impact of the applied thermodynamic model.
We inject and extract hydrogen in a radially symmetric, dome-shaped aquifer under a caprock. The aquifer is characterized by an isotropic, intrinsic permeability of k = 1×10 −12 m 2 and a porosity of Φ = 0.2. Only a slice of the total storage is considered, which has the shape of a piece of cake and is shown in Figure 16. The slice of the domain covers an angle of 4.5°, has a radius of 1,200 m, an arch-length of 100 m and a height of 50 m. The temperature and pressure conditions are based on previous work of Pfeiffer (Pfeiffer et al., 2017) with about 25°C and 20-70 bar depending on injection or extraction process. On the tip of the slice a Neumann boundary condition is defined to inject and extract pure hydrogen into the storage. The injection/extraction area of the well is represented by the green rectangular on the left side in Figure 16. The other boundaries of the domain are implemented as no-flow boundary conditions, except for the right side of the domain, where the pressure is constant set at 40 bar through a Dirichlet boundary condition to model a horizontal open storage. The Brooks and Corey parameters are assumed as λ = 0.6 and entry pressure p d = 1,189 Pa. At the start of the simulation, the domain is already filled with the gas phase (hydrogen and evaporated water) at equilibrium in most parts of the domain to simulate a hydrogen storage in idle state. This starting condition at t = 0 s can also be seen in Figure 16. The horizontal gas/water interface is located at a depth of 50 m below the highest point of the well. Hydrogen is injected with a mass flow of 3 kg/s until a total injected mass of 30 t of hydrogen is reached and afterward extracted with a rate of 1.25 kg/s. Between the injection and extraction process a relaxation time of 4.5 hr is implemented. Under these thermodynamic conditions the liquid and vapor phase are assumed to consist mainly of pure water and hydrogen, respectively. Therefore, the transport properties of the binary mixture in each phase are assumed to be the pure component properties. We perform and evaluate simulations of the simplified subsurface hydrogen storage using the correlated transport properties from the PC-SAFT EoS and entropy scaling approach. For that purpose, the fluid and transport properties determined by our thermodynamic model were tabulated by pressure and temperature, read in by DUMUX and were linearly interpolated during the simulations. We note, that the Figure 15. Viscosity of water from entropy scaling and the IAPWS relation (Wagner & Kretzschmar, 2008). resolution of the tables was high enough, so that further refinement of the tables lead to no significant change in the simulation results.
Throughout injection, extraction and idle phases, a significant change in saturation could not be noticed. Instead, hydrogen was mostly stored by compressing the hydrogen rich gas phase. Therefore, in the following, the differences in pressure of the gas phase will be used to point out deviations between the simulation scenarios with different combinations of the thermodynamic relationships. Simulations using the newly developed thermodynamic properties of hydrogen, water and their mixture are summarized under the term "PC-SAFT," whereas simulations using the relationships of comparison (Table 1) are summarized under the term "DUMUX." The differences between simulations using PC-SAFT versus the thermodynamic models implemented in DUMUX becomes most noticeable close to the injection/extraction well, where the highest flow rates are realized (see Figure 17). Further away from the well, the influence of the thermodynamic model becomes less significant.
The pressure of the gas phase in the well over time during injection phase, idle phase, and extraction phase is shown in Figure 18 for two temperatures. Simulations using properties from PC-SAFT show a larger change in pressure over time than simulations using properties from literature, both during injection and extraction phase. We attribute this to the higher flow resistance due to a larger hydrogen viscosity according to PC-SAFT. Steady state is not reached during injection or extraction phase, which is why the difference between the simulations using different thermodynamic models grows with time during these phases. Consequently, the largest pressure difference between the simulations can be observed at the end of the extraction phase. The difference during idle phase is negligible, since the pressure in the well equals the equilibrium pressure close to the well location inside the domain, which is given by the boundary conditions. The absolute gas pressure in the well during injection or extraction phase increases with increasing temperature, and so does the difference in the pressures predicted by the simulations. Since the density of hydrogen is generally smaller for higher temperatures in both thermodynamic models, and the viscosity of hydrogen increases with temperature in both thermodynamic models, higher temperatures mean that larger gas volumes are injected or extracted against a higher resistance. As a result, accurate transport properties are of highest importance here, and underestimating the hydrogen viscosity causes larger discrepancies.
We investigate the maximum pressure difference in the well during injection and extraction more closely. Depending on the thermodynamic model used to determine the properties viscosity of water and hydrogen, density of water and hydrogen, as well as phase equilibrium and diffusion coefficients, the maximum pressure difference changes, as can be seen in Figure 19. At maximum, the differences in pressure between the models amount to over 5 bar during the injection phase and 8 bar during the extraction phase. In an advection-driven process such as the storage of hydrogen, we expect viscosity and density to be the decisive parameters. Figure 19 confirms that the viscosity of hydrogen is a key quantity, followed by the density of hydrogen, both of which show a major influence on the pressure in the well. Viscosity and density of water play a minor role, even though the thermodynamic models determine different values for those properties. However, since the water phase barely moves in our simulations, accurate transport properties for water are not required. As expected for such an advection-dominated system, phase equilibrium and diffusion coefficients have a negligible impact on simulation results as well.
In summary our results underline the importance of reliable thermodynamic models to determine transport properties like viscosity and density, in particular for hydrogen. For higher temperatures and closer to the well, accurate transport properties of hydrogen are most significant. We note that for  higher flow rates than used in this study, the interface between gas and water in the storage may move (Hagemann et al., 2015). In this case, viscosity and density of water may also play a role.

Analytical Investigation of Caprock Diffusion
In the following, we study diffusion of hydrogen through a caprock that is fully saturated with pure water, as described by Fick's Law In such a diffusion-dominated problem, differences in diffusive fluxes as calculated by different models directly depend on the determination of diffusion coefficients and phase equilibria. For the evaluation of Equation 24, we used a 0D-model, where the gradient of molar fractions is computed by a difference quotient of the molar fractions directly above and below the caprock. The molar fraction of hydrogen above the caprock was assumed to be zero, pertaining to an aquifer with horizontal flow transporting leaking hydrogen away and keeping the gradient at a maximum. The molar fraction within our storage aquifer is calculated in dependence of the temperature T Figure 19. Pressure differences of the gas phase during injection and extraction of the simulations using different combinations of transport properties from DUMUX and Perturbed-Chain Statistical Associating Fluid Theory (PC-SAFT). The pressure differences are calculated in reference to the simulations using properties from entropy scaling and PC-SAFT equation of state of this work. and the pressure p of the storage aquifer. Further, a porosity of ϕ = 0.01, a water saturation of H 2 O = 1 and a thickness of the caprock (diffusion length) of Δz = 50 m are assumed. Based on literature, this represents a relatively thin caprock. Figure 20 shows the resulting diffusive flux through the caprock, using the correlated PC-SAFT EoS and the models implemented in DUMUX, over temperature and for different pressures. The pressure dependence is mainly dominated by the pressure dependence of the phase equilibrium (see Figure 6 for reference), whereas the diffusion coefficients are mostly dominated by temperature. The different slopes of the diffusion coefficients in the liquid phase with regard to the temperature dependence (see Figure 13 for reference) clearly determine the slope of the diffusive fluxes, which leads to lower fluxes at low temperatures and higher fluxes at high temperatures using the correlations from PC-SAFT compared to the alternative models from literature.
These results show the relevance of robust and accurate thermodynamic models in diffusion-dominated scenarios. Especially diffusion coefficients with their strong temperature dependence have a major influence on the diffusive fluxes. The PC-SAFT EoS in combination with the entropy scaling approach deliver a clear improvement over the empirical models here.

Conclusion
Thermodynamic models for transport coefficients (shear viscosity and diffusion coefficients) as required for the storage of hydrogen in the subsurface were developed and compared to empirical approaches. For this purpose, the PC-SAFT EoS was correlated to experimental data to derive thermodynamic properties like densities, entropies, as well as phase equilibria of hydrogen and water. The pVT behavior of hydrogen modeled by PC-SAFT was compared to the highly parametrized empirical GERG2008 EoS and reference data from NIST. We observe excellent agreement between the three. Viscosity and diffusion coefficients were determined using an entropy scaling approach. The storage of hydrogen in a porous aquifer was studied and results obtained from the proposed model for thermodynamic and transport properties (PC-SAFT) were compared to results of a model that employs an empirical description of thermodynamic and transport coefficients. Second, we studied the diffusive fluxes of hydrogen through a caprock.
The PC-SAFT EoS in combination with the entropy scaling approach has proven to be a powerful tool for predicting thermodynamic and transport properties in a large pressure and temperature range. A significant advantage is in the consistency of the model and in a low number of adjustable parameters, so that many fluid properties can be derived from a single model. In contrast to this are the numerous empirical relationships that can be found in the literature, each with multiple fitting parameters, which are only able to describe a single thermodynamic property, resulting in a multitude of equations and model assumptions. Only the density anomaly of water cannot be represented satisfactorily by the PC-SAFT EoS. Here it is recommended to use a more comprehensively parametrized and specialized model like the one presented by the IAPWS.
According to the CFD simulation scenarios, density and viscosity of hydrogen have the biggest impact on the pressure field in our advection-dominated simulation scenario. Diffusion coefficients and phase equilibrium play a negligible role in this case. A sophisticated thermodynamic model to determine density and viscosity of hydrogen in particular, as the one presented in this work, is therefore central to an accurate prediction of the hydrodynamic behavior of the storage. In contrast to that, diffusion coefficients play a vital role with diffusion through the caprock. Our newly developed thermodynamic model leads to a much improved approximation of the diffusive flux, especially with changes in temperature and pressure.
Considering the model for the representation of the dependence of the diffusion coefficients on the composition, no verification of the model with experimental data was possible here, because data was not available. In order to verify the model, molecular simulations could be performed, which is a need for further research.