Large-scale geological storages of hydrogen (H₂) and carbon dioxide (CO₂) in saline aquifers present feasible options for a sustainable energy future. We compared the plume migration of CO₂ and H₂ in aquifers using the FluidFlower benchmark, incorporating the state-of-the-art thermophysical and petrophysical properties. The H₂ plume, with its higher buoyancy and mobility compared to CO₂, remains predominantly in the gas phase due to its lower solubility, increasing the chances of escaping through fractures or migration to distant regions. This additionally leads to a higher pressurized reservoir, which, along with higher buoyancy, increases the chance of caprock penetration. Dissolution trapping of CO₂ into brine increases over time due to its fingering, while H₂ does not show fingering. Our findings show that while geological carbon storage (GCS) benefits significantly from all structural, dissolution, and residual trapping, underground hydrogen storage (UHS) relies mainly on structural trapping, making the integrity of sealing elements of the system a key factor in its performance.

Thermal-Hydro-Mechanical-Compositional analysis is crucial for addressing challenges like wellbore stability, land subsidence, and induced seismicity in the geo-energy applications. Numerical simulations of coupled thermo-poromechanical processes provide a general-purpose tool for evaluating these phenomena across laboratory and field scales. However, efficient integration of the coupled equations for fluid mass, energy and momentum poses multiple numerical and implementation difficulties, such as combining different numerical methods on staggered grids and associated limitations on admissible grids. This paper introduces a novel fully-implicit Finite Volume Method (FVM) for modeling thermal compositional flow in thermo-poroelastic rocks. The scheme employs gradient-based, coupled multi-point approximations of fluid mass, momentum and heat fluxes.

The novelty of the scheme lies in its integration of temperature as a parameter in the flux approximation process. The scheme supports a wide range of cell topologies, arbitrary heterogeneity and anisotropy as well as various boundary conditions, while respecting local flux balance under temperature gradients. Overall, the scheme represents a unified FVM-based approach for the integration of all conservation laws relevant to geo-energy applications on a cell-centered collocated grid. Additionally, the implemented two-stage block-partitioned preconditioning strategy enables the efficient solution of obtained linear systems.

The framework, implemented in the open-source Delft Advanced Research Terra Simulator (open-DARTS), leverages the Operator-Based Linearization (OBL) technique for flexibility in compositional fluid properties. Rigorous validation demonstrates the framework’s capabilities in capturing advanced phenomena, including thermal expansion, thermo-poroelastic effect and compositional flow with phase transitions. The performance of preconditioning strategy is assessed using the mechanical extension of the SPE10 benchmark model.

The activation of fracture networks poses significant risks and raises safety concerns for projects involving such geological structures. Consequently, an accurate and efficient simulation strategy is essential for modeling highly fractured subsurface formations. While the extended finite element method (XFEM), coupled with the penalty method, effectively models slip-stick conditions along fracture surfaces and fracture propagation under compression, its efficiency declines when handling dense fracture networks. To address this challenge, a multiscale XFEM (MS-XFEM) approach is developed and presented. MS-XFEM approximates fine-scale displacement field by interpolating solutions on a coarser-scale mesh using algebraically constructed basis functions. All extra degrees of freedom (DOFs) are incorporated within the basis functions matrix, rendering the coarse-scale system a standard finite element-based system. In each propagation step, basis functions are algebraically and locally updated to capture fracture propagation. Through four proof-of-concept test cases, the accuracy and efficiency of MS-XFEM in simulating fractured geological formations are demonstrated, emphasizing its potential for real-world applications.

To safely and efficiently utilize porous reservoirs for underground hydrogen storage (UHS), it is essential to characterize hydrogen transport properties at multiple scales. In this study, hydrogen/brine multiphase flow at 50 bar and 25 °C in a 17 cm Berea sandstone rock core was characterized and visualized at the pore and core scales using micro X-ray CT. The experiment included a single drainage and imbibition cycle during which relative permeability hysteresis was measured, and two no-flow periods to study the redistribution of hydrogen in the pore space during storage periods. An end-point relative permeability of 0.043 was found at Sw = 0.56, and the residual gas saturation was measured to be 0.32. Despite extensive pre-equilibration, significant dissolution of hydrogen into brine occurred near the core inlet due to elevated pressures and the corresponding increase in hydrogen solubility. During drainage, many disconnected hydrogen ganglia were observed further down the core which could be explained by the exsolution of the dissolved hydrogen. During imbibition, the dissolution of hydrogen led to the formation of preferential flow paths near the inlet, and eventually removed most of the trapped hydrogen in the final stage of the experiment. The two no-flow periods were characterized by the fragmentation of medium-sized hydrogen ganglia and the growth of a few larger ganglia, providing evidence for hydrogen re-connection through the dissolution-driven process of Ostwald ripening. These results demonstrate that despite the low solubility of hydrogen in brine, hydrogen dissolution can significantly influence the observed multiphase flow and trapping behavior in the reservoir and should be considered in UHS modeling.

Hydrogen, regarded as a clean energy carrier, holds potential for subsurface porous media storage to balance energy supply and demand. Investigating the interactions between hydrogen and porous media, especially the potential residual trapping during cyclic injection and production, becomes imperative. Experimental studies were conducted to examine the influence of flow rate and cyclic injection/production on hydrogen storage and recovery efficiency. Findings reveal generally low hydrogen saturation in brine-saturated porous media during storage process, due to notable flow instability. However, higher injection rates contribute to increased hydrogen saturation. Moreover, observations indicate that cyclic hydrogen storage and drainage may lead to hysteresis in hydrogen storage process. That is, cyclic hydrogen storage and drainage may result in undesired residual trapping of hydrogen in porous media. This study underscores the need for a comprehensive investigation to validate these observations, shedding light on the complexities of hydrogen storage dynamics in subsurface porous media.

H₂-CO₂ mixtures find wide-ranging applications, including their growing significance as synthetic fuels in the transportation industry, relevance in capture technologies for carbon capture and storage, occurrence in subsurface storage of hydrogen, and hydrogenation of carbon dioxide to form hydrocarbons and alcohols. Here, we focus on the thermodynamic properties of H₂-CO₂ mixtures pertinent to underground hydrogen storage in depleted gas reservoirs. Molecular dynamics simulations are used to compute mutual (Fick) diffusivities for a wide range of pressures (5 to 50 MPa), temperatures (323.15 to 423.15 K), and mixture compositions (hydrogen mole fraction from 0 to 1). At 5 MPa, the computed mutual diffusivities agree within 5% with the kinetic theory of Chapman and Enskog at 423.15 K, albeit exhibiting deviations of up to 25% between 323.15 and 373.15 K. Even at 50 MPa, kinetic theory predictions match computed diffusivities within 15% for mixtures comprising over 80% H₂ due to the ideal-gas-like behavior. In mixtures with higher concentrations of CO₂, the Moggridge correlation emerges as a dependable substitute for the kinetic theory. Specifically, when the CO₂ content reaches 50%, the Moggridge correlation achieves predictions within 10% of the computed Fick diffusivities. Phase equilibria of ternary mixtures involving CO₂-H₂-NaCl were explored using Gibbs Ensemble (GE) simulations with the Continuous Fractional Component Monte Carlo (CFCMC) technique. The computed solubilities of CO₂ and H₂ in NaCl brine increased with the fugacity of the respective component but decreased with NaCl concentration (salting out effect). While the solubility of CO₂ in NaCl brine decreased in the ternary system compared to the binary CO₂-NaCl brine system, the solubility of H₂ in NaCl brine increased less in the ternary system compared to the binary H₂-NaCl brine system. The cooperative effect of H₂-CO₂ enhances the H₂ solubility while suppressing the CO₂ solubility. The water content in the gas phase was found to be intermediate between H₂-NaCl brine and CO₂-NaCl brine systems. Our findings have implications for hydrogen storage and chemical technologies dealing with CO₂-H₂ mixtures, particularly where experimental data are lacking, emphasizing the need for reliable thermodynamic data on H₂-CO₂ mixtures.

Characterization of the microbial activity impacts on transport and storage of hydrogen is a crucial aspect of successful Underground Hydrogen Storage (UHS). Microbes can use hydrogen for their metabolism, which can then lead to formation of biofilms. Biofilms can potentially alter the wettability of the system and, consequently, impact the flow dynamics and trapping mechanisms in the reservoir. In this study, we investigate the impact of microbial activity on wettability of the hydrogen/brine/rock system, using the captive-bubble cell experimental approach. Apparent contact angles are measured for bubbles of pure hydrogen in contact with a solid surface inside a cell filled with living brine which contains sulphate reducing microbes. To investigate the impact of surface roughness, two different solid samples are used: a "rough" Bentheimer Sandstone sample and a "smooth" pure Quartz sample. It is found that, in systems where buoyancy and interfacial forces are the main acting forces, the impact of biofilm formation on the apparent contact angle highly depends on the surface roughness. For the "rough" Bentheimer sandstone, the apparent contact angle was unchanged by biofilm formation, while for the smooth pure Quartz sample the apparent contact angle decreased significantly, making the system more water-wet. This decrease in apparent contact angle is in contrast with an earlier study present in the literature where a significant increase in contact angle due to microbial activity was reported. The wettability of the biofilm is mainly determined by the consistency of the Extracellular Polymeric Substances (EPS) which depends on the growth conditions in the system. Therefore, to determine the impact of microbial activity on the wettability during UHS will require accurate replication of the reservoir conditions including surface roughness, chemical composition of the brine, the microbial community, as well as temperature, pressure and pH-value conditions.

CO₂ capture and storage is a viable solution in the effort to mitigate global climate change. Deep saline aquifers, in particular, have emerged as promising storage options, owing to their vast capacity and widespread distribution. However, the task of proficiently monitoring and simulating CO₂ behavior within these formations poses significant challenges. To address this, we introduce the physics-constraint neural network for CO₂ storage (CO₂PCNet), a model specifically designed for simulating and monitoring CO₂ storage in deep saline aquifers during injection and post-injection periods. Recognizing the significant challenges in accurately modeling the distribution and movement of CO₂ under varying permeability conditions, the CO₂PCNet integrates the principles of physics with the robustness of deep learning, serving as a powerful surrogate model. The architecture of CO₂PCNet starts with an encoder that adeptly processes spatial features from overall mole fraction (z_CO2) and pressure fields (P_l), capturing the complex dynamics of a CO₂ trajectory. By incorporating permeability information through a conditioning step, the network ensures a faithful representation of the influences on CO₂ behavior in subsurface conditions. A ConvLSTM module subsequently discerns temporal evolutions, reflecting the real-world progression of CO₂ plumes within the reservoir. Lastly, the decoder precisely reconstructs the predictive spatial profile of CO₂ distribution. CO₂PCNet, with its integration of convolutional layers, recurrent mechanisms, and physics-informed constraints, offers a refined approach to CO₂ storage simulation. This model offers the potential of utilizing advanced computational methods in advancing CCS practices.

The modeling of fluid flow in heterogeneous, anisotropic and fractured porous media is relevant in many applications, including hydrocarbon and groundwater extraction, dispersion of contaminants, hydrogen or carbon dioxide (CO₂) storage. Thus, accurate and scalable simulation of fluid flow through these formations continues to be a great challenge. The presence of fractures that are explicitly modeled, ranging from flow barriers to highly conductive channels, significantly increases the complexity of the numerical simulation. The Embedded Discrete Fracture Model (EDFM) produces results that can be as accurate as those obtained using equidimensional Discrete Fracture Models (DFM) at a much lower computational cost for highly conductive fractures, but it is not adequate for impermeable barriers. In contrast, the pEDFM (projection-based Embedded Discrete Fracture Model) produces accurate results for both, channels and barriers by using additional matrix-fracture connectivities. In its current versions, it is restricted to cartesian and corner-point grid geometries and to the classical Two Point Flux Approximation (TPFA) scheme. In this work, for the first time, the pEDFM is extended to handle unstructured tetrahedral meshes (pEDFM-U). In addition, interface fluxes are approximated by using the Multipoint Flux Approximation method with a Diamond stencil (MPFA-D). This allows for simulation of highly heterogeneous and anisotropic geo-models. The developed method is robust and can deal with full permeability tensors on arbitrary tetrahedral meshes. The advective terms are discretized considering an implicit First Order Upwind (FOU) method. Our method is implemented within the DARSim (Delft Advanced Reservoir Simulation) open-source simulator framework. Through several test cases, the proposed method showed to be robust and capable to accurately capture the effects of both high and low permeability fractures for general tetrahedral meshes, under arbitrary heterogeneous and anisotropic permeability tensors for the rock matrix.

The role of Thermal-Hydro-Mechanical-Compositional analysis in the development of geo-energy resources has been amplified in recent years. As an example, challenges such as wellbore stability, land subsidence and induced seismicity highlight the necessity for comprehensive geomechanical evaluations which are then coupled with thermo-hydrodynamical processes within the reservoir. Numerical simulations of the coupled thermo-poromechanical processes provide a general-purpose tool capable of performing these evaluations at both continuum laboratory and field scales. However, efficient integration of the coupled system of fluid mass, energy and momentum conservation equations poses multiple numerical and implementation difficulties, such as combining different numerical methods on staggered grids and associated limitations on admissible grids. This paper introduces a new fully-implicit scheme of the Finite Volume Method (FVM) for modeling thermal compositional flow in thermo-poroelastic rocks. The scheme uses the gradient-based variant of coupled multi-point approximations of fluid mass, momentum, heat convection and conduction fluxes, which are derived from their respective local balances. The novelty of the scheme is that it incorporates temperature into the approximation of these fluxes. Consequently, the approximation of displacement gradients depends on temperatures, while the approximation of temperature itself is derived from the balance of heat conduction fluxes. At the same time, we utilize a single-point upstream weighting for the temperature-dependent terms in heat convection fluxes. The resulting scheme respects the local balance of fluxes in the presence of temperature gradients. Besides, it also supports star-shaped and various boundary conditions. Overall, the scheme represents a unified FVM-based approach for the integration of all conservation laws relevant to geo-energy applications on a cell-centered collocated grid. Furthermore, the implemented two-stage block-partitioned preconditioning strategy enables the efficient solution of obtained linear systems. The proposed modeling framework has been implemented in an open-source Delft Advanced Research Terra Simulator (DARTS). Moreover, the flexibility regarding compositional fluid properties is reinforced by the Operator-Based Linearization (OBL) technique incorporated into DARTS. The proposed modeling framework has undergone rigorous validation in convergence study, and comparisons against established analytical and numerical solutions. The framework covers advanced physical phenomena including thermal expansion and contraction, porosity dependent on pressure, temperature and strain, and multiphase flow with phase changes and chemical alterations. The framework capabilities and the performance of the preconditioning strategy have been assessed in the mechanical extension of the 10th SPE Comparative study (SPE10) model.

Thermodynamic factors for diffusion connect the Fick and Maxwell-Stefan diffusion coefficients used to quantify mass transfer. Activity coefficient models or equations of state can be fitted to experimental or simulation data, from which thermodynamic factors can be obtained by differentiation. The accuracy of thermodynamic factors determined using indirect routes is dictated by the specific choice of an activity coefficient model or an equation of state. The Permuted Widom’s Test Particle Insertion (PWTPI) method developed by Balaji et al. enables direct determination of thermodynamic factors in binary and multicomponent systems. For highly dense systems, for example, typical liquids, it is well known that molecular test insertion methods fail. In this article, we use the Continuous Fractional Component Monte Carlo (CFCMC) method to directly calculate thermodynamic factors by adopting the PWTPI method. The CFCMC method uses fractional molecules whose interactions with their surrounding molecules are modulated by a coupling parameter. Even in highly dense systems, the CFCMC method efficiently handles molecule insertions and removals, overcoming the limitations of the PWTPI method. We show excellent agreement between the results of the PWTPI and CFCMC methods for the calculation of thermodynamic factors in binary systems of Lennard-Jones molecules and ternary systems of Weeks-Chandler-Andersen molecules. The CFCMC method applied to calculate the thermodynamic factors of realistic molecular systems consisting of binary mixtures of carbon dioxide and hydrogen agrees well with the NIST REFPROP database. Our study highlights the effectiveness of the CFCMC method in determining thermodynamic factors for diffusion, even in densely packed systems, using relatively small numbers of molecules.

Pore pressure fluctuation in subsurface reservoirs and its resulting mechanical response can cause fault reactivation. Numerical simulation of such induced seismicity is important to develop reliable seismic hazard and risk assessments. However, modeling of fault reactivation is quite challenging, especially in the case of displaced faults, i.e., faults with non-zero offset. In this paper, we perform a systematic benchmarking study to validate two recently developed numerical methods for fault slip simulation. Reference solutions are based on a semi-analytical approach that makes use of inclusion theory and Cauchy-type singular integral equations. The two numerical methods both use finite volume discretizations, but they employ different approaches to represent faults. One of them uses a conformal discrete fault model (DFM) while the other employs an embedded (non-conformal) fault model. The semi-analytical test cases cover a vertical frictionless fault, and inclined displaced faults with constant friction and slip-weakening friction. It was found that both numerical methods accurately represent pre-slip stress fields caused by pore pressure changes. Moreover, they also successfully cope with a vertical frictionless fault. However, for the case with an inclined displaced fault with a constant friction coefficient, the embedded method can not converge for the post-slip phase, whereas the DFM successfully coped with both constant and slip-weakening friction coefficients. In its current implementation, the DFM is therefore the model of choice when accurate simulation of local faulted systems is required.

Large-scale storage technologies are crucial to balance consumption and intermittent production of renewable energy systems. One of these technologies can be developed by converting the excess energy into compressed air or hydrogen, i.e., compressed gas, and storing it in underground solution-mined salt caverns. Salt caverns are proven seals towards compressed air and hydrogen. However, several challenges, including fast injection/production cycles and operation of systems of caverns, are yet to be resolved to allow for a safe scale-up of energy storage in salt caverns. To address these challenges, it is important to identify key parameters that impact both the safety and efficiency of the operations. For this purpose, the present study conducts sensitivity analyses to show the importance of different parameters on the time-dependent mechanical behavior of salt caverns, individually and in a multi-cavern system. The impact of different deformation mechanisms (e.g. transient and reverse creep), model calibration, cavern shape, presence of interlayers and multi-cavern interactions are investigated in this study. The constitutive model adopted in this work and the mathematical formulation are presented in detail. Additionally, an open-source three-dimensional simulator, named “SafeInCave”, is developed for the numerical solution of the non-linear governing equations. The findings provide insights into improving the reliability of numerical simulations for the safe and efficient operation of salt caverns in energy storage applications.

The storage of renewable hydrogen in salt caverns requires fast injection and production rates to cope with the imbalance between energy production and consumption. This raises concerns about the mechanical stability of salt caverns under such operational conditions. The use of appropriate constitutive models for salt mechanics is an important step in investigating this issue, therefore many constitutive models with several parameters have been presented in the literature. However, a robust calibration strategy to reliably determine which model and parameter set represents the given rock, based on stress–strain data sets, remains an unsolved challenge. In this paper, for the first time in the community, we present a multi-step strategy to determine a single parameter set based on many deformation data sets for salt rocks. Towards this end, we first develop a comprehensive constitutive model able to capture all relevant nonlinear deformation physics of transient, reverse, and steady-state creep. The determination of the single set of representative material parameters is then achieved by framing the calibration process as an optimization problem, for which the global Particle Swarm Optimization algorithm is employed. To allow for dynamic data integration, a multi-step calibration strategy is developed for a situation where experiments are included one at a time, as they become available. Additionally, due to the existing mild heterogeneity in the experimental rock samples, our optimization strategy is made flexible to allow for this slight heterogeneity. The devised optimization strategy, based on the multi-physics comprehensive constitutive modeling framework, results in a single set of representative material properties of all the deformation data sets. As a rigorous mathematical analysis for the presented method and the lack of relevant experimental data sets, we consider a wide range of synthetic experimental data sets, inspired by the existing sparse relevant data in the literature. The results of our performance analyses show that the proposed calibration strategy is robust. Moreover, the results become increasingly more accurate as more data sets become available.

Underground hydrogen (H₂) storage in saline aquifers is a viable solution for large-scale H₂ storage. Due to its remarkably low viscosity and density, the flow of H₂ within saline aquifers exhibits strong instability, which needs to be thoroughly investigated to ensure safe operations at the storage site. For the first time, we develop a lattice Boltzmann model tailored for pore-scale simulations of the H₂-brine system under typical subsurface storage conditions. The model captures the significant contrast of fluid properties between H₂ and brine, and it offers the flexibility to adjust the contact angle to suit varying wetting conditions. We show that the snap-off is enhanced in a system with a high capillary number and a small contact angle. These conditions lead to a low recovery factor, which is unfavorable for H₂ production from the aquifer. Moreover, the relative permeability curves, computed from the simulation results, exhibit distinct behaviors for H₂ and brine. In the case of the wetting phase, the relative permeability can be quantified using the quadratic expression, whereas for the non-wetting phase, the relative permeability exhibits a nearly linear behavior, and saturation alone appears insufficient to characterize the relative permeability at large saturations of non-wetting phase. This implies that different formula for liquid and gas phases may be employed for continuum-scale simulations.

Quantification of the poromechanical response of subsurface formations due to human-induced pore pressure fluctuations is critical for the performance and stability assessment of many geo-energy systems. In particular, natural faults in the subsurface introduce the hazard of induced seismicity. Numerical modeling of fault reactivation is challenging, while the specific details of induced stresses and fault slip in reservoirs with displaced (i.e. non-zero offset) faults may cause additional challenges depending on the type of numerical formulation employed. To facilitate the systematic development and testing of numerical tools for the simulation of induced seismicity in faulted reservoirs we developed a set of semi-analytical test problems of increasing complexity, based on inclusion theory and Cauchy singular integral equations. With these we investigate the accuracy of two recently developed Finite Volume (FV) schemes with collocated and staggered arrangements of unknowns. One of them employs a conformal discrete fault model (DFM) which can guarantee sufficient accuracy at the cost of adaptive mesh refinement but may suffer from modelling and computational challenges when addressing large-scale realistic geological configurations. The second one employs an embedded (or non-conformal) discrete fault model (EDFM) which avoids the need for excessive mesh refinement, but of which the accuracy and the range of applicability are still to be investigated. We found that both numerical schemes accurately represent the pre-slip Coulomb stresses, but show different degrees of accuracy in representing the resulting depletion-induced fault slip. The semi-analytical benchmark data are available via DOI 10.4121/22240309.

In fractured geological formations, as a result of the in-situ stress changes, fractures can propagate or slide. This phenomenon can be efficiently modeled by the extended finite element method (XFEM) when there are only a few fractures present. However, geological reservoirs contain many fractures which can also cross and are densely populated. Therefore, the classical XFEM is too expensive to be applied for the simulation of propagating fractures in geological formations. To reduce the costs, typically, homogenization or upscaling is used. However, they result in inaccurate solutions, since no separation of scales exists in this process. To resolve this challenge, in this work, a multiscale XFEM (MS-XFEM) for propagating fractures is developed and presented. In each time step, given the current geometries of the fractures, local XFEM-based basis functions are constructed or adaptively updated. The adaptive update takes place in certain regions where fracture geometries are changed due to propagation. Using these basis functions, a very efficient FEM-based coarse-scale system is developed since it has no extra degrees of freedom (DOFs). Once the coarse-scale system is solved, its solution is prolonged to the fine-scale original resolution using the basis functions. This approximate fine-scale solution is then used to estimate the group of growing fractures tips and their growing angles. This allows for exploiting the locality of the propagation process fully while solving a global system. To control the error, an iterative procedure is also developed. Proof-of-concept test cases are presented to study the developed MS-XFEM algorithm. It is shown that MS-XFEM results are capable of predicting the propagating paths for complicated fracture patterns. As such, MS-XFEM casts a promising method for field-scale applications.

We present an efficient compositional framework for simulation of CO2 storage in saline aquifers with complex geological geometries during a lifelong injection and migration process. To improve the computation efficiency, the general framework considers the essential hydrodynamic physics, including hysteresis, dissolution and capillarity, by means of parameterized space. The parameterization method translates physical models into parameterized spaces during an offline stage before simulation starts. Among them, the hysteresis behavior of constitutive relations is captured by the surfaces created from bounding and scanning curves, on which relative permeability and capillarity pressure are determined directly with a pair of saturation and turning point values. The new development also allows for simulation of realistic reservoir models with complex geological features. The numerical framework is validated by comparing simulation results obtained from the Cartesian-box and the converted corner-point grids of the same geometry, and it is applied to a field-scale reservoir eventually. For the benchmark problem, the CO2 is injected into a layered formation. Key processes such as accumulation of CO2 under capillarity barriers, gas breakthrough and dissolution, are well captured and agree with the results reported in literature. The roles of various physical effects and their interactions in CO2 trapping are investigated in a realistic reservoir model using the corner-point grid. It is found that dissolution of CO2 in brine occurs when CO2 and brine are in contact. The effect of residual saturation and hysteresis behavior can be captured by the proposed scanning curve surface in a robust way. The existence of capillarity causes less sharp CO2-brine interfaces by enhancing the imbibition of the brine behind the CO2 plume, which also increases the residual trapping. Moreover, the time-dependent characteristics of the trapping amount reveals the different time scales on which various trapping mechanisms (dissolution and residual) operate and the interplay. The novelty of the development is that essential physics for CO2 trapping are considered by the means of parameterized space. As it is implemented on corner-point grid geometries, it casts a promising approach to predict the migration of CO2 plume, and to assess the amount of CO2 trapped by different trapping mechanisms in realistic field-scale reservoirs.

Hydrogen is a promising energy carrier for a low-carbon future energy system, as it can be stored on a megaton scale (equivalent to TWh of energy) in subsurface reservoirs. However, safe and efficient underground hydrogen storage requires a thorough understanding of the geomechanics of the host rock under fluid pressure fluctuations. In this context, we summarize the current state of knowledge regarding geomechanics relevant to carbon dioxide and natural gas storage in salt caverns and depleted reservoirs. We further elaborate on how this knowledge can be applied to underground hydrogen storage. The primary focus lies on the mechanical response of rocks under cyclic hydrogen injection and production, fault reactivation, the impact of hydrogen on rock properties, and other associated risks and challenges. In addition, we discuss wellbore integrity from the perspective of underground hydrogen storage. The paper provides insights into the history of energy storage, laboratory scale experiments, and analytical and simulation studies at the field scale. We also emphasize the current knowledge gaps and the necessity to enhance our understanding of the geomechanical aspects of hydrogen storage. This involves developing predictive models coupled with laboratory scale and field-scale testing, along with benchmarking methodologies.

Underground Hydrogen Storage (UHS) is an emerging large-scale energy storage technology. Researchers are investigating its feasibility and performance, including its injectivity, productivity, and storage capacity through numerical simulations. However, several ad-hoc relative permeability and capillary pressure functions have been used in the literature, with no direct link to the underlying physics of the hydrogen storage and production process. Recent relative permeability measurements for the hydrogen-brine system show very low hydrogen relative permeability and strong liquid phase hysteresis, very different to what has been observed for other fluid systems for the same rock type. This raises the concern as to what extend the existing studies in the literature are able to reliably quantify the feasibility of the potential storage projects. In this study, we investigate how experimentally measured hydrogen-brine relative permeability hysteresis affects the performance of UHS projects through numerical reservoir simulations. Relative permeability data measured during a hydrogen-water core-flooding experiment within ADMIRE project is used to design a relative permeability hysteresis model. Next, numerical simulation for a UHS project in a generic braided-fluvial water-gas reservoir is performed using this hysteresis model. A performance assessment is carried out for several UHS scenarios with different drainage relative permeability curves, hysteresis model coefficients, and injection/production rates. Our results show that both gas and liquid relative permeability hysteresis play an important role in UHS irrespective of injection/production rate. Ignoring gas hysteresis may cause up to 338% of uncertainty on cumulative hydrogen production, as it has negative effects on injectivity and productivity due to the resulting limited variation range of gas saturation and pressure during cyclic operations. In contrast, hysteresis in the liquid phase relative permeability resolves this issue to some extent by improving the displacement of the liquid phase. Finally, implementing relative permeability curves from other fluid systems during UHS performance assessment will cause uncertainty in terms of gas saturation and up to 141% underestimation on cumulative hydrogen production. These observations illustrate the importance of using relative permeability curves characteristic of hydrogen-brine system for assessing the UHS performances.