Thermal compositional simulators perform multiphase flash calculations to determine the number of phases as well as mole fractions of components in different phases. This information is used for property calculations and as input to mass and energy balance residuals. Multiphase flash calculations are challenging and could fail in situations where the phase envelope is extremely narrow. For example, this often occurs around steam injection wells in thermal EOR recovery processes.
Recently, Heidari and Stone (2021) proposed a thermal stability check that replaces these multiphase flashes in thermal compositional simulations. The model was designed for cases where components are not soluble in the aqueous phase and water is the only component in the aqueous phase. Although this is a valid assumption for many thermal processes such as Steam-Assisted Gravity Drainage (SAGD), it is not accurate for operations with significant solubility of non-aqueous components in the aqueous phase such as CO2 injection. CO2 is highly soluble in water and in order to accurately simulate CO2 injection into depleted reservoirs or saline aquifers, CO2 aqueous solubility must be accounted for.
In this work, an extension of the Heidari and Stone (2021) stability check is presented to account for solubility of other components in the aqueous phase. It requires two sets of equilibrium ratios (k-values) between vapour-liquid and aqueous-vapor which are only functions of pressure and temperature. The method is used for phase appearance while disappearance is derived from a phase volume fraction (saturation) update after a linear solve; i.e. if a phase saturation goes below zero during the non-linear solver cycle then that phase disappears.
The new thermal stability check has been implemented in a reservoir simulator and tested rigorously for both synthetic and real field cases with varying degrees of component solubility in the aqueous phase. Physical results of simulations using this method are either identical to those with multiphase flashes or more accurate if flash failures occur. The new method is robust, easy to implement, computationally inexpensive and doesn’t have the complexity of flashes that account for component solubility in the aqueous phase.

Linear solvers account for among the largest portion of runtime in most kinds of numerical simulation applications. This also holds for reservoir simulations. Thus, efficient linear solvers methods such as System-AMG are a key to the successful application of reservoir simulations. It has been demonstrated earlier how different building blocks of multigrid can be combined to solver strategies that are bespoke to the application with certain kinds of simulations, such as Black-Oil, thermal or coupled geomechanics.
These methods, however, still comprise a lot of options for a fine-grained control. An optimal setting of detailed parameters is a rather volatile trade-off between robustness and computational efficiency. It is determined by individual properties of a particular simulation, computing environment and accuracy requirements. That is, while an efficiently designed overall solver strategy already significantly accelerates a simulation, quite some further potential for optimization remains in many cases. This, however, is hardly exploitable manually in all details.
Instead, we are proposing an autonomous control mechanism that can select parameters and methods individually. We use methods of machine learning via genetic algorithms, as it has been done earlier also in other applications. In our control mechanism, however, we combine this with a tree-based approach to limit the parameter optimization to search spaces that are considered reasonable in advance. This reduces the learning efforts. Surrogate-based techniques allow for transferring results from previous, possibly different runs to further guide the learning process.
A deep integration in the solver method allows for accessing all relevant data for decision and learning processes and helps to reduce overhead costs. It also allows for reducing the number of solver setups within a simulation run and guarantees robustness by quickly reacting to convergence break-downs.
We will demonstrate the benefits that such a control mechanism can provide in reservoir simulations and beyond. And we will discuss aspects of reproducibility of results.

This work proposes an original preconditioner coupling the Constrained Pressure Residual (CPR) method with block preconditioning for the efficient solution of the linearized systems of equations arising in fully implicit reservoir simulation flow models. This preconditioner has been specifically designed for Lagrange multipliers-based discretizations, like the Mixed Hybrid Finite Element (MHFE) or the Mimetic Finite Difference (MFD) scheme. Here, we focus on a MHFE-based discretization of the two-phase flow model in porous media, where the phase flux continuity is strongly enforced in the mass balance equation to stabilize the nonlinear solver convergence. While the system unknowns consist of cell and face pressures and cell saturations, the Jacobian matrix, which is non-symmetric, usually large and ill-conditioned, exhibits a 3x3 block structure. The simulator performance in transient simulation essentially pivots on the acceleration capability offered by the preconditioner and the interplay with the linear solver, usually a Krylov subspace method. This is a well-known problem in literature, where the CPR has been established as a robust and efficient tool, but most research focused on solutions for Two-Point Flux Approximation (TPFA)-based discretizations that do not immediately and readily extend to our problem formulation. Therefore, we designed a dedicated multi-stage strategy, inspired by the CPR algorithm, where a block preconditioner is used for the pressure part with the aim at exploiting the inner 2x2 block structure. Approximated decoupling factors of the Jacobian are used to recast the Schur complement. Their computation is performed in parallel, and, at an algebraic level, they are obtained through restriction operators, where the size of the restricted subspaces is either statically or dynamically adapted. A closer inspection reveals that recomputing the decoupling factors from scratch at each nonlinear iteration is not necessary, but they can be updated following the waterfront advancement, thus improving the preconditioner efficiency. The proposed preconditioning framework has been tested with broad experimentation, comprising both synthetic and realistic applications in Cartesian and non-Cartesian domains, highlighting its advantages and weaknesses.