Phase equilibrium calculations play an important role in a wide variety of applications in chemical and petroleum engineering. In this work, we focus on CO₂-hydrocarbon mixtures, with applications ranging from enhanced oil recovery processes to CO₂ storage. In compositional reservoir simulation, both robustness and efficiency are of utmost importance. The conventional approach for multiphase equilibrium consists of a sequence of phase stability and flash calculations. At each level of the stepwise process, stability testing is performed starting from several initial guesses; therefore, reducing the number of stability calls and using judiciously the information from stability to initialize a phase split are key points in developing an efficient stability-flash algorithm. Two new initialization strategies for multiphase flash calculations are proposed. The first one (improved stepwise initialization) follows the conventional procedure, but uses additional initial guesses. In the second one (improved multiple initialization), a three-phase split is initiated if at least three minima of the tangent plane distance (TPD) function are detected by stability analysis of feed composition. Both proposed methods are using all information from phase stability testing at each stage. Unlike in previous formulations, compositions at all minima of the TPD function, including trivial and positive TPDs are used to generate initial equilibrium constants. Highly robust routines are used, based on successive substitution iterations (SSI) in early iteration stages, followed by Newton iterations with modified Cholesky factorization and line search, in both stability and flash calculations. The proposed methods are tested and compared with the conventional procedure for several benchmark mixtures from the literature, containing hydrocarbon components and CO₂. Phase diagrams are constructed in the P-Z plane, focusing on the number of stationary points of the TPD functions found in each step of the multiphase stability-flash algorithm and on how they must be efficiently used in initialization. For all the test mixtures, in the proposed stability-flash strategy, the number of calls of the stability and flash routines and the number of iterations in flash calculations are significantly reduced as compared to previous approaches, recommending the new approach as a useful tool in compositional simulation.

In this work, we present a kinetic simulation model for gas hydrates in porous media using the Operator-Based Linearization (OBL) technique. The OBL approach introduces algebraic operators that represent the physical terms in the mass and energy balance equations. Operators are calculated only in supporting points comprising the discretized parameter space, and operator values and partial derivatives for linear system assembly are readily obtained through (multi-)linear interpolation. Taking advantage of this setup, the implementation of advanced thermodynamic models for hydrate formation and dissociation under kinetic assumptions is simplified. We test the assumptions for thermodynamic modelling by analysing the Gibbs energy surfaces of the fluid and hydrate phases and demonstrate that, in the limit, the thermodynamic equilibrium for both kinetic and equilibrium reaction models is equivalent. We compare the simulation results with the published experimental results for CH₄-hydrates and extend the assessment to a CO₂-hydrate formation experiment in a semi-batch, constant-pressure configuration. The model reproduces the main pressure–temperature transients and hydrate evolution for both CH₄- and CO₂-systems. We demonstrate applicability at core scale for hydrate formation and, at field scale, for gas production from CH₄-hydrates by thermal stimulation and depressurization. The interaction of thermal-compositional phenomena (phase changes, adiabatic expansion, kinetic rates, and reaction enthalpy) gives rise to highly nonlinear physics that an appropriate OBL discretization resolves. Overall, the patterns of hydrate formation and dissociation are highly sensitive to the kinetic-rate inputs; hence, the appropriate choice of the reaction model remains a key consideration from both physical and numerical perspectives.

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.

In this work, we present a thermal-compositional simulation framework for modelling of CO₂ sequestration in de- pleted hydrocarbon reservoirs. The parametrization technique utilizes thermodynamic state-dependent operators expressing the governing equations for the thermal-compositional system to solve the nonlinear problem. This approach provides flexibility in the assembly of the Jacobian, which allows straightforward implementation of advanced thermodynamics. Taking advantage of the flexibility of operator-based linearization (OBL), multiphase thermodynamic modelling at arbitrary state specifications is implemented. The use of a hybrid-EoS approach to combine equations of state for aqueous and hydrocarbon phases and advanced initialization schemes for multi- phase equilibrium calculations improves the accuracy and efficiency of the simulation. Careful phase identifica- tion is required for the simulation of multiphase flow, in particular with the potential occurrence of multiple liquid phases in CO₂-hydrocarbon mixtures. We apply the simulation framework to model a set of CO₂ injection cases at conditions typical for depleted hydrocarbon fields. We demonstrate that important thermophysical phenomena resulting from the interaction of CO₂ and impurities with reservoir fluids can be accurately captured using the OBL approach. The consistency of compositional simulation is supported by robust and efficient modelling of multiphase equilibria between brines, hydrocarbons and CO₂. The method is shown to be robust for capturing the thermal effects related to expansion, mixing and phase transitions. This work presents a highly flexible and efficient framework for modelling of multiphase flow and transport in CCUS-related subsurface applications. Ro- bust modelling of thermodynamic equilibria at arbitrary state specification captures the complex thermophysical interactions between CO₂ and reservoir fluids.

This study uses the concept of exergy-return on exergy-investment (ERoEI) to evaluate the life-cycle exergetic efficiency and CO₂ intensity (grams CO₂ per MJ of electricity) of (diabatic and adiabatic) compressed air energy storage (CAES) systems. Several CAES configurations are assessed under defined system boundaries, including diabatic systems powered by methane (CH₄) or hydrogen (H2), and adiabatic system with a thermal energy storage (TES) facility.

The results show that conventional (diabatic) CAES system powered by natural gas has the lower exergetic efficiency and higher CO2 intensity compared to adiabatic CAES due to the heat dissipation during compression stage and additional fuel requirements for reheating the air during expansion. Integrating carbon capture and storage (CCS) plant with conventional diabatic CAES can nearly halve the CO₂ intensity for electricity generation although the additional exergy investment for the CCS process reduces the exergetic efficiency of the system. Transitioning to green H2 (produced from low-carbon electricity) as the primary turbine fuel in the diabatic CAES results in a 65–76 % reduction in CO₂ intensity. However, the average exergetic efficiency of system decreases by around 10 %, mainly due to the substantial exergy investment associated with hydrogen production. It is also found that the adiabatic CAES system integrated with TES demonstrates the highest thermodynamic and environmental performance. When 100 % of compression heat is captured and reused during discharge phase, the system reaches ERoEI values up to 61 % with CO2 intensity of 12–26 g CO₂ per MJe.

Disclaimer: The results and performance metrics presented in this study are based on modelled scenarios and literature-derived parameters under defined system boundaries. Actual performance of CAES systems may vary depending on site-specific conditions, technology maturity, and operational configurations. All efficiency values, CO₂ intensity estimates, and comparative assessments should be interpreted within the context of the assumptions and limitations described herein. This study does not constitute a commercial endorsement or performance guarantee. The authors have made every effort to ensure accuracy but accept no liability for decisions made based on this analysis.

The SPE11 comparative solution project presents a benchmark for geological carbon storage in an aquifer, as the development of sufficiently accurate CO₂ sequestration models is critical for predicting the distribution of CO₂ during and after injection. In this paper we present a convergence analysis of the SPE11 benchmark simulation using the Delft Advanced Research Terra Simulator (open-DARTS). Open-DARTS, an open-source simulation framework designed for forward and inverse modeling, as well as uncertainty quantification, employs a unified thermal-compositional formulation and operator-based linearization. In our convergence analysis the SPE11b (2D - reservoir conditions) starts to converge at a grid resolution of 1340 × 240, after which added resolution provides diminishing returns. In addition the three-dimensional SPE11c benchmark is simulated with 8M grid blocks. However, 2D results from SPE11b suggest that a greater resolution is required for a truly converged solution. Furthermore, we extend the SPE11b benchmark to include H₂S as a trace impurity in the injection stream.

Geothermal energy offers a sustainable source of heat and electricity but alters reservoir pressure and temperature, affecting in-situ stress and potentially triggering fault reactivation and induced seismicity. Deep geothermal reservoirs are valuable for their high temperatures but pose challenges like low permeability and fracture-dominated flow, increasing the risk of fault instability.

This study explores two approaches to assess stress changes: a semi-analytical geomechanical proxy and a fully-coupled Thermo-Hydro-Mechanical (THM) model using open-DARTS. The THM model simulates coupled thermal, hydraulic, and mechanical processes in complex rock formations, while the proxy method approximates displacements and stress changes using reservoir simulation outputs and homogeneous geomechanical rock properties assumptions.

The proxy model has been applied to matrix- and fault-dominated systems, including the Brugge dataset. Results include pressure, temperature, displacements, stress changes predictions over 30 years. Fault stability is evaluated using Mohr-Coulomb criteria with a constant friction coefficient.

In fracture-dominated systems, faults often control flow but. Discrete Fracture Model (DFM) has been used for flow modelling.

Combining proxy and THM models can optimize the balance between accuracy and computational cost. The study emphasizes the differing impacts of pressure and temperature on fault stability during geothermal operations.

Three-phase equilibrium calculations for water-CO2-hydrocarbon mixtures are required in the compositional simulation of various applications in CO2 storage, geothermal systems, and enhanced oil recovery. The very low solubility of hydrocarbon components in water leads to a special mathematical structure of the problem. Several techniques were suggested, such as the free-water flash (FWF) and the augmented free-water flash (AFWF); in the former, the aqueous phase is pure water, while in the latter only certain components, CO2 or methane for example, are dissolved in the aqueous phase. However, only the first-order successive substitution method was used in the previous published approaches, making them unattractive for compositional simulations in which a significant number of phase equilibrium calculations are performed. In this work, a robust and efficient AFWF method is proposed, using combined successive substitutions-modified Newton iterations. The new method is general, allowing partial solubility of any selected component in the water-rich phase, depending on the specific compositions and operating conditions. A detailed description of second-order methods in a Gibbs energy minimization framework for the general AFWF is presented. In the AFWF, the dimension of the problem and the number of function evaluations (thus the computation time) are significantly reduced. Moreover, it is shown that the augmented method always has better convergence properties than its conventional multiphase flash counterpart, in both first- and second-order methods. The new AFWF method is tested for various hydrocarbon-water-CO2 mixtures and proved to be robust and efficient, systematically outperforming the conventional approach. Unlike in previous AFWF formulations, the number of components soluble in water is not limited, leading to a controlled accuracy with respect to a full three-phase equilibrium, even at high pressures and/or large amounts of CO2.

Subsurface CO₂ sequestration is a promising method to advance carbon neutrality and support the shift toward sustainable energy. However, the unique behavior of CO₂ in these operations, particularly for cold CO₂ injection in depleted hydrocarbon reservoirs, poses challenges to wellbore injectivity, reservoir containment, and reservoir capacity. These challenges necessitate the development of a numerical model to better understand and optimize the interplay between wellbore dynamics and reservoir processes. In this work, we present the development of an open-source coupled wellbore-reservoir numerical model, named DARTS-well, which is tailored to CO₂ disposal in subsurface reservoirs. To this end, a multi-segment, multi-phase, non-isothermal wellbore model is first developed using the Drift-Flux Model (DFM), and its results for selected CO2 injection scenarios are validated against the commercial transient wellbore simulator OLGA. The multi-segment wellbore model is then coupled with the Delft Advanced Research Terra Simulator (DARTS) which is used in this study as the reservoir simulator. DARTS is widely used and validated for energy transition applications. The coupled model utilizes the Operator-Based Linearization (OBL) technique, employing state-dependent operators for thermodynamic properties interpolated from predefined tables or generated on the fly. This OBL parametrization approach addresses challenges associated with complex physics and reduces computational time, making it well-suited for modeling subsurface CO₂ sequestration.

Physics-Informed Neural Networks (PINNs) gains attentions as a promising approach for applying deep neural networks to the numerical solution of nonlinear partial differential equations (PDEs). However, due to the challenging regions within the solutions of 'stiff' PDEs, e.g., shock front of CO₂ immiscible flooding, adaptive methods are essential to ensure the neural network accurately addresses these issues. In this work, we introduce a novel method for adaptively training PINNs, named Self-Adaptive PINNs (SA-PINNs). This approach employs fully trainable adaptation weights that are applied individually to each training point. Consequently, the neural network autonomously identifies challenging regions of the solution space and focuses its learning efforts on these areas. This method is hereby used to simulate a two-phase immiscible flooding in a low-permeability oil reservoir, with considering gas dissolution and the threshold pressure gradient of oil phase in low-permeability oil reservoirs, i.e., modified Buckley-Leverett (B-L) problem. The model is capable of generating a precise physical solution, accurately capturing both shock and rarefaction waves under the specified initial and boundary conditions, though the introduction of complicated physics increases the nonlinearity of the governing PDEs. The self-adaptive mechanism modifies the behavior of the deep neural network by simultaneously minimizing the losses and maximizing the weights. It, thus, can effectively capture the non-linear characteristics of the solution, thereby overcoming the existing limitations of PINNs. In these numerical experiments, the SA-PINNs demonstrated superior performance compared to other state-of-the-art PINN algorithms in terms of L2 error. Moreover, it was also achieved with a reduced number of training epochs. SA-PINNs can effectively model the dynamics of complex physical systems by optimizing network parameters to minimize the residuals of the PDEs.

Geothermal energy has the potential to decarbonize heating, cooling, and power production. However, efficient and sustainable exploitation is challenging due to the limited data, which restricts our ability to characterize and quantify the multi-scale, hierarchical geological structures of geothermal reservoirs. This study proposes a scenario-based data assimilation framework that enables efficient modelling of multiple geological scenarios and is coupled with flow and heat transfer simulations for uncertainty analysis. We demonstrate the framework on a synthetic but geologically consistent low-enthalpy geothermal case, where heat is produced from a doublet in a channelized fluvial sandstone reservoir. Using the open-source Rapid Reservoir Modelling (RRM) tool, we efficiently create multiple geological scenarios with different reservoir architectures constrained by wellpath facies data. Reservoir properties for each scenario are modelled using geostatistics to generate a geologically plausible and sufficiently diverse ensemble of realizations. These are subjected to heat and flow simulations using the open-source Delft Advanced Research Terra Simulator (open-DARTS) to quantify uncertainties in temperatures and pressures. Finally, ensemble smoother with multiple data assimilation (ESMDA) is employed to assimilate temperature and pressure profiles at the injection well, monitoring borehole, and production well across all ensemble realizations for the different geological scenarios. Based on the deviation of modelled and observed well temperature and pressure profiles, the framework falsifies geological scenarios with poor data assimilation performance that is unlikely to reflect the actual reservoir architecture, and identifies plausible scenarios that yield more reliable production forecasts. The workflow effectively constrains geological uncertainties and improves forecast reliability in geothermal systems.

Efficient geothermal resource development remains challenging due to inherent geological uncertainty and limited subsurface data. A proof-of-concept for a digital twin for a fluvial geothermal reservoir, similar to the Delft campus geothermal project, is presented. This digital twin has the aim to integrate geological scenario modeling, production simulation, uncertainty analysis, and data assimilation to mitigate operational risks, reduce maintenance costs, extend reservoir longevity, and enhance the overall sustainability of this project. In this contribution, we assess the efficiency of the ensemble smoother with multiple data assimilation (ESMDA) for subsurface property inversion of a fluvial geothermal system. First, we developed an efficient method that allows for the swift creation of multiple geological scenarios of channelized reservoir geometries, fully constrained to well information, using Rapid Reservoir Modeling (RRM). Next, we generated an ensemble containing multiple geological realizations for a given scenario representing the geothermal system using stochastic reservoir modelling. For a single scenario and its ensemble of stochastically generated property distributions, heat flow and production rates were simulated using the Delft Advanced Research Terra Simulator (DARTS). One of the ensemble members and its simulated production data were taken as the “truth” (or reference) case. ESMDA was then employed to invert the property distribution within the fluvial channels of all other ensemble members, using the “observed” temperature and pressure data along the injection and production well from the “truth” case. We also consider the presence of a monitoring borehole to analyze how additional monitoring data impacts the convergence of ESMDA. The simulation results of the posterior models demonstrated a significant reduction in root mean square error for temperature and pressure data which align more closely with the “observations” compared to the prior models. This outcome confirms the feasibility of applying ESMDA for data assimilation in fluvial geothermal systems, such as the Delft campus geothermal project.

Halite precipitation during CO2 injection can significantly reduce injectivity and impact long-term storage in saline aquifers and depleted reservoirs. However, the impact of geological porous media on salt precipitation and brine movement is not fully undersood. This study explores salt precipitation dynamics from the pore to core scale using microfluidic experiments and core-flooding tests. Microfluidic results reveal three distinct phases of salt deposition: slow evaporation, rapid evaporation, and complete dry-out. Heterogeneous pore structures retain more initial water, leading to localized salt accumulation due to capillary effects. Core-scale experiments show that permeability strongly influences salt penetration and porosity reduction. In high-permeability cores, salt fronts extend deeper into the rock, causing up to 70% porosity reduction. In contrast, heterogeneous cores experience limited salt penetration but increased surface accumulation, suggesting that capillary pressure and brine redistribution control final deposition patterns. These findings highlight the complex interactions between fluid flow, rock properties, and salt crystallization, which provides valuable insights for predicting injectivity loss and optimizing CO2 storage strategies. Understanding these mechanisms is essential for improving reservoir management and ensuring the long-term stability of geological CO2 sequestration projects.

The effective management of geo-energy systems heavily relies on robust modeling frameworks that integrate diverse simulation capabilities, including flow and transport, phase equilibrium, geochemistry and geomechanics. While a multiphysics simulation engine within a unified framework has its advantages, integrating specialized modeling packages often enhances viability. Efficient and seamless communication between these engines be- comes crucial for improving the performance and scalability of the integration. Advanced parametrization tech- niques can facilitate this integration by efficiently approximating and interpolating coupling data, ensuring both speed and accuracy. In this study, we compare the efficiency of different interpolation techniques used for the parametrization of complex many-component fluid systems in compositional simulation. We employ an Operator- Based Linearization (OBL) framework that leverages the general formulation of corresponding conservation laws. OBL effectively learns the operators required for assembly of the laws while interpolation delivers fast evalua- tion of operators and their derivatives for all physical states in a simulation domain. Multilinear interpolation is a simple and robust approach, yet it has poor scaling properties with respect to the dimension of the physical state. To alleviate interpolation costs in multiple dimensions, we study the performance and accuracy of other interpolation techniques, including linear interpolation with standard and Delaunay triangulation. Overall, this approach provides great flexibility, saves development costs and simplifies the incorporation of thermodynamics and geochemistry engines for precise modeling of phase equilibrium, reactive transport, dissolution-precipitation and kinetics of chemical reactions. This research extends the scalability of the OBL framework and addresses the challenges of high dimensionality in compositional modeling. Consequently, this approach holds significant potential for integrating various complex multiphysics problems, enabling the creation of more comprehensive digital twins for geo-energy systems management.

Long-term geothermal production is subject to considerable uncertainty due to limited data availability and inherent geological heterogeneity. While observation and data acquisition improve our understanding of the reservoir, they also contribute significantly to project costs. It is essential to identify the most informative observation strategy. In this study, we apply a previously developed scenario-based data assimilation framework that integrates rapid geological modelling, efficient numerical simulation, and Ensemble Smoother with Multiple Data Assimilation (ESMDA) to constrain uncertainties in reservoir properties and production forecasts to a synthetic but geologically realistic fluvial geothermal system and conduct a data worth analysis to evaluate the impact of different observations (production temperature and injection pressure, well temperature and pressure profiles, etc.) on uncertainty reduction. Results show that production temperature and injection pressure alone, though cost-effective, are insufficient to significantly reduce uncertainties in reservoir performance forecasts. In contrast, well temperature and pressure profiles exhibit substantially higher data worth, leading to much better-constrained predictions. Moreover, incorporating a monitoring borehole further constrains uncertainty by capturing subsurface dynamics between the injector and producer. These findings underscore the importance of monitoring pressure and temperature profiles in the wells of a geothermal doublet.

Direct Use Geothermal Systems (DUGS) are rapidly and densely deployed to meet the growing demand for renewable energy with less carbon emissions globally. The simulation of DUGS can provide a reservoir-scale understanding of geothermal resource assessment, where the geothermal system's lifetime and the injection well Bottom Hole Pressure (BHP) are used as performance indicators. However, there are inherent errors from numerical simulations of any engineering problems, due to approximating continuous partial differential equations by their discretized approximation in time and space. In this work, we establish an optimal numerical setup with reduced errors across the homogeneous, stratified and heterogeneous models for the simulation of a geothermal system. Next, we develop a standardized method for calculating recoverable Heat In Place (HIP) and an analytical solution for evaluating the HIP recovery factor across various geological models using a single forward simulation. We present reference examples on the design of DUGS simulations using the open-source software Delft Advanced Research Terra Simulator (open-DARTS). The open-DARTS platform enables accurate and efficient sensitivity and uncertainty analysis. Using Distance-Based Generalized Sensitivity Analysis (DGSA), we identify reservoir depth and discharge rate as the most influential parameters for geothermal projects across all three types of geological models.

Subsurface CO2 sequestration plays a crucial role in advancing carbon neutrality and supporting the transition to sustainable energy. However, the unique behavior of CO2, particularly during cold CO2 injection into depleted hydrocarbon reservoirs, poses challenges to wellbore injectivity. Addressing these challenges requires a numerical model that captures the complex interplay between wellbore dynamics and reservoir processes. In this work, we present DARTS-well, an open-source, fully coupled wellbore-reservoir model developed using the Operator-Based Linearization (OBL) technique. To this end, a transient, multi-segment, two-phase, non-isothermal wellbore model based on the Drift-Flux Model (DFM) is first developed and then coupled with the Delft Advanced Research Terra Simulator (DARTS) as the reservoir simulator. The model is demonstrated through test cases involving CO2 injection into a depleted reservoir, illustrating its potential to enhance the design and optimization of subsurface CO2 disposal systems.

Capillary pinning refers to the immobilization of CO₂ at capillary barriers when the uprising CO₂ pressure is lower than the capillary entry pressure of the overlaying pore throats. Also known as local capillary trapping, it has been proposed as a fifth geologic CO₂ storage mechanism, alongside structural, solubility, residual, and mineral trapping. Despite extensive research, the fragmented terminology surrounding capillary pinning has led to confusion, making it challenging to synthesize findings effectively. Often conflated with mechanisms such as residual and hysteresis trapping, capillary pinning is commonly underestimated or completely overlooked in reservoir-scale models. Furthermore, difficulties in characterizing and upscaling small-scale geologic heterogeneities that influence capillary pinning contribute to significant uncertainties, with estimates of CO₂ trapped via this mechanism ranging from 3 % to 100 % of total CO₂ trapped via capillary actions. This review explores the fundamental mechanisms, experimental findings, and modeling approaches for assessing CO₂ capillary pinning in carbon capture and storage (CCS). It seeks to bridge the gap between the reservoir engineering community, with its extensive expertise in hydrocarbon recovery but that needs adjustments for CCS applications, and the subsurface storage community, which stands to benefit from this knowledge but often lacks access to relevant literature. Additionally, the study identifies key research opportunities to advance the understanding of capillary pinning in sedimentary rocks, ultimately enhancing the efficacy and reliability of CCS operations.

The 11th Society of Petroleum Engineers Comparative Solution Project (shortened SPE11 herein) benchmarked simulation tools for geological carbon dioxide (CO2) storage. A total of 45 groups from leading research institutions and industry across the globe signed up to participate, with 18 ultimately contributing valid results that were included in the comparative study reported here. This paper summarizes the SPE11 results. A comprehensive introduction and qualitative discussion of the submitted data are provided, together with an overview of online resources for accessing the full depth of data. A global metric for analyzing the relative distance between submissions is proposed and used to conduct a quantitative analysis of the submissions. This analysis attempts to statistically resolve the key aspects influencing the variability between submissions. The study shows that the major qualitative variation between the submitted results is related to thermal effects, dissolution-driven convective mixing, and resolution of facies discontinuities. Moreover, a strong dependence on grid resolution is observed across all three versions of the SPE11. However, our quantitative analysis suggests that the observed variations are predominantly influenced by factors not documented in the technical responses provided by the participants. We therefore identify that unreported variations due to human choices within the process of setting up, conducting, and reporting on the simulations underlying each SPE11 submission are at least as impactful as the computational choices reported.

Modeling of multiphase flow and reactive mass transport in porous media remains a pivotal challenge in the realm of subsurface energy storage, demanding a nuanced understanding across varying scales. This review paper presents a comprehensive overview of the latest advancements in multiscale modeling techniques that address the inherent complexity of these processes. Three cutting-edge approaches are presented: hybrid multiscale simulation, which leverages both continuum and discrete modeling frameworks to enhance model fidelity; approximated physics, which simplifies complex reactions and interactions to expedite computations without significantly sacrificing accuracy; and machine-learning-assisted multiscale simulation, which integrates predictive analytics to refine simulation outputs. Each method presents distinct advantages and hurdles, collectively advancing the precision and computational efficiency of subsurface modeling. Despite the substantial progress, we recognize the persistent challenges, such as the need for more robust coupling techniques, the balance between model complexity and computational feasibility, and effectively combining machine learning with traditional physical models. Promising directions for future work are discussed to address these challenges, aiming to push the boundaries of current multiscale modeling capabilities.