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.

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.

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.

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.

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 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.

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.