# New initialization procedures from phase stability testing in three-phase flash calculations for CO2-hydrocarbon mixtures

Juan Heringer, Michiel Wapperom, Catinca Secuianu, Denis Voskov, Dan

Vladimir Nichita

► To cite this version:

Juan Heringer, Michiel Wapperom, Catinca Secuianu, Denis Voskov, Dan Vladimir Nichita. New initialization procedures from phase stability testing in three-phase flash calculations for CO2-hydrocarbon mixtures. Fluid Phase Equilibria, 2026, 604, pp.114653. 10.1016/j.fluid.2025.114653. hal-05422454

HAL Id: hal-05422454

https://hal.science/hal-05422454v1

Submitted on 18 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# New initialization procedures from phase stability testing in three-phase flash calculations for CO₂-hydrocarbon mixtures

Juan Heringer<sup>a</sup>, Michiel Wapperom<sup>b</sup>, Catinca Secuianu<sup>c</sup>, Denis Voskov<sup>b,d</sup>,
Dan Vladimir Nichita<sup>a, e</sup>

## ARTICLE INFO

Keywords:

Phase stability

Multiphase flash

Stepwise strategy

Initialization

Iterations

## ABSTRACT

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.

## 1. Introduction

Phase equilibrium calculations play an important role in a wide variety of applications in chemical and petroleum engineering. Beyond the traditional CO₂ injection used in oil recovery processes, in recent years, a great attention has been paid to applications in carbon capture, utilization and storage (CCUS), processes that are intended to reduce the amount of CO₂ in the atmosphere. In those applications, complex phase

behavior occurs, and three or more phases might be present [1], making the compositional simulation more challenging. In reservoir compositional simulators, a failure in correctly describing the phase behavior will result in either an erroneous physical modeling of the problem or in significant error propagation that may compromise a simulation run. Moreover, in field-scale simulations, the phase equilibrium routine is called a significant number of times, consuming an important part of the computational time. Therefore, both robustness and efficiency are of

<footnote><sup>a</sup> CNRS UMR 5150, Laboratoire des Fluides Complexes et leurs Réservoirs, Université de Pau et des Pays de l'Adour, Pau, France</footnote>

<footnote><sup>b</sup> Department of Geoscience and Engineering, TU Delft, the Netherlands</footnote>

<footnote><sup>c</sup> Department of Inorganic Chemistry, Physical Chemistry and Electrochemistry, Faculty of Chemical Engineering and Biotechnologies, National University of Science and Technology "Politehnica" Bucharest, Romania</footnote>

<footnote><sup>d</sup> Department of Energy Science & Engineering, Stanford University, Stanford, CA, USA</footnote>

<footnote>* Corresponding author.
E-mail address: dnichita@univ-pau.fr (D.V. Nichita).</footnote>
utmost importance.

Phase stability testing aims to ascertain whether a chemical system maintains equilibrium as a single phase or splits into two or more phases. It uses the tangent plane distance (TPD) function [2–4], and the stationary points are calculated in order to verify the stability of the mixture and to provide initial guesses for flash calculations. Michelsen’ [5] “two-sided” initialization (using ideal equilibrium constants) works well in two-phase vapor-liquid flash calculations. However, for phase equilibria involving three or more equilibrium phases, or for liquid-liquid two-phase equilibrium, stability testing is required for initialization of phase splitting calculations [5]. The traditional approach for multiphase equilibrium consists of combining phase stability analysis and multiphase flash calculations in a stepwise manner [5]. First, a stability test is carried out for the feed and the composition of the global minimum of the TPD function, and the feed composition is used to generate initial K-values for the two-phase flash. After performing the two-phase split, the stability of one of the equilibrium phases is tested, and the compositions of the global minimum of the TPD function, and of the two equilibrium phases are used to generate two sets of K-values to initialize the three-phase split. The sequence is repeated until a stable state of the mixture is detected. Detailed worked examples for three-phase flash examples can be found in textbooks [6,7].

Even if a negative TPD is found from the first tried initial guess (thus it is certain that at least an additional phase will form), this minimum may be a local one, thus the search must continue with stability runs from all the initial guesses, to enhance the chances of finding the global minimum of the TPD function. Therefore, reducing the number of stability testing calls is a key point in developing an efficient stability-flash algorithm, which can be used in compositional reservoir simulation. Various sets of initial guesses and strategies for stepwise calculations and for addressing the phase stability and phase split problems were presented, with applications for CO₂-hydrocarbon mixtures [8–18].

The phase stability testing problem, consisting of a constrained minimization of the TPD function with respect to mole fractions, was formulated by Michelsen as an unconstrained minimization of a modified TPD function (having the same stationary points with the same sign as the original function) with respect to formal mole numbers [4]. Multiphase split calculations consist of an unconstrained minimization of the Gibbs free energy with respect to phase mole numbers [5]. Alternatively, natural logarithms of the equilibrium ratios can be used as independent variables [5]; in this case, the negative flash [19,20] option is available. In this work, highly robust and efficient phase equilibrium routines are used, based on successive substitution iterations (SSI) in early iteration stages, followed by Newton iterations with modified Cholesky factorization [21,22] and line search [23], in both stability and two- and three-phase flash calculations.

Recently, Huang et al [14] presented an initialization strategy for CO₂-hydrocarbon mixtures, called Multiple Identification (MI), in which all compositions of minima with negative TPDs in stability testing for feed composition are used to identify a possible three-phase equilibrium, leading to the elimination of a number of two-phase splits and stability analysis.

In this work, two new initialization strategies for two- and three-phase flashes are proposed and analyzed in detail for several CO₂-hydrocarbon mixtures. The first method, called Improved Stepwise Initialization (ISI), uses the same stepwise procedure as in the conventional multiphase flash but with improved construction of initial guesses for the K-values. In the second one, called Improved Multiple Initialization (IMI), a direct three-phase flash is performed using information from stability testing on feed composition if the TPD function has three minima (with negative, positive or zero TPDs) and two sets of initial equilibrium constants can be formed. An important element of novelty consists of using also the composition of minima with a strict positive values of the TPD function. The IMI method avoids a significant number of stability testing calls and two-phase flashes. Moreover, unlike in the MI method from Huang et al [14], the feed composition is used to

construct K-values only if the trivial solution corresponds to a minimum of the TPD function. Another difference with respect to Ref [14] is that in this work, all initial guesses honor the component material balance, and some unnecessary three-phase calculations performed in the MI method are avoided.

The paper is structured as follows. First, we give a brief presentation of robust and efficient calculations used in phase stability testing and two-phase and multiphase flash calculations, then the conventional method and the two proposed initialization procedures from phase stability testing methods (which use information from all minima of the TPD function) are presented. The new procedures are validated on four CO₂-hydrocarbon mixtures from the literature, with a detailed analysis in three- and two-phase domains before concluding. Two Appendices present the flowcharts of the conventional and proposed methods and the input data for the test examples.

## 2. Stability testing and phase split

In this section, the governing equations and the robust and efficient calculation framework used in phase stability testing and two-phase and multiphase flash calculations are presented.

### 2.1. Phase stability testing

The modified dimensionless tangent plane distance (TPD) function (Michelsen [4]) is

$$D(\mathbf{Y}) = \left(1 - \sum_{i=1}^{nc} Y_i\right) + \sum_{i=1}^{nc} Y_i(\ln Y_i + \ln \varphi_i(\mathbf{Y}) - h_i) \quad (1)$$

where $Y_i$ can be formally viewed as mole numbers and $h_i = \ln z_i + \ln \varphi_i(z_i)$.

This function has the same sign and same stationary points as the TPD function. An unconstrained minimization problem is solved, with only bound $Y_i > 0$. The stability analysis problem consists of minimizing the TPD functions, starting from a number of initial guesses.

The gradient of D is

$$g_i = \frac{\partial D}{\partial Y_i} = \ln Y_i + \ln \varphi_i - h_i; i = 1, nc \quad (2)$$

The Newton iteration equation is

$$H \Delta Y = -g \quad (3)$$

where

$$H_{ij} = \frac{\partial^2 D}{\partial Y_i \partial Y_j} = \frac{\partial g_i}{\partial Y_j} = \frac{\delta_{ij}}{Y_i} + \frac{\partial \ln \varphi_i}{\partial Y_j}; i, j = 1, nc \quad (4)$$

Michelsen [4] also proposed the independent variables $\alpha_i = 2\sqrt{Y_i}$ (which corresponds to an optimum scaling [24]) In this case, the gradient is

$$g_i^* = \frac{\partial D}{\partial \alpha_i} = g_i \sqrt{Y_i} \quad (5)$$

and the Hessian matrix is

$$H_{ij}^* = \frac{\partial^2 D}{\partial \alpha_i \partial \alpha_j} = \delta_{ij} + \sqrt{Y_i} \sqrt{Y_j} \frac{\partial \ln \varphi_i}{\partial Y_j} + \frac{1}{2} \delta_{ij} g_i^* \quad (6)$$

where the second term corresponds to a generally low effective rank matrix and the last term vanishes at the solution.

The iteration equation of the first-order SSI method (Michelsen [4]) is obtained directly from the zero-gradient equation, $g(\mathbf{Y}) = 0$, that is

$$\ln Y_i = h_i - \ln \varphi_i(\mathbf{Y}) \quad (7)$$

The SSI method is extremely robust but can be extremely slow and a
switch is performed to the second-order Newton method; this takes advantage of the strengths and avoids disadvantages of each method.

## 2.2. Multiphase equilibrium

The objective function for multiphase flash calculations is the dimensionless Gibbs free energy, given by

$$G = \sum_{i=1}^{nc} \sum_{k=1}^{np} n_{ik} \ln f_{ik} \quad (8)$$

An unconstrained minimization of the function G is performed with respect to the mole numbers of np-1 phases, within their physical bounds $0 \le n_{ik} \le z_i$, $i = 1, nc$; $k = 1, np$; $k \ne R$ (where R indicates the reference phase). For a component i, the reference phase R is taken to be the one with the highest molar amount [5,9]

The elements of the gradient vector g are:

$$\frac{\partial G}{\partial n_{ik}} = \ln K_{ik} + \ln \phi_{ik} - \ln \phi_{iR}; \quad i = 1, nc; k = 1, np; k \neq R \quad (9)$$

where the equilibrium constants are

$$K_{ik} = x_{ik} / x_{iR} = \phi_{iR} / \phi_{ik}; \quad i = 1, nc; k = 1, np; k \neq R \quad (10)$$

The SSI iteration updates the equilibrium constants from

$$\ln K_{ik} = \ln \phi_{iR} - \ln \phi_{ik}; \quad i = 1, nc; k = 1, np; k \neq R \quad (11)$$

The iteration equation for the Newton method is

$$H \Delta n = -g \quad (12)$$

where the vector n collects the mole numbers of the independent phases and H is the Hessian matrix, with elements

$$H_{ij} = \frac{\partial^2 G}{\partial n_{im} \partial n_{jk}} = U_{ij} + \Phi_{ij}; \quad i, j = 1, nc; k, m = 1, np; k, m \neq R \quad (13)$$

where

$$U_{ij} = \frac{\delta_{mk}}{\theta_m} \left( \frac{\delta_{ij}}{x_{im}} - 1 \right) + \frac{1}{\theta_R} \left( \frac{\delta_{ij}}{x_{iR}} - 1 \right) \quad (14a)$$

and

$$\Phi_{ij} = \delta_{mk} \frac{\partial \ln \phi_{im}}{\partial n_{jm}} + \frac{\partial \ln \phi_{iR}}{\partial n_{jR}} \quad (14b)$$

If the natural logarithms of equilibrium constants, ln K<sub>ik</sub>, are taken as independent variables [5,25] the nonlinear system in ln K is

$$J \Delta \ln K = -g \quad (15)$$

$$\text{where } J_{ij} = \frac{\partial^2 G}{\partial n_{im} \partial \ln K_{ik}}, \quad i, j = 1, nc; k, m = 1, np; k, m \neq R.$$

In this case, the Rachford-Rice system of nonlinear equations must be solved

$$R_j = \sum_{i=1}^{nc} \frac{z_i (K_{ij} - 1)}{E_i} = 0; \quad j = 1, np - 1; j \neq R \quad (16)$$

where

$$E_i = 1 + \sum_{\substack{k=1 \\ k \neq R}}^{np} \beta_k (K_{ik} - 1); \quad i = 1, nc \quad (17)$$

in which β<sub>k</sub> are the phase mole fractions.

Michelsen [26] proposed to replace the resolution of the nonlinear system by the minimization of a convex function. Written in terms of constant equilibrium ratios as proposed by Leibovici and Nichita [27], Michelsen's function is

$$Q(\beta) = - \sum_{i=1}^{nc} z_i \ln E_i(\beta) \quad (18)$$

The constrained minimization of Q in np-1 variables, subject to the linear inequality constraints E<sub>i</sub> > 0 allows the negative flash option [19, 20, 28]. The Rachford-Rice equations correspond to the zero-gradient conditions of the function Q. In this work, we use a Newton method with line search to minimize Q, using the restricted feasible domain proposed by Okuno [29]; an inexact line search procedure is used, as described in Ref [24].

However, in two-phase flashes, it is not convenient to minimize the function Q. In this case, the Rachford-Rice equation is efficiently solved using convex transformations [30] and restricted solution windows [31], with a guaranteed convergence of second-and third-order Newton methods without any interval control.

The Euclidian norm of the gradient vector

$$S_g = \left( \sum_{i=1}^{niv} g_i \right)^{0.5} \quad (19)$$

is used in the stopping criteria, where niv is the number of independent variables, that is, nc in phase stability analysis and two-phase flashes and 2xnc in three-phase flashes. Iterations are stopped when S<sub>g</sub> < e.

## 2.3. Calculation framework for phase stability and multiphase flash

In this work, highly robust calculation procedures are used, based on successive substitution iterations (SSI) in early iteration stages, followed by Newton iterations with modified Cholesky factorization (to ensure a descent direction) and line search (to ensure a decrease of the objective function), in both stability and flash calculations. A similar calculation framework (except SSI, which cannot be used in volume-based equilibrium calculations), was used for volume-based phase stability testing at temperature and volume specifications [32] and at pressure and temperature specifications [33] and for two-phase volume-based flash at pressure and temperature specifications [34] (except SSI, which cannot be used in volume-based equilibrium calculations). Here, the algorithms are adapted for three-phase equilibrium calculations. We consider a generic function F (can be the TPD function or the Gibbs free energy) of independent variables ξ (can be Y<sub>i</sub> or α<sub>i</sub> for stability testing or mole numbers for flash calculations). The Newton iteration is

$$H \Delta \xi = -g \quad (20)$$

### 2.3.1. Modified Cholesky factorization

During iterations, the Newton method may fail to produce a descent direction (as revealed during a Cholesky factorization if negative diagonal elements are obtained). In the modified Cholesky factorization [21, 22], if some eigenvalues become negative, the Hessian matrix is modified using a diagonal correction, that is, by adding some diagonal elements such that becomes safely positive definite. The Newton iteration equation becomes

$$(H + E) \Delta \xi = -g \quad (21)$$

where E is a diagonal matrix.

In this work, the modified Cholesky factorization version of Schnabel and Eskow [22,35,36] is used.

### 2.3.2. Line search

The independent variables ξ = (ξ<sub>1</sub>, ..., ξ<sub>niv</sub>)<sup>T</sup> are updated in Newton iteration by

$$\xi^{(k+1)} = \xi^{(k)} + \lambda_k \Delta \xi^{(k)} \quad (22)$$

where Δξ<sup>(k)</sup> = −H<sup>−1</sup>g<sup>(k)</sup> is the search direction and the scalar λ<sub>k</sub> > 0 is the step length.
Given a descent direction, a two-stage line search procedure is used
(see Refs [32–34]):

i) In a first stage, it is checked whether the updated variables lie within the feasible region for the full Newton step ($\lambda_k = 1$). If some of the bounds are violated, the step length must be reduced, using $\lambda_{max} = \min_i (\lambda_i)$, that is, the maximum step length keeping the variables within their physical bounds, where $\lambda_i$ are the distance on the search direction at which a bound is violated. At the end of this step $\lambda_k = \lambda_{max} \le 1$ and all required inequalities are satisfied, namely $Y_i > 0, i = 1, nc$ (or $\alpha_i > 0, i = 1, nc$) for stability testing and $0 \le n_{ik} \le z_i; i = 1, nc; k = 1, np; k \ne R$ for multiphase flash calculations. This step is computationally inexpensive (no function evaluations are required);

ii) A second stage is performed if $F(\xi^{(k)} + \lambda_k \Delta\xi^{(k)}) > F(\xi^{(k)})$, that is, the objective function is not decreased. A one dimensional minimization of $F$ (the objective function is considered as a function depending only of $\lambda_k$) is performed. An efficient inexact line search using quadratic and cubic backtracking [37] is used. A detailed description can be found in Ref [24].

The procedure described above ensures a sequence of decreasing
objective functions during iterations, eventually converging to a local
minimum.

3. The multiphase strategy

The classical multiphase equilibrium calculations are based on
alternating stability testing and phase split calculations. The stability
test predicts whether or not the mixture is stable, and it also provides
initial guesses for flash calculations.

3.1. Stepwise identification method (SI)

The Stepwise Identification (SI) method, that is, the conventional
strategy, consists of the following steps:

1. Perform one phase stability test on the feed composition, starting with $nc_p+4$ different initial guesses, where $nc_p$ is the number of pure components used to obtain initial guesses of K-values. Store the different minima found in the stability test. In case at least one minimum with negative TPD is found, proceed to **step 2**, otherwise, finish the algorithm with a one-phase solution with feed composition.

2. Select the lowest minimum found to form a set of K-values, using the composition related to the minima and the feed. These equilibrium ratios will serve as initial guesses for the two-phase split calculation. Next, follow **step 3**.

3. Initialize two-phase split calculation with the set of equilibrium ratios constructed in **step 2**. Phase fractions and compositions of the two phases are saved. Proceed to **step 4**.

4. Start a phase stability test on one of the phases using again $nc_p+4$ different initial guesses. Save the distinct minima obtained from stability tests. In case a minimum with negative TPD is noticed, go to **step 5**; if not, **finish** with a two-phase solution with compositions already calculated.

5. Construct two sets of K-values via the combination of the composition related to the minimum which has the lowest TPD, and the compositions of the two phases calculated in **step 3**. Go to **step 6**.

6. Three-phase split calculation is performed and information on phase fractions and mole fractions is stored. **Finish** with a three-phase solution.

The flowchart describing the conventional strategy is given in
Fig. A1.

3.2. Improved stepwise identification (ISI)

An Improved Stepwise Identification procedure (ISI) is proposed. Similar to the conventional method, the ISI procedure goes through all the calculation levels until reaching the stable state of the mixture, but in this case, improved initialization procedures are carried out before each phase split calculations. The key difference with respect to the SI method is that initial equilibrium ratios for two-phase split calculation are built using the two minima that have the lowest TPD, if two or more different TPDs are obtained. In this case, trivial solutions and minima with positive TPDs are allowed to be used for construction the K-values. Following the two-phase split calculation, a stability test is performed on one of the phases and two sets of K-values are created based on the number of minima found in the second stability test. In case three minima are found, all of them are used to initialize the three-phase split. The main point of this method is to show that even though the procedure is similar to the conventional one, it has improved efficiency just by using all information from stability tests to build initial guesses for phase split calculations.

The ISI method consists of the following steps:

1. Perform one phase stability test on the feed composition, starting with $nc_p+4$ different initial guesses. Store the different minima found in the stability test. In case at least one minimum with negative TPD is found, proceed to **step 2**; otherwise, finish the algorithm with a one-phase solution with feed composition.

2. Use the composition of the two lowest minima encountered, which include ones that are the trivial solution and ones with positive or negative TPD, to form a set of K-values that will serve as initial guesses for the two-phase split calculation. Such a procedure is consistent with the geometrical interpretation of the flash problem, where the group of minima obtained in the stability test has compositions that are relatively close to the final flash solution compositions, which yields a better initial guess and thus a smaller number of iterations in phase split calculations. Next, follow **step 3**.

3. Start a two-phase split calculation with the set of equilibrium ratios constructed in **step 2**. Phase fractions and compositions of the two phases are saved. Proceed to **step 4**.

4. A phase stability test is performed on one of the phases using again $nc_p+4$ different initial guesses. Record the distinct minima obtained from stability tests. In case instability was detected, go to **step 5**; if not, **finish** with a two-phase solution with compositions already calculated.

5. Construct two sets of K-values by using all the minima that have the lowest TPD. If three different minima are present, use all of them to construct the two pairs of equilibrium ratios. In case only two are noticed, form the two groups of equilibrium ratios considering the two minima and the composition of the trial phase. Finally, if only one minimum is encountered, construct the K-values with its composition and the compositions of the two phases obtained in **step 3**. Continue in **step 6**.

6. A three-phase split calculation is performed. **Finish** with a three-phase solution.

The flowchart of ISI method presented in Fig. A2 gives more details
about the implementation.

Concerning the choice of compositions used in Step 2, for a mixture
in a two-phase state, if there are precisely two negative minima, the
conditions are inside the spinodal in the P-T plane and the trivial solu-
tion corresponds to a maximum. The initial K-values are generated using
both trial-phase compositions corresponding to the minima, instead of
using the global minimum and the feed composition. This was also
suggested earlier by Whitson and Brule [38]. If there are three negative
minima, the possible choices for generating initial K-values for a
two-phase flash are to either use the global minimum (or the minimum
with the second lowest TPD) and the feed composition or to use two of
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Sour gas mixture results at T = 178.8 K (a) phase boundaries; (b) total number of minima found in one-phase stability test; (c) number of minima with negative TPD found in one-phase stability test; (d) number of minima with zero TPD found in one-phase stability test.</figcaption></figure>

the negative minima (the ones that have the lowest TPD). Extensive testing indicated that the latter choice is globally better.

### 3.3. Improved multiple identification (IMI)

The Improved Multiple Identification method (IMI) has the same principle as in ISI, by using all information from the stability test to reduce the number of iterations in phase split calculations. For this method, a kind of "bypass" procedure is used if at least three minima are found in stability testing on feed and a direct three-phase split calculation is initiated, as in Ref [14], by constructing two sets of K-values from the compositions of the three minima.

The main steps of the IMI procedure are as follows:

1. Perform one phase stability test on the feed composition starting with $nc_p+4$ different initial guesses. Store the different minima found in the stability test. In case at least one minimum with negative TPD is found, proceed to **step 2**, otherwise finish the algorithm with a one-phase solution with feed composition.

2. In case three distinct minima are found, use the composition of the three minima encountered (including the trivial solution positive

TPDs) to form two sets of K-values which will serve as initial guesses for the three-phase split calculation, and proceed to **step 6**. Otherwise, create one set of K-values with the two minima and initiate **step 3**.

3. Start a two-phase split calculation. Phase fractions and compositions of the two phases are saved. If the last step was **step 2**, proceed to **step 4**. Otherwise, if the last step was **step 6**, **finish** with a two-phase solution.

4. Phase stability test is performed on one of the phases using $nc_p+4$ different initial guesses. Save the distinct minima obtained from stability tests. In case a minimum with negative TPD is noticed, go to the **step 5**, if not, **finish** with a two-phase solution with compositions already calculated.

5. Construct two sets of K-values by using the minima that have the lowest TPD. If three different minima are present, use all of them to construct the two pairs of equilibrium ratios. In case only two are noticed, form the two groups of equilibrium ratios considering the two minima and the composition of the trial phase. Finally, if only one minimum is encountered, construct the K-values with its composition and the compositions of the two phases obtained in **step 3**. Continue in **step 7**.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Total number of multiphase split iterations for sour gas mixture results in terms of SSI phase split iterations at T = 178.8 K (a) conventional; (b) ISI; (c) IMI.</figcaption></figure>

6. Three iterations of the three-phase flash are performed. In case an unphysical composition is identified, use the composition of the two phases that have positive phase fractions to build a set of K-values and proceed to step 3. Otherwise, save the compositions of all phases and continue calculations in step 7.

7. Perform three-phase split calculation. Finish with a three-phase solution.

The flowchart of the IMI method is given in Fig. A3.

Note that if we start in the three-phase domain with three nontrivial TPDs, the method converged to a three-phase solution in all test examples. In this case, the two-phase split calculation and the two-phase stability calls are avoided, which leads to a more efficient flash algorithm. Moreover, for the conditions detailed in this paper (that is, stand-alone calculations for CO₂-hydrocarbon mixtures with up to three phases), the K-values built from the compositions of the three minima reduced the number of iterations in the three-phase split (as will be detailed in the next section) and no false solutions were encountered.

### 3.4. Initial guesses for phase stability testing

The initialization strategy follows the recommendation in Ref [5]

(the “two-sided” methodology), as well as the additional initial guesses recommended in Ref [8,9,12,14,15]. Therefore, there are nc+4 initial guesses; four of them are defined as: $Y_i^{(1)} = z_i K_{iW}$, $Y_i^{(2)} = z_i/K_{iW}$, $Y_i^{(3)} = z_i K_{iW}^{1/3}$ and $Y_i^{(4)} = z_i/K_{iW}^{1/3}$, where $K_{iW}$ is the Wilson relation [39] for “ideal” phase equilibrium ratios

$$K_{iW} = \frac{P_{ci}}{p} \exp \left[ 5.37(1 + \omega_i) \left( 1 - \frac{T_{ci}}{T} \right) \right] \quad (23)$$

For the rest of the initial guesses, nearly pure composition for each nc components present in the mixture are defined as

$$\left\{ \begin{array}{l} K_i = c/z_i; c > 0.9 \\ K_j = 0.1/(nc_p - 1)/z_i; j \neq i \end{array} \right. \quad (24)$$

In practice, however, it is not necessary to obtain initial guesses from all the pure components; instead, a reduced number of selected components, $nc_p \le nc$, can be used. At least the key components must be selected, in this case the CO₂.

### 4. Results and discussion

The proposed methods are tested on four mixtures from the literature [8]: a sour gas at T = 178.8K, Maljamar reservoir oil at T = 305.35K,
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Sour gas mixture results at *P* = 20 bar and *T* = 178.8 K (a) phase fraction distribution; (b) number of minima found in stability test on feed composition; (c) number of phase split iterations; (d) number of stability calls.</figcaption></figure>

Maljamar separator oil at *T* = 305.35K, and North Ward Estes (NWE) oil at *T* = 301.48K, all of them combined with different amounts of CO₂. The two proposed methods, ISI and IMI are compared with the conventional procedure (SI). Initialization in stability test follows Li and Firoozabadi [8], in which multiple initial guesses are tried in order to increase the chances of finding all minima of the TPD in the compositional space. The Peng-Robinson equation of state [40,41] is used in all examples. In the stopping criterion, the tolerance for the gradient norm in both stability testing and flash calculations is ε = 10⁻¹⁰. The switch from SSI to Newton iterations is performed if the gradient norm falls below εsw = 0.1.

### 4.1. Problem 1 - sour gas - CO₂ mixture

Composition, component properties, and non-zero BIPs are listed in Table B1. In Fig. 1a, the P-Z phase envelope discriminating the number of phases is depicted. The same phase distribution was found in Li and Firoozabadi [8], Petitfrere and Nichita [9], and Huang et al [14]. This mixture has a complex phase distribution with three regions with one phase, three areas with two phases, and a wide region with three phases, as presented in Fig. 1a. The number of minima found in stability test on

feed composition is drawn in Fig. 1b (total), Fig. 1c (with negative TPD), and Fig. 1d (trivial). In red, three distinct TPDs are present, and since enough information is obtained, IMI tries automatically to generate K-values from the compositions of three points and initiate a direct three-phase split calculation. The two proposed methods presented are using information not only from minima with negative TPD, but also from all the minima found, which can be negative, positive or trivial. Therefore, as compared to the methods that only use the negative TPDs, the proposed methods perform better in a wider domain.

Proceeding with convergence analysis for the sour gas mixture, Fig. 2 shows the total number of multiphase split iterations for different pressures and CO₂ concentrations. Fig. 2a depicts the number of iterations for the conventional method, which is similar to what was presented by Petitfrere and Nichita [9]. The color maps for the ISI method in Fig. 2b, and IMI method, in Fig. 2c, show a significant reduction in the number of iterations, in both two- and three-phase regions, as compared to the conventional method. In SI, a ridge seems to be forming in the middle of the two-phase L1-L2 and in the three-phase region. On the other hand, in both proposed methods, this elevation is flattened and reaches values even lower than the regions near the phase boundary; in IMI, this decay is even more pronounced in the three-phase region. In
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Comparison of initial and converged K-values for the CO₂-sour gas mixture at *P* = 20 bar and *T* = 178.8 K, (a) for SI in the first set of ratios; (b) for IMI in the first set of ratios; (c) for SI in the second set of ratios; (d) for IMI in the second set of ratios.</figcaption></figure>

the two-phase liquid-liquid region, IMI requires three more iterations than ISI, due to the three-phase flashes performed to test if a three-phase state will actually occur. In this case, an unphysical result is found, and a two-phase split calculation is started with K-values formed from the three-phase split test. By doing this, not only a better initial guess is obtained for the two-phase split, but also information from the three-phase split test is not lost (as in the MI method from Ref [14]), and a number of stability tests is avoided.

In the following, six different pressures are selected for the sour gas mixture at 178.8 K for detailed study in different regions for different pressures and CO₂ concentrations.

### 4.1.1. Phase equilibrium at 20 bar and 178.8 K

At a pressure of 20 bar, the phase fraction distribution calculated using the three methods is presented in Fig. 3a. At lower concentrations of CO₂, a vapor (V)-liquid (L1) region is formed, and as the carbon dioxide fraction increases, a third phase appears (a second liquid phase L2). As the concentration further increases, the liquid phase L1 is reduced until it disappears at a CO₂ concentration of about 85 %, starting a two-phase region. A further increase in the CO₂ amount forces the disappearance of the vapor phase, the mixture remaining only in the

liquid L2 phase. All methods give the same results, which signify that the two proposed methods are consistent with the SI procedure.

Details on the number and nature of minima found in phase stability analysis (PSA) on feed composition are depicted in Fig. 3b. In the middle of the domain, marked in dashed red lines, the trivial solution is not a minimum, and three negative TPDs are identified. Note that whenever three minima are found, IMI is initiated. If this method was based only on the number of minima with negative TPDs, IMI would be used in a smaller region. In this case, it would be losing the chance to enhance the flash calculation efficiency, since the three minima area, bordered by a blue dashed line, is wider than that of the three minima with negative TPDs. When a total of two minima are found and one of them is the trivial solution, ISI, IMI and SI have equivalent procedures.

The total number of phase split iterations as a function of CO₂ concentration is given in Fig. 3c. The SI has the highest number of phase split iterations in a significant portion of the domain, where no trivial solution is found from PSA. At the same time, IMI generally has the best performance among the three methods inside the three minima domain in the blue dashed line. The exception is the region at around 50 % of CO₂, where ISI and IMI reach the same number of iterations. The flattening of the triangular shape of the number of iterations in SI by ISI and
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Sour gas mixture results at *P* = 25 bar and *T* = 178.8 K (a) phase fraction distribution; (b) number of minima found in stability test on feed composition; (c) number of phase split iterations; (d) number of stability calls.</figcaption></figure>

IMI clearly appears in Fig. 3c. The advantage of IMI is not only the improved initial guesses for the equilibrium ratios but also the fact that two-phase PSA and two-phase flashes are avoided. The number of stability calls is presented in Fig. 3d In a significant portion of the domain, where there are three minima from the PSA on feed, the number of phase stability calls in IMI is the lowest. This happens due to a direct initialization of the three-phase split.

To better understand why the performance of the SI is considerably poorer than IMI in the three-phase domain, a study on the K-values distribution is performed. In the SI method, initial K-values are built with the composition of the two phases calculated from the MS and one composition related to the lowest TPD calculated in the PSA performed on one of the phases. In IMI, the three minima found in PSA are utilized to create the equilibrium ratios necessary to start the three-phase split calculation. In Fig. 4, two sets of equilibrium ratios obtained after flash convergence, in blue, are compared with the initial two sets of K-values obtained with ISI and IMI, in black. In Fig. 4a, initial K-values for the SI method have a good agreement with the K-values at the solution only near the phase boundary. Moreover, they reach the highest deviation at around 55 % of CO₂, which is exactly the concentration at which the SI has the highest number of iterations in Fig. 3c. In Fig. 4b, the same comparison is made between initial and final K-values, but the method

used is the IMI. Initial K-values in this case are much closer to the solution than the ones from SI. In the same figure, a line representing the IMI limit is presented. This line separates the part where a direct three-phase split is performed, on the right, and the part where no direct three-phase split can be launched, on the left. For the second set of K-values, presented in Fig. 4c and Fig. 4d, the deviation in the SI is also more pronounced than in the IMI. The second set, however, has a lower deviation in general in both SI and IMI.

### 4.1.2 Phase equilibrium at 25 bar and 178.8 K

At a higher pressure, 25 bar, Fig. 5 presents the phase fraction distribution, information about the minima found in PSA on feed composition, the number of phase split iterations necessary to reach the solution for various concentrations of CO₂, and the total number of stability calls. The results are similar to the ones at 20 bar, but with a slightly smaller three-phase region, as presented in Fig. 5a, and similar behavior of the TPDs, as displayed in Fig. 5b The variation in terms of number of iterations is also minor, but a higher value is observed in SI procedure, and an increased difference between ISI and IMI, as drawn in Fig. 5c. IMI requires less iterations as compared with SI and ISI, except in the case that not enough information is available to build the two sets of K-values, or when enough information is available but the solution is
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Sour gas mixture results at *P* = 30.20 bar and *T* = 178.8 K (a) phase fraction distribution; (b) number of minima found after performing stability test on feed composition; (c) number of phase split iterations; (d) number of stability calls.</figcaption></figure>

actually a two-phase. In the latter case, from 81.2 % to 84.8 % of CO₂, three additional three-phase SSI iterations are carried out, and a two-phase solution is identified. One set of equilibrium ratios is built using the compositions found in the test. As a consequence, not only are initial guesses available for the two-phase split by not losing the information obtained in the three-phase test, but also two-phase stability tests are avoided. As noticed for this two-phase region, ISI has a lower number of iterations as compared to IMI, but the latter has a lower number of stability calls, as shown in Fig. 5d

### 4.1.3. Phase equilibrium at 30.2 bar and 178.8 K

In this case, the values of pressure and temperature are selected so that the phase fraction distribution for the sour gas, see Fig. 6a, can be compared with the one from Ref [14]. Qualitatively, they have analogous behavior. Otherwise, results from the previous examples, with lower pressure, and this case are alike. The exceptions are: a smaller three-phase region, the number of minima found in PSA and their quality and the number of iterations, in which the SI actually has an increase and ISI and IMI have the opposite trend. In Fig. 6b, the region with two minima with negative TPDs is reduced, as compared with the examples at lower pressures. The number of iterations for the three methods, SI, ISI and IMI, is drawn in Fig. 6c. As expected, for a large

domain, the IMI performs better than its counterparts. It is clear, however, that in the two-phase liquid-liquid interval, from 73.6–84.8 % of CO₂, the ISI has a lower number of iterations than IMI, but the number of PSA calls for the latter is lower than for ISI and SI, as shown in Fig. 6d

### 4.1.4. Phase equilibrium at 35 bar and 178.8 K

The phase distribution at 35 bar and 178.8 K for different CO₂ concentrations is presented in Fig. 7a. As compared with previous examples at lower pressures, there is a wider two-phase L1-L2, a reduced three-phase region and a lower amount of vapor phase. The area with three negative TPDs is considerably smaller than for lower pressures, as can be seen in Fig. 7b. In terms of total numbers of phase split iterations, see Fig. 7c, for the narrower three-phase region, ISI and IMI perform significantly better than SI, except near the two-phase L1-L2 region where the difference between methods is reduced. In a large part of this two-phase L1-L2 domain, that is from 72 % to 85 % of CO₂, three minima are identified, one with negative TPD, one with positive TPD and one trivial. In this case, the same initial K-values are constructed for ISI and SI, so they have the same number of iterations. At the same time, IMI has a slightly higher number of iterations (from one to three additional iterations), as a result of the three-phase test, but it avoids two-phase stability tests. In this case, the compositions of the identified present
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Sour gas mixture results at *P* = 35 bar and *T* = 178.8 K (a) phase fraction distribution; (b) number of minima found in stability test on feed composition; (c) number of phase split iterations; (d) number of stability calls.</figcaption></figure>

phases are used to form a set of K-values to initialize the two-phase split.

### 4.1.5. Phase equilibrium at 40 bar and 178.8 K

In this case, the three methods considered in this work are tested and compared in two-phase conditions. As verified previously, in the three-phase domain, IMI has the best performance when compared to the other algorithms, since two-phase stability tests are avoided and good initial guesses are obtained directly from the one-phase stability test. Here, two-phase liquid-vapor and liquid-liquid cases are analyzed. In Fig. 8a, identical phase fraction distribution resulted from the three methods, SI, ISI and IMI are drawn. Moreover, the phase distribution is consistent with the phase diagram built in Fig. 1a, showing a small one-phase zone at around 45 % of CO₂, where only liquid L1 is present.

The number of minima found after performing a one-phase stability test is presented in Fig. 8b. The total number of minima found after the stability test is three in almost the integrity of the liquid-liquid region. With three minima, IMI tests a three-phase split calculation, and since an unphysical response is noticed, the composition of the two phases identified as present are used to build a set of equilibrium ratios. The tiny L1 single-phase region around 45 % of CO₂ and the L2 region for concentrations higher than 85 % are correctly identified (only non-negative TPDs found). Only the trivial solution is found after

performing PSA for CO₂ concentrations higher than 95 %.

In terms of the number of phase split iterations, ISI generally outperforms IMI and SI, as seem in Fig. 8c. The regions with zero iterations are the ones with only one phase present, requiring only stability tests. The best performance of both ISI and IMI is in the region in which two negative TPDs are present, and no trivial solution is identified. The region where only two TPDs are found, with one negative and one zero, all the methods require the same number of split iterations. In a large part of the domain, three minima are identified. ISI uses two minima and provides a relevant reduction in the number of iterations compared to the SI method, flattening the triangular shape. The reduction in IMI is also relevant, but it requires three iterations in testing the potential presence of three phases; therefore, it has more iterations than ISI. On the other hand, IMI also flattens the iteration curve obtained in SI, and this improved method reduces the number of stability calls, as seen in Fig. 8d

### 4.1.6. Phase equilibrium at 45 bar and 178.8 K

The final test for the sour gas mixture is performed at 45 bar. At this pressure, by changing the CO₂ concentration, there is a liquid-vapor (L1-V), liquid-liquid (L1-L2) and two single phase regions, L1, at lower CO₂ concentrations, and L2, at higher concentrations, as it can be seen in
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Sour gas mixture results at *P* = 40 bar and *T* = 178.8 K (a) phase fraction distribution; (b) number of minima found in stability test on feed composition; (c) number of phase split iterations; (d) number of stability calls.</figcaption></figure>

Fig. 9a. Information on the minima found in stability test on feed composition is given in Fig. 9b, which has not a significant change from the previous lower pressure results; it shows, however, an increase in the intervals in which the number of negative TPD minima is zero, therefore wider one-phase regions. In Fig. 9c, two sets of K-values are compared, one derived from initial guesses of the SI procedure and the other coming from the equilibrium ratios obtained after convergence. Around 59 % of CO₂, there is a discontinuity in the equilibrium ratio distribution, which is the place where the deviation is the highest. From Fig. 9d, it is evident that the initial K-values from IMI are closer to the correct equilibrium ratio in blue when compared to the values estimated by the SI method.

The proximity of the initial equilibrium ratios derived from SI and IMI to the ones at the solution gives some clues on the efficiency of each method. In Fig. 9e, the number of phase split iterations for SI, ISI and IMI is presented. The number of iterations in the SI method forms a triangular shape, in which the maximum is at the same CO₂ concentration as in the discontinuity identified in Fig. 9c. The number of iterations for SI decreases as it approaches the two-phase L1-L2 boundary. In the same figure, a wide region with three TPD minima is found, and a direct three-phase split is tried. Since only two phases are present, the IMI performs a negative flash test, which explains the higher number of iterations with

respect to ISI, which has, in general, the lowest number of iterations. Even though IMI has more iterations than ISI, the former not only reduces the number of iterations with respect to SI in the triangular area, but it also has a significant reduction in the number of stability calls, which is seen in Fig. 9f. In this figure, in general the IMI has a lower number of stability calls, which leads to a relevant increase in efficiency, since for each call the number of stability iterations range in this case from 6 iterations to 25.

## 4.2. Problem 2 - Majamar reservoir - CO₂ mixture

Composition, component properties, and non-zero BIPs are listed in Table B2. In Fig. 10a, the phase envelope discriminating the number of equilibrium phases at *T* = 305.35 K is depicted. The results obtained are close to the ones presented in Refs [8,14]. In the previous problem, in a significant portion of the phase diagram, it was possible to start a three-phase split calculation by only using information from three stationary points acquired after a one-phase stability test. On the other hand, for the Majamar reservoir mixture, not enough information was obtained to start directly a three-phase split calculation, as shown in Fig. 10b, with a maximum of two minima of the TPD.

The number of iterations in the *P-Z* plane is plotted in Fig. 11a(SI),
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Sour gas mixture results at *P* = 45 bar and *T* = 178.8 K (a) phase fraction distribution; (b) number of minima found in stability test on feed composition; (c) initial and converged K-values in SI; (d) initial and converged K-values in IMI; (e) number of phase split iterations; (f) number of stability test calls.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 10. Maljamar reservoir mixture results at T = 305.35 K (a) phase boundaries; (b) number of minima with negative TPDs in one-phase stability.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Total number of multiphase split iterations for Maljamar reservoir mixture results at T = 305.35 K: (a) conventional; (b) ISI; (c) IMI.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Maljamar reservoir mixture results at *P* = 80 bar and *T* = 305.35 K: (a) phase fraction distribution; (b) number of minima found in stability test on feed composition; (c) number of phase split iterations; (d) number of stability calls.</figcaption></figure>

Fig. 11b (ISI), and Fig. 11c (IMI). Despite the lack of enough information from a one-phase stability test to start a direct three-phase split calculation, ISI and IMI methods can reduce the number of total phase split iterations by building phase equilibrium ratios that give better initial guesses using information of two stationary points obtained in the stability test. A significant reduction of the number of iterations in both two-and three-phase domains is achieved by ISI and IMI.

The results of phase equilibrium calculations at 80 bar and 305.35 K are presented in Fig. 12a–12d. Fig. 12a depicts the phase fraction distribution using the three different methods. At a lower CO₂ concentration, a vapor-liquid (V-L1) region is formed. As the amount of carbon dioxide increases, a three-phase region forms, in which the vapor phase has a drastic decrease, while the L2 phase has the opposite behavior and the L1 phase lowers its phase fraction. In this case, the number of minima found in the stability test on the feed composition is two, as presented in Fig. 12b. This means that not enough information is obtained to generate initial guesses for a direct three-phase split calculation in IMI. Therefore, ISI and IMI are equivalent in terms of efficiency, and are more efficient than the SI method, as can be seen from the number of total phase split iterations given in Fig. 12c; the two peaks in the number of iterations for SI are avoided by ISI and IMI. At the same time, the number of stability calls is the same for the three methods, as shown in Fig. 12d. In this

particular case, the number of minima identified in the one-phase stability test is, in general, lower than the previous problem (sour gas mixture). It explains why the reduction of the number of iterations in phase split calculations is not as significant as in the previous problem.

### 4.3. Problem 3 - Maljamar separator - CO₂ mixture

Another mixture considered for study is the Maljamar separator oil mixed with different amounts of CO₂ at *T* = 305.35 K. Composition, component properties, and non-zero BIPs are listed in Table B3. As compared to Refs [9,14], the same phase envelopes were obtained with all the methods studied in this paper, SI, ISI, and IMI.

The number of phases as a function of CO₂ concentration and pressure is shown in Fig. 13a. Fig. 13b presents the number of minima found in the one-phase stability test. Two sets of K-values are built, and a direct initialization of a three-phase split calculation is performed. In these regions, IMI will not only have a good initial guess for the three-phase split, but it will also avoid a number of stability tests. In the one-phase region, the number of minima with negative TPDs is zero, as expected, as depicted in Fig. 13c. In the same figure, there is a wide region with two negative TPDs and a smaller region with only one negative TPD. The number of minima with trivial solution is provided in Fig. 13d; there is a
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Maljamar separator mixture results at $T = 305.35$ K: (a) phase boundaries; (b) total number of minima in one-phase stability; (c) number of minima with negative TPDs in one-phase stability; (d) number of minima with trivial solution in one-phase stability.</figcaption></figure>

region with one minimum (in dark blue) and a wide region (in light blue) in which no trivial solution is found.

By comparing the total number of phase split iterations required by the three methods in Fig. 14, it is clear that in the three-phase region, the IMI method is the most efficient, and ISI is more efficient than SI. In the two-phase parts of the phase diagram, both ISI and IMI clearly outperform SI.

### 4.3.1. Phase equilibrium at 69.4 bar and 305.35 K

Fig. 15a presents the phase fraction distribution with the three methods (SI, ISI, and IMI) for the Maljamar reservoir mixture at 69.4 bar and 305.35 K. Our results match those presented in Ref [14]. The isobar passes through a liquid-liquid region, then to a three-phase region, and after that, the liquid-vapor phases are identified. In Fig. 15b, the number of minima found in the one-phase stability test is, in the majority of the cases, three, even though in most of the domain there are only two phases. That means that enough information is obtained to generate initial guesses for K-values to perform the three-phase split, but actually, only two phases appear. Fig. 15c presents the number of total phase split iterations; in the three-phase area, IMI has the lowest number of iterations, followed by ISI and SI. On the other hand, in the two-phase L1-L2

domain, where three negative TPD minima are found, the ISI method requires fewer phase split iterations as compared to SI and IMI. The latter method, however, requires less stability calls as compared to the ISI method, as can be seen in Fig. 15d

### 4.3.2. Phase equilibrium at 72.0 bar and 305.35 K

For this case, a wide two-phase liquid-liquid domain (for 68.2–97.6 % of CO₂) is identified, as well as a narrow three-phase region at higher carbon dioxide concentration (for 97.65–99.35 % of CO₂), as can be seen from the phase fraction distribution depicted in Fig. 16a; the results are the same for the three methods. In terms of the number of minima identified in one-phase stability test, a wider region is covered with two different TPD minima and a smaller region with three TPD minima, as presented in Fig. 16bFig. 16c gives the number of phase split iterations using the three methods. The region with three TPD minima does not cover a three-phase part, but instead a two-phase one; this is the reason why IMI has a higher number of iterations than ISI in that region, since the former tests whether or not to proceed with the three-phase split calculation. However, IMI is more computationally efficient than ISI, since it reduces the number of stability calls with respect to ISI, as can be seen from Fig. 16d. Both ISI and IMI methods decrease considerably the
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Total number of multiphase split iterations for Maljamar separator mixture results at *T* = 305.35 K for the following methods: (a) conventional; (b) ISI; (c) IMI.</figcaption></figure>

number of iterations as compared to the SI method, which presents a triangular shape with a peak in number of iterations; this format then flattened out when using ISI and IMI.

### 4.4. Problem 4 - NWE oil + impure CO₂ mixture

The fourth mixture is a North Ward Estes oil with impure CO₂ (95 % CO₂ and 5 % C₁) at *T* = 301.48 K. Composition, component properties, and non-zero BIPs are listed in Table B4. The phase boundaries shown in Fig. 17a match those obtained in Refs [9,14]. Fig. 17b depicts the total number of TPD local minima found in one-phase stability test; it is clear that a significant portion of the P-Z plane is filled with two different TPDs. In this case, no direct three-phase split calculation is started, but K-values can be built with the compositions related to the two TPDs to improve the initial guess for the two-phase split calculation. Part of the single-phase region has two TPDs, but evidently none of them are negative, as identified in Fig. 17c. After verifying Fig. 17d, the two TPDs found in the one-phase region are actually one positive (the conditions correspond to a kind of two-phase shadow region [20]) and one trivial solution. Two negative TPDs are encountered in a large portion of the plane, including two-and three-phase regions; this means that enough information is available to improve the initial guess for the two-phase split. The number of multiphase split iterations for the three methods

in the P-Z plane is given in Fig. 18a(SI), Fig. 18b (ISI), and Fig. 18c (IMI). Both ISI and IMI achieve a spectacular reduction of the number of iterations with respect to the conventional method.

#### 4.4.1. Phase equilibrium at 84 bar and 301.48 K

The phase fraction distribution on the isobar of *P* = 84 bar is drawn in Fig. 19a and matches the results from Ref [14]. A narrow three-phase region is identified in the CO₂ concentration range of 67.3 % to 77.5 %, in which two liquid and one vapor phase are present. From lower to higher CO₂ concentration, one phase, L1, is present in the mixture, and its percentage in the mixture starts to be lowered when a vapor phase, V, starts to form. The vapor phase has a sharp decrease in its concentration when a third phase, L2, appears, showing a retrograde behavior. As the amount of carbon dioxide increases, the vapor phase disappears while the L2 phase continues leveling up and eventually the L1 phase is no longer present and the mixture turns into a single L2 phase. In terms of number of TPD minima shown in Fig. 19b, in the entire domain, there are only two minima; thus, no direct initialization of a three-phase split is possible, since not enough information is obtained. In the part of the domain where the trivial solution is identified as a minimum, ISI and IMI require the same number of iterations, smaller than for the conventional method, as depicted in Fig. 19c. Similar to what has been noticed in other cases, there is a significant increase in number of iterations in the
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Maljamar separator mixture results at *P* = 69.40 bar and *T* = 305.35 K: (a) phase fraction distribution; (b) number of minima found after performing stability test on feed composition; (c) number of phase split iterations; (d) number of stability calls.</figcaption></figure>

SI method, creating a triangular shape; the number of iterations is significantly reduced by ISI and IMI.

### 4.4.2. Phase equilibrium at 100 bar and 301.48 K

In Fig. 20a, the phase fraction distribution at *P* = 100 bar is drawn for the three methods. A two-phase liquid-liquid region is identified. The total number of TPD minima and their nature are presented in Fig. 20bA maximum of two minima is obtained from the stability test for the entire concentration domain. In part of it, the trivial solution is identified as a local minimum, and the K-values are generated in a similar way as in the conventional method. On the other hand, a wider region has no trivial solution, and the initialization of the two-phase split in ISI and IMI is different from the conventional method. The initial guesses for the K-values in the SI method and the final K-values after convergence are drawn in Fig. 20c. The same comparison is made for the ISI and IMI methods (which has the same results in this case), as illustrated in Fig. 20d ISI and IMI give equilibrium ratios that are closer to the solution than the SI method. At around 66.5 % of CO₂, there is a discontinuity in the equilibrium ratio distribution for the conventional method, as can be seen in Fig. 20c. At this point, the initial guess is poor, with a relatively high deviation from the K-values from the solution, which increases the number of iterations, as presented in Fig. 20e. Since ISI and IMI produce

better initial guesses, closer to the converged K-values, their number of iterations are lower than in the conventional method and the peak in Fig. 20e is flattened by ISI and IMI. When the trivial solution is a minimum in stability test, the three methods have the same initial procedure and the same initial guesses for the K-values, as seen in Fig. 20c and in Fig. 20d, thus the same number of iterations. In this case, since the number of minima found in stability test is lower than three, no direct initialization is possible and the same number of stability calls is seen in Fig. 20f for the three methods.

The two proposed procedures, ISI and IMI, take full advantage of the information from stability tests to build K-values and initialize either two- or three-phase split calculation. There are several important differences in IMI with respect to the multiple identification method of Huang et al [14]. Firstly, only compositions of negative and zero (trivial) TPDs are used in Ref [14]. The novelty is that we also use compositions of minima with a strict positive value of the TPD function. Even though this may seem counterintuitive, since a positive TPD indicates stability of a phase, there is a continuity in compositions beyond the domain with only one or two negative minima; using also the positive TPDs leads to improved estimates of K-values, as shown graphically in two examples, for two-and three-phase cases.

Secondly, unlike in Ref [14], the trivial solution is used in forming
<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Maljamar separator mixture at *P* = 72 bar and *T* = 305.35 K: (a) phase fraction distribution; (b) number of minima found in stability test on feed composition; (c) number of phase split iterations; (d) number of stability calls.</figcaption></figure>

K-values only if it corresponds to a minimum. Thirdly, the compositions of initial guesses are always satisfying the component material balance constraints by solving the Rachford-Rice equations. In Ref [14], a least square estimation is used to initialize phase mole fractions, using a merit function based on the residual of the component material balance constraints. If the minimum of the objective function is zero, this is equivalent to solving the Rachford-Rice equations; otherwise, if this function is positive, compositions form an infeasible set of K-values for a three-phase flash. Failures are reported in Ref [14], due to an attempt to solve a three-phase Rachford-Rice with infeasible starting variables. Additional iterations in MI are simply spared by not launching a three-phase flash in these cases (we do not attempt to generate initial guesses for a three-phase split if not enough information is available to build two sets of K-values).

ISI is also an attractive method, since better initial guesses make this method more efficient than SI. For instance, in a two-phase split, SI will use the compositions of the global minimum of the TPD function and the feed composition, while ISI will use the compositions of the minima with the lowest TPD, leading to the correct type of two-phase equilibrium. For example, if the real two-phase state is LL, the flash may converge to a VL state in SI and to an LL state in ISI. In this case, instability is detected in the two-phase stability test by SI (triggering an unnecessary three-phase

flash), and a stable LL state is ascertained by stability testing in ISI.

The idea of using the second-lowest minimum of the TPD function for the construction of an alternative initial guess for a multiphase-flash calculation has appeared already in Pan et al [12] and Jex et al [42]. In Ref [42], if the three-phase flash from the lowest TPD-value in two-phase stability fails, another two-phase flash is performed using the second lowest TPD minimum from two-phase stability. Note that the initial K-values for the two-phase flashes are generated using the composition of a negative TPD (global minimum first and second lowest minimum if the first attempt fails) and the feed composition. The point in our methods is to use all available information provided by the stability test, that is, all compositions of the minimum TPDs to generate initial K-values; thus, if two minima with negative TPDs are found in two-phase stability testing, their compositions and the composition of the tested phase (from the two-phase flash) are used to generate initial K-values for the three-phase flash in both ISI and IMI methods (see flowcharts in figures A2 and A3). If three minima are found in two-phase stability testing, all three compositions are used to generate the two sets of K-values. Our methods are avoiding the additional two-phase flash in Ref [42].

The four test examples used in this work are benchmark ones, used by several authors before to test three-phase equilibrium methods.
<figure><img src="image_18.png" /><figcaption>**Fig. 17.** NWE mixture with impure CO₂ (95 % CO₂+5 % C1) results at T = 301.48 K: (a) phase boundaries; (b) total number of minima found in one-phase stability; (c) number of minima with negative TPDs found in one-phase stability; (d) number of minima with trivial solution found in one-phase stability.</figcaption></figure>

However, in a general sense, the converged solution can be a false solution, not reaching the correct set of three phases, for example if the three-phase flash is performed in a four-phase region, or if a minimum of the TPD function corresponds to an incipient fourth phase. This was not the case in this work, in which all test examples involve only three-phase systems.

If a fourth phase exists, the solution of the three-phase flash will be an intermediary unstable solution in any version of a sequential stability-flash approach [4,5]; by performing a stability test on one of the three equilibrium phases, the existence of a fourth phase can be detected. At this stage, if a four-phase flash is required, the conventional (SI) method can be applied; we are currently investigating the extension of the ISI and IMI procedures to more than three equilibrium phases using the same principles. Limited preliminary tests are encouraging, showing that if four minima of the TPD function are found in the four phase region, IMI will skip two- and three-phase flashes and related stability tests, converging directly to the correct four-phase solution. Further tests are required to extend the proposed methods for handling the four possible three-phase states and six possible two-phase states. However, when applying the ISI and IMI methods as presented in this work, we must be certain that at most three equilibrium phases can exist.

If a non-hydrocarbon component other than CO₂ that may trigger the formation of a fourth phase (a third liquid phase rich in that component) is present in the mixture, a careful screening must be performed for the expected composition, pressure and temperature ranges before using the methods within a compositional reservoir simulator, to check if a fourth phase can or cannot form.

Moreover, the nonlinearity of the Gibbs free energy depends at first on the thermodynamic model used, not the computational algorithm. Therefore, the algorithm's robustness cannot be generalized and require further research for its extensions to systems with more than three-phases or using different thermodynamic models. The case of hydrocarbon-water mixtures is special, because of the very low solubility of most hydrocarbon components in the aqueous phase. In some cases, the solubility of all or most of the hydrocarbon components in the aqueous phase can be neglected. Whether a simplified model is used or not, the hydrocarbon-water equilibrium has a special mathematical structure that can be exploited. Even though the proposed methods can be used for hydrocarbon-water mixtures, we found that some modifications taking into account the particularity of such mixtures would improve the efficiency of these methods and that a customized version of the proposed algorithms must be designed. This remarks hold for any
<figure><img src="image_19.png" /><figcaption>Fig. 18. Total number of multiphase split iterations for NWE mixture with impure CO₂ (95 % CO₂+5 % C1) at T = 301.48 K: (a) SI; (b) ISI; (c) IMI.</figcaption></figure>

component that can form an almost pure third phase.

In a simulation scenario where most regions are in two phases, but three phases can exist at some conditions, testing everywhere the stability from all initial guesses would not be an extra cost instead of a benefit, as it may seem intuitively. The number of phases is not known a priori and testing stability using all the initial guesses from the given set (see Section 3.4) must be performed (this is the standard practice for multiphase flashes [5]). Whatever the final state of the mixture, even if a negative TPD is detected, already asserting instability, we don't know if it is a local or the global minimum and the search must continue to increase the probability that the global minimum is found. In the simulation case mentioned above, for most conditions a maximum of two minima will be found. If there are three minima, IMI will solve directly a three-phase flash when this is possible and it can also handle the situation when three minima exist in the two-phase region (indicating an incipient third phase) and converge to the correct two-phase solution. If however we certainly know in advance that only a VL behavior is expected (as in standard petroleum reservoir engineering applications), using only Michelsen's two-sided initialization [4] for stability suffices. In single-phase stability testing, an important number of stability calls can be further avoided using the stability analysis bypassing method [43], designed to be used in compositional flow simulation; its

implementation is out of the scope of this work.

It must be noted that even when a three-phase flash algorithm fails to converge to the correct solution or completely diverges in stand-alone calculations, there are numerous ways to restart without following the same iteration route for failure. However, the idea in this work is to reduce to a minimum the situations in which calculations need to be restarted when a false solution is obtained or divergence occurs. No failures were detected in the test examples using the proposed methods. In these test cases, P-Z domains covers a variety of situations, some conditions being very difficult ones, such as the vicinity of bi-critical points [8].

## 5. Conclusions

In this work, two new initialization strategies for multiphase flash calculations are proposed. In the first one (improved stepwise initialization), the conventional stepwise identification procedure is improved by using additional initial guesses. In the second one (improved multiple initialization), a three-phase split is initiated if the stability testing of the feed composition detects at least three minima of the tangent plane distance function. Both proposed methods use all information (number of minima and their composition) from phase stability testing at each
<figure><img src="image_20.png" /><figcaption>**Fig. 19.** NWE mixture with impure CO₂ (95 % CO₂+ 5 % C1) results at *P* = 84 bar and *T* = 301.48 K: (a) phase fraction distribution; (b) number of minima found in stability test on feed composition; (c) number of phase split iterations; (d) number of stability calls.</figcaption></figure>

phase split 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. In both phase stability and flash calculations, highly robust calculation procedures are used, based on successive substitution iterations (SSI) in early iteration stages, followed by Newton iterations with modified Cholesky factorization and line search. The proposed methods are tested and compared with the conventional procedure for several benchmark mixtures from the literature, containing hydrocarbon components and mixed in different proportions with CO₂. Phase diagrams are constructed in the pressure-composition plane, focusing on the number of stationary points of the TPD functions found in stability testing at each step of the multiphase stability-flash algorithm and on how they must be efficiently used in initialization of flash calculations. 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 which both efficiency and robustness are of utmost importance. The proposed methods were designed for, but are not limited to three-phase systems; both ISI and IMI procedures can be extended to more equilibrium phases using the same principles; the extension to the four-phase equilibrium case is currently being investigated. The proposed algorithms were tested for stand-alone calculations for CO₂-hydrocarbon mixtures with up to three phases. The results recommend the proposed methods as strong candidates to be included in compositional reservoir simulators. This would represent a more severe test of the robustness, since in numerical reservoir simulations calculations are performed a significant number of times for each time step in every grid block of the discretized reservoir model, and a larger variety of PTZ conditions may challenge the phase equilibrium routine. The code is currently being implemented in a compositional reservoir simulator.
<figure><img src="image_21.png" /><figcaption>**Fig. 20.** NWE mixture with impure CO₂ (95 % CO₂+5 % C1) at *P* = 100 bar and *T* = 301.48 K: (a) phase fraction distribution; (b) number of minima found in stability test on feed composition; (c) initial and converged K-values in SI; (d) initial and converged K-values in IMI; (e) number of phase split iterations; (f) total number of stability test calls.</figcaption></figure>
CRediT authorship contribution statement

**Juan Heringer:** Writing – original draft, Software, Methodology, Investigation, Formal analysis, Conceptualization. **Michiel Wapperom:** Software, Investigation, Conceptualization. **Catinca Secuianu:** Writing – review & editing, Validation, Methodology, Investigation, Conceptualization. **Denis Voskov:** Writing – review & editing, Validation, Supervision, Methodology, Investigation, Formal analysis, Conceptualization. **Dan Vladimir Nichita:** Writing – review & editing, Writing – original draft, Validation, Supervision, Resources, Project administration, Methodology, Investigation, Formal analysis, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

JH and MW thank Total Energies for financial support. CS and DVN thank the ISIFOR Carnot Institute for a mobility grant (MODCO project). This work was first presented at the 2024 AIChE Annual Meeting, San Diego, CA, October 26–31, 2024.

Appendix A. Flowcharts of the initialization procedures

This appendix presents the flowcharts describing the steps of the conventional (SI), ISI and IMI methods.

<figure><img src="image_22.png" /><figcaption>Fig. A1. Flowchart of the conventional (SI) procedure.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. A2. Flowchart of the ISI method.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Fig. A3. Flowchart of the IMI method.</figcaption></figure>

## Appendix B Composition, component properties and non-zero BIPs

<table><caption>**Table B1**

Properties for Sour gas mixture.</caption><thead><tr><th>Component</th><th>T<sub>c</sub> (K)</th><th>P<sub>c</sub> (bar)</th><th>ω</th><th>Mol%</th><th>K<sub>I,CO2</sub></th><th>K<sub>I,N2</sub></th><th>K<sub>I,H2S</sub></th></tr></thead><tbody><tr><td>CO<sub>2</sub></td><td>304.20</td><td>73.90</td><td>0.225</td><td>70.59</td><td>-</td><td>-</td><td>-</td></tr><tr><td>N<sub>2</sub></td><td>126.20</td><td>33.50</td><td>0.040</td><td>7.03</td><td>-0.020</td><td>-</td><td>-</td></tr><tr><td>H<sub>2</sub>S</td><td>373.20</td><td>89.40</td><td>0.081</td><td>1.97</td><td>0.120</td><td>0.200</td><td>-</td></tr><tr><td>C<sub>1</sub></td><td>190.60</td><td>45.40</td><td>0.008</td><td>6.86</td><td>0.125</td><td>0.031</td><td>0.100</td></tr><tr><td>C<sub>2</sub></td><td>305.40</td><td>48.20</td><td>0.098</td><td>10.56</td><td>0.135</td><td>0.042</td><td>0.080</td></tr><tr><td>C<sub>3</sub></td><td>369.80</td><td>41.90</td><td>0.152</td><td>2.97</td><td>0.150</td><td>0.091</td><td>0.080</td></tr></tbody></table>
<table><caption>**Table B2**
Properties for the Maljamar Reservoir Oil.</caption><thead><tr><th>Component</th><th>T<sub>c</sub> (K)</th><th>P<sub>c</sub> (bar)</th><th>ω</th><th>Mol%</th><th>K<sub>I,CO2</sub></th><th>K<sub>I,C1</sub></th></tr></thead><tbody><tr><td>CO<sub>2</sub></td><td>304.20</td><td>73.80</td><td>0.225</td><td>-</td><td>-</td><td>-</td></tr><tr><td>C<sub>1</sub></td><td>190.60</td><td>45.40</td><td>0.008</td><td>29.32</td><td>0.115</td><td>-</td></tr><tr><td>C<sub>2</sub></td><td>305.40</td><td>48.20</td><td>0.098</td><td>10.19</td><td>0.115</td><td>-</td></tr><tr><td>C<sub>3</sub></td><td>369.8</td><td>41.90</td><td>0.152</td><td>8.35</td><td>0.115</td><td>-</td></tr><tr><td>nC<sub>4</sub></td><td>425.2</td><td>37.5</td><td>0.193</td><td>3.31</td><td>0.115</td><td>-</td></tr><tr><td>C<sub>5-7</sub></td><td>516.7</td><td>28.8</td><td>0.265</td><td>12.04</td><td>0.115</td><td>0.045</td></tr><tr><td>C<sub>8-10</sub></td><td>590.0</td><td>23.7</td><td>0.364</td><td>15.81</td><td>0.115</td><td>0.055</td></tr><tr><td>C<sub>11-14</sub></td><td>668.6</td><td>18.6</td><td>0.499</td><td>8.23</td><td>0.115</td><td>0.055</td></tr><tr><td>C<sub>15-20</sub></td><td>745.8</td><td>14.8</td><td>0.661</td><td>5.28</td><td>0.115</td><td>0.060</td></tr><tr><td>C<sub>21-28</sub></td><td>812.7</td><td>12.0</td><td>0.877</td><td>2.76</td><td>0.115</td><td>0.080</td></tr><tr><td>C<sub>29+</sub></td><td>914.9</td><td>8.5</td><td>1.279</td><td>4.64</td><td>0.115</td><td>0.280</td></tr></tbody></table>

<table><caption>**Table B3**
Component data for the Maljamar Separator Oil.</caption><thead><tr><th>Component</th><th>T<sub>c</sub> (K)</th><th>P<sub>c</sub> (bar)</th><th>ω</th><th>Mol%</th><th>K<sub>I,CO2</sub></th></tr></thead><tbody><tr><td>CO<sub>2</sub></td><td>304.20</td><td>73.90</td><td>0.225</td><td>-</td><td>-</td></tr><tr><td>C<sub>5-7</sub></td><td>516.7</td><td>28.8</td><td>0.265</td><td>23.54</td><td>0.115</td></tr><tr><td>C<sub>8-10</sub></td><td>590.0</td><td>23.7</td><td>0.364</td><td>32.95</td><td>0.115</td></tr><tr><td>C<sub>11-14</sub></td><td>668.6</td><td>18.6</td><td>0.499</td><td>17.13</td><td>0.115</td></tr><tr><td>C<sub>15-20</sub></td><td>745.8</td><td>14.8</td><td>0.661</td><td>10.99</td><td>0.115</td></tr><tr><td>C<sub>21-28</sub></td><td>812.7</td><td>12.0</td><td>0.877</td><td>5.74</td><td>0.115</td></tr><tr><td>C<sub>29+</sub></td><td>914.9</td><td>8.5</td><td>1.279</td><td>9.65</td><td>0.115</td></tr></tbody></table>

<table><caption>**Table B4**
Component data for the NWE with impure CO2.</caption><thead><tr><th>Component</th><th>T<sub>c</sub> (K)</th><th>P<sub>c</sub> (bar)</th><th>ω</th><th>Mol%(Initial)</th><th>Mol%(Secondary)</th><th>K<sub>I,CO2</sub></th></tr></thead><tbody><tr><td>CO<sub>2</sub></td><td>304.20</td><td>73.90</td><td>0.225</td><td>0.77</td><td>95</td><td>-</td></tr><tr><td>C1</td><td>190.6</td><td>28.8</td><td>0.265</td><td>20.25</td><td>5</td><td>0.120</td></tr><tr><td>PC1</td><td>343.64</td><td>23.7</td><td>0.364</td><td>11.8</td><td>-</td><td>0.120</td></tr><tr><td>PC2</td><td>466.41</td><td>18.6</td><td>0.499</td><td>14.84</td><td>-</td><td>0.120</td></tr><tr><td>PC3</td><td>603.07</td><td>14.8</td><td>0.661</td><td>28.63</td><td>-</td><td>0.090</td></tr><tr><td>PC4</td><td>733.79</td><td>12.0</td><td>0.877</td><td>14.9</td><td>-</td><td>0.090</td></tr><tr><td>PC5</td><td>923.2</td><td>8.5</td><td>1.279</td><td>8.81</td><td>-</td><td>0.090</td></tr></tbody></table>

## Data availability

No data was used for the research described in the article.

## References

[1] F.M. Orr, A.D. Yu, C.L. Lien, Phase behavior of CO₂ and crude oil in low-temperature reservoirs, SPE J. 21 (1981) 480-492.

[2] J.W. Gibbs, A method of geometrical representation of the thermodynamic properties of substances by mean of surfaces, Trans. Conn. Acad. 2 (1873) 382-404.

[3] L.E. Baker, A.C. Pierce, K.D. Luks, Gibbs energy analysis of phase equilibria, Soc. Pet. Eng. J. 22 (1982) 731-742.

[4] M.L. Michelsen, The isothermal flash problem. Part I. Stability, Fluid Phase Equilib. 9 (1982) 1-20.

[5] M.L. Michelsen, The isothermal problem. Part II. Phase split calculation, Fluid Phase Equilib. 9 (1982) 21-40.

[6] A. Firoozabadi, Thermodynamics of Hydrocarbon Reservoirs, McGraw Hill, 1999.

[7] H.A. Li, Multiphase Equilibria of Complex Reservoir Fluids, Springer, 2022.

[8] Z. Li, A. Firoozabadi, General strategy for stability testing and phase-split calculations in two and three phases, SPE J. 17 (2012) 1096-1107.

[9] M. Petifrere, D.V. Nichita, Robust and efficient Trust-region based stability analysis and multiphase flash calculations, Fluid Phase Equilib. 362 (2014) 51-68.

[10] H. Pan, Y. Chen, J. Sheffield, Y.B. Chang, D. Zhou, Phase-behavior modeling and flow simulation for low-temperature CO₂ injection, SPE Reservoir Eval. Eng. 18 (2015) 250-263.

[11] W. Pang, H. Li, An augmented free-water three-phase Rachford-Rice algorithm for CO₂/hydrocarbons/water mixtures, Fluid Phase Equilib. 450 (2017) 86-98.

[12] H. Pan, M. Connolly, H. Tchelepi, Multiphase equilibrium calculation framework for compositional simulation of CO₂ injection in low-temperature reservoirs, Ind Eng Chem Res 58 (2019) 2052-2070.

[13] M. Imai, H. Pan, M. Connolly, H. Tchelepi, M. Kurihara, Reduced variables method for four-phase equilibrium calculations of hydrocarbon-water-CO₂ mixtures at a low temperature, Fluid Phase Equilib. 497 (2019) 151-163.

[14] G. Huang, B. Yuan, W. Zhang, X. Zhu, B. Dindoruk, B. Wang, A novel method to accelerate the three-phase flash calculation of the hydrocarbon-CO₂ system, Fluid Phase Equilib. 585 (2024) 114127.

[15] S. Xu, G. Wang, B. Gao, J. Tian, A line-search-based algorithm for multiphase flash calculations with CO₂-hydrocarbon system, ACS Omega 9 (2024) 34370-34379.

[16] L. Li, H. Jing, J. Liu, H. Pan, Z. Fang, T. Kuang, Y. Lan, J. Guo, The artificial neural network-based two-phase equilibrium calculation framework for fast compositional reservoir simulation of CO₂ EOR, Fluid Phase Equilib. 585 (2024) 114151.

[17] R. Sun, H. Pan, H. Tchelepi, New algorithm of three-phase equilibrium calculations for CO₂-hydrocarbon-water systems, Geoenergy Sci. Eng. 244 (2025) 213426.

[18] S. Li, L. Xu, H. Li, A modified trust-region-based algorithm for robust four-phase equilibrium calculations, in press soc, Petrol. Eng. J. (2025), https://doi.org/10.2118/228395-PA.

[19] C.H. Whitson, M.L. Michelsen, The negative flash, Fluid Phase Equilib. 53 (1990) 51-72.

[20] D.V. Nichita, D. Broseta, F. Montel, Calculation of convergence pressure/temperature and stability test limit loci of mixtures with cubic equations of state, Fluid Phase Equilib. 261 (2007) 176-184.

[21] P.E. Gill, W. Murray, Newton type methods for unconstrained and linearly constrained optimization, Math. Program. 7 (1974) 311-350.

[22] R.B. Schnabel, E. Eskow, A new modified cholesky factorization, SIAM J. Sci. Stat. Comput. 11 (1990) 1136-1158.

[23] J. Nocedal, S.J. Wright, Numerical Optimization, 2nd edition, Springer, 2006.
[24] D.V. Nichita, M. Petitfrere, Phase equilibrium calculations with quasi-Newton methods, Fluid Phase Equilib. 406 (2015) 194–208.

[25] M. Petitfrere, D.V. Nichita, Multiphase equilibrium calculations using a reduction method, Fluid Phase Equilib. 401 (2015) 110–126.

[26] M.L. Michelsen, Calculation of multiphase equilibrium, Comput. Chem. Eng. 18 (1994) 545–550.

[27] C.F. Leibovici, D.V. Nichita, A new look at multiphase rachford-rice equations for negative flashes, Fluid Phase Equilib. 267 (2008) 127–132.

[28] C.F. Leibovici, J. Neoschil, A solution of Rachford-Rice equations for multiphase systems, Fluid Phase Equilib. 112 (1995) 217–221.

[29] R. Okuno, R.T. Johns, K. Seehrnoori, A new algorithm for Rachford-Rice for multiphase compositional simulation, SPE J. 15 (2010) 313–325.

[30] D.V. Nichita, C.F. Leibovici, A rapid and robust method for solving the Rachford-Rice equation using convex transformations, Fluid Phase Equilib. 353 (2013) 38–49.

[31] D.V. Nichita, C.F. Leibovici, Improved solution windows for the resolution of the Rachford-Rice equation, Fluid Phase Equilib. 452 (2017) 69–73.

[32] D.V. Nichita, Fast and robust phase stability testing at isothermal-isochoric conditions, Fluid Phase Equilib. 447 (2017) 107–124.

[33] D.V. Nichita, Volume-based phase stability testing at pressure and temperature specifications, Fluid Phase Equilib. 458 (2018) 123–141.

[34] D.V. Nichita, A volume-based approach to phase equilibrium calculations at pressure and temperature specifications, Fluid Phase Equilib. 461 (2018) 70–83.

[35] R.B. Schnabel, E. Eskow, Software For a New Modified Cholesky Factorization; CUCS-443-89, Colorado U. at Boulder, Dept. of Computer Science, 1989.

[36] E. Eskow, R.B. Schnabel, Algorithm 695: software for a new modified Cholesky factorization, ACM Trans. Math. Software 17 (1991) 306–312.

[37] J.E., Jr. R.B. Dennis, Schnabel, Numerical Methods For Unconstrained Optimization and Nonlinear Equations, Classics in Applied Mathematics, Society for Industrial and Applied Mathematics, Philadelphia, 1996.

[38] C.H. Whitson, M.R. Brule, Phase Behavior, Society of Petroleum Engineers, Richardson, Texas, 2000.

[39] G. Wilson, in: A Modified Redlich-Kwong Equation of State, Application to General Physical Data Calculations, paper no. 15C presented at the AIChE 65th National Meeting, Cleveland, Ohio, May 4-7, 1969.

[40] D.Y. Peng, D.B. Robinson, A new two-constant equation of State, Ind. Eng. Chem. Fundam. 15 (1976) 59–64.

[41] D.B. Robinson, D.Y. Peng, 1978. The characterization of the heptanes and heavier fractions for the GPA Peng-Robinson programs. Gas Processors Association, Research Report RR-28.

[42] M. Jex, J. Mikyska, A. Firoozabadi, Robust multiphase-split calculations based on improved successive substitution schemes, SPE J. 29 (2024) 3702–3717.

[43] H. Pan, H.A. Tchelepi, Compositional flow simulation using reduced-variables and stability-analysis bypassing, in: SPE Reservoir Simulation Conference, The Woodlands, Texas, USA, 2011, https://doi.org/10.2118/142189-MS. February.