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.

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.

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.

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.

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.