Enhancing Hydrogen Recovery from Saline Aquifers: Quantifying Wettability and Hysteresis Influence and Minimizing Losses with a Cushion Gas

: This study aims to numerically assess the impact of wettability and relative permeability hysteresis on hydrogen losses during underground hydrogen storage (UHS) and explore strategies to minimize them by using an appropriate cushion gas. The research utilizes the Carlson model to calculate the saturation of trapped gas and the Killough model to account for water hysteresis. By incorporating the Land coefficient based on laboratory-measured data for a hydrogen/brine system, our findings demonstrate a significant influence of gas hysteresis on the hydrogen recovery factor when H 2 is used as a cushion gas. The base model, which neglects the hysteresis effect, indicates a recovery factor of 78% by the fourth cycle, which can be improved. In contrast, the modified model, which considers hysteresis and results in a trapped gas saturation of approximately 17%, shows a hydrogen recovery factor of 45% by the fourth cycle. Additionally, gas hysteresis has a notable impact on water production, with an observed 12.5% increase in volume in the model that incorporates gas hysteresis. Furthermore, optimization of the recovery process was conducted by evaluating different cushion gases such as CO 2 , N 2 , and CH 4 , with the latter proving to be the optimal choice. These findings enhance the accuracy of estimating the H 2 recovery factor, which is crucial for assessing the feasibility of storage projects.


Introduction
The anticipated global H 2 demand is projected to surpass 600 million metric tons within the next few decades [1].Hence, the need for feasible and sufficient H 2 storage is required.Hydrogen storage is a critical concept associated with the broader goals of energy transition, energy security and sustainability.Some storage technologies have been developed to contain hydrogen and convert it to electricity or utilize it as a fuel when needed.Those technologies were limited by the size of the vessels, or tanks used to store the hydrogen, and required high-energy cost and expensive materials.Therefore, researchers considered underground hydrogen storage to overcome these challenges by offering a large storage capacity and a cost-effective solution while ensuring high safety and energy security [2].Depleted gas reservoirs and deep saline aquifers are among the appealing options for undergoing hydrogen storage.Table 1 provides a summary of the surface and subsurface storage technologies developed or currently being researched for storing hydrogen.Poor geological characterization [4,5] However, H 2 has unique physical and chemical properties that are unlike any other geologically stored fluid.It is a highly reactive, highly buoyant, and very flammable substance [3].Dedicated research initiatives are applied to unravel the complex storage mechanisms inherent in the underground hydrogen storage (UHS) process.Primarily, the mitigation of excessive hydrogen loss emerges as a critical research domain, necessitating the systematic evaluation of geologically pertinent conditions to ensure heightened storage efficiency and economic viability.Herein, we shed light on the influence of wettability and hysteresis phenomena on the trapping mechanism during UHS, which can significantly impact the recovery of H 2 .
Wettability, an important petrophysical parameter, is the tendency of a surface to adhere to one fluid in the presence of another fluid governed by intermolecular interactions [4], carrying significant implications for the mobility and storage efficacy of hydrogen within these geological formations.The mobility of H 2 , on the other hand, is impacted by the significant phenomenon of relative permeability hysteresis.Alterations in hydrogen flow direction can impact the flow pathways of hydrogen and other fluids, such as water, owing to the hysteresis effect of relative permeability, thereby indirectly influencing storage efficiency and hydrogen recovery.Cyclic injection and production may induce the detachment and entrapment of hydrogen bubbles, impeding their subsequent flow or recovery.Consequently, the fluid saturation and the fluids' flow behavior are influenced by the history of the drainage cycles (the increase in the saturation of the nonwetting phase) and the imbibition cycles (the increase in the saturation of the wetting phase), which comprises the hysteresis effect [6].Therefore, this study investigates the influence of wettability on hydrogen storage performance.
Additionally, we analyze the impact of the hysteresis effect on hydrogen loss and assess the feasibility of the storage project in sandstone saline aquifers.Numerous experimental studies have investigated the wettability of sandstone in H 2 /water systems.Wettability characterization is predominantly conducted through contact angle measurements.Contact angles falling within the 0-60 • range typically denote strongly water-wet conditions, and the wetness to water decreases when contact angles exceed 70 • .In this study, the term "neutrally wet" is used to describe a contact angle that is greater than 70 • but less than 120 • .This characterizes a surface that is neutral and has no preference for imbibing one fluid over the other in a hydrogen and brine system [7].
Iglauer et al. [8] and Hashemi et al. [9] employed the tilted plate and captive bubble methods on both pure quartz and aged quartz under reservoir conditions.Their findings revealed contact angles in the range of 25-70 • , signifying that sandstone exhibits characteristics ranging from strongly water-wet to intermediate-wet in the presence of H 2 /water.Other researchers performed surface complexation modeling and 3D microcomputed tomography (CT) and agreed that sandstone would remain water-wet during hydrogen storage [10].However, studies focusing on the influence of bacterial growth on wettability observed that sulfate-reducing bacteria can induce a substantial increase in contact angle measurements.Consequently, this alteration shifts the wettability from a water-wet condition to a neutrally wet condition [11,12].These observations substantiate the variability of sandstone wettability, indicating potential states of either water-wet or neutrally wet conditions.Nonetheless, determining the optimal wettability condition for hydrogen storage and assessing its impact on recovery enhancement or reduction remains controversial in the literature.
In the context of hydrogen loss in underground storage, the cushion gas plays a crucial role in preventing the escape of injected gas and maintaining the necessary operational pressure, mitigating a portion of the inherent loss due to residual entrapment.The cushion gas is therefore intended to be kept in the reservoir and not to be recovered until the storage operation is concluded [13].Therefore, the inclusion of cushion gas substantially impacts the overall capital expenditure.H 2 is known to be a valuable commodity, and to use a vast volume of it as cushion gas could increase the total cost and create a limitation for feasible and economical storage [14].To address this issue, researchers recommended using alternative gases to be injected as cushion gas to decrease H 2 loss in the formation, increasing H 2 recovery, and reducing the project's budget [15,16].In alignment with the operational frameworks observed in European town gas storage projects, where H 2 constitutes up to 60% of the storage volume along with a blend of other gases including CH 4 , N 2 , and CO 2 , researchers have proposed employing these gases as an alternative cushion gas in underground H 2 storage [17].It is crucial to emphasize that commercialscale underground H 2 projects in depleted gas reservoirs or saline aquifers have not been carried out yet.It is also important to note that numerical simulation studies play a vital role in understanding, analyzing, and optimizing underground H 2 storage systems to establish successful and feasible projects.
This area of study has been relatively unexplored, with only a few researchers taking an interest [18].In 2021, Lysyy et al. [19] studied the effect of injecting formation gas as a cushion gas on H 2 recovery in depleted reservoirs using a black-oil simulation model (ECLIPSE-100) [19].Their fluid model consisted of a PVT black-oil table for oil, water, and generic gas phases and used default values for natural gas/H 2 relative permeability as a function of total gas saturation.Another study by Zamehrian and Sedaee [20] investigated the role of cushion gas but targeted the condensate gas reservoir.They examined different cushion gases and studied their effect on H 2 recovery and purity in a partially depleted condensate gas reservoir located in the Middle East.Their results suggested using N 2 as the gas cushion for maximizing H 2 recovery.However, CO 2 was effective in terms of maximizing condensate production.Multiple researchers agreed on having the highest recovery of H 2 when using CH 4 and N 2 as the cushion gas.Kanaani et al. [16] carried out a similar study to test different cushion gases (CH 4 , N 2 , and CO 2 ) on H 2 production in depleted oil reservoirs.They concluded that using N 2 as the cushion gas resulted in the highest average reservoir pressure.However, the maximum H 2 recovery was achieved by injecting CH 4 as the cushion gas.The results highlighted the significance of the cushion gas type, surpassing other factors that could influence the cumulative H 2 recovery, including the volume of the cushion gas.Similarly, Heinemann et al., 2021 [21], focused their study on the role of cushion gas for H 2 storage in an open saline aquifer.They investigated how reservoir depth, the shape of the trap, and reservoir permeability would impact the ratio of cushion gas to working gas, using H 2 itself as the cushion gas.No other types of gases were investigated.Even though using a cushion gas for underground H 2 storage seems critical, researchers are still trying to understand the multiphase flow that is associated with this process.Wang et al. [22] performed a modeling study to investigate the flow behavior of H 2 with CO 2 as cushion gas and analyze the hydrodynamic mechanisms of H 2 storage.They considered important factors such as viscous fingering, capillary force, and gravity segregation.However, they made major assumptions that affected their results such as neglecting gas dissolution in brine.Neglecting gas dissolution, particularly, may result in an overestimation of the produced H 2 impurity, because the dissolution of CO 2 in brine is considered one of the main trapping mechanisms that will minimize the amount of mobile cushion gas.The main observation of their study was that H 2 recovery performance is better when the gravity force is larger than the capillary force and they achieved the maximum purity of H 2 in this case study.This is because the density between H 2 and CO 2 is different, and the mixing zone is mainly controlled by gravity.On the other hand, when viscous force is dominant, the purity of the produced H 2 is decreased.
Although all these publications bring important insights and add knowledge to our understanding of the rule of cushion gas and its impact on the H 2 storage feasibility, the approaches used to run these simulations are rather simple and will not lead to the most accurate results.
The previous experience of modeling methane reservoirs and CCS projects in numerical simulation is critical to develop and model hydrogen storage operations through advanced simulators.However, it is imperative to recognize that each type of operation presents distinct considerations owing to the differing physical and chemical properties of gases, as well as the specific objectives of the projects.
In this paper, we propose an innovative methodology to precisely model fluid behavior and gas mixing in a fully composition approach to study the impact of using different cushion gases and analyze their impact on H 2 recovery, H 2 purity, water production, and the overall performance of the storage facility.We take into consideration modeling the behavior of a specific fluid or system accurately by tuning the equation of state, viscosity, and solubility models to match the laboratory data.Our approach involves employing the accurate relative permeability, which characterizes the flow behavior of gas mixtures, thereby ensuring a heightened level of simulation accuracy.There exists a research gap in the examination of H 2 storage within saline aquifers, and therefore, we are performing our study in a saline aquifer reservoir.We present three alternative options to be used as the cushion gas, besides the H 2 itself; CO 2, which has been intensively studied to be sequestrated permanently in saline aquifers [23][24][25][26], CH 4 , which is also a potent greenhouse gas that we are aiming to limit its emission [27][28][29], and N 2 , which was historically used as a mixture with H 2 and other gases, and is known to be a cheaper gas, and non-reactive [30,31].
Besides that, this study also aims to enhance the knowledge of the feasibility of hydrogen storage processes in saline aquifers by investigating the influence of wettability and hysteresis through numerical simulations.Previous studies on gas storage, such as carbon capture and storage (CCS), have underscored the significance of incorporating the hysteresis effect in numerical simulations [24,32,33].Including hysteresis is crucial for accurately accounting for trapped CO 2 and significantly influencing storage security [34,35].In contrast to CCS, the inherent trapping mechanisms in hydrogen underground storage may adversely affect the feasibility of the project [36,37].Neglecting the hysteresis effect can lead to a substantial underestimation of residual trapped gas volume.Conversely, overestimations are likely for mobile hydrogen and extracted water quantities.This underscores the imperative need to consider relative permeability hysteresis to ensure precise evaluations in subsurface hydrogen storage.
Finally, we have made some major assumptions in this work.Since the study is mainly focused on the petrophysical aspects and understanding the fluids' behavior, we assumed no constraint is imposed on water production as it may obscure the clarity and precision of the research results.Studying the geochemical, or bacterial activity in UHS is beyond the scope of our study, so we did not account for these phenomena, however, some studies on this topic can be found in the literature [38,39].

Geological Modeling
This section will describe the methodological steps we will consider in modeling and examining the different wettability conditions, and hysteresis in UHS.A full description of the geological model and the operational strategy are described in detail below.
A synthetic anticline aquifer model was simulated using CMG-GEM, version 2022. 10 [40], to perform multiple studies during UHS.CMG-GEM is known for its versatile capability and adaptability in predicting and analyzing various thermodynamic properties and fluid behavior within underground hydrogen storage systems.It employs finite volume and finite difference methods to discretize material balance and energy balance equations in space and time.
The three-dimensional model is symmetrical, with an average reservoir thickness of 90 m.We intentionally chose a simple structure to avoid the impact of heterogeneous geological complexity.The model consists of 150 × 50 × 50 grid blocks, and each grid is 100 × 100 × 2 m 3 in size.The horizontal permeability and the porosity have a heterogeneous distribution with an average of 20% and 200 mD, respectively.Figure 1 shows the horizontal permeability profile.

PÚBLICA
Finally, we have made some major assumptions in this work.Since the study is mainly focused on the petrophysical aspects and understanding the fluids' behavior, we assumed no constraint is imposed on water production as it may obscure the clarity and precision of the research results.Studying the geochemical, or bacterial activity in UHS is beyond the scope of our study, so we did not account for these phenomena, however, some studies on this topic can be found in the literature [38,39].

Geological Modeling
This section will describe the methodological steps we will consider in modeling and examining the different wettability conditions, and hysteresis in UHS.A full description of the geological model and the operational strategy are described in detail below.
A synthetic anticline aquifer model was simulated using CMG-GEM, version 2022.10 [40], to perform multiple studies during UHS.CMG-GEM is known for its versatile capability and adaptability in predicting and analyzing various thermodynamic properties and fluid behavior within underground hydrogen storage systems.It employs finite volume and finite difference methods to discretize material balance and energy balance equations in space and time.
The three-dimensional model is symmetrical, with an average reservoir thickness of 90 m.We intentionally chose a simple structure to avoid the impact of heterogeneous geological complexity.The model consists of 150 × 50 × 50 grid blocks, and each grid is 100 × 100 × 2 m 3 in size.The horizontal permeability and the porosity have a heterogeneous distribution with an average of 20% and 200 mD, respectively.Figure 1 shows the horizontal permeability profile.For simplicity, we assumed a single well that operates as the injector and the producer, placed at the top center position of the anticlinal.To model the phase behavior and evaluate the hydrogen gas properties in saline aquifers, the Peng-Robinson cubic equation-of-state (EoS) is used [41].Henry's Law [42] is considered to account for gas solubility.Table 2 summarizes the reservoir properties.Here, we present a synthetic model that does not correspond to the properties of an actual aquifer.However, to maintain realistic and representative values, the selected properties fall within the average range of geological characteristics typical of deep saline aquifers [6,7].For simplicity, we assumed a single well that operates as the injector and the producer, placed at the top center position of the anticlinal.To model the phase behavior and evaluate the hydrogen gas properties in saline aquifers, the Peng-Robinson cubic equation-of-state (EoS) is used [41].Henry's Law [42] is considered to account for gas solubility.Table 2 summarizes the reservoir properties.Here, we present a synthetic model that does not correspond to the properties of an actual aquifer.However, to maintain realistic and representative values, the selected properties fall within the average range of geological characteristics typical of deep saline aquifers [6,7].By delving into a more advanced research sphere such as CO 2 sequestration simulations, we noticed the emphasis researchers had shown for the role of grid sizing and its impact on the accuracy of gas sweep efficiency, solubility, and plume front [43,44].The existing literature exhibits a perceptible gap in studies concerning grid refinement within the domain of H 2 storage.Researchers have advocated for further investigations in this area to develop the required knowledge in understanding numerical dispersion impact in H 2 storage simulations, and to address the magnitude of the accuracy of the results [44][45][46].We therefore took the initiative to study this aspect and applied a local grid refinement to the original model we presented previously.The initial grids were divided into 5 parts in the x, y, and z direction resulting in refined grids of 20 × 10 × 0.4 m 3 .Grid refinement beyond this size resulted in increasing simulation time and showed negligible differences in results.
The hydrogen is injected for one year as a cushion gas at a rate of 250,000 m 3 /day.Following that, we set the injection constraint to be 250,000 m 3 /day over six 6-month periods (spring and summer) for the working gas.We shut down the well for three months and started producing at a rate of 500,000 m 3 /day for three months (winter season), as depicted in Figure 2. Our model included four complete cycles of hydrogen injection/production over four years, and we performed prolonged production for some analysis.The rate of H 2 injection is within the average values we've observed in the literature [47,48].Considering the reservoir depth, and to ensure no fracturing of the reservoir rock, we have set a maximum injection pressure constraint of 27,000 kPa, which is below the fracturing pressure, and a minimum bottom-hole pressure of 3500 kPa.In our model, we included 4 complete cycles of H 2 injection/production over 4 years.Our study aims to understand the petrophysics properties and flow behavior phenomena, therefore, we had no constraint imposed on water production and assumed the facility could handle the amount of produced water.By delving into a more advanced research sphere such as CO2 sequestration simulations, we noticed the emphasis researchers had shown for the role of grid sizing and its impact on the accuracy of gas sweep efficiency, solubility, and plume front [43,44].The existing literature exhibits a perceptible gap in studies concerning grid refinement within the domain of H2 storage.Researchers have advocated for further investigations in this area to develop the required knowledge in understanding numerical dispersion impact in H2 storage simulations, and to address the magnitude of the accuracy of the results [44][45][46].We therefore took the initiative to study this aspect and applied a local grid refinement to the original model we presented previously.The initial grids were divided into 5 parts in the x, y, and z direction resulting in refined grids of 20 × 10 × 0.4 m³.Grid refinement beyond this size resulted in increasing simulation time and showed negligible differences in results.
The hydrogen is injected for one year as a cushion gas at a rate of 250,000 m³/day.Following that, we set the injection constraint to be 250,000 m 3 /day over six 6-month periods (spring and summer) for the working gas.We shut down the well for three months and started producing at a rate of 500,000 m 3 /day for three months (winter season), as depicted in Figure 2. Our model included four complete cycles of hydrogen injection/production over four years, and we performed prolonged production for some analysis.The rate of H2 injection is within the average values we've observed in the literature [47,48].Considering the reservoir depth, and to ensure no fracturing of the reservoir rock, we have set a maximum injection pressure constraint of 27,000 kPa, which is below the fracturing pressure, and a minimum bottom-hole pressure of 3500 kPa.In our model, we included 4 complete cycles of H2 injection/production over 4 years.Our study aims to understand the petrophysics properties and flow behavior phenomena, therefore, we had no constraint imposed on water production and assumed the facility could handle the amount of produced water.

Relative Permeability Curves at Different Wettability Conditions
Underground hydrogen storage in aquifers is a relatively new research area, and therefore there are limited experimental studies investigating the petrophysical parameters and fluid behavior for H2, such as relative permeability measurements.In 2018, Yekta

Relative Permeability Curves at Different Wettability Conditions
Underground hydrogen storage in aquifers is a relatively new research area, and therefore there are limited experimental studies investigating the petrophysical parameters and fluid behavior for H 2 , such as relative permeability measurements.In 2018, Yekta et al. [49] carried out the first measurement of hydrogen and brine relative permeability in a water-wet sandstone.Following this study, few other researchers have experimentally studied the hydrogen-brine relative permeability [50][51][52].All the experimental work that measured the relative permeability for hydrogen-brine performed their experiment in water-wet conditions.However, using pore network modeling, Hashemi et al. [53] generated the relative permeability for different wettability conditions (brine contact angles of 51 • and 83 • ).In our model, we input the data from Yekta et al. [49] for the water-wet condition (brine contact angle of 35 • ) and the data from Hashemi et al. [53] using the modeled relative permeability for a contact angle of 83 • to simulate a neutrally wet condition.These curves were previously modified by Pan et al. [54] to smooth the experimental data, and the modified curves were implemented in this study.Figure 3 represents the relative permeability data used for this case study.
Hydrogen 2024, 5, FOR PEER REVIEW 7 PÚBLICA et al. [49] carried out the first measurement of hydrogen and brine relative permeability in a water-wet sandstone.Following this study, few other researchers have experimentally studied the hydrogen-brine relative permeability [50][51][52].All the experimental work that measured the relative permeability for hydrogen-brine performed their experiment in water-wet conditions.However, using pore network modeling, Hashemi et al. [53] generated the relative permeability for different wettability conditions (brine contact angles of 51° and 83°).In our model, we input the data from Yekta et al. [49] for the water-wet condition (brine contact angle of 35°) and the data from Hashemi et al. [53] using the modeled relative permeability for a contact angle of 83° to simulate a neutrally wet condition.These curves were previously modified by Pan et al. [54] to smooth the experimental data, and the modified curves were implemented in this study.Figure 3 represents the relative permeability data used for this case study.

Gas Relative Permeability Hysteresis
Previous studies on the influence of hysteresis on gas storage performance have highlighted its significant role [43].However, there is a limited number of experimental studies in the literature investigating this phenomenon for UHS.Furthermore, this parameter has been neglected in most published simulation studies.In this study, the data of drainage relative permeability curves and the maximum residual gas saturation by the end of a complete drainage cycle are used to calculate the residual gas saturation (end-point gas saturation) of the scanning curves through the Land gas trapping and Carlson models [55,56].The relative permeability curves for the hydrogen-brine system were analyzed [54] to calculate the trapped gas saturation ( ) used on the following equation: where  is the initial saturation of gas at the beginning of the imbibition and C is the Land coefficient calculated as follows:

Gas Relative Permeability Hysteresis
Previous studies on the influence of hysteresis on gas storage performance have highlighted its significant role [43].However, there is a limited number of experimental studies in the literature investigating this phenomenon for UHS.Furthermore, this parameter has been neglected in most published simulation studies.In this study, the data of drainage relative permeability curves and the maximum residual gas saturation by the end of a complete drainage cycle are used to calculate the residual gas saturation (end-point gas saturation) of the scanning curves through the Land gas trapping and Carlson models [55,56].The relative permeability curves for the hydrogen-brine system were analyzed [54] to calculate the trapped gas saturation (S gt ) used on the following equation: where S gi is the initial saturation of gas at the beginning of the imbibition and C is the Land coefficient calculated as follows: where S grmax is the maximum residual gas saturation and S gmax is the maximum historical gas saturation.
The Land coefficient (C) was calculated from Equation (2) to be 2.5, which is within the average range of what has been reported previously by Bo et al. [50], who stated that the values of the Land coefficient are between 2 and 3.2.The trapped gas saturation (S gt ), based on the calculated Land coefficient, is determined to be 17% roughly.Notably, the Land coefficient values for H 2 are relatively higher than those observed for CO 2 hysteresis studies, where the Land coefficient falls between 1 and 1.7 [57].

Water Relative Permeability Hysteresis
We reviewed the relevant literature for CO 2 storage that investigated the water relative permeability hysteresis, and most of the studies agreed on the same conclusion: water relative permeability exhibited non-hysteretic behavior, unlike gas [6].However, some of the experimental studies for hydrogen drainage and imbibition relative permeability curves demonstrated an obvious water hysteresis [4].Henceforth, we employed the modified Killough model [58] to examine the influence of water hysteresis.Among alternative hysteresis empirical models such as the Beattie et al. model [59], Kjosavik et al. model [60], and Spiteri et al. model [32], the Killough model demonstrated superior accuracy in predicting the hysteresis effects on relative permeability within gas/water systems [61].However, the scanning water relative permeability curve measured from the original Killough equation can intersect with the primary drainage curve, resulting in cross-curve points with inadequate K rw values.Therefore, the modified Killough model was used here to overcome this obstacle.To model the water hysteresis effect, the Killough model requires the bounding curves (both drainage and imbibition curves) to calculate the scanning curves mathematically by the following equation: The maximum value of K r at a specific scanning curve K Im rw (S Nr ) is calculated as follows: where λ is a given parameter, here we assumed λ to be 0.5 (default value).The modified Killough model is calculated as the following: where Nr , maximum possible residual or trapped nonwetting saturation; K Im rw , imbibition water-phase relative permeability; K Dr rw , drainage water-phase relative permeability; K * Im rw , experimental or analytical curve for imbibition water-phase relative permeability, which falls between the maximum possible water saturation and the maximum possible trapped water saturation.
The relative permeability data presented previously [46,49,50] are utilized for the water hysteresis analysis.However, the experiment by Yekta et al. [49] did not include the imbibition curve; hence, the imbibition curve modeled by Hashemi et al. [53] will be implemented.Although the water relative permeability is modeled at a contact angle of 83 • , we assumed that its applicability extends to a contact angle of 35 • as well.
Finally, the simulation does not incorporate the phenomenon of capillary pressure hysteresis because no data on capillary pressure hysteresis were available in the literature.However, Nazari [62] conducted an investigation utilizing nitrogen capillary pressure data as an alternative and concluded that the impact of capillary hysteresis is negligible.The chemical/biochemical trapping is beyond the scope of this study; therefore, we assumed no hydrogen loss due to chemical trapping.Other studies on hydrogen consumption by bacteria or geochemical reactions can be found in the literature [39,63].
Our investigation comprises three models for comparative analysis: the initial model excludes any hysteresis, the second model incorporates only hydrogen hysteresis effects, and the third model integrates both gas and water hysteresis impacts.The primary objective is to assess the influence of hysteresis on hydrogen plume distribution, hydrogen recovery efficiency, and water production rates.

Rock-Fluid Modeling for Cushion Gas Assessment
Accurate application of the relative permeability curves ensures a meaningful prediction of fluid displacement, fluids' mobility, and H 2 recovery.Furthermore, relative permeability plays a crucial role in reservoir simulation for underground H 2 storage.These curves represent the variation in fluid flow through porous media under different saturation conditions and describe the displacement behavior of fluids within a reservoir.However, when considering an alternative cushion gas other than the H 2 itself, it is important to account for the relative permeability curves that are associated with that gas, and the mixture of gases.Therefore, we considered utilizing the interpolation function to enable a better estimation of relative permeability values at different saturations and various gas compositions.Measuring relative permeability curves in the laboratory is demanding in terms of time and cost, and prone to inaccuracy [64].The existing literature lacks experimental data on relative permeabilities for a mixture of gases with brine.The novelty of the technique we are presenting in this work serves as a bridge between the existing gaps in experimental work, and the heightened accuracy required for numerical simulations related to underground H 2 storage.
In this specific context, we concentrate on the interpolation of multiple relative permeability sets using a two-step method of interpolation that is based on the phase molar fraction of a component as the first step; the interpolation is applied to the saturation (horizontal axis) endpoints, and the normalized curves will be modified accordingly.Then, the relative permeability values (vertical axes) are interpolated at the required saturation.Considering two relative permeability sets, A and B, the residual gas saturation for each set (S A grw , S B grw ) are interpolated to calculate S x i grw , where x i is the phase molar fraction of the specified component [40].
At a specific phase molar fraction (x i ), Equation ( 1) is used to compute S x i grw : where ω i is the weighting factor calculated from Equation (7) as follows: where x A , x B represent the mole fraction of the component specified for each relative permeability set A and B. And for the second stage, Kr x i g is similarly computed using Equation (8): Likewise, the values for water relative permeability are calculated.Relative permeability sets are collected from the literature for H 2 /water, and similarly, CH 4 /water, CO 2 /water and N 2 /water are measured for sandstone formation [49,51,57,65].The relative permeability curves are depicted in Figures 4 and 5.

Impact of Wettability on Underground Hydrogen Storage
Previous research on wettability had conflicting findings regarding its impact on hydrogen production, as reported in various experimental studies [66].However, the outcomes of our study demonstrated that hydrogen recovery increased by 8.5% in the water-  [51], CH 4 in blue [65], H 2 in gray [49], and CO 2 in yellow [57]) in gas/brine system.Gas relative permeability curves for different cushion gases (N2 in orange [51], CH4 in blue [65], H2 in gray [49], and CO2 in yellow [57]) in gas/brine system.

Impact of Wettability on Underground Hydrogen Storage
Previous research on wettability had conflicting findings regarding its impact on hydrogen production, as reported in various experimental studies [66].However, the outcomes of our study demonstrated that hydrogen recovery increased by 8.5% in the water-  [51], CH 4 in blue [65], H 2 in gray [49], and CO 2 in yellow [57]) in gas/brine system.

Impact of Wettability on Underground Hydrogen Storage
Previous research on wettability had conflicting findings regarding its impact on hydrogen production, as reported in various experimental studies [66].However, the outcomes of our study demonstrated that hydrogen recovery increased by 8.5% in the water-wet case compared to the neutrally wet case by the end of the fourth cycle.Of the entire working gas mass, 76.5% was recovered under strongly water-wet conditions and 68% under neutrally wet conditions.Figure 6 below exhibits the cumulative hydrogen mass injection results versus the cumulative mass production for the two studied cases.
Hydrogen 2024, 5, FOR PEER REVIEW 11 ÚBLICA wet case compared to the neutrally wet case by the end of the fourth cycle.Of the entire working gas mass, 76.5% was recovered under strongly water-wet conditions and 68% under neutrally wet conditions.Figure 6 below exhibits the cumulative hydrogen mass injection results versus the cumulative mass production for the two studied cases.Similarly, our results showed a 10% increase in water production in the water-wet condition compared to the neutrally wet, as illustrated in Figure 7.This result is directly attributed to the higher water relative permeability values and greater water mobility in the water-wet condition.Moreover, when water is in the wetting phase, it will occupy the corners of the pores, which facilitates the connectivity of the flow and thus increases water production.On the other hand, in a naturally wet condition, neither the water nor the hydrogen will exhibit a preference for large versus small pores.This will still lead to the entrapment of some water in the pores.As a result, the flow connectivity of water decreases, leading to a substantial decline in water production.This phenomenon can be attributed to the lack of effective connectivity between the water molecules.
To better understand the impact of the wettability condition on the system, we investigated the bottom-hole pressure (BHP) in the injector well for both water-wet and neutrally wet scenarios.This investigation was crucial in determining the stability of the system and assessing its reliability for long-term operations.Figure 8 demonstrates the steady BHP for the water-wet scenario, contrasting with the oscillations observed for the neutrally wet scenario.Similarly, our results showed a 10% increase in water production in the water-wet condition compared to the neutrally wet, as illustrated in Figure 7.
Hydrogen 2024, 5, FOR PEER REVIEW 11 PÚBLICA wet case compared to the neutrally wet case by the end of the fourth cycle.Of the entire working gas mass, 76.5% was recovered under strongly water-wet conditions and 68% under neutrally wet conditions.Figure 6 below exhibits the cumulative hydrogen mass injection results versus the cumulative mass production for the two studied cases.Similarly, our results showed a 10% increase in water production in the water-wet condition compared to the neutrally wet, as illustrated in Figure 7.This result is directly attributed to the higher water relative permeability values and greater water mobility in the water-wet condition.Moreover, when water is in the wetting phase, it will occupy the corners of the pores, which facilitates the connectivity of the flow and thus increases water production.On the other hand, in a naturally wet condition, neither the water nor the hydrogen will exhibit a preference for large versus small pores.This will still lead to the entrapment of some water in the pores.As a result, the flow connectivity of water decreases, leading to a substantial decline in water production.This phenomenon can be attributed to the lack of effective connectivity between the water molecules.
To better understand the impact of the wettability condition on the system, we investigated the bottom-hole pressure (BHP) in the injector well for both water-wet and neutrally wet scenarios.This investigation was crucial in determining the stability of the system and assessing its reliability for long-term operations.Figure 8 demonstrates the steady BHP for the water-wet scenario, contrasting with the oscillations observed for the neutrally wet scenario.This result is directly attributed to the higher water relative permeability values and greater water mobility in the water-wet condition.Moreover, when water is in the wetting phase, it will occupy the corners of the pores, which facilitates the connectivity of the flow and thus increases water production.On the other hand, in a naturally wet condition, neither the water nor the hydrogen will exhibit a preference for large versus small pores.This will still lead to the entrapment of some water in the pores.As a result, the flow connectivity of water decreases, leading to a substantial decline in water production.This phenomenon can be attributed to the lack of effective connectivity between the water molecules.
To better understand the impact of the wettability condition on the system, we investigated the bottom-hole pressure (BHP) in the injector well for both water-wet and neutrally wet scenarios.This investigation was crucial in determining the stability of the system and assessing its reliability for long-term operations.Figure 8 demonstrates the steady BHP for the water-wet scenario, contrasting with the oscillations observed for the neutrally wet scenario.This finding implies that displacing water in a neutrally wet environment poses significant challenges.However, further analysis requires additional data on relative permeability curves for contact angles exceeding 90°.

Impact of Gas and Water Hysteresis on Underground Hydrogen Storage
The comprehensive investigation into the influence of relative permeability hysteresis on the injected H2 is pivotal for achieving accurate predictions concerning the fate of a storage project.This encompasses the precise estimation of lost and recovered H2.Here, we compare Case 1 (no hysteresis), Case 2 (gas hysteresis only), and Case 3 (water and gas hysteresis) to illustrate the impact of this phenomenon.Firstly, Case 1 and Case 2 are compared to study the impact of gas hysteresis.Secondly, Case 2 and Case 3 are compared to investigate the importance of including water hysteresis.

The Impact of Gas Relative Permeability Hysteresis
Forced drainage and spontaneous imbibition are the two primary mechanisms governing saline aquifers' hydrogen injection and withdrawal cycles.The alternation of flow directions in the repetitive injection and withdrawal cycles induces a hysteresis effect, necessitating thorough investigation [32].Here, we investigated the impact of gas hysteresis and its influence on hydrogen recovery.The operation comprises four complete cycles conducted over four years, followed by an extended production phase until depletion.The results indicate that the hysteresis effect exerts a substantial impact on recovery, leading to a decrease in the H2 recovery factor from 76.5% by the fourth cycle for Case 1 to a hydrogen recovery factor of 45% when accounting for hysteresis in Case 2, as illustrated in Figure 9.This finding implies that displacing water in a neutrally wet environment poses significant challenges.However, further analysis requires additional data on relative permeability curves for contact angles exceeding 90 • .

Impact of Gas and Water Hysteresis on Underground Hydrogen Storage
The comprehensive investigation into the influence of relative permeability hysteresis on the injected H 2 is pivotal for achieving accurate predictions concerning the fate of a storage project.This encompasses the precise estimation of lost and recovered H 2 .Here, we compare Case 1 (no hysteresis), Case 2 (gas hysteresis only), and Case 3 (water and gas hysteresis) to illustrate the impact of this phenomenon.Firstly, Case 1 and Case 2 are compared to study the impact of gas hysteresis.Secondly, Case 2 and Case 3 are compared to investigate the importance of including water hysteresis.

The Impact of Gas Relative Permeability Hysteresis
Forced drainage and spontaneous imbibition are the two primary mechanisms governing saline aquifers' hydrogen injection and withdrawal cycles.The alternation of flow directions in the repetitive injection and withdrawal cycles induces a hysteresis effect, necessitating thorough investigation [32].Here, we investigated the impact of gas hysteresis and its influence on hydrogen recovery.The operation comprises four complete cycles conducted over four years, followed by an extended production phase until depletion.The results indicate that the hysteresis effect exerts a substantial impact on recovery, leading to a decrease in the H 2 recovery factor from 76.5% by the fourth cycle for Case 1 to a hydrogen recovery factor of 45% when accounting for hysteresis in Case 2, as illustrated in Figure 9.This finding implies that displacing water in a neutrally wet environment poses significant challenges.However, further analysis requires additional data on relative permeability curves for contact angles exceeding 90°.

Impact of Gas and Water Hysteresis on Underground Hydrogen Storage
The comprehensive investigation into the influence of relative permeability hysteresis on the injected H2 is pivotal for achieving accurate predictions concerning the fate of a storage project.This encompasses the precise estimation of lost and recovered H2.Here, we compare Case 1 (no hysteresis), Case 2 (gas hysteresis only), and Case 3 (water and gas hysteresis) to illustrate the impact of this phenomenon.Firstly, Case 1 and Case 2 are compared to study the impact of gas hysteresis.Secondly, Case 2 and Case 3 are compared to investigate the importance of including water hysteresis.

The Impact of Gas Relative Permeability Hysteresis
Forced drainage and spontaneous imbibition are the two primary mechanisms governing saline aquifers' hydrogen injection and withdrawal cycles.The alternation of flow directions in the repetitive injection and withdrawal cycles induces a hysteresis effect, necessitating thorough investigation [32].Here, we investigated the impact of gas hysteresis and its influence on hydrogen recovery.The operation comprises four complete cycles conducted over four years, followed by an extended production phase until depletion.The results indicate that the hysteresis effect exerts a substantial impact on recovery, leading to a decrease in the H2 recovery factor from 76.5% by the fourth cycle for Case 1 to a hydrogen recovery factor of 45% when accounting for hysteresis in Case 2, as illustrated in Figure 9.The gas hysteresis effect led to an extended duration required to produce the majority of the mobile gas, which was two years more than the duration required to extract the majority of H 2 in Case 1. Figure 10 depicts the cumulative mass of injected H 2 over the initial four years with a black line, while the green and dashed blue lines represent the cumulative production for Case 1 and Case 2, respectively.Furthermore, the cumulative production curve for Case 1 demonstrates a plateau by 2031, indicating nearly five years of prolonged production before depletion is reached.In contrast, Case 2 only achieves a depletion state by 2033, with a 5% decrease in cumulative H 2 production relative to Case 1.

PÚBLICA
The gas hysteresis effect led to an extended duration required to produce the majority of the mobile gas, which was two years more than the duration required to extract the majority of H2 in Case 1. Figure 10 depicts the cumulative mass of injected H2 over the initial four years with a black line, while the green and dashed blue lines represent the cumulative production for Case 1 and Case 2, respectively.Furthermore, the cumulative production curve for Case 1 demonstrates a plateau by 2031, indicating nearly five years of prolonged production before depletion is reached.In contrast, Case 2 only achieves a depletion state by 2033, with a 5% decrease in cumulative H2 production relative to Case 1.This phenomenon is attributed to the hysteresis effect, resulting in a higher residual gas per grid block and a wider distribution of hydrogen around the wellbore, thereby increasing the difficulty of subsequent recovery operations.This is exemplified through a comparative analysis of gas saturation near the wellbore between Cases 1 and 2, as shown in Figure 11.This phenomenon is attributed to the hysteresis effect, resulting in a higher residual gas per grid block and a wider distribution of hydrogen around the wellbore, thereby increasing the difficulty of subsequent recovery operations.This is exemplified through a comparative analysis of gas saturation near the wellbore between Cases 1 and 2, as shown in Figure 11.

PÚBLICA
The gas hysteresis effect led to an extended duration required to produce the majority of the mobile gas, which was two years more than the duration required to extract the majority of H2 in Case 1. Figure 10 depicts the cumulative mass of injected H2 over the initial four years with a black line, while the green and dashed blue lines represent the cumulative production for Case 1 and Case 2, respectively.Furthermore, the cumulative production curve for Case 1 demonstrates a plateau by 2031, indicating nearly five years of prolonged production before depletion is reached.In contrast, Case 2 only achieves a depletion state by 2033, with a 5% decrease in cumulative H2 production relative to Case 1.This phenomenon is attributed to the hysteresis effect, resulting in a higher residual gas per grid block and a wider distribution of hydrogen around the wellbore, thereby increasing the difficulty of subsequent recovery operations.This is exemplified through a comparative analysis of gas saturation near the wellbore between Cases 1 and 2, as shown in Figure 11.Consequently, the entrapment of a portion of the H 2 gas within smaller pores results in more space for the water to occupy within the larger pores.This configuration facilitates a more effortless flow of the wetting phase, leading to an increased relative permeability curve for water during the imbibition process, which is in agreement with what has been seen in the literature for CO 2 /brine systems [5].Therefore, larger volumes of water are produced during withdrawal cycles.Figure 12 demonstrates the 12.5% increase in water volume by including the hysteresis effect compared to its neglect, providing empirical evidence that supports our elucidation.

PÚBLICA
Consequently, the entrapment of a portion of the H2 gas within smaller pores results in more space for the water to occupy within the larger pores.This configuration facilitates a more effortless flow of the wetting phase, leading to an increased relative permeability curve for water during the imbibition process, which is in agreement with what has been seen in the literature for CO2/brine systems [5].Therefore, larger volumes of water are produced during withdrawal cycles.Figure 12 demonstrates the 12.5% increase in water volume by including the hysteresis effect compared to its neglect, providing empirical evidence that supports our elucidation.

The Impact of Water Relative Permeability Hysteresis
The water relative permeability hysteresis has been neglected in most CCS numerical simulation studies because water exhibited non-hysteretic behavior [33,43,66,67,68,69].Contrary to the previous observations, recent experiments and pore-scale studies investigating the fluid behavior of hydrogen-brine systems yielded divergent results.The residual gas saturation of hydrogen has been observed to be relatively higher compared to other studied gases.This has a significant effect on the endpoint of the water relative permeability curve during the imbibition process, which in turn affects the hysteresis behavior.Figure 13 exhibits various H2/brine relative permeability for gas/water systems including drainage and imbibition curves [49][50][51][52]54,70].
As part of our investigation, we undertook a comparative study to assess the influence of water hysteresis by analyzing Cases 2 and 3. Case 2, which was adapted from the previous section, accounted for gas hysteresis only.Case 3, however, incorporates both gas and water hysteresis.We aimed to assess the impact of water hysteresis on the overall storage performance and determine whether it is essential or can be safely neglected for UHS numerical simulation assessments.Our results indicated that the inclusion of water hysteresis resulted in a 21% reduction in H2 recovery by the end of the fourth cycle compared to Case 2. As a result, the recovery factor will decrease from 45% to 35% by the year 2027.Figure 14 depicts the difference in cumulative H2 recovery as a result of including water and gas hysteresis versus including gas hysteresis only.The decline in the recovery factor can be primarily attributed to the reduction in bottom-hole pressure, which eventually reached the constraint limit and, consequently, we shut-in the producer.
On the other hand, a notable impact on cumulative water production was observed.Figure 15 displays a comparison of the cumulative water production, which revealed a 37% decrease in water volume with the inclusion of water and gas hysteresis as opposed to including the gas hysteresis solely.
Reaching the constraint of the bottom-hole pressure earlier in the case that included water and gas hysteresis (Case 3) indicates that water hysteresis introduces some difficulties in the fluid flow that may substantially impact the storage performance.The

The Impact of Water Relative Permeability Hysteresis
The water relative permeability hysteresis has been neglected in most CCS numerical simulation studies because water exhibited non-hysteretic behavior [33,43,[66][67][68][69]. Contrary to the previous observations, recent experiments and pore-scale studies investigating the fluid behavior of hydrogen-brine systems yielded divergent results.The residual gas saturation of hydrogen has been observed to be relatively higher compared to other studied gases.This has a significant effect on the endpoint of the water relative permeability curve during the imbibition process, which in turn affects the hysteresis behavior.Figure 13 exhibits various H 2 /brine relative permeability for gas/water systems including drainage and imbibition curves [49][50][51][52]54,70].
Hydrogen 2024, 5, FOR PEER REVIEW 15 simulation was then extended with prolonged production until depletion was reached, and the results were compared with Case 1 and Case 2, as depicted in Figure 16.The black line represents the cumulative injection, the dashed green line represents the cumulative production from Case 1, the dashed blue line represents the cumulative production from Case 2, and the solid red line is the cumulative production from Case 3. The results illustrated in Figure 16 indicate that the inclusion of water and gas hysteresis leads to a rise in the amount of trapped H2, consequently lowering the recovery factor of the stored gas.The recovery percentages, calculated at the time of depletion for Case 1, Case 2, and Case 3, were found to be 94%, 90%, and 88%, respectively.It is therefore evident that water hysteresis should be considered when assessing hydrogen storage feasibility in saline aquifers, and neglecting this parameter can result in an overestimation of H2 recovery.As part of our investigation, we undertook a comparative study to assess the influence of water hysteresis by analyzing Cases 2 and 3. Case 2, which was adapted from the previous section, accounted for gas hysteresis only.Case 3, however, incorporates both gas and water hysteresis.We aimed to assess the impact of water hysteresis on the overall storage performance and determine whether it is essential or can be safely neglected for UHS numerical simulation assessments.Our results indicated that the inclusion of water hysteresis resulted in a 21% reduction in H 2 recovery by the end of the fourth cycle compared to Case 2. As a result, the recovery factor will decrease from 45% to 35% by the year 2027.Figure 14 depicts the difference in cumulative H 2 recovery as a result of including water and gas hysteresis versus including gas hysteresis only.The decline in the recovery factor can be primarily attributed to the reduction in bottom-hole pressure, which eventually reached the constraint limit and, consequently, we shut-in the producer.On the other hand, a notable impact on cumulative water production was observed.Figure 15 displays a comparison of the cumulative water production, which revealed a 37% decrease in water volume with the inclusion of water and gas hysteresis as opposed to including the gas hysteresis solely.

Impact of Relative Permeability Interpolation
The reported relative permeability curves for various gases and brine are distinct for Reaching the constraint of the bottom-hole pressure earlier in the case that included water and gas hysteresis (Case 3) indicates that water hysteresis introduces some difficulties in the fluid flow that may substantially impact the storage performance.The simulation was then extended with prolonged production until depletion was reached, and the results were compared with Case 1 and Case 2, as depicted in Figure 16.The black line represents the cumulative injection, the dashed green line represents the cumulative production from Case 1, the dashed blue line represents the cumulative production from Case 2, and the solid red line is the cumulative production from Case 3.
The results illustrated in Figure 16 indicate that the inclusion of water and gas hysteresis leads to a rise in the amount of trapped H 2 , consequently lowering the recovery factor of the stored gas.The recovery percentages, calculated at the time of depletion for Case 1, Case 2, and Case 3, were found to be 94%, 90%, and 88%, respectively.It is therefore evident that water hysteresis should be considered when assessing hydrogen storage feasibility in saline aquifers, and neglecting this parameter can result in an overestimation of H 2 recovery.

Impact of Relative Permeability Interpolation
The reported relative permeability curves for various gases and brine are distinct for each gas and exhibit differences from one another.Therefore, accounting for the mobility of each gas and interpolating the relative permeability for the mixture of gases is important for UHS analysis.We employed the interpolation function with two sets of relative permeability curves that are associated with the utilized cushion gas.After a comprehensive study of the three suggested cushion gases, a consistent pattern emerged.Accordingly, for the sake of brevity and clarity, we will provide a detailed analysis of the case scenario utilizing CH4 as the cushion gas, since the findings in this case are representative across all three cases.To assess the significance of incorporating the interpolation function, we compared two cases; the first case included only the hydrogen-brine relative permeability, and the second case included two sets of relative permeability curves, and a linear interpolation method was applied.Investigating the volume of the recovered hydrogen, the recovered cushion gas, H2 impurity, and water production were the main criteria to evaluate the potential impact of having multiple relative permeability curves.The volume of water production and the volume of the recovered cushion gas exhibited the most substantial impact by employing the interpolation function, as seen in Figures 17  and 18.The implementation of the interpolation function yielded a 10% increase in water production volume compared to the base scenario.This indicates that the relative permeability to water is subject to differ when we consider the mixture of gases in the reservoir.

Results of Cushion Gas Assessment 6.1. Impact of Relative Permeability Interpolation
The reported relative permeability curves for various gases and brine are distinct for each gas and exhibit differences from one another.Therefore, accounting for the mobility of each gas and interpolating the relative permeability for the mixture of gases is important for UHS analysis.We employed the interpolation function with two sets of relative permeability curves that are associated with the utilized cushion gas.After a comprehensive study of the three suggested cushion gases, a consistent pattern emerged.Accordingly, for the sake of brevity and clarity, we will provide a detailed analysis of the case scenario utilizing CH 4 as the cushion gas, since the findings in this case are representative across all three cases.To assess the significance of incorporating the interpolation function, we compared two cases; the first case included only the hydrogen-brine relative permeability, and the second case included two sets of relative permeability curves, and a linear interpolation method was applied.Investigating the volume of the recovered hydrogen, the recovered cushion gas, H 2 impurity, and water production were the main criteria to evaluate the potential impact of having multiple relative permeability curves.The volume of water production and the volume of the recovered cushion gas exhibited the most substantial impact by employing the interpolation function, as seen in Figures 17 and 18.The implementation of the interpolation function yielded a 10% increase in water production volume compared to the base scenario.This indicates that the relative permeability to water is subject to differ when we consider the mixture of gases in the reservoir.Relative permeability influences the distribution of fluids within the reservoir and reflects how easily each phase can flow.Therefore, our results reported a more accurate estimation of H2 and cushion gas recovery when we accounted for the distribution and movement of the mixture of gases, rather than assuming a single relative permeability curve.When we compared the amount of cumulative CH4 mass that we have produced over the years, we observed a 10% decrease in the modified case that takes into account the interpolation of relative permeability curves compared to the base case.This result indicates that the mobility of the cushion gas was better controlled when we assigned the correct relative permeability curves, and consequently, producing less cushion gas will lead to a higher purity of hydrogen, which is a positive outcome.
The cumulative H2 mass recovered over five years was the least impacted by employing the interpolation function.However, an increment of 5% was observed in the modified case.This indicates that accounting for the mixture of CH4-H2 properly improved the displacement efficiency.Another piece of evidence that supports this finding is the lower bottom-hole pressure observed during injection and production in the modified case scenario, as seen in Figures 19 and 20.This is attributed to the increase in water displacement efficiency.The fact that we produced larger volumes of water led to increasing the storage capacity and improved the overall operation.This implies that using an alternative cushion gas not only serves as a denser layer of gas that prevents excessive H2 spreading and loss but also acts as a displacement agent that enhances the efficiency of water Relative permeability influences the distribution of fluids within the reservoir and reflects how easily each phase can flow.Therefore, our results reported a more accurate estimation of H 2 and cushion gas recovery when we accounted for the distribution and movement of the mixture of gases, rather than assuming a single relative permeability curve.When we compared the amount of cumulative CH 4 mass that we have produced over the years, we observed a 10% decrease in the modified case that takes into account the interpolation of relative permeability curves compared to the base case.This result indicates that the mobility of the cushion gas was better controlled when we assigned the correct relative permeability curves, and consequently, producing less cushion gas will lead to a higher purity of hydrogen, which is a positive outcome.
The cumulative H 2 mass recovered over five years was the least impacted by employing the interpolation function.However, an increment of 5% was observed in the modified case.This indicates that accounting for the mixture of CH 4 -H 2 properly improved the displacement efficiency.Another piece of evidence that supports this finding is the lower bottom-hole pressure observed during injection and production in the modified case scenario, as seen in Figures 19 and 20.This is attributed to the increase in water displacement efficiency.The fact that we produced larger volumes of water led to increasing the storage capacity and improved the overall operation.This implies that using an alternative cushion gas not only serves as a denser layer of gas that prevents excessive H 2 spreading and loss but also acts as a displacement agent that enhances the efficiency of water displacement and optimizes the storage capacity, leading to improved performance and operational outcomes.
Hydrogen 2024, 5, FOR PEER REVIEW 18 displacement and optimizes the storage capacity, leading to improved performance and operational outcomes.Another result of implementing the interpolation of relative permeability curves was observed in the lateral spreading of the cushion gas and the H2 across the reservoir.Our results indicated that the distribution of the cushion gas and the working gas was better controlled and therefore resulted in the mitigation of H2 spreading, leading to the reduction in H2 plume size.As seen in Figure 21, the size of the H2 plume was reduced by 300 m laterally, demonstrating that ignoring the application of the interpolation function will overestimate the size of the H2 plume.Another result of implementing the interpolation of relative permeability curves was observed in the lateral spreading of the cushion gas and the H 2 across the reservoir.Our results indicated that the distribution of the cushion gas and the working gas was better controlled and therefore resulted in the mitigation of H 2 spreading, leading to the reduction in H 2 plume size.As seen in Figure 21, the size of the H 2 plume was reduced by 300 m laterally, demonstrating that ignoring the application of the interpolation function will overestimate the size of the H 2 plume.The reduction in H2 plume size is a favorable outcome in modeling UHS operations, and minimizing the H2 lateral spreading reduces the chance of losing H2 into the reservoir, as reported previously in the literature [71].

Impact of Cushion Gas Type
Underground H2 storage in saline aquifers would require a large volume of cushion gas to prevent H2 loss.Considering the additional cost this will add to the project budget, we, therefore, tested different cushion gases to be used as an alternative to H2 and assessed The reduction in H 2 plume size is a favorable outcome in modeling UHS operations, and minimizing the H 2 lateral spreading reduces the chance of losing H 2 into the reservoir, as reported previously in the literature [71].

Impact of Cushion Gas Type
Underground H 2 storage in saline aquifers would require a large volume of cushion gas to prevent H 2 loss.Considering the additional cost this will add to the project budget, we, therefore, tested different cushion gases to be used as an alternative to H 2 and assessed their impact on H 2 recovery, purity, and storage performance.The analysis revealed that CH 4 and N 2 served as a better cushion gas, leading to the highest H 2 recovery (77%) compared to CO 2 (74%), as shown in Figure 22.Although the three cushion gases performed similarly in terms of H 2 recovery, we observed notable differences in water production, H 2 purity, H 2 distribution, and bottom-hole pressure.We thus provide a detailed screening to compare each cushion gas and the overall impact they have on storage efficiency.Although the gases behave similarly, each gas has different characteristics regarding molecular weight, density, and solubility in water.In our study, we injected the cushion gases at an exact rate (volume/day), yet the volume of the cushion gas in reservoir conditions would differ based on each gas density and solubility.CO2 for example, is known to be highly soluble in water compared to CH 4 , H2, and N2, and exhibits the highest volumetric mass.In turn, this results in a reduction in the cushion gas plume of CO2 at reservoir conditions compared to the other gases, causing more loss of H2 into the formation and allowing more water to flow into the production system.Within a comparable framework, N2 has a higher density compared to CH4 and it is heavier in mass, and therefore, it exhibits less dispersion under the same condition, leading to a reduction in the cushion gas plume and potentially allowing water to flow more readily than if CH4 were used as the cushion gas.As seen in Figure 23, the cumulative water production is the highest in the case of utilizing CO2 as the cushion gas, followed by N2, and CH4, and lastly in the case of H2 used as the cushion gas.The observation of varying dispersion behavior for each cushion gas also explains the distribution behavior of H2 (the working gas) in each case.When H2 itself is used as the cushion gas, the low mass of H2 and its tendency to migrate upward will lead to a massive spreading of the injected gas, and therefore, a larger plume size.On the other hand, between the three other gases (N2, CH4, and CO2), CH4 is likely to exhibit greater Although the gases behave similarly, each gas has different characteristics regarding molecular weight, density, and solubility in water.In our study, we injected the cushion gases at an exact rate (volume/day), yet the volume of the cushion gas in reservoir conditions would differ based on each gas density and solubility.CO 2 for example, is known to be highly soluble in water compared to CH 4, H 2, and N 2, and exhibits the highest volumetric mass.In turn, this results in a reduction in the cushion gas plume of CO 2 at reservoir conditions compared to the other gases, causing more loss of H 2 into the formation and allowing more water to flow into the production system.Within a comparable framework, N 2 has a higher density compared to CH 4 and it is heavier in mass, and therefore, it exhibits less dispersion under the same condition, leading to a reduction in the cushion gas plume and potentially allowing water to flow more readily than if CH 4 were used as the cushion gas.As seen in Figure 23, the cumulative water production is the highest in the case of utilizing CO 2 as the cushion gas, followed by N 2, and CH 4, and lastly in the case of H 2 used as the cushion gas.Although the gases behave similarly, each gas has different characteristics regarding molecular weight, density, and solubility in water.In our study, we injected the cushion gases at an exact rate (volume/day), yet the volume of the cushion gas in reservoir conditions would differ based on each gas density and solubility.CO2 for example, is known to be highly soluble in water compared to CH 4 , H2, and N2, and exhibits the highest volumetric mass.In turn, this results in a reduction in the cushion gas plume of CO2 at reservoir conditions compared to the other gases, causing more loss of H2 into the formation and allowing more water to flow into the production system.Within a comparable framework, N2 has a higher density compared to CH4 and it is heavier in mass, and therefore, it exhibits less dispersion under the same condition, leading to a reduction in the cushion gas plume and potentially allowing water to flow more readily than if CH4 were used as the cushion gas.As seen in Figure 23, the cumulative water production is the highest in the case of utilizing CO2 as the cushion gas, followed by N2, and CH4, and lastly in the case of H2 used as the cushion gas.The observation of varying dispersion behavior for each cushion gas also explains the distribution behavior of H2 (the working gas) in each case.When H2 itself is used as the cushion gas, the low mass of H2 and its tendency to migrate upward will lead to a massive spreading of the injected gas, and therefore, a larger plume size.On the other hand, between the three other gases (N2, CH4, and CO2), CH4 is likely to exhibit greater The observation of varying dispersion behavior for each cushion gas also explains the distribution behavior of H 2 (the working gas) in each case.When H 2 itself is used as the cushion gas, the low mass of H 2 and its tendency to migrate upward will lead to a massive spreading of the injected gas, and therefore, a larger plume size.On the other hand, between the three other gases (N 2, CH 4, and CO 2 ), CH 4 is likely to exhibit greater dispersion, owing to its lower molecular mass.But taking into consideration its relatively higher density and weight compared to hydrogen, CH 4 acted as a barrier that limited the excessive spreading of hydrogen, yielding to the smallest H 2 plume, as seen in Figure 24.Such a circumstance is deemed favorable in underground H 2 storage.One of the important aspects to highlight when comparing the cushion gas performance is the impact on H2 purity.Using an alternative gas to increase the H2 recovery might come at the price of lowering H2 purity.Therefore, targeting the maximum recovery and maintaining an elevated level of H2 purity is the optimum case.Hence, we compared the H2 purity in terms of gas mole fraction, as seen in Figure 25, and our results indicated that among the gases considered, H2 exhibited the highest purity, followed by CH4 and N2.CO2, however, reduced the purity of the recovered H2 during the first four cycles, but its performance converged with that of the other gases later.Multiple factors have influenced these outcomes.The relative permeability of CO2 surpasses that of CH4 and N2, indicating an anticipation of higher CO2 flow than CH4 or N2 under equivalent water saturation conditions.The other reason is related to the size of the cushion gas plume in reservoir conditions.As previously noted, CO2 has a higher density than the rest of the proposed gases, mitigating dispersion effects and constraining its distribution.Consequently, under an equivalent rate of injected cushion gas, the thickness of the cushion gas plume reduces when CO2 is employed, leading to an increase in cushion gas concentration per grid block around the wellbore.However, by increasing the volume of the working gas over the cycles, the displacement of the cushion gas away from the wellbore takes place, hereby enhancing H2 purity across all scenarios.One of the important aspects to highlight when comparing the cushion gas performance is the impact on H 2 purity.Using an alternative gas to increase the H 2 recovery might come at the price of lowering H 2 purity.Therefore, targeting the maximum recovery and maintaining an elevated level of H 2 purity is the optimum case.Hence, we compared the H 2 purity in terms of gas mole fraction, as seen in Figure 25, and our results indicated that among the gases considered, H 2 exhibited the highest purity, followed by CH 4 and N 2 .CO 2 , however, reduced the purity of the recovered H 2 during the first four cycles, but its performance converged with that of the other gases later.Multiple factors have influenced these outcomes.The relative permeability of CO 2 surpasses that of CH 4 and N 2 , indicating an anticipation of higher CO 2 flow than CH 4 or N 2 under equivalent water saturation conditions.The other reason is related to the size of the cushion gas plume in reservoir conditions.As previously noted, CO 2 has a higher density than the rest of the proposed gases, mitigating dispersion effects and constraining its distribution.Consequently, under an equivalent rate of injected cushion gas, the thickness of the cushion gas plume reduces when CO 2 is employed, leading to an increase in cushion gas concentration per grid block around the wellbore.However, by increasing the volume of the working gas over the cycles, the displacement of the cushion gas away from the wellbore takes place, hereby enhancing H 2 purity across all scenarios.
quently, under an equivalent rate of injected cushion gas, the thickness of the cushion gas plume reduces when CO2 is employed, leading to an increase in cushion gas concentration per grid block around the wellbore.However, by increasing the volume of the working gas over the cycles, the displacement of the cushion gas away from the wellbore takes place, hereby enhancing H2 purity across all scenarios.Thus, the notable increase in H2 recovery when using an alternative cushion gas is counterbalanced by a compromise in H2 purity.However, as seen in Figure 26, targeting the maximum purity of H2 by utilizing H2 itself as the cushion gas resulted in a 63% Thus, the notable increase in H 2 recovery when using an alternative cushion gas is counterbalanced by a compromise in H 2 purity.However, as seen in Figure 26, targeting the maximum purity of H 2 by utilizing H 2 itself as the cushion gas resulted in a 63% recovery factor of the working gas compared to 74-77% when other alternative cushion gases are used.The relevant literature suggested that alternative cushion gases such as CH 4 and CO 2 have experimentally shown a decreasing effect on H 2 adsorption in water, therefore preventing H 2 loss [72].
Hydrogen 2024, 5, FOR PEER REVIEW 22 recovery factor of the working gas compared to 74-77% when other alternative cushion gases are used.The relevant literature suggested that alternative cushion gases such as CH4 and CO2 have experimentally shown a decreasing effect on H2 adsorption in water, therefore preventing H2 loss [72].Figure 27 shows the bottom-hole pressure profile for the producer during the cyclic injection and production period.It is noteworthy that the utilization of N2 and CO2 as cushion gases resulted in the most significant pressure decline.This phenomenon can be attributed to the maximum production of water associated with these two scenarios, as previously mentioned.Thus, we specified the main criteria to assess and identify the optimum cushion gas.The criteria include maximizing H2 production, minimizing water production, maintaining operational pressure, and sustaining elevated levels of H2 purity.Based on the out- Figure 27 shows the bottom-hole pressure profile for the producer during the cyclic injection and production period.It is noteworthy that the utilization of N 2 and CO 2 as cushion gases resulted in the most significant pressure decline.This phenomenon can be attributed to the maximum production of water associated with these two scenarios, as previously mentioned.
Hydrogen 2024, 5, FOR PEER REVIEW 22 recovery factor of the working gas compared to 74-77% when other alternative cushion gases are used.The relevant literature suggested that alternative cushion gases such as CH4 and CO2 have experimentally shown a decreasing effect on H2 adsorption in water, therefore preventing H2 loss [72].Figure 27 shows the bottom-hole pressure profile for the producer during the cyclic injection and production period.It is noteworthy that the utilization of N2 and CO2 as cushion gases resulted in the most significant pressure decline.This phenomenon can be attributed to the maximum production of water associated with these two scenarios, as previously mentioned.Thus, we specified the main criteria to assess and identify the optimum cushion gas.The criteria include maximizing H2 production, minimizing water production, maintaining operational pressure, and sustaining elevated levels of H2 purity.Based on the outcomes, CH4 emerges as the optimal choice, followed by N2 and subsequently CO2.Thus, we specified the main criteria to assess and identify the optimum cushion gas.The criteria include maximizing H 2 production, minimizing water production, maintaining operational pressure, and sustaining elevated levels of H 2 purity.Based on the outcomes, CH 4 emerges as the optimal choice, followed by N 2 and subsequently CO 2 .

Conclusions
This study presented some numerical outcomes to quantify the hydrogen losses due to wettability and relative permeability hysteresis for aqueous and gaseous phases during underground storage in saline aquifers.Regarding the mechanism of H 2 loss due to wettability and hysteresis, we have the following conclusions:

•
Our investigation into underground hydrogen storage revealed that the variation in wettability has a limited impact on the efficiency of the storage system.The water-wet condition resulted in a recovery factor that is 8% greater than that achieved under the neutrally wet condition.

•
Furthermore, the results indicated that the hysteresis effect exerted the most substantial influence on the recovery of hydrogen.Our results suggested that accounting for the gas hysteresis effect resulted in a substantial drop in recovery from 78% to 45% by the fourth cycle.The inclusion of water relative permeability hysteresis, however, resulted in a notable reduction in both hydrogen recovery and water production for each cycle, which suggests that accounting for water hysteresis in UHS simulation studies is important to have a precise estimation of the recovery factor and to predict the operational statue accurately.

•
Our analysis indicates the necessity for further experimental studies measuring the relative permeability of the H 2 /brine system during drainage and imbibition processes.
It is imperative to acknowledge that the residual gas saturation data reported from the laboratory experiments usually exhibit a substantial margin of error and that consequently influences the accuracy of the utilized data.Additional experimental investigations will enhance our comprehensive knowledge of hydrogen behavior and maximize the precision of the simulation results to assess the viability of the storage projects.Furthermore, no published data related to hysteresis in a weak waterwet/neutrally wet formation was available.However, based on our results, a marginal increase in the trapped gas saturation is anticipated, with no significant impact on hydrogen recovery.
To overcome part of the hydrogen losses, findings highlight the significance of injecting cushion gas during H 2 storage in saline aquifers.The results not only underscore the significance of this practice but also advance our understanding of the diverse impacts associated with the use of various cushion gases on H 2 recovery and storage performance.It is worth pointing out the following:

•
Based on our research outcomes, CH 4 emerges as the optimal choice, followed by N 2, and subsequently CO 2 .The results indicated that the relative permeabilities for different gases and their mixture have a significant impact on both H 2 recovery and water production, therefore we emphasize the importance of considering the associated relative permeability curves for all the injected gases in case of utilizing a different cushion gas than H 2 .

•
Furthermore, the utilization of CH 4 as a cushion gas not only enhances H 2 recovery but also mitigates H 2 plume spreading, thereby creating an optimal condition for UHS.Lastly, we found that increasing the cushion gas volume, beyond the required amount to maintain the aquifer pressure, improves the recovery in the first cycle but its impact tapers off over a longer duration.
This study represents a significant advancement in the field by being the first to incorporate a precise fluid model that accounts for the fluid behavior of the mixture of gases and addresses the fluid dynamics with high precision through the interpolation of relative permeability.Additionally, we conduct a thorough investigation of the impact of hysteresis and wettability, providing comprehensive insights that have not been previously explored.Our novel approach to model underground hydrogen storage simulations in saline aquifers enhances the reliability of our findings and establishes a benchmark for subsequent studies in the field.

Figure 1 .
Figure 1.The horizontal permeability distribution for the model used in this study in mD.

Figure 1 .
Figure 1.The horizontal permeability distribution for the model used in this study in mD.
N , nonwetting phase saturation; S Nr , residual or trapped water saturation;

Figure 6 .
Figure 6.Cumulative hydrogen injection in brown and production based on the wettability condition: water-wet (contact angle = 35°) in red and neutrally wet (contact angle = 83°) in dashed blue.

Figure 7 .
Figure 7. Cumulative water production in strongly water-wet (in red) and neutrally wet (in blue) conditions.

Figure 6 .
Figure 6.Cumulative hydrogen injection in brown and production based on the wettability condition: water-wet (contact angle = 35 • ) in red and neutrally wet (contact angle = 83 • ) in dashed blue.

Figure 6 .
Figure 6.Cumulative hydrogen injection in brown and production based on the wettability condition: water-wet (contact angle = 35°) in red and neutrally wet (contact angle = 83°) in dashed blue.

Figure 7 .
Figure 7. Cumulative water production in strongly water-wet (in red) and neutrally wet (in blue) conditions.

Figure 7 .
Figure 7. Cumulative water production in strongly water-wet (in red) and neutrally wet (in blue) conditions.

Figure 8 .
Figure 8.Comparison of injector BHP at different wettability conditions: water-wet in red and neutrally wet in blue.

Figure 9 .
Figure 9. Cumulative hydrogen mass for injection (in brown) and production cycles considering hysteresis (in orange) and with no hysteresis (in green).

Figure 8 .
Figure 8.Comparison of injector BHP at different wettability conditions: water-wet in red and neutrally wet in blue.

Hydrogen 2024, 5 , 12 ÚBLICAFigure 8 .
Figure 8.Comparison of injector BHP at different wettability conditions: water-wet in red and neutrally wet in blue.

Figure 9 .
Figure 9. Cumulative hydrogen mass for injection (in brown) and production cycles considering hysteresis (in orange) and with no hysteresis (in green).

Figure 9 .
Figure 9. Cumulative hydrogen mass for injection (in brown) and production cycles considering hysteresis (in orange) and with no hysteresis (in green).

Figure 10 .
Figure 10.Comparison of the cumulative gas mass.The black line represents the cumulative injection, the green line represents the cumulative production from Case 1, and the dashed blue line is the cumulative production from Case 2.

Figure 11 .
Figure 11.Comparison of H2 saturation profiles with the inclusion (B) and neglect of gas hysteresis (A).

Figure 10 .
Figure 10.Comparison of the cumulative gas mass.The black line represents the cumulative injection, the green line represents the cumulative production from Case 1, and the dashed blue line is the cumulative production from Case 2.

Figure 10 .
Figure 10.Comparison of the cumulative gas mass.The black line represents the cumulative injection, the green line represents the cumulative production from Case 1, and the dashed blue line is the cumulative production from Case 2.

Figure 11 .
Figure 11.Comparison of H2 saturation profiles with the inclusion (B) and neglect of gas hysteresis (A).

Figure 11 .
Figure 11.Comparison of H 2 saturation profiles with the inclusion (B) and neglect of gas hysteresis (A).

Figure 12 .
Figure 12.Comparison of cumulative water production influenced by the neglect (in green) and in the presence of hysteresis (in orange).

Figure 12 .
Figure 12.Comparison of cumulative water production influenced by the neglect (in green) and in the presence of hysteresis (in orange).

Figure 14 .
Figure 14.Comparison of cumulative hydrogen production with (in blue) and without water hysteresis (in red).

Figure 14 .
Figure 14.Comparison of cumulative hydrogen production with (in blue) and without water hysteresis (in red).

Hydrogen 2024, 5 , 16 Figure 15 .
Figure 15.Comparison of cumulative water production with (in blue) and without water hysteresis (in red).

Figure 16 .
Figure 16.Comparison of cumulative gas mass as an influence of the inclusion of water and gas hysterias, and the neglect of both.

Figure 15 .
Figure 15.Comparison of cumulative water production with (in blue) and without water hysteresis (in red).

Figure 15 .
Figure 15.Comparison of cumulative water production with (in blue) and without water hysteresis (in red).

Figure 16 .
Figure 16.Comparison of cumulative gas mass as an influence of the inclusion of water and gas hysterias, and the neglect of both.

Figure 16 .
Figure 16.Comparison of cumulative gas mass as an influence of the inclusion of water and gas hysterias, and the neglect of both.

Figure 19 .
Figure 19.Bottom-hole pressure for the injector well for the modified case with Kr curve interpolation (dashed red line) and the base case (blue solid lines).

Figure 19 .Figure 19 .
Figure 19.Bottom-hole pressure for the injector well for the modified case with Kr curve interpolation (dashed red line) and the base case (blue solid lines).

Figure 20 .
Figure 20.Bottom-hole pressure for the injector and producer well for the modified case with Kr curve interpolation (dashed red line) and the base case (blue solid lines).

Figure 20 .
Figure 20.Bottom-hole pressure for the injector and producer well for the modified case with Kr curve interpolation (dashed red line) and the base case (blue solid lines).

Hydrogen 2024, 5 ,Figure 21 .
Figure 21.H2 saturation profiles displaying the lateral spreading of H2 plume in (a) base case and (b) modified case.

Figure 21 .
Figure 21.H 2 saturation profiles displaying the lateral spreading of H 2 plume in (a) base case and (b) modified case.

Figure 23 .
Figure 23.The cumulative volume of water production when using various types of cushion gas (namely, CO2 in green, H2 in blue, N2 in red, and CH 4 in black ).

Figure 22 .
Figure 22.Cumulative mass of H 2 injection (in dark red) and production while utilizing different cushion gases (CO 2 in green, CH 4 in black, and N 2 in red).

Figure 23 .
Figure 23.The cumulative volume of water production when using various types of cushion gas (namely, CO2 in green, H2 in blue, N2 in red, and CH 4 in black ).

Figure 23 .
Figure 23.The cumulative volume of water production when using various types of cushion gas (namely, CO 2 in green, H 2 in blue, N 2 in red, and CH 4 in black).

Hydrogen 2024, 5 , 21 Figure 24 .
Figure 24.Top view of the H2 plume in different scenarios using different types of cushion gas, reported in H2 mole fraction.

Figure 25 .
Figure 25.H2 purity in terms of gas mole fraction during the first five cycles using different cushion

Figure 24 .
Figure 24.Top view of the H 2 plume in different scenarios using different types of cushion gas, reported in H 2 mole fraction.

Figure 25 .
Figure 25.H2 purity in terms of gas mole fraction during the first five cycles using different cushion gases (CO2 in green, H2 in blue, N2 in red, and CH4 in black).

Figure 25 .
Figure 25.H 2 purity in terms of gas mole fraction during the first five cycles using different cushion gases (CO 2 in green, H 2 in blue, N 2 in red, and CH 4 in black).

Figure 26 .
Figure 26.Base case scenario of hydrogen production (in blue) utilizing H2 as the cushion gas.Results indicated a 63% recovery factor of the working gas mass.

Figure 27 .
Figure 27.Well bottom-hole pressure at the producer for various scenarios using different cushion gases (CO2 in green, H2 in blue, N2 in red, and CH4 in black ).

Figure 26 .
Figure 26.Base case scenario of hydrogen production (in blue) utilizing H 2 as the cushion gas.Results indicated a 63% recovery factor of the working gas mass.

Figure 26 .
Figure 26.Base case scenario of hydrogen production (in blue) utilizing H2 as the cushion gas.Results indicated a 63% recovery factor of the working gas mass.

Figure 27 .
Figure 27.Well bottom-hole pressure at the producer for various scenarios using different cushion gases (CO2 in green, H2 in blue, N2 in red, and CH4 in black ).

Figure 27 .
Figure 27.Well bottom-hole pressure at the producer for various scenarios using different cushion gases (CO 2 in green, H 2 in blue, N 2 in red, and CH 4 in black).

Table 1 .
Summary of hydrogen storage technologies including surface and underground methods.