Hydrogen and air storage in salt caverns: a thermodynamic model for phase equilibrium calculations

. When storing gas in a salt cavern, it occupies most of the excavated volume, but the lower part of the cavern inevitably contains residual brine, in contact with the gas. The design of hydrogen and compressed air storage in salt caverns requires to have a thermodynamic model able to accurately predict both phase properties such as densities, and phase equilibrium (gas solubility and water content of the vapour phase). This work proposes a parameterization of the e-PPC-SAFT equation of state in this context. Experimental data of pure components and mixtures of light gas + pure water and light gas + salted water are reviewed and used to ﬁ t pure component parameters for hydrogen, nitrogen, oxygen, and the brine, and binary interaction parameters between H 2 , O 2 , N 2 + water and H 2 , O 2 , N 2 + ions (Na þ and Cl (cid:2) ), for temperature ranging from 273 to 473 K and salinities up to NaCl saturation (6 mol/kg). The model developed delivers good accuracy in reproducing data: the average deviation between experiments and calculated data is between 3% and 9% for gas solubility in saturated brine. More interestingly, the model has been validated on its capability to predict data not included in the parameterization database, including the composition of the vapor phase, and its extension to a mixture, such as air. Finally, it has been used in a case study of Compressed Air Energy Storage (CAES) to evaluate the water content of the gas produced during injection-withdrawing cycles.


Introduction
In the context of the energy transition, renewable energies like solar and wind are promising solutions to replace fossil fuels [1]. However, such kind of energy production is intermittent. To overcome this issue, these renewable energies can be stored during high-energy production phases, and used when needed [2]. Among all storage technologies, underground storage is well known and provides several advantages such as safety, longevity, and high storage capacity [3]. More specifically, storage in salt caverns is very promising for light gas storage due to very low permeability. The electricity produced from renewables must first be turned into an energy carrier that can be stored, like reactive hydrogen [4], or be used for Compressed Air Energy Storage (CAES) [5].
Salt caverns are developed by injecting water into a geological salt deposit, thus leaching the salt and then withdrawing the brine. Once the desired geometrical volume is obtained, the first gas fill is made by injecting gas and withdrawing the brine. At the end of this so-called "debrining" stage, an amount of brine typically between 1.000 m 3 and 10.000 m 3 stay at the bottom of the cavern. The cavern is then put into operation. This presence of the residual brine impacts the gas storage operations. First, it can solubilize a part of the injected gas. Second, while it is injected dry, the gas produced contains a higher amount of water, making it necessary to use dehydration surface processes. In the case of compressed air, the withdrawn humid air causes aggressive corrosion conditions for the well and surface piping and equipment. It is thus of primary importance to have a thermodynamic model able to predict phase equilibrium between brine and gas in large conditions of temperature and pressure, and up to saturation in terms of salinity.
Concerning the phase equilibrium model between hydrogen and brines, only a few thermodynamic models have been published. Lopez-Lazaro et al. [6] proposed an approach based on the Soreide and Whitson equation of state [7] calibrated on a few existing experimental and molecular simulation data. Recently, Chabab et al. [8] published new experimental data in temperature and salinity ranges in better accordance with underground storage, and they used both Soreide and Whitson and an electrolyte version of the Cubic Plus Association equation of state [9] to successfully model their data. Li et al. [10] proposed a heterogeneous approach based on the Pitzer activity model and Henry constants for the aqueous phase and a purehydrogen equation of state for the vapour phase. This model assumes that the vapour phase behaves ideally, which leads to non-negligible deviations at high pressure, as demonstrated in a previous study [11]. Recently, Zhu et al. [12] also proposed a heterogeneous approach using a Pitzer model for the liquid phase. The fugacity of hydrogen in the vapour phase is calculated with a pure-hydrogen model. A correction is brought to evaluate the water content in this phase. Although accurate on data used for model calibration, these heterogeneous approaches of Li et al. and Zhu et al. are rather devoted to low-pressure conditions and do not allow an accurate calculation of the vapour phase composition. Another heterogeneous approach was recently proposed by Torín-Ollarves and Trusler, based on the Setchenow coefficient and Henry constant correlations [13]. This model fits well with experimental data but remains quite empirical with a large number of parameters. Again, this model does not allow an accurate calculation of the vapour phase composition. Roa Pinto et al. proposed to use an electrolyte version of the PPC-SAFT equation of state and show good agreement with previous models, but also a good capacity in extrapolation when compared to data not used in the calibration database [11]. However, this model exhibits non-negligible deviations in the compressibility factor of pure hydrogen at high pressure, opening the way to improvements to get better accuracy of gas volume produced during storage withdrawal.
Concerning phase equilibrium between air and saturated brines, no specific model has been found in the literature probably due to the lack of experimental data, although molecular simulation is more and more used as a complementary tool to experiments in this field [14]. Nevertheless, we can mention the modeling work recently proposed by Chabab et al. on N 2 + brine and O 2 + brine systems with the Soreide and Whitson equation of state [15]. This model has been validated on both gas solubility and water content in the vapor phase at a large range of temperatures and salinities. Coupled with their previous parameterization for other gases like H 2 , CO 2, and CH 4 , this is probably today the most complete thermodynamic model to deal with gas storage in salt caverns.
In this work, we propose a thermodynamic model able to predict both phase equilibrium: H 2 + brine and Air + brine. Air is considered a mixture of 80% mol. of nitrogen and 20% mol. of oxygen. Thus, binary systems N 2 + brine and O 2 + brine must be studied first, before evaluating the model in prediction on the multi-component mixture air + brine. Brine is assumed to be made of water and NaCl since this salt is the most abundant in salt caverns. The model proposed is a revision of the initial model of Roa Pinto et al. [11]. It is based on the e-PPC-SAFT equation of state: due to its strong physical background, this model has a good capability of extrapolation, which is an important requirement in the present case of scarce experimental data. This model will be parameterized to reproduce not only gas solubility in brine but also the composition of the vapor phase. The capability of the model to predict accurately the water content of the vapor phase is a key benefit compared to existing modeling works based on a heterogeneous approach [10,12,13]. Furthermore, the accuracy of a SAFT-based thermodynamic model, which includes specific terms to deal with association and electrolyte interactions, is expected to be at least the same as a cubic-based model such as the Soreide and Whitson equation of state, and to have better predictive behavior.
Another key benefit of the model proposed in this work is to be based on an ion-parameterization approach, and not a salt-parameterization approach like with the Soreide and Whitson approach. Consequently, such a model allows calculating specifically the liquid fugacity of each ion individually, which can be of prior importance if this model has to be further used in a reactive equilibria scheme where ions are involved, for instance in typical fluid-rock interaction calculations.
In the first part of this article, the e-PPC-SAFT framework will be recalled, with a specific focus on the model parameters. In the second part, a review of available experimental data is proposed. The systems investigated cover pure components (hydrogen, nitrogen, oxygen), brine, and mixtures of light gases with pure water and brine. Then, in the third part, the model parameterization is carried out for both pure components and mixtures. In the fourth part, the model is validated on its capability to predict data and properties not used in the parameterization database, and on its capability to compute multicomponent mixtures such as air only from binary considerations. Finally, an application of this model will be presented to calculate gas moisture during typical Compressed Air Energy Storage injection-withdrawal cycles.

The e-PPC-SAFT model
The e-PPC-SAFT equation of state [16] is used to model the different forces acting in an electrolytic multiphasic system. It gives the reduced residual Helmholtz energy A res as a sum of different contributions, each contribution representing a specific interacting force between particles: where R is the ideal gas constant and T the temperature. The first four contributions (hard-sphere, chain, dispersion, and association) derive directly from the original PC-SAFT version of Gross and Sadowski [17]. The hard sphere term describes the repulsive interactions between hard spheres. Its expression is the following: The Author(s): Science and Technology for Energy Transition 78, 10 (2023) with where, q is the molar density and total volume of the fluid, N A the Avogadro's number (6.02214 Â 10 23 mol À1 ), r HS i is the hard sphere diameter, e i is the dispersive energy parameter, k is the softness parameter which is equal to 0.12 (except for water in the model used in this work: k i = 0.203) and x i is the molar fraction of component i.
The chain term is used to model the chain formation between segments. This contribution is given by: where m i is the chain length parameter and g hs ij is the radial distribution function given by Boublik [18] and Mansoori et al. [19]: and The dispersive contribution describes the short-range Van der Waals attraction between segments. The corresponding expression is: with where x ¼ r r is the reduced radial distance around a segment and Z hc is the compressibility factor of hard chain.
Mixing parameters m 2 e r 3 and m 2 e 2 r 3 are given by: where e ij and r HS ij are evaluated using the Lorentz-Berthelot mixing rules: where k ij is an empirical interaction binary parameter. The association term describes the strong short-range attractions due to association forces between association sites. It derives from the Wertheim theory [20] and is given by [21]: with and X Ai is the fraction of association sites non-bonded to site A of molecule i. M i is the number of association sites in molecule i, and Á Ai Bj the association force between site A of molecule i and site B of molecule j. It is characterized by two binary parameters: the association energy (e Ai Bj ) and is the association volume (k A i B j ), determined from the CR-1 mixing rule: where w ij and u ij are empirical interaction binary parameters. From this original version of the PC-SAFT equation of state, Nguyen-Huynh et al. [22] added a specific polar term from the Jog and Chapman theory [23] to include interactions between dipoles and quadrupoles. The expression of this term is the following: where l i is the dipole moment of the molecule and x pi its polar fraction. I 2;ij and I 3;ijk are correlations whose expressions can be found in [23]. The quadripolar moment Q i of the molecule and its quadrupolar fraction x Qi can also be taken into account, following equations given in [23]. The addition of a Non-Additive Hard-Sphere (NAHS) term to the PPC-SAFT equation of state was proposed by Trinh et al. to improve prediction with hydrogen-containing mixtures [24] by correcting the additivity assumption of the hard spheres. The expression of this term is: where g SWS ij is the radial distribution function of squarewell potential [25] and d ij the cross-diameter calculated by: where l ij is an empirical binary interaction parameter. The last two terms of equation (1) are specific to electrolytic systems. For long-range electrolytic interactions, the use of MSA (Mean Spherical Approximation) is used as suggested by many authors [11,16,[26][27][28][29][30][31][32], although an approach based on Debye-Hückel theory is expected to give a similar performance [33]. Its expression is given by: where e r is the dielectric constant of the solution, e o the dielectric constant in the vacuum (8.85418 Â 10 À12 C mol À1 ), n i the number of moles of ion i, e is the electron charge (1.60218 Â 10 À19 C), Z i the ion charge, r MSA i the solvated ion diameter and C the shielding parameter calculated as follows: In equation (1), the Born term is also classically used to model short-range ion solvation forces [11,16,27]. It is defined by: where r Born i is the Born diameter for each ion. This diameter is determined from experimental Gibbs energies of solvation Á s G i [34], through the equation resulting from equation (25): Both MSA and Born terms involve the dielectric constant of the aqueous solution. In this work, it is calculated using the model of the dielectric constant of pure water proposed by Schreckenberg et al. [27] to which the salt correction of Pottel is added [35]: with where N T is the total number of moles, V is the volume of the system, d v and d T are constants with a value for the water of 0.3777 dm 3 /mol and 1403 K respectively, and Finally, Table 1 provides a summary of what parameters are used for each term of the e-PPC-SAFT model. Knowing the residual Helmholtz energy from equation (1), the fugacity coefficient u p i of each component i in a given phase p can be calculated from: where V p is the volume of the phase p, n i is the number of mole of component i, P is the pressure and n is the total mole number in phase p. For all component i, phase equilibrium is calculated by equalizing fugacities of i in liquid and vapor phases, under mass balance constraint: where x and y stand for the liquid and vapor molar fraction, respectively, and " x and " y the vectors of compositions of these phases.
For light gases showing a very low solubility in the aqueous phase, it is common to express solubility in terms of the Henry constant. It is defined as: where f i is the fugacity of i in the liquid or vapor phase. From an equation of state, it can be calculated by: where P r S is the saturation pressure of the solvent, and the superscript 1 stands for infinite dilution of i in the solvent. It can be noticed that Henry's constants computed with equation (32) are for real mixtures and not salt-free mixtures.
The computation of the Henry constant is carried out in two steps. First, a "bubble point" flash of the real solvent (i.e. with mole fraction including ions) is carried out. Then, the fugacity coefficient at infinite dilution for the gas is calculated in this real solvent. Other equilibria computed in this work (more specifically to plot pressure-composition and temperature-composition diagrams) are performed using a diphasic TP flash, always using the mole fractions of the real mixture, including ions. This flash algorithm equals the fugacity of all species (including ions) between liquid and vapor phases with respect to mass balance. It is worth noticing that the Born term used in the SAFT model makes the liquid fugacity coefficient of ions extremely low, resulting thus in a vapor composition very close to zero. In the context of underground gas storage, it is of primary importance for economic evaluations to calculate with very good accuracy the volume of the gas produced. Consequently, the parameterization of the light gases involved in this study (hydrogen, nitrogen, and oxygen) is carried out to reproduce experimental compressibility factors Z in temperature and pressure ranges of industrial interest: The compressibility factors of all these pure components are collected from the NIST database [36] for temperatures ranging from 273 to 473 K and pressures ranging from 1 to 300 bar.

Brine
The brine is considered to be an aqueous NaCl solution since halite is the most abundant salt found in salt caverns. NaCl is assumed fully dissociated in Na + and Cl À ions. One of the most important properties to be correctly modeled is the Mean Ionic Activity Coefficient (MIAC) of NaCl salt since it describes the deviation from the ideality of the aqueous solution. MIAC is defined from a single ion activity coefficient by: The e-PPC-SAFT equation of state computes fugacity coefficients u i , and not directly activity coefficients. The relationship between both for component i is: where the exponent * is related to the reference state (infinite dilution for ions). The MIAC experimental data covering a temperature range from 273 K to 473 K and salinities up to saturation are taken from [37]. Note that experimental ion activity coefficients are very often given in the molality scale, while the equation of state computes MIAC from equation (34) in the molar fraction scale. A conversion can be done using the molar fraction of the solvent: where the exponents (m) and (x) stand for molality and molar fraction scale, respectively. The composition of the gas produced during storage cycles is also of primary importance to designing surface facilities. More specifically, the water content of the vapor phase must be accurately known. Thus, special attention to our modeling is also paid to correctly reproduce the vapor pressure of the brine. Experimental data for this property on the same temperature and salinity ranges are taken from [38,39]. Finally, a set of experimental brine densities [40] are also used to parameterize the model.

Light gas + pure water systems
Concerning hydrogen, Lopez-Lazaro et al. published recently a review of experimental hydrogen solubility in pure water data [6]: around 300 data from 273 to 623 K were collected and reconciliation in terms of Henry constants. These data were fitted following the empirical correlation suggested by Trinh et al. [41]: A assoc e AB (association energy) w ij (correction in the association energy combining rule) K A i (association volume) u ij (correction in the association volume combining rule) The Author(s): Science and Technology for Energy Transition 78, 10 (2023) where T r is the reduced temperature of pure water (T r = T/T C where T C is the critical temperature of pure water: 647 K). P r s T ð Þ is the vapor pressure of the solvent, taken from the DIPPR database [42]. a, b, c are three empirical parameters given in Table 2.
A similar review is proposed in this work for the solubility of nitrogen and oxygen in pure water. For both nitrogen and oxygen, experimental data are collected and converted in terms of the Henry constant. In the literature, the solubility data can be directly reported as Henry constant, or also in terms of Kuenen coefficient S i , Ostwald coefficient L i , Bunsen coefficient a i , and molar fraction as a function of temperature and pressure ("TPx" data). For this last kind of data, we retain here only data at low pressure, since conversion into Henry constant requires to be close to infinite dilution. The reference of experimental data collected is given in Table 3. The conversion rules in Henry constant recommended by IUPAC [43] are: where T h and P h are standard temperature (273.15 K) and pressure (1.013 bar), v l;Ã S is the liquid molar volume of pure water at T, taken from the DIPPR databank, and M S the molar mass of pure water. From this molar fraction and equation (32), Henry constants are then estimated by: where it is assumed that component i is pure in the vapor phase (composition of water is neglected), and u vap i is the fugacity coefficient of pure component i at temperature T and pressure P calculated with e-PPC-SAFT model and parameters of Table 4. To check the assumption that water in the vapor phase is negligible, the pure water vapor pressure is calculated for each experimental data and compared to the experimental total pressure. It is found that on average, the pure water vapor pressure represents around 3% of the total pressure, which is less of the experimental uncertainty. This suggests that this assumption is valid.
The reconciliation experimental data are plotted in Figure 1. An average experimental uncertainty is estimated to be ±10%. These data are fitted using empirical equation (38) (bold line in Fig. 1). The parameters of this correlation are given in Table 2. Note that for O 2 , it can be observed that a few experimental data published are not consistent with most of the data. These inconsistent data have not been considered to establish the correlation.

Light gas + brine systems
Experimental hydrogen solubility data in saline water are relatively scarce in literature and cover only a limited range of temperature (up to 333 K) and salinities (2 mol/kg) [44][45][46][47]. Lopez-Lazaro et al. reviewed these data and proposed additional ones at higher temperatures (up to 500 K) from molecular simulations [6]. This set of solubility data converted into the Henry constant is used to parameterize our model. Very recently, Chabab et al. [8] and Torín-Ollarves and Trusler [13] published new solubility data for temperatures up to 423 K and salinities up to 5 mol/kg. These new data are used in this work as a validation dataset for the model.
In contrast to hydrogen, the experimental data of oxygen solubility in salted water cover a wider range of temperature (up to 473 K) and salinities (5.8 mol/kg) [48][49][50][51]. As indicated in Table 3, the recent experimental oxygen solubility data of Chabab et al. [52], not used in the parameterization database, will be used to validate the model.
Concerning nitrogen, experimental solubility data in saline water are very scarce in the literature. The published Henry constants (or data which can be converted in Henry constant following Eq. (39)) concern a very restricted temperature range (up to 298 K) [53,54], and cannot be used alone to parameterize our model in an underground gas storage context. Thus, in contrast to other systems, the binary interaction parameter between N 2 and ions is directly fitted to match "PTx" experimental data at elevated pressure, covering a wider range of temperature (303 to 398 K) and salinities (up to 5.5 mol/kg) [55,56]. For hydrogen, nitrogen, and oxygen, the experimental compressibility factors are used to adjust the dispersion energy and segment diameter of each of these light gas, while segment numbers are taken from previous studies [92,93]. Table 4 gives the optimized parameters, and Figure 2 compares calculated values and reference data.
The average deviations between model results and the reference data are equal to 0.1%, 0.2%, and 0.1% for hydrogen, nitrogen, and oxygen, respectively. These deviations are less than the uncertainties reported by the NIST for the reference data (0.2%, 0.5%, and 1.2%, respectively). The parameters fitted are consistent with other PC-SAFT parameter sets available in the literature [94][95][96][97].

Brine
For the brine (mixture of H 2 O, Na +, and Cl À ), special attention should be paid to accurately reproduce MIAC close to NaCl saturation, since it corresponds to conditions encountered in salt caverns. It requires tuning not only pure ion parameters but also binary parameters between cation and anion since at high ionic concentration, ion pairing can be observed [98]. For pure ion parameters, we consider the hard sphere diameter equal to the Pauling diameter, and a segment number equal to 1.
In the e-PPC-SAFT framework, ions are assumed to interact with other particles only through association forces, and not through dispersion forces. Consequently, the energy dispersion parameter of ions is equal to 0. Finally, there are three parameters per ion to be adjusted: association energy, association volume, and MSA diameter. To better mimic ion-pairing formation at high salt concentration, the cross-binary parameters on association energy w ij and association volume u ij (Eq. (16)) between Na + and Cl À are also adjusted. A total of seven parameters is then used for brine parameterization. The e-PPC-SAFT parameters of pure water are directly taken from a previous study [16]. Finally, the optimized parameters are reported in Table 4 (for pure component parameters), and on Table 5 (for binary parameters). Figure 3 compares experimental MIAC [37], vapor pressures [38,39], and densities [40] with calculated values.
A very good agreement is found in the whole range of temperatures and salinities. The average deviations between experimental and calculated MIAC, vapor pressures, and densities are equal to 2.0%, 1.3%, and 1.4%, respectively.

Light gas + pure water systems
Correlation 37 is used as reference data to adjust interaction binary parameters of the e-PPC-SAFT model between H 2 and H 2 O, N 2 and H 2 O, and O 2 and H 2 O. This optimization is done for temperatures ranging from 273 to 473 K in order to focus on conditions of interest for gas underground  [48][49][50][51] Oxygen solubility data, from [52] Air + H 2 O + NaCl -None (no experimental data available) Note: the segment diameter of pure water is given by: r HS = 2.2423 + 0.51212 exp(0.001126 T) + 9904.13/T 2 .
The Author(s): Science and Technology for Energy Transition 78, 10 (2023) storage. As mentioned in Table 1, two binary parameters can be tuned: the k ij and l ij parameters, involved in the cross energy dispersion and the cross diameter of the NAHS term, respectively. A previous study based on molecular simulations of the H 2 + H 2 O system showed that the l ij parameter has a stronger influence on Henry's constant calculations [6]. Consequently, only this binary parameter is adjusted, and the k ij parameter is set to 0. A preliminary optimization using a constant value of the l ij parameter led to important deviations between experimental and calculated values. Furthermore, the model was not able to reproduce the maximum in the Henry constant curves observed experimentally. Thus, it was decided to introduce a temperature dependency for this binary parameter: The adjusted l ij parameters are given in Table 5. It is highly recommended to use this correlation only in the subcritical domain (temperature less than 647 K) to avoid the too-large value of l ij . Figure 3 presents the experimental data calculated data. This optimization leads to a good agreement, with an average deviation between 1.8%, 1.1%, and 1.3% for systems H 2 + H 2 O, N 2 + H 2 O, and O 2 + H 2 O, respectively.   [36]. Lines: this model.

Light gas + brine systems
In the e-PPC-SAFT model, ions are assumed to interact with solvent and neutral solutes only through repulsive and associative forces (no dispersion interaction). Thus, the l ij interaction parameter of the NAHS is the only binary parameter available to tune the model. This binary parameter between Na + and X (X = H 2 , O 2 or N 2 ) and between Cl À and X are taken equal since it is not possible to distinguish these two interactions independently in the brine. For H 2 + salted water and O 2 + salted water systems, the optimization of this parameter shows that a constant value (given in Tab. 5) is sufficient to reach a good accuracy as shown in Figure 5 and Figure 6 with an average deviation between experimental and calculated Henry constant of 3% and 8%, respectively. Note that for oxygen, experiments are not reported for fixed salinity values. Thus, the graphs in Figure 6 show experimental data gathered in specific salinity ranges.
For the N 2 + saline water system, a temperature dependency is found necessary to reach an acceptable average deviation (9%). The optimized parameters are reported in Table 5, and Figure 7 shows a comparison between experimental and calculated data.

Model validation
For each system, a set of experimental data is excluded from the model parameterization database and is used to validate our model. The model is validated not only on its capability to correctly predict data not used in the parameterization database, but also to extrapolate to mixtures. This last validation is of primary importance to accurately evaluates phase equilibrium between brine and air for CAES units. Air is considered in this work as a mixture of N 2 (80%mol.) and O 2 (20%mol.). Table 3 provides a synthetic view of the references of experimental data used for both model parameterization and model validation.
For H 2 + H 2 O, N 2 + H 2 O, and O 2 + H 2 O systems, the model predicts very satisfactory solubility data not used in the parameterization database, with an average deviation of 5%. The comparisons between experimental and pre-dicted bubble curves are reported in the Supplementary Material (Fig. S1). More interestingly, the model predicts very well the water content in the vapor phase, as illustrated in Figure 8, with an average deviation of 3.5%. Finally, the extension of the model to multicomponent systems is also validated by its capacity to accurately predict air solubility in water in a large range of temperatures and pressure, as illustrated in Figure 9, with an average deviation of 7%.
As for gas + pure water systems, the model proposed for gas + salted water systems is validated on its capability to predict data not included in the parameterization database as indicated in Table 3. For the H 2 + H 2 O + NaCl system, recent hydrogen solubility data are used for validation [8]. As illustrated in Figure 10, our model predicts with good accuracy these data in the whole temperature range (323 to 373 K) and salinity range (1-5 mol/kg). The average deviation between experiments and prediction is about 5%, while the experimental uncertainty estimated by the authors is very close (4%).
For O 2 + H 2 O + NaCl system, the recent experimental solubility data of Chabab et al. [52], not used in the parameterization database, are used to validate the model. As illustrated in Figure 11, the model predictions are found in good agreement with these experiments, with an average deviation of 7%. It can be nevertheless noticed that predictions are less accurate at the highest salinity (deviation of 15% at 4 mol/kg).
For N 2 + H 2 O + NaCl system, due to the lack of experimental data, all the existing data have been used for the model parameterization. Finally, for Air + H 2 O + NaCl, no data have been found in the literature to validate model predictions. It could be interesting for further studies to acquire new experimental data on such a system and compare them to our predictions.

Application to air storage in salt caverns
The e-PPC-SAFT equation of state parameterized in this work is used to evaluate some thermodynamic properties in the case of Compressed Air Energy Storage (CAES) in the salt cavern. Air is assumed to be a mixture of 80%   [37][38][39][40]. Lines: this model.   (Table 3). Blue diamonds: 0.5 mol/kg. Red squares: 1 mol/kg. Green triangles: 2 mol/kg.     For each temperature/pressure couple of the cycle, the e-PPC-SAFT model is used to evaluate the compressibility factor Z of the vapor phase. The evolution of this factor over the cycles is illustrated in Figure S3 of the Supplementary Material. This factor enables computing the total quantity of gas mixture in this phase. The e-PPC-SAFT is also used to compute phase equilibrium for each point of the cycle in order to determine more specifically the water content of the vapor phase. The knowledge of this gas water content during CAES operation is of primary importance, since water may increase corrosion issues in surface facilities and promote hydrate formation in pipes. Hence, Figure 12 illustrates the evolution of the water content in the gas produced during the cycles, as well as the corresponding mass of gaseous water present in the cavern. In order to highlight the effect of salinity, two calculations are carried out: in a hypothetical salt-free system, and a saturated halite brine, representative of salt caverns currently operated. The quantity of air solubilized in the brine is also plotted in Figure 13.  11. Predictions of oxygen solubility in salted water at 323 K (left) and 373 K (right). Lines: this model. Symbols: experimental data not used in the parameterization database [52]. Purple squares: 0.5 mol/kg. Orange diamonds: 1 mol/kg. Gray triangles: 4 mol/kg.  These figures show that in a salt cavern, with residual saturated brine at its bottom, the quantity of air solubilized is lower than what would be found in a hypothetical cavern with residual pure water, as a result of the "salting-out" effect. In presence of salt, water molecules mainly solvate ions and are less available to solubilize gas, leading also to lower volatility. As a result of the mass balance, the water content in the vapor phase is 30% less in salt cavern than in pure water. Note that the calculated quantities result only from thermodynamic considerations and neglect kinetic aspects, and more specifically the time needed to reach the thermodynamic equilibrium. Temperature gradients in the cavern are also neglected, although due to vaporization heat being absorbed in the brine-gas interface only, at the bottom of the cavern, while condensation heat release is distributed on the whole cavern volume, some conditions can create a cold gas layer just above the brine, preventing natural convection and may increase the time to reach this equilibrium in the whole cavern [101]. These assumptions necessarily cause an overestimation of the water content for short injection/withdrawing cycles. Nevertheless, this approach provides a first conservative estimation of gas moisture, useful to design surface dehydration units [102].

Conclusion
In this work, a parameterization of the e-PPC-SAFT equation of state has been carried out to predict phase equilibrium in salt caverns for hydrogen or compressed air storage. The amount of gas stored in the cavern being a key economical factor, the pure component parameters of the light gases (H 2 , N 2 , O 2 ) have been fitted to accurately match reference compressibility factors. Deviations are in the same range of experimental uncertainties. The brine properties have also to be carefully modeled, to accurately predict gas solubilities, phase density, and water volatility. Brine is assumed to be a mixture of water, Na + , and Cl À , and model parameters have been fitted to match experimental mean ionic activity coefficients, brine vapor pressure, and brine densities. The average deviation between experimental and calculated data are 2.0, 1.3, and 1.4%, respectively, for temperatures ranging from 273 to 473 K, and salinities up to saturation (6 mol/kg). The solubility of the light gases is first investigated in pure water. For N 2 and O 2 , a new Henry constant correlation is proposed based on a literature review of existing data, while for H 2 an existing recent correlation is used. These correlations are used to fit the binary interaction parameters of the e-PPC-SAFT model between H 2 + H 2 O, N 2 + H 2 O, and O 2 + H 2 O. Calculated Henry constants in pure water reproduce reference experimental correlation with good accuracy (average deviation less than 2%).
A similar approach is carried out for systems involving saline water. A literature review of the Henry constant of light gases in saline water is carried out, and binary interaction parameters of the e-PPC-SAFT model between H 2 + ion, N 2 + ion, and O 2 + ion have been fitted to match experimental data. Average deviations between 3% and 9% are obtained, in the same range of experimental uncertainties.
This parameterized model is validated on its capability to accurately predict data or properties not used in the parameterization database, but also on its capability of predicting mixture properties. Thus, it has been used to predict Air + pure water phase equilibrium with good accuracy (deviation less than 7% on air molar fraction solubilized), the air being considered as a mixture of 80%mol. N 2 + 20%mol. O 2 . The model has also been used to predict air + saline water phase equilibrium to evaluate the water content of the vapor phase for an application of CAES in a salt cavern. It can be noticed that there are no experimental data available for air + saline water phase equilibrium. Thus, it could be useful in further studies to acquire new experimental data and to compare the predictions generated by this new model. It is also worth noticing that the case study investigated in this work only deals with a thermodynamic approach. Kinetic and diffusive effects are neglected, but this approach can nevertheless be used as a first conservative evaluation to design surface dehydration units.

Supplementary material
The supplementary material of this article is available at https://stet-review.org/10.2516/stet/2023004/olm. Figure S1.  Figure S2. Temperature and pressure variations during 8 cycles of CAES unit. Solid line: temperature. Dashed lines: pressure. Figure S3. Evolution of the vapour phase compressibility factor during the cycles, calculated from the e-PPC-SAFT model.