Underground hydrogen storage (UHS) is a potential technology that can resolve renewable energy supply-demand challenge at seasonal (terawatt-hours) scales. Enabling this technology and optimizing its performance require a wide range of analyses from hydrodynamics to geomechanics and biogeochemistry, among which understanding the transport (and trapping) of hydrogen in porous rocks stands out. A key parameter in quantification of hydrogen transport in partially brine-saturated geological formations is its relative permeability (k_rh). In this study, we develop a theoretical k_rh model using upscaling concepts from effective medium approximation and percolation theory. Our theoretical model, developed based on pore-scale characteristics, estimates k_rh from pore size distribution, capillary pressure curve or mercury intrusion capillary pressure curve, and critical hydrogen saturation, S_hc, at which k_rh approaches zero. We evaluate the proposed model using eight experimental datasets and eleven pore-network simulations. Discrepancies are observed for some of the carbonate samples, likely due to secondary porosity effects (e.g., presence of vugs and/or fractures), and in some of the sandstone rocks, possibly due to imprecise S_hc estimation. These observations highlight the importance of improving pore structure characterization to better account for such heterogeneities and enhance model accuracy for reliable quantification of the k_rh relevant to UHS applications. These findings also highlight the critical role of accurate parameter estimation, such as determining the S_hc in estimating k_rh. Overall, the study demonstrates that the proposed approach provides a cost-effective and practical alternative to extensive experiments and simulations, offering a promising tool for quantifying k_rh relevant to UHS applications.

In the past three decades, a wide array of computational methodologies and simulation frameworks have emerged to address the complexities of modeling flow and transport processes in fractured porous media. The conformal mesh approaches which explicitly align the computational grid with fracture surfaces are considered by many to be the most accurate. However, such methods require excessive fine-scale meshing, rendering them impractical for large or complex fracture networks. The Embedded Discrete Fracture Model (EDFM) offers a good balance between accuracy and efficiency and has gained a lot of traction in recent years. Nonetheless, it is not free of drawbacks; EDFM can, and often will, generate fracture cells that have orders of magnitudes smaller volumes than the matrix cells, which significantly impacts the convergence performance of nonlinear solvers. In this work, we propose to learn the complex flow and transport dynamics in fractured porous media with graph neural networks (GNN). GNNs are well suited for this task due to the unstructured topology of the computation grid resulting from EDFM discretization. We propose two deep learning architectures, a GNN and a recurrent GNN. Both networks follow a two-stage training strategy: an autoregressive one step roll-out, followed by a fine-tuning step where the model is supervised using the whole ground-truth sequence. We demonstrate that the two-stage training approach is effective in mitigating error accumulation during autoregressive model rollouts in the testing phase. Our findings indicate that both GNNs generalize well to unseen fracture realizations. While the second stage of training proved to be beneficial for the GNN model, its impact on the recurrent GNN model was less pronounced. Finally, the performance of both GNNs for temporal extrapolation is tested. The recurrent GNN significantly outperformed the GNN in terms of accuracy, thereby underscoring its superior capability in predicting long sequences.

This study extends the Algebraic Dynamic Multilevel (ADM) method for simulating contaminant transport and retention in vadose zones. Building upon a fully implicit scheme that couples variably saturated flow and contaminant transport, the developed ADM framework effectively predicts contaminant plume migration across both unsaturated and saturated media under heterogeneous conditions. During the simulation, ADM dynamically adjusts grid resolution based on the spatial gradients of primary variables, applying fine-scale grids in regions with steep gradients and coarsening the mesh where fields remain smooth. These dynamic adjustments are achieved through prolongation and restriction operators that transfer solutions across multilevel grid systems. As both water content and contaminant concentration evolve spatiotemporally, dual coarsening criteria are introduced to simultaneously capture flow and transport dynamics. Results show that the developed model reproduces the contaminant migration obtained from the fully resolved solution using substantially fewer grids. Moreover, it offers the flexibility to trade off numerical accuracy against computational cost by selecting an appropriate coarsening criterion.

This study presents a numerical investigation of pressure solution creep and its influence on the mechanical behavior of salt caverns for underground hydrogen storage. A 3D modeling framework, implemented in the open-source simulator SafeInCave, incorporates both dislocation and pressure solution creep mechanisms and is applied to caverns with varying geometries, depths, temperatures, and interlayer positions under realistic conditions. The creep models are appropriately calibrated against experimental results from the literature to account for both stress and temperature effects. Results show that pressure solution creep becomes increasingly significant over time, particularly in shallow and cold formations, where it dominates deformation. It is more active away from cavern walls, where stresses and temperatures are low, while dislocation creep concentrates near the cavern walls and governs behavior at greater depths and higher temperatures. Overall, the study demonstrates that accurately capturing the effect of pressure solution creep is essential for reliable prediction of deformation and structural integrity in underground hydrogen storage caverns.

This study introduces a multiscale simulation framework, termed Projection-based Embedded Discrete Fracture Modeling with Algebraic Dynamic Multilevel method (pEDFM-ADM), which integrates an embedded discrete fracture network representation with a fully algebraic, front-tracking-based mesh adaptation strategy. Incorporating a fully implicit scheme, compositional thermodynamics, and algebraic multilevel operators, the framework captures essential subsurface processes such as buoyancy-driven migration, convective dissolution, phase partitioning, and fracture-matrix interactions under geologically realistic conditions. The method constructs a hierarchy of multilevel grids and localized multiscale basis functions that introduce fine-scale heterogeneities at each coarse level. Adaptive mesh refinement and coarsening are driven by local variations in CO₂ mass fraction and executed through algebraic prolongation and restriction operators, enabling efficient projection between grid levels. The framework is systematically evaluated across a sequence of test cases with increasing complexity, including systems with low-permeability flow barriers, highly conductive fractures, striking a trade-off between computational resource and detailed simulation accuracy. Overall, the pEDFM-ADM framework provides a scalable, fully algebraic, and physically adaptive modeling tool for large-scale CO₂ storage simulations in fractured porous media, supporting predictive simulation and risk assessment for long-term carbon sequestration.

Multiscale simulation frameworks are essential to quantify the CO2 trapping and migration in large-scale saline aquifers, which entail highly-resolved fine-scale heterogeneous properties. However, classical upscaling approaches which aim to define effective properties on larger grid sizes can lead to significant and systematic overestimation of the solubility and residual trapping mechanisms. Reliable assessment of these two trapping mechanisms is crucial to ensure the integrity of the storage process and properly mitigate the leakage risks. Therefore, it is essential to develop advanced simulation technologies that are both accurate and efficient (i.e., scalable) for simulation of complex CO2 plume dynamics within large-scale heterogeneous reservoir models. To overcome this challenge, in this work three advanced strategies are developed and investigated: Effective Values (EV) for parameters, Local Grid Refinement (LGR) and Algebraic Dynamic Multilevel (ADM). The numerical investigations specially include a set of consistent models in the Ponta Aguda saline aquifer, with a total area of 40,000 km2[jls-end-space/], located offshore the Brazilian coast. The results indicate that the ADM is a promising method, delivering stable and robust results in a representative section of the field. This encourages further extensions of this method for real-field deployment. Specially, LGR and EV are found to be limited in their scopes for field simulations, since they depend on a matching pre-procedure (against a reference solution) for their upscaled parameters before any new simulations can be run. In addition, their tuned parameters cannot be transferred from one model to another. ADM, on the other hand, does not require any upscaling procedure, as the multiscale basis functions allow for consistent mapping across resolutions.

Underground hydrogen storage (UHS) in underground geological reservoirs is a promising solution for large-scale energy storage. However, several challenges, particularly geomechanical ones, must be resolved before UHS can be widely and safely deployed. The interactions between hydrogen, brine, and reservoir rock, combined with the cyclic stresses resulting from hydrogen injection and withdrawal may affect the mechanical integrity of the reservoir, the caprock, as well as its surrounding formations. This is an experimental investigation into the geomechanical impact of a 6 month exposure of clay-rich sandstone (Yellow Felser) rocks to hydrogen and/or brine. Cm-scale samples were exposed to hydrogen-saturated brine at 150 bar and in an autoclave for the period of six months. Afterwards, triaxial cyclic loading experiments were conducted on the samples under confining pressures of 10, 20, and 30 MPa. The results are compared with those from the reference samples, which have been exposed to brine only, for the same time period. Each mechanical test included eight stress cycles in the linear stress regime (below the brittle yield point), followed by loading to failure. The frequency, amplitude, and stress conditions were tailored to each confining pressure. The results showed that six months of hydrogen-saturated brine exposure had no noticeable effect on the failure envelope, elastic properties, inelastic strain, and acoustic properties of the Yellow Felser sandstone compared to exposure to brine alone. Internal friction, P-wave velocity, and Young’s modulus each showed a change of around 3%, which is on the same order as the repeatability and therefore indicating minimal geomechanical alteration. Complementary qualitative and quantitative scanning electron microscopy (SEM) analyses revealed negligible microstructural changes. When eight stress cycles were applied within the linear stress regime, the majority of inelastic strain occurred during the first cycle, with no progressive accumulation thereafter. A comparison with samples tested under monotonic loading to failure confirmed that cyclic loading under these conditions does not affect the rock strength of Yellow Felser sandstone. These findings provide new insights into the combined effects of cyclic stress and hydrogen/brine/rock interactions on the geomechanical behavior of clay-rich sandstones under reservoir-relevant pressure and temperature conditions.

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.

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.

This work addresses numerical instabilities that can appear when computing the mean stress in linear elasticity and coupled poroelasticity problems discretized with low-order finite elements. The linear elasticity and coupled poroelasticity models are solved using both primal and mixed finite element formulations. Stabilization is obtained by enriching the finite element approximation with an approximation of the Laplacian of displacements. This Laplacian is then evaluated with the Physical Influence Scheme (PIS) by leveraging the underlying governing equation. A key step in the proposed stabilization is the calculation of a parameter h, often computed in the literature as a characteristic length of the element. In this work, we calculate h by solving an optimization problem at the element level. To avoid the high computational cost associated with this procedure, a machine learning model is proposed to predict the optimal h. The benefit of combining PIS with an appropriate computation of h is that the resulting stabilization scheme does not rely on any type of heuristic or user-specified tuning parameter, as often required in other stabilization methods. The results show that the proposed stabilization strategy can effectively remove both saddle-point and Gibbs mean stress oscillations in linear elasticity. We also report, for the first time, that mean stress oscillations can also appear when solving coupled poroelasticity problems, and, differently from pore pressure oscillations (which naturally vanish with time), mean stress instabilities are persistent throughout the whole simulation time, unless deliberately removed. The proposed stabilized mixed formulation is able to remove both pore pressure and mean stress oscillations in coupled poroelasticity problems. Finally, the calculation of h is shown to be critical for the quality of the stabilization, with the machine learning-based approach providing the best compromise between numerical diffusion and accuracy.

This extensive review highlights the central role of classical molecular simulation in advancing hydrogen (H2) technologies. As the transition to a sustainable energy landscape is urgently needed, the optimization of H2 processes, spanning production, purification, transportation, storage, safety, and utilization is essential. To this end, accurate prediction of thermodynamic, transport, structural, and interfacial properties is important for overcoming engineering challenges across the entire H2 value chain. Experimental measurements, despite being the traditional way of obtaining these properties, can be limited by the distinctive nature of H2, harsh operating conditions, safety constraints, and extensive parameter spaces. Free from such limitations, classical molecular simulations, in the general frameworks of Monte Carlo and Molecular Dynamics, provide an optimal balance between computational efficiency and accuracy, bridging the gap between quantum mechanical calculations and macro-scale modeling. This review also systematically covers molecular simulation methods and force fields for computing key properties of H2 systems, such as phase and adsorption equilibria and transport coefficients. Beyond property prediction, we explore how molecular simulation reveals fundamental mechanisms governing hydrate formation and dissociation, membrane permeations, and H2 embrittlement. When possible, data from multiple sources are compared and critically assessed, while effort is put on evaluating the force fields used and methodological approaches followed in the literature. Finally, this review aims at identifying research gaps and future opportunities, emphasizing emerging approaches, such as molecular simulation in the era of artificial intelligence.

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 storage in porous media is promising for large-scale energy storage. However, its technical and financial effectiveness is heavily dependent on a reliable site selection strategy. In this review, we critically assess the available literature across disciplines to identify the most influential criteria for reliable site selection. Drawing from this evaluation, we propose a systematic, multidisciplinary framework for early stage reservoir screening, integrating key criteria from reservoir performance, geomechanics and containment, location and techno-economics, and biogeochemistry. Our framework allows for rapid identification and ranking of the most suitable reservoirs by proposing 11 elimination criteria and 15 screening criteria. The presented framework consists of practically applicable and scientifically grounded criteria to support consistent, early stage decisions based on readily available data while allowing for detailed site-specific analysis in later project development phases. By unifying diverse disciplinary insights into a structured methodology, this study contributes to more informed, inclusive, and effective site selection.

As underground hydrogen storage (UHS) is expected to play a key role in future renewable energy systems, understanding the potential geomechanical risks, such as induced seismicity, is essential. Therefore, this study aims to assess the probability of induced seismicity associated with the prospect of large-scale UHS plans. We commence by developing simulation models with increasing complexity, starting from the basic characteristics of the salt formation, salt cavern, and operational conditions, and progressing to the inclusion of structural features within the salt formation as well as in the overburden and sideburden. A 2D finite element simulator is used to incorporate deformation and simulate creep behaviour, which is subsequently coupled with a rate-and-state Coulomb threshold model to compute the seismicity rate from stress changes. The developed framework accounts for the geological and mechanical characteristics of the heterogeneities that influence local stress fields, allowing us to identify conditions that may increase seismic risk or enhance stability.

Regional-scale saline aquifers are promising candidates for geological CO2 storage but present significant modeling challenges due to their vast extent, heterogeneity, and limited subsurface data. This study introduces a multiscale modeling framework that was applied to assess CO2 storage in the Ponta Aguda saline aquifer (Santos Basin, Brazil, 40,000 km2 area). Consistency of the multiscale models is checked by combining boundaries conditions for pressure match and verification of trapping mechanisms representativity. Four different methods were evaluated regarding the trapping mechanisms accuracy in coarse models: Local Grid Refinement, Effective Values, and Algebraic Dynamic Multilevel. Compositional simulations conducted with CMG-GEM and DARSim2 demonstrate that coarse-scale models systematically overestimate CO2 trapping due to numerical artifacts, particularly in solubility and hysteresis behavior. These artifacts arise from mismatched CO2/brine volumes in large cells, leading to artificially enhanced trapping efficiency. Among the evaluated methods, Algebraic Dynamic Multilevel delivers the most reliable predictions, providing a general solution that aligns closely with fine-scale reference simulations while remaining computationally feasible. The results show the importance of scale-consistent modeling approaches for accurate CO2 storage assessment and highlight the risks of relying on overly simplified coarse models in the design and optimization of carbon storage projects in large saline aquifers.

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.

The North Sea’s potential as a Green Energy Hub depends on large-scale CO2 storage in shallow-marine sandstones, but the effects of geologic heterogeneity, such as permeability barriers and capillary entry pressure contrasts, remain underexplored. This study uses multiphase flow simulations on geologically realistic, surface-based reservoir models informed by outcrop analogue data from wave-dominated shoreface sandstones. We investigate how sedimentological heterogeneity influences CO2 plume migration, pressure evolution, and storage capacity.

Preliminary results show that capillary barriers tied to facies architecture and early cementation, conditioned to clinoform geometries, significantly control plume movement. These barriers promote lateral spreading and residual trapping, representing a potential upper limit on long-term CO2 storage when stable. Clinoform-related heterogeneity also induces flow compartmentalization, limiting pressure dissipation and enhancing anisotropy, which may reduce injectivity and cause spatially variable pressure buildup.

Comparisons with waterflood simulations reveal contrasting dynamics: water advances more uniformly, while CO2 migration is more sensitive to fine-scale architecture due to its lower interfacial tension and capillary entry pressures. These findings underscore the need to incorporate realistic sedimentological heterogeneity in dynamic models to avoid misestimating injectivity, pressure behavior, and storage security. This approach offers a robust framework for early-stage screening and risk assessment in complex storage settings.

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.

Finite-size effects of transport properties computed from molecular dynamics simulations are investigated for Weeks-Chandler-Andersen systems at reduced densities of 0.05 (dilute gas), 0.45 (dense gas), and 0.85 (fluid close to the solid-liquid transition). Viscosities, self-diffusivities, Onsager coefficients, and electrical conductivities are computed for various system sizes ranging from 64 to 8192 WCA particles at each density. At dilute and intermediate densities, finite-size corrections to the transport properties significantly deviate from the widely used Yeh–Hummer correction, which was originally developed for the liquid phase.