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.

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.

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.

In the original published article, one of the affiliation details for corresponding author "Denis Voskov" were incorrectly given as "Department of Energy Resources Engineering, Stanford University, CA, USA" but should have been "Department of Energy Science and Engineering, Stanford University, CA, USA". It is updated in this correction. The original article has been corrected.

The open Delft Advanced Research Terra Simulator (open-DARTS) framework is an open-source reservoir simulation software. The open-DARTS focused on energy transition applications, such as geothermal energy production and carbon sequestration. It enables the modeling of compositional thermal flow, coupled with a geomechanical solver based on the Finite Volume discretization and adjoints method for inverse modeling. The open-DARTS supports different grid types (structured, corner-point geometry, and unstructured), discrete fracture networks, contact mechanics, and various thermal-chemical interactions. The recently proposed generic nonlinear formulation supports the most general nonlinear PDEs designed for various energy transition applications. The open-DARTS has been implemented in C++ and Python to optimize hardware utilization while ensuring flexibility. The most computationally expensive part is written in C++ and compiled into libraries, which are subsequently exposed to Python using pybind11. This allows the extension and overriding of C++ functions by user-defined Python code. For example, using only a Python interface, one can adjust a timestep strategy, nonlinear solver, or properties output. Besides, the Python interface of open-DARTS provides straightforward coupling with other Python-based numerical modeling packages, including the meshing, file storage, caching, and visualization modules. The open-DARTS core uses the advantages of C++ language, such as efficient low-level memory management, object-oriented programming, compile-time polymorphism, and parallelization with OpenMP. One of the advantages of open-DARTS is the Operator-Based Linearization (OBL) technique, which can resolve challenges associated with complex physics and reduce the computation time, especially for ensemble-based simulations. We would also like to share our experience on the project, repository, and the development workflow configuration using gitlab.com, including the build system (cmake), handling merge requests, automated testing in CI/CD pipelines, documentation management (gitlab.io), wiki utilization, and release publishing. Additionally, Python’s integration into open-DARTS offers the advantage of straightforward installation via PyPI and simplifies defining requirements for users who prefer to avoid compiling code from source files.

In this study, we conduct a comprehensive history matching study for the FluidFlower benchmark model. This benchmark was prepared and organized by the University of Bergen, the University of Stuttgart, and Massachusetts Institute of Technology, for promoting understanding of the complex physics of geological carbon storage (GCS) through in-house experiments and numerical simulations. This paper synthesizes the experiences of history matching the benchmark data encountered by the Delft-DARTS and CSIRO participants. History matching is first performed based on a low-dimensional-zonated structured model using a simple Poisson-like solver. The permeability of six facies and two faults is inferred in this stage to match the digitized concentration data. The history matching is then further enhanced to richer spatial and physical models to capture the spatial variation of permeability and buoyancy effects, using an unstructured grid. Efficient adjoint methods are used to evaluate the gradient used in the optimization of data misfits or equivalent Bayesian log-likelihoods. With efficient optimization methods available for both maximum a posteriori model inference and Randomized Maximum Likelihood methods for model uncertainty, we perform history matching using both binary and continuous concentration observations. The results show that the tracer plumes in the enriched model match the experimental plumes more accurately compared with the results from the parsimonious-zonated model. The history matching results based on the concentration observations provide more similar plume shapes compared with the case based on the binary observations. The permeability difference between the model before and after history matching reveals that the tracer plume zone and the high permeable zone are the regions of high sensitivity in terms of data misfit between the model response and observations. Surprisingly, CO ₂ concentration plume forecasts based on these history-matched models were not especially sensitive to the improvements observed in the enhanced model.

In this work, we present a unified framework for the simulation of CO₂ sequestration problems at various time and space scales. The parametrization technique utilizes thermodynamic state-dependent operators expressing the governing equations for the thermal-compositional-reactive system to solve the nonlinear problem. This approach provides flexibility in the assembly of the Jacobian, which allows straightforward implementation of advanced thermodynamics. We validate our simulation framework through several simulation studies including complex physical phenomena relevant to CCUS. The proposed simulation framework is validated against a set of numerical and experimental benchmark tests, demonstrating the efficiency and accuracy of the modeling framework for CCUS-related subsurface applications. Important physical phenomena resulting from the complex thermodynamic interactions of CO₂ and impurities with reservoir fluids can be accurately captured now in detailed dynamic simulation. The investigated simulation scenarios include a reproduction of lab experiments at the core scale, investigation of macro-scale analog model and simulation of large-scale industrial application. The simulation time can also span from hours to years among various applications. Complex thermal-compositional-reactive phenomena can be addressed at each of these space and time scales. The unified thermodynamic description allows us to perform all these simulations for a reasonable CPU time due to advanced parametrization techniques and efficient GPU capabilities in our in-house reservoir simulator DARTS.

In this work, we describe our decisions made to perform the FluidFlower simulation study and discuss various aspects of the benchmark that are different from our usual subsurface simulation practice. We will discuss the impact of various modeling choices on the outcomes of the simulation models, such as gridding, discretization, and solver strategies, and the lessons learned, taking into account the different conditions of the FluidFlower study compared to conditions commonly dealt with in subsurface simulation. We will start with a brief description of the DARTS framework utilized for compositional simulation, the thermodynamic and physical modeling related to the atmospheric CO2 -brine system, and the modeling workflow used in our benchmark submission. Additionally, we describe a custom nonlinear solver developed for the atmospheric benchmark conditions to improve convergence including the linear solver strategy since our default two-stage preconditioner does not perform effectively. To make meaningful comparisons between each of the modeling choices, we define a baseline model which is a simplified version of our setup in the main FluidFlower benchmark. The baseline model is then used to study the effect of Cartesian and unstructured meshes and a two-point flux approximation compared with a multi-point flux approximation for capturing the physics at play. We conclude our work with lessons learned and future recommendations.

Successful deployment of geological carbon storage (GCS) requires an extensive use of reservoir simulators for screening, ranking and optimization of storage sites. However, the time scales of GCS are such that no sufficient long-term data is available yet to validate the simulators against. As a consequence, there is currently no solid basis for assessing the quality with which the dynamics of large-scale GCS operations can be forecasted. To meet this knowledge gap, we have conducted a major GCS validation benchmark study. To achieve reasonable time scales, a laboratory-size geological storage formation was constructed (the “FluidFlower”), forming the basis for both the experimental and computational work. A validation experiment consisting of repeated GCS operations was conducted in the FluidFlower, providing what we define as the true physical dynamics for this system. Nine different research groups from around the world provided forecasts, both individually and collaboratively, based on a detailed physical and petrophysical characterization of the FluidFlower sands. The major contribution of this paper is a report and discussion of the results of the validation benchmark study, complemented by a description of the benchmarking process and the participating computational models. The forecasts from the participating groups are compared to each other and to the experimental data by means of various indicative qualitative and quantitative measures. By this, we provide a detailed assessment of the capabilities of reservoir simulators and their users to capture both the injection and post-injection dynamics of the GCS operations.