The idea of a vast, untapped reservoir of natural hydrogen, (1) one that could transform global energy systems, is as enticing as it is elusive. If a substantial, economically viable hydrogen field exists and can be exploited with minimal leakage to the atmosphere, it would mark a paradigm shift, offering a relatively clean energy source with low operational emissions. Unlike fossil fuels, which require carbon-intensive extraction and combustion, or green hydrogen, (2) which depends on large-scale renewable energy input, natural hydrogen could provide an energy supply that does not require conversion from electricity or fossil fuels, making it a potentially simpler and more continuous resource. (3) However, because hydrogen is an indirect greenhouse gas, any future exploitation must include robust monitoring strategies to quantify and minimize atmospheric losses in order to preserve its environmental advantages.

But there is a fundamental challenge: despite numerous documented hydrogen occurrences worldwide, no truly large-scale, commercially viable reserve has been found. This reality forces us to confront a crucial question, are we on the verge of discovering a new energy frontier, or will natural hydrogen remain a scientific curiosity (4) with limited economic impact?

To answer this, we must address two key uncertainties: Where should we look for a major natural hydrogen deposit? And what steps are needed to confirm its economic viability? The answer requires an exploration strategy that combines geological understanding, advanced detection techniques, and a willingness to take financial risks. If we fail to find a large reserve, alternative strategies, such as engineered subsurface hydrogen production (hydrogen farming (5)/stimulated hydrogen (6)), will become increasingly necessary, albeit at a higher cost.

Residual trapping is a critical mechanism influencing the efficiency of Underground Hydrogen Storage (UHS). This study investigates the underlying processes of residual trapping by bypassing, through bifurcating geometries, focusing on how geometrical parameters and flow characteristics affect the trapping process. We develop a dynamic simulation framework based on the lattice Boltzmann method (LBM) to simulate full drainage/imbibition cycles. Various geometries, based on the pore doublet model, were investigated and supported by theoretical analysis. In addition, trapping behavior of hydrogen was compared to that of CO₂ and CH₄. It is found that the channel width ratio, specially across the local bifurcating geometries, and the roundness of the grains, are among the key factors which control hydrogen trapping. Results indicate that the suited reservoirs for underground hydrogen storage have narrower channel-size ratios and smoother edges at micro-scale. Operational conditions also play a significant role. Lower flow rates enhance bypassing, which increases trapping.

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.

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.

This work introduces a novel application of the Algebraic Dynamic Multilevel (ADM) method for simulating CO₂ storage in deep saline aquifers. By integrating a fully implicit coupling strategy, fully compositional thermodynamics, and adaptive mesh refinement, the ADM framework effectively models phenomena such as buoyancy-driven migration, convective dissolution, and phase partitioning under various subsurface conditions. The method starts with the construction of a hierarchy of multilevel grids and the generation of localized multiscale basis functions, which account for heterogeneities at each coarse level. During the simulation, the ADM method dynamically refines areas with significant overall CO₂ mass fraction gradients while coarsening smooth regions, thus optimizing computational resources without compromising the accuracy required to capture essential flow and transport characteristics. This dynamic grid adjustment is facilitated by algebraic prolongation and restriction operators, which map the fine-scale system onto a coarser grid suited to the evolving distribution of the CO₂ plume. This feature allows the ADM to navigate various coarsening thresholds efficiently, striking a trade-off between computational economy and detailed simulation accuracy. Even at relatively higher thresholds, key trapping mechanisms are captured with sufficient detail for quantification. These capabilities make the ADM framework well suited for long-term CO₂ sequestration in highly heterogeneous reservoirs, where large-scale models may otherwise become impractically expensive, offering a practical balance between the need for detailed simulations and manageable computational requirements. Overall, the ADM framework proves to be a robust tool for designing, monitoring, and analyzing large-scale CO₂ storage operations, supporting reliable and cost-effective solutions in carbon management.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

Considering the storage capacity and already existing infrastructures, underground porous reservoirs are highly suitable to store green energy, for example, in the form of green gases such as hydrogen and compressed air. Depending on the energy demand and supply, the energy-rich fluids are injected and produced, which induces cyclic change of state-of-the-stress in the reservoir and its surrounding. Detailed analyses of the geo-mechanical deformations under variable storage conditions i.e., storage frequency and fluid fluctuating pressures, are crucially important for safe and efficient operations. The present work presents an integrated analysis, based on experimental and constitutive modeling aspects, to investigate sandstones’ geomechanical response to cyclic loading relevant to underground energy storage (UES). To this end, sandstone rock samples were subjected to cyclic loading above and below the onset of dilatant cracking under different frequencies and loading amplitudes. Axial strains and Acoustic Emissions (AE) were measured in both regimes to quantify the total deformation (strain) of the rock and its AE characteristics. It is found that the inelastic strain and number of AE events is the highest in the first cycle and reduce subsequently cycle after cycle. Moreover, cyclic inelastic deformations are affected by the mean stress, amplitude, and frequency of the stress waveform. On the one hand, the higher the mean stress and the amplitude, the higher the total inelastic strains. On the other hand, the lower the frequency, the higher the total inelastic strain. From the modeling perspectives, five types of deformation mechanisms were identified based on the governing physics: elastic, viscoelastic, compaction-based cyclic inelastic, inelastic brittle creep, and dilatation-based inelastic deformation. To model elastic, viscoelastic, and brittle creep, the Nishihara model was used. A cyclic modified cam clay model (MCC) and hardening–softening model were applied to capture plastic deformation. The results show a very good fit of the constitutive model with the experimental results, which could help in studying the response of reservoirs to injection and production.

Reservoir simulations for subsurface processes play an important role in successful deployment of geoscience applications such as geothermal energy extraction and geo-storage of fluids. These simulations provide time-lapse dynamics of the coupled poromechanical processes within the reservoir and its over-, under-, and side-burden environments. For more reliable operations, it is crucial to connect these reservoir simulation results with the seismic surveys (i.e., observation data). However, despite being crucial, such integration is challenging due to the fact that the reservoir dynamics alters the seismic parameters. In this work, a coupled reservoir simulation and time-lapse seismic methodology is developed for multiphase flow operations in subsurface reservoirs. To this end, a poromechanical simulator is designed for multiphase flow and connected to a forward seismic modeller. This simulator is then used to assess a novel methodology of seismic monitoring by isolating the reservoir signal from the entire reflection response. This methodology is shown to be able to track the development of the fluid front over time, even in the presence of a highly reflective overburden with strong time-lapse variations. These results suggest that the proposed methodology can contribute to a better understanding of fluid flow in the subsurface. Ultimately, this will lead to improved monitoring of reservoirs for underground energy storage or production.