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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.