Full-EoS based thermal multiphase compositional simulation of CO2 and steam injection processes

Abstract Several compositional reservoir simulators based on equations of state (EoS) have been designed. Yet, few of them can deal with thermal processes such as steam injection and no fully compositional thermal simulation of the steam injection process has been proposed with extra-heavy oils yet. Those simulations appear essential and will offer better tools to decide whether to carry out the exploitation of a heavy oil field or not. In those processes, the accurate modelling of the water/steam phase plays an important role and accurate multiphase equilibrium calculations are necessary. Thermodynamics generate highly nonlinear problems. Besides, in reservoir simulation a huge amount of phase equilibrium calculations is required, and a single failure may cause significant error propagations leading to false solutions. This study presents several improvements leading to a more robust and efficient phase equilibrium calculation (stability and flash) program. A general multiphase flash implementing all the developed algorithms is presented and tested against experimental and literature data. It can handle any number of phases; a numerical example consists in a four-phase simulation of CO2 injection. Simulations of steam flooding are performed with highly heterogeneous reservoirs. Besides, a fully compositional simulation of the SAGD process for an extra heavy bitumen is performed, which appears to be the one of the first simulation of the kind in the literature.


Introduction
The first reservoir simulation techniques dealing with thermal processes were designed in the late sixties (Spillette and Nielsen, 1968;Baker, 1969;Shutler, 1969Shutler, , 1970. These models were based on a three-phase linear model (Shutler, 1969), then extended to two dimensions (Shutler, 1970). Later on, Vinsome (1974) introduced an IMPES method to simulate the steam drive process. Coats (1976) proposed the first simulator to couple thermal and compositional modelling. Coats (1978) presented an extension of Coats (1976) on the bare of a three-dimensional, compositional model to simulate steam injection. Later on, Ishimoto (1985) proposed to compute the oil and hydrocarbon properties by means of an equation of state. They used Raoult's law to compute the solubility of water in the oil. Then, Chien et al. (1989) proposed a general compositional simulator which could deal with thermal options. The simulator could handle both K-value and EoS-based formulation. Once again, the water phase was assumed to be an ideal solution.
More recently, Cicek and Ertekin (1996), and later Cicek (2005) developed a fully implicit compositional multiphase simulator to model steam injection problems. In their formulation, water is treated within the equilibrium calculations. The simulation technique is based on K-value models which can be tuned to match phase-diagrams. However, the tables are generally functions of one oil component, the pressure and the temperature. The K-values approach generally does not give correct solubility of hydrocarbon components in the water phase and the solubility of water in the hydrocarbon liquid phase.
It has been observed (Fig. 1, McKetta and Katz, 1948) however, that for high temperatures, the solubility of the water component in the hydrocarbon rich phase was not negligible. The current commercial simulation software makes some approximations to model steam injection problems. Those approximations are even more important in case of heavy oils. Most of the modern thermal-compositional simulation formulations are working with K-Values, or treat water separately from the equilibrium flash calculation which can lead to significant inaccuracies. For this reason, some authors started to simulate steam injection with EoS-based full compositional reservoir simulators. Using the fact that the solubility of hydrocarbons within the water phase is negligible for some ranges of pressures-temperatures, some authors dealing with steam flood problems have developed a free water flash based on this assumption Luo and Barrufet (2005), and more recently applied to heavy oil SAGD problems (Heidari, 2014). The same methodology has been applied to deal with in-situ upgrading by Lapene (2010). Several very recent studies focused on improved multiphase equilibrium calculations for hydrocarbon-water mixtures (Connolly et al., 2019;Li, 2018, 2019).
For cold CO 2 injection problems, no assumption can be made, and some authors have extended simulator models to deal with full fourphase multiphase flash problems. Varavei and Sepehrnoori (2009) and Okuno (2009) developed four-phase EoS-based simulators. As the solubility of water within hydrocarbons is negligible for those small temperatures, they did not include water within the equilibrium calculations. Brantferger (1991) developed a simulator using an equation of state to calculate the thermodynamic properties of each phase (including the water phase). They proposed an isenthalpic flash, taking enthalpy as primary variable instead of temperature. Rasmussen et al. (2006) implemented a criterion for deciding when it is justified to bypass the stability analysis in compositional simulations, leading to significant reductions of the computation time spent on phase equilibrium calculations. Voskov et al. (2009) proposed a general-purpose reservoir simulator and developed an efficient way to calculate the phase equilibrium based on the parametrization of the compositional space using tie-lines. An adaptive tabulation scheme was used to replace the phase stability test in compositional simulators, leading to very large improvements (at least an order of magnitude) in the performance of phase behavior computations. Varavei and Sepehrnoori (2009) and more recently Zaydullin et al. (2014) developed full three phase flash simulations for steam injection problems and Feizabadi (2013) extended compositional simulations to simulate isothermal solvent injection problems with a full four-phase EoS-based compositional simulator. Zaydullin et al. (2014) also proposed a strategy to bypass the phase identification of mixtures that can form three or more phases, as well as some expensive flash calculations, by using information from a limited number of key parameterized tie-simplexes.
However, to the best of our knowledge, no case of fully EoS-based compositional simulations of steam injection (like the SAGD process) with extra-heavy oil has been proposed in the literature yet, treating the water phase within the equilibrium calculations (and without any assumption). Those simulations seem to be important and will offer improved tools for predictive studies of the heavy oil fields.
Phase equilibrium problems (phase stability test and multiphase flash calculations) have to be solved repeatedly, often a very large number of times, in conventional compositional reservoir simulation and account for a large part of the computational effort; this becomes even more intensive if more than two phases can co-exist at equilibrium, which is usually relevant to complex tertiary recovery processes, such as gas steam-solvent co-injection. The computational effort can be significant, and it is increasing with the number of components in the mixture. Moreover, a single failure in thermodynamic calculations may cause significant error propagations leading to false solutions or failures of the simulation. Thus, it is imperative that phase equilibrium calculation algorithms are efficient, highly robust and balanced. Stability testing and flash calculations require advanced techniques, such as trust-region methods to ensure a descent direction of the objective function (together with line-search procedures to ensure a decrease of the function at each iteration) to deal with a variety of potential difficult conditions encountered in a compositional simulation run. This method may be required when initial guesses are poor or the successive substitution iteration (SSI) method switches to Newton iterations too early. It is especially pronounced when conditions are close to singularities (critical point for phase splitting and stability test limit locus, STLL, for stability testing), or when iterates are crossing a region of negative curvature of the objective function which always occurs within a certain pressure range, more or less extended, above the STLL at given temperature (Nichita, 2016).
Another approach, given by the so-called reduction methods, represent an attractive technique in the attempt to reduce the computational time, by significantly reducing the dimensionality of the problem, The basic idea behind the reduction method is to express the fugacity coefficients in terms of a reduced number of variables, instead of expressing them as a function of compositions. The number of independent variables does not depend on the number of components in the mixture but only on the number of components with non-zero binary interaction parameters (BIPs) in the equation of state (EoS), and the computational time increases linearly with the number of components (in conventional methods this dependence is at least quadratic). Extensive numerical experiments (Petitfr� ere and Nichita, 2015a;2015b) have shown that reduction methods seem not to be suitable for compositional simulation, where the number of components is limited, typically to less than a dozen. On the contrary, for phase stability testing, it was shown that reduction methods are more efficient than conventional ones even for mixtures with few components, which makes them suitable for compositional reservoir simulation. The reduced stability testing method of Nichita and Petitfrere (2013) was found to be the most effective formulation. This observation is very important since, in some compositional formulations (natural variables family), no explicit phase split is required. Instead, flash calculations become a part of a larger nonlinear problem solution where local phase equilibrium constraints are coupled with global flow equations and solved only in those control volumes where stability test yields two-phase state.
In this paper, the developed multiphase equilibrium algorithm is presented in a first part. It is tested in standalone (i.e. without any reservoir simulator) for various mixtures, and the results are compared with literature and experimental data. Our multiphase equilibrium code (denoted MFlash) has been integrated into two reservoir simulators to perform compositional simulations: IHRRS platform (internal reservoir simulator from Total S.A.) and AD-GPRS from the Stanford University. Different processes are simulated. First, a two-phase simulation in presence of water treated separately from equilibrium calculations is tested and benchmarked against a commercial simulator.
CO 2 injection at low temperature can exhibit a CO 2 -reach phase stripping most of the oil. This process is highly difficult to simulate and requires robust phase equilibrium algorithms. To test our algorithm, we therefore present tests of three and four-phase fully compositional simulations of CO 2 injection problems. To the best of our knowledge, this is the first full EoS based four phase reservoir simulation performed in the literature. We then present the results of fully compositional simulations of the steam-flooding process in a heterogeneous reservoir and conclude with a fully compositional simulation of the SAGD process performed on an Athabasca bitumen.
In this last part, different methodologies to improve computational times in the reservoir simulations are compared for all the presented simulations.
� The reduced variables for the stability analysis against conventional variables � The bypass approach proposed by Rasmussen et al. (2006), which avoids stability computations (called B-R) in this work � The procedure from Voskov and Tchelepi (2007), Voskov and Tchelepi (2009a), Voskov and Tchelepi (2009b), Iranshahr et al. (2010a) and more recently Zaydullin et al. (2016) based on the tie-simplex parameterization of the space will be used together with MFlash.

Phase equilibrium calculations
In this work, the multiphase-equilibrium procedure is based on a sequential two-step process, as first suggested by Michelsen (1982a), alternating phase-split calculations and stability analysis algorithms based on the tangent plane distance (TPD) function, with multiple initial guess, as proposed by Li and Firoozabadi (2012). A general form of two-parameter cubic EoS, incorporating the Soave-Redlich-Kwong (SRK) EoS (Soave, 1972) and the Peng-Robinson (PR) EoS (Peng and Robinson, 1976) is used.

Phase stability testing
The phase stability testing procedure consists in the unconstrained minimization of the modified tangent plane distance function (Michelsen, 1982a) where Y i are (formally) mole numbers, ϕ i is the fugacity coefficient of component i and d i ¼ ln z i þ ln ϕ i ðzÞ. Different methods can be used to minimize this function. The most common methods are the successive substitution (SSI) and the Newton-Raphson methods with Y i as primary variables. The gradient of the function with respect to Y i is given by and the Hessian matrix is A Newton step can be carried out to update Y i To improve scaling, Michelsen (1982a) proposed the change of variables α i ¼ 2 ffi ffi ffi ffi ffi Y i p for the minimization of the modified TPD function, giving where g ¼ T g, H ¼ THT T , and T is a diagonal matrix, with diagonal elements T i;i ¼ ffi ffi ffi ffi ffi Y i p . If lnY i are taken as the independent variables, the Newton iterations are The SSI update is given by and it corresponds to J ¼ I in Eq. (6). First, a number of SSI iterations (which are very robust but may be very slow) is carried out, then a switch is performed to the second-order Newton method.

Phase split calculations
The objective function to be minimized is the (dimensionless) Gibbs free energy of the mixture (Michelsen, 1982b), given by where n ij are the mole numbers of the component i in the phase j; one of the phases in dependent, thus there are nc � (np-1) independent variables. Differentiation with respect to n ij gives the gradient where the index F designates the reference (dependent) phase and K ij ¼ x ij =x iF are the equilibrium constants of the component i.

The Hessian matrix is
with i; j ¼ 1; nc; k; m ¼ 1; np; k; m 6 ¼ F. Again, different methods and independent variable choices can be used to minimize this function. The most common methods are the SSI method and the Newton-Raphson with n ij or ln K ij ¼ ln x ij À ln x iF as primary variables.
The SSI update gives and the Newton iteration equation is and a switch from the former to the latter method is performed according to predetermined criteria. It must be noted that in the simulator used in this work Zaydullin et al., 2014), the flash calculation routine is called only: i) in the case of new phase appearance, to initialize Newton iterations for solving the global system of equations and ii) for grid blocks where the mixture is near-critical and the global Newton method had not converged in a predetermined number of iterations (here 10 iterations); in the latter case, flash calculations are forced in order to help the convergence for the global system.

Trust-region methods
For multiphase-split calculations or for stability analysis in difficult conditions (that is to say close to the STLL), we use the Trust-Region method with Y i as primary variables, developed in Petitfr� ere and Nichita (2014), that offers robustness and high efficiency, especially for ill-conditioned systems (Conn et al., 2000).
A Taylor development of the objective function f (which can be either D M for stability or G for flash) around x k gives Taylor's development of the gradient gives Minimizing the function f means getting the gradient g ¼ rf to be zero. Therefore, trying to put gðx k þsÞ ¼ 0 we obtain the Newton method, which is quadratic and consists in solving For solving the system by Newton iterations, the Hessian must be positive definite, otherwise the method cannot be applied. In the Trust-Region method, the Hessian is corrected to ensure positive definiteness. If B k is not positive definite, by adding a diagonal element to the Hessian, H k ¼ B k þ λI, it becomes positive definite, for a value λ > maxð0; À λ 1 þ εÞ, where λ 1 is the smallest eigenvalue of the Hessian matrix.
The Trust-Region algorithm consists in minimizing the approximation of fðx þsÞ at the second order According to Nocedal and Wright (2006), solving this subproblem is equivalent to solving the following: for λ � 0 and s feasible. The Trust Region method is optimizing the order of convergence, performing steps between the gradient descent (first-order when the system is ill-conditioned) and the Newton (second order when the Hessian is positive definite) steps which are generally super-linear. However, for stability analysis, for most of the points, the calculations can be shortened using reduction variables as developed in Nichita and Petitfrere (2013); both Trust-Region and reduced phase stability methods are evaluated in this work.

Reduction methods
Calling C ij the binary interaction coefficient between component i and j in the mixture, the spectral decomposition (Hendriks and van Bergen, 1992) where m is the rank of C, λ α are the non-zero eigenvalues of C, q' αi ; α ¼ 1; m; i ¼ 1; nc are the elements of the corresponding eigenvectors, Λ ¼ diagðλ 1 ; :::; λ m Þ and C 0 is an nc � m matrix (of elements q' αi ). The matrix C exhibits a double bordered structure and is rank deficient in most cases; the rank is m � 2c þ 1 (Petitfr� ere and Nichita, 2015a;2015b), where c is the number of components having non-zero BIPs with the remaining components in the mixture. For many mixtures of interest, the rank of C is m << nc, for instance in mixtures containing many components belonging to homologous series, such as hydrocarbon mixtures.
The reduction parameters (Hendriks, 1988), Q ¼ ðQ 1 ; Q 2 ; :::; Q M Þ T , are defined as where M ¼ m þ 1 and q αi are elements of the reduction matrix (Hendriks, 1988); in the first m reduction parameters The fugacity coefficients can be expressed as functions of the reduction parameters (Nichita and Graciaa, 2011) and the coefficients h kα are functions only of reduction parameters, h kα ¼ h kα ðQ k Þ: and δ 1;2 ¼ 1 � ffi ffi ffi 2 p for the Peng-Robinson EoS. The error equations in the reduction method (Nichita and Petitfrere, 2013) are and the elements of the Jacobian matrix are The resulting linear system of equations in the Newton method is where e ¼ ðe 1 ; :::; e Mþ1 Þ T . It was shown later that this formulation of reduced phase stability testing corresponds to a constrained minimization of the TPD function with respect to a specific set of variables and constraints and an improved formulation (taking advantage of symmetry) was given (Petitfr� ere and Nichita, 2015a; 2015b).

Flow equations
As soon as fluid is injected and/or produced, the forces are not in equilibrium and fluid advection and diffusion appears. Theoretically, the fluid is not in thermodynamic equilibrium either. However, a generally accepted assumption in the reservoir simulation framework is to consider the fluid in local thermodynamic equilibrium. The continuum model can be represented as a algebraic system of equations after applied space and time discretization (Islam et al., 2010).
The equation solved in reservoir simulations is the mass conservation equation for each component i and for non-reactive flows where ρ j is the molar density of the phase j (mol/m 3 ) (all the units will be given in S.I.), x ij is the mole fraction of component i in phase j, S j is the saturation of phase j, φ is the porosity and q is the source term (mol.m 3 . m.s À 1 ). At macroscopic (Darcy) scale, for a 3D flow system under the gravitational force, the generalized Darcy's law can be applied: where p j is the phase pressure (Pa), p ij ¼ p i À p j is the capillary pressure between phases i and j, k rj is the relative permeability of the phase j, ρ m;j is the mass density of the phase j (kg m À 3 ). g is gravity acceleration (m s À 2 ), K is the permeability tensor (m 2 ) and μ j is the viscosity of the phase j (Pa s).
The energy balance may be given in the form: where P np j¼1 h j ρ j q j accounts for the source term, and q j is the same as for Eq. (27).

Fluid properties
Concerning phase viscosities, for light oils, the LBC (Lohrentz-Bray-Clark) viscosity model (Lohrenz et al., 1964), usually incorporated in standard commercial reservoir simulators, is used. For heavy oil applications, other models have been integrated. For the oil phase, an exponential law is used for each component in which the parameters a i and b i was computed by solving a leastsquare problem to fit the experimental data. For the water phase, the viscosity was computed using the STARS For the gas phase, the viscosity model is where the parameters c i and d i was computed by regression against experimental data. Finally, the viscosity of the fluid can be computed for each phase based on the relation In all test examples, phase densities were computed based on the EoS model (the Peng-Robinson EoS with volume translation was used), except the extra heavy oil case (SAGD), in which the liquid densities were computed using the STARS formulation with input parameters fitted to the experimental data.

Relative permeability
For three phase systems, the Stone I model (Stone, 1973) is used to interpolate between two-phase relative permeabilities. For a four phase flow, we choose a Brooks and Corey (1964) model as suggested in Varavei (2009) andFeizabadi (2013).

Three-phase bypass
For a given pressure and temperature, if for two different feed compositions z 1 and z 2 the same compositions are obtained at the minimum of the Gibbs energy, z 1 and z 2 belongs to the same tie-simplex. Voskov and Tchelepi (2009a) noticed that for gas injection problems, the solution path involves a limited number of tie-simplex in thermodynamics. Thus, based on a parameterization procedure, they developed a CSAT algorithm (Compositional Space Adaptive Tabulation) which re-uses the information from previous equilibrium calculation results to bypass the stability for redundant compositions. The procedure has been developed in Voskov and Tchelepi (2009b), Voskov and Tchelepi (2009c), Iranshahr (2012), Iranshahr et al. (2012). More recently, Zaydullin et al. (2013) and Zaydullin et al. (2016) extended the procedure to make a more robust algorithm: the three-phase bypass method (called 3 PB afterwards). They reported great improvements in computational times. A presentation of the method is given in the following, starting with the mathematical background.

Working space
The system of equations to solve when performing multiphase flash calculations is given by 8 > > > > > > > > < > > > > > > > > : The overall mole fractions z 2 R nc are usually used as independent variables in reservoir simulations. One of them can be expressed as a linear combination of the others, z nc ¼ 1 À P ncÀ 1 i¼1 z i . Therefore, the working compositional space is given by: The compositional space can be expressed as a linear combination of the mole fractions x This means that for a given tie-simplex Δ defined by its composition at the equilibrium x each feed composition which belongs to Δ is only a function of the phase mole fractions. Therefore, for each tie-simplex, each feed composition (vector of nc-1 dimensions) can be parameterized with only np dimensions. The np dimensional tie-simplex is defined by: When the dimension np of the tie-simplex is np ¼ 2, the tie-simplex is called a tie-line, when np ¼ 3, it is called a tie-triangle. Each tie-simplex Δ np represents a convex geometrical form in the compositional space (such as line, triangle) where each vertex corresponds to the composition of one of the phases. In Fig. 2 is represented a ternary diagram (for a given pressure and temperature). A tie-triangle (corresponds to the three-phase region of composition [x i , y i , w i ]) is represented in grey, and different sets of tie-lines are represented in red, blue and green (which correspond to different two-phase regions). Voskov and Tchelepi (2009c) showed that for a strictly np-phase system, there is a unique tie-simplex which intersects a given feed composition. Voskov and Tchelepi (2007) proved that a tie-line parameterization of the space is possible, for two phases. Voskov and Tchelepi (2009a) extended the idea for any number of phases. As long as the tie-simplex space is not degenerated (i.e. the tie-simplex of dimension N which parameterizes the space does not becomes a tie-simplex of dimension N-1), and the phase mole fractions θ can be greater or smaller than 0, the tie-simplex space can be used to parameterize the equilibrium thermodynamics problem. Based on Voskov and Tchelepi (2009a) observation in which few tie-simplexes are accessed during a reservoir simulation, this parameterization was proposed in order to skip most of the stability analysis during a simulation.
To illustrate this idea, an example is given in Fig. 2. To simplify, let's consider the temperature and pressure as fixed. Supposing the tiesimplexes represented in Fig. 2 have been stored and supposing a stability analysis is required for a composition z i . as z i lies within the parameterized tie-triangle, it can be known directly that at equilibrium, three phases are present for this composition and stability analysis can be avoided. In practice the method is more complex and a description of the method is proposed in the next subsections.

Continuity of a simplex
Tie-simplexes can be used to parameterize the whole system because Iranshahr et al. (2013) proved that the tie-simplex parameterization in the compositional space was a continuous function of pressure, temperature and composition. The continuity property makes a discretization possible. Not only is the tie-simplex continuous, but a tie-simplex which parameterizes the thermodynamics in the tie-simplex space has also been shown to minimize the Gibbs energy . The Gibbs energy surface of a mixture is given by (Baker et al., 1982) FðzÞ ¼ Iranshahr et al. (2012) showed that the Gibbs free energy G is convex:

Parameterization
In this section, the parameterization of the tie-simplex space is presented, first keeping the pressure and temperature constant and then with the variation of both parameters.

Parameterization of tie-simplex giving the maximum number of phases (Δ np )
Starting from a composition z ¼ P nc j¼1 θ j =np, in the middle of Δ np , a negative flash is performed to locate the first tie-simplex leading to the maximum number of phases. An example is given here for three components three phases (Fig. 5a). In this figure, the tie-triangle is represented in red. The composition x gives the np vertices of the tie-simplex. This composition is stored.

Parameterization of tie-simplex plane
Starting from the edge of the tie-simplex, in the same plane, a mesh is adaptively created in the space ½z; θ� (see Fig. 5. a.). This mesh can be more or less accurate depending on its size. The status (number of phases   obtained with the equilibrium program) is saved for each vertex of the mesh. The 3 PB algorithm simplifies the problem approximating continuous phase envelopes by means of discrete envelopes (due to the construction of the mesh). In Fig. 3, a comparison from Iranshahr et al. (2012) between the phase envelopes computed with CSAT (previous version of 3 PB) is carried out with the real one (the same applies for 3 PB).
Once the plane of the tie-simplex has been meshed, a new composition (taken to be at a distance d from the previous plane) is computed. As long as the new composition lies within the tie-simplex in np dimensions, the whole procedure is repeated (see Fig. 4). Only a small number of tie-triangles are stored to keep an efficient procedure. The general space is meshed and interpolations are made. In most of the cases, the composition does not lie exactly within the plane of a stored tie-triangle. This is the same for pressures and temperatures.

Parameterization in temperatures and pressures
A good estimation of the range of pressures and temperatures accessed during the injection process in the reservoir simulator can be known before the beginning of the simulation. Iranshahr et al. (2010a) developed a parameterization procedure for the pressure and temperature for CSAT (the same procedure is used in 3 PB). Here is a description of the methodology. For each variable, a discrete grid is created giving regular values of T and p between the a-priori Tmin and Tmax (pmin and pmax respectively) that will be accessed during the simulation.
For a given number of parameterized temperatures (N T ) and pressures (N P ), a grid is computed with a step: During the process, the interpolation is performed between two adjacent grids, p i < p < p iþ1 and T j < T < T jþ1 . However, just below the MCP (Minimal critical pressure), the parameterization is more difficult. The MCP is the minimal pressure for which the composition is intersected by a critical tie-line (a critical tie-line is given by x i ¼ y i ; 8i ¼ 1; nc). Two nearby discrete temperature values may have quite distinct MCPs. This is why a refinement must be applied to avoid bad interpolations. A refinement procedure is used close to the MCP is used (Iranshahr et al., 2010a). If some pressures and temperatures appear outside of the bounds, the parameterization is extended to include the new equilibrium conditions.

Projection in a tie-simplex space
In the parameterization, only few planes of the tie-simplex space are parameterized for memory and CPU time purposes. Most of the time the compositions z in the tie-simplex space do not belong to any parameterized plane. To obtain the status for a given composition z, a projection of z is made to the closest parameterized tie-simplex. Any point that belongs to the subspace of the tie-simplex is defined by From a given composition z, a projection is carried out in each tiesimplex to find the closest one (in Δ np ). The solution of this projection gives the phase mole fractions θ i . The distance from a feed composition   to a tie-simplex is given by From the tie-simplex giving the smallest distance, if the phase mole fractions θ i ; i ¼ 1; np are all within ½0; 1�, the solution is found. If not, depending on the positive θ i , it is easy to locate the face of the tiesimplex (and then the np-1 region in which the projection belongs to).
Within this framework, if F > ε, it means that new tie-simplexes have to be computed and stored, otherwise, the solution is found.

Status identification
Once the closest tie-simplex from a given composition has been localized, the projection lies within one element of the parameterized mesh presented earlier. Giving the status of each vertex, the stability analysis can be passed by or not: � if all the vertices give the same status, the status is assumed to be the same and the global Newton can be performed directly (Fig. 6b). � if one status is different from the others, the equilibrium flash calculation is computed to determine the real status (Fig. 5b).

Pressure and temperature parameterization and interpolation
The interpolation is first carried out in pressure: Then in temperature: where ψ represents tie-simplex equilibrium compositions and The interpolated tie-simplex is not identical to the actual tie-simplex through the composition and is acceptable only if they are close to each other within a tolerance. If the composition is far from the interpolated tie-line, a new table of tie-lines in the P-T grid is generated.
In each point that is interpolated, one needs to pre-compute an estimate of the MCP to be sure that a parameterization is possible. Moreover, one problem which can occur is when for a given pressure, the two adjacent interpolations give different number of phases. In this case, the multiphase flash is performed to secure the computation.

Two phase case: match with a commercial reservoir simulator
To test MFlash in the reservoir simulation framework, a two-phase case has first been selected. It corresponds to a gas lift problem which is the third comparative SPE test problem (SPE3, Kenyon and Behie, 1987). The simulation domain measures 4022 m � 1609 m � 48.768 m, the porosity is 13%, K v ¼ 10K z , and the medium is heterogeneous. The grid is shown in Fig. 8. The Peng-Robinson equation of state is used, with component properties as given in Table 1 and the non-zero BIPs from Table 2.
In this case, MFlash was integrated within IHRRS (reservoir simulator from Total S.A.). A comparison is proposed between MFlash þ IHRRS and the ECLIPSE commercial simulator by Schlumberger.
The injection wells are shown in yellow in Fig. 8. The production wells are shown in blue in the same figure. Water is injected at 100 STB/ day, gas is produced at a control gas rate of 500 STB/day.
In Fig. 7 and Fig. 8, the water saturation (respectively the oil saturation) at the end of the simulations for MFlash þ IHRRS and ECLIPSE is represented (both simulations lead to the same results). In Fig. 9 the total production rates of water, gas and oil are plotted, for the solution generated with ECLIPSE and the solution generated with MFlash þ IHRRS. The total production rates are identical in the two simulations, validating the two-phase flash developed in this work.
Simulations presented in the following subsections were performed by coupling the reservoir simulator AD-GPRS from Stanford University Zaydullin et al., 2014;Volkov and Voskov, 2016) with the code MFlash used as an external library.

Full isothermal three and four phase compositional simulations of CO 2 gas injection
In this section, three and four phase simulations of CO 2 injection are presented. In our simulations, it has been noticed that a solvent-rich phase could appear when co-injecting solvent with the steam in the reservoir. To develop a robust package, the presence of a fourth phase was necessary to avoid discontinuities which could lead to possible crash of the simulation. So, we investigated the presence of a fourth phase, and developed a four phase CO 2 injection case.
The mixture is the Bob Slaughter Block (BSB) West Texas oil. The characterization of the mixture is taken from Khan et al. (1992). From the experimental data, they provided BIPs and critical properties to fit experimental data; these parameters are given in Table 3. At the reservoir conditions, T ¼ 313.71 K, the P-z phase diagram of the BSB oil is  . 7. Water saturation at the end of the simulation, SPE3 case.
represented in Fig. 10 for pressures varying between 35 and 135 bar and compositions of CO 2 between 10% mol and 99.9% mol. The experimental values of the liquid-vapor/liquid transition and the three phase boundaries, obtained by Khan et al. (1992) are represented in red. The results obtained with MFlash are very similar to the experimental data. The BSB test case parameters (as given in Okuno, 2009) are shown in Table 4. A gas mixture of 5% C 1 and 95% CO 2 is injected at a BHP of 89.63 bar. The producer produces at a BHP of 62.05 bar. The initial reservoir pressure is 75.84 bar. The results obtained with MFlash are represented In Fig. 11 (the injection well is located in the left-hand side of the reservoir and the producer is located in the right hand side). Finally, except for a more important diffusion effect that can be noticed with MFlash, the saturation profiles are similar to those obtained by Okuno (2009).
By injecting gas into the reservoir, close to the injection well, the gas sweeps the oil in place. However, close to the front of saturations both phases are present. At this stage, the oil is not well displaced by the gas, because at these conditions both phases does not reach miscibility. The reservoir pressure is progressively increasing during the injection. Gas and oil become miscible at these conditions and a CO 2 -rich phase appears. The saturation of the oil phase becomes very low, while the CO 2 rich phase saturation becomes significantly increased. This indicates an extraction of the oil components to the CO 2 rich phase. Some authors showed that the CO 2 injection increase the production due to the extraction of medium and heavy oil components to the CO 2 rich phase (Creek and Sheffield, 1993).The density of the CO 2 rich phase is close to the density of the oil phase, indicating a near miscibility of both phases associated with a higher extraction of the heavy components. In the same time, the viscosity of the CO 2 -rich phase remains generally lower than the viscosity of the oil phase.
A four-phase system is obtained by adding water to the BOB fluid. In this simulation, the assigned water composition is z water ¼ 0.1, and the composition of the BOB oil is scaled accordingly. The BIPs between water and the remaining components are set to 0.5. The water properties are given in Table 5. The P-z phase diagram of the BOB oil in presence of water, at T ¼ 313.71 K is given in Fig. 12.
The same simulation test case as for the three-phase CO 2 injection problem is used here in presence of water.
In Fig. 13 a,b,c,d, the pressure and saturations of gas, CO 2 and oil phases are shown, respectively, after 90 days of simulation. The saturations are comparable to those obtained with the previous three-phase CO 2 injection case, yet, stronger diffusion effects take place this time. One of the possible explanations comes from the relative permeability model difference for three and four phases (see the relative permeability section). In this case, at low temperature, the water does not mix with the other phases significantly (water saturation changes only close to the well). These results show that a four-phase simulation is possible with the developed equilibrium code. The flash successfully handles the four-   phase case, providing a good regularity in the saturations (see Fig. 13. b, Fig. 13. c and Fig. 13d).

Stand-alone simulations
In stand-alone, water-hydrocarbon mixtures were tested at different conditions. The first mixture contains three components: C 1 , H 2 O, nC 10 , whose properties are given in Table 6 (fluid property, light oil). Fig. 14a, b shows ternary diagrams. Each system is represented with a different color. The phase transitions obtained by Iranshahr et al. (2010b) are represented in black squares. Good agreements are found with MFlash results. The phase transitions are also well modelled.

Simulation with a medium oil in a highly heterogeneous reservoir
A simulation with a complex steam injection problem (the highly heterogeneous SPE10 reservoir test case) was performed. The mixture is medium oil whose properties are defined in Table 7 (Fluid property, medium steam injection case). The reservoir parameters are given in Table 8. The logarithm of permeability fields of the first and third layers are plotted in Fig. 15a,b.
The initial reservoir parameters are given below: � reservoir pressure: 30 bar at a depth of 2750 m, then the remaining pressure are initialized based on hydro-static equilibrium initialization of the reservoir. � reservoir temperature: 290 K All the wells are vertical. The production wells are located in the four corners of the domain and the injection well is located in the center of the reservoir.   Fig. 11. Three-phase simulation of the CO 2 injection process. Injection parameters: � Injection of steam and CO 2 : 90% of water and 10% of CO 2 . � Temperature of injection: 500 K. � Injection pressure: 50 bar.
Only one time step cut occurred during the entire simulation, revealing a robust flash, with a good phase identification. In Fig. 16. a,b, the gas saturation for the first and third layers of the SPE10 cases are plotted after 600 days of simulation. The saturations reflect the heterogeneities of the porous media. The oil saturations for the first and third layers are shown in Fig. 17. a,b. Finally, the water saturations are plotted in Fig. 18a,b.
The high permeability zone on the right-hand side of the field acts as a barrier and prevents the expansion of the steam there. The steam and water flushes the oil to the two producers located in the left-hand side.
The simulation is stopped with the water breakthrough. During the steam injection in the reservoir, the high injection temperature decreases the viscosity of the oil and increases its mobility. The oil is swept more easily leading to the enhanced oil recovery. Fig. 16. b shows the high gas saturation around the injection well. Then, away from the injection well, the water and light components condense at the contact with the cold region. Fig. 18. b shows the higher saturation of water when the gas disappears. Finally, the oil saturation field analysis shows that the major volume of oil has been produced from the left part of the reservoir. It must be noted that with three dimensional heterogeneous reservoirs, more thermodynamic paths are accessed during a simulation. Once again, this example reveals the robustness of the developed algorithm.

Simulation of SAGD process with a synthetic heavy oil and an extraheavy oil
Two horizontal wells are located one under the other (the injection well above). As steam is injected, a steam chamber is formed over the injection well. At the chamber edge, the water and the light oil components condense. Besides, the elevated temperature decreases the oil viscosity and the oil becomes mobile. The condensed water and the extra-heavy oil flows to the production well under the gravity.
Two simulations of the SAGD process have been performed. The first considers a synthetic heavy oil, initially at T ¼ 285 K. The oil is not heated before the steam injection. In the second case, a real Athabasca bitumen has been modelled. In this case, the oil is pre-heated before the steam injection.  Fig. 13. Four-phase simulation of CO 2 injection, t ¼ 50 days. The composition of the synthetic heavy oil is given in Table 9. The initial parameters are: � reservoir pressure: 8.11 bar at the depth 180 m, then the remaining pressures are initialized based on hydrostatic equilibrium initialization of the reservoir. � reservoir temperature: 285 K � the initial oil viscosity oil viscosity at T ¼ 285 K is 8.86eþ7 cP The injection well is located at the depth 117.5 m and the production well, 5 m below. Injection parameters are: � fraction of steam (99.4 mol %) and CO 2 -C 1 (0.6 mol %). � injection temperature: 500 K � injection pressure: 25 bar.
Only half of the steam chamber (SC) cross section is shown in Fig. 19 and Fig. 20. Each cell is a cube of length 0.5 m, the grid dimension is made by 101 � 1 � 52. Fig. 19. a,b and Fig. 20. a represent the saturations of water, gas and oil, respectively, after 150 days of simulation. The temperature is given    in Fig. 20. b for the same time. The SC develops with the characteristic triangular shape of the SAGD process. In Fig. 19. a, an accumulation of the water phase can be seen along the SC. Along the SC, the high temperature decreases the liquid viscosities, enabling the oil to flow to the production well under gravity drainage. The second oil used for SAGD modelling, is an Athabasca oil containing extra-heavy components. The fluid properties are given in Table 10. The initial parameters and the grid are the same as for the synthetic heavy oil. At T ¼ 283K, the initial oil viscosity is μ ¼ 6.64eþ11 cP. Steam is injected (99.5 mol %) with solvent (0.5 mol %). The solvent is the M0 pseudo-component made of C 3 -C 10 components. Fig. 21a,b and Fig. 22. a show the saturation fields (of water/gas and oil respectively) at t ¼ 360 days. At this time, the SC has been developed, growing mainly vertically. Water and oil accumulation can be seen along the steam chamber edge. Both heated liquids are flowing to the production well by gravity drainage. Moreover, high concentration of gas can be seen at the top of the SC. The light components with a lower density move to the top of the SC and encounter the liquid phases (oil and water) at the edge of the chamber. Fig. 22. b shows the temperature profile in the reservoir for the same time. Due to the gravity effects, the gas mixture moves to the top of the reservoir. Once the steam chamber reaches the top, it starts growing horizontally. Finally, at t ¼ 410 days,        the steam chamber reaches the top of the reservoir; Fig. 23 shows the saturation of the gas phase at this moment. Such thermal compositional simulations are quite complex because important changes in the conditions are observed (composition, pressure and temperature). Indeed, for the steam injection process, compositional effects occur such as steam distillation and water/light oil components condensation. With heavy oils, the gravity creates concentration gradients which imply many changes in the oil composition. When injecting steam, it flushes progressively all the other components. The phase envelope of the resulting composition becomes narrow and the multiphase problem approaches ill-conditioning where a small variation in the conditions can create an important change in the equilibrium state. Simulation of the SAGD process up to 410 days revealed the capability of our algorithms to overcome these difficulties and to converge even for extremely difficult (nearly ill-conditioned) problems.

Comparisons between different compositional acceleration procedures
Different acceleration techniques were also tested with MFlash: � the reduction variables (RV) for the stability analysis. It was noticed that under certain circumstances, the reduction was decreasing the computation time of the stability analysis. We propose here to test the variables in the reservoir simulation framework. � the bypass method from Rasmussen et al. (2006) (labelled B-R). � the three-phase bypass method (3 PB) presented in section 4 (labelled 3 PB).
These algorithms were compared with simulations using MFlash with conventional variables (labelled Conventional).
Computation times were tested for three different cases: the threephase CO 2 injection case, the light oil steam injection case and the SPE10 steam injection case. Table 11 show the total time of simulation for the different algorithms. Table 12 gives the percentage of the equilibrium calculations in the total CPU time of the simulation. Tables 11 and 12 reveal a different CPU time for the CO 2 injection and the steam injection cases. When dealing with CO 2 and liquid-liquid phase behavior, more initial guesses are needed than for steam injection problems. This is why the equilibrium computation represents 70% of the total simulation for the CO 2 injection case, whereas it is around 40% for the steam injection test cases. Tables 11 and 12 also show that the total simulation time is generally decreased when the reduction method is used. It is particularly true for the SPE10 test case (8166 s for reduced variables vs 11,865 s for conventional variables). However, in this study, the use of the B-R method as developed in Rasmussen et al. (2006), in presence of water does not improve the computational time for the tested cases. Indeed, the water initial guess for the stability analysis often converges to a minimum which is different from the trivial solution and no phase equilibrium can be bypassed.
The use of the 3 PB coupled with MFlash is the option leading to the smallest computation times for all the cases. The 3 PB method uses a negative flash procedure to compute supporting tie-simplex and uses multiphase calculations to adaptively parameterize the tie-simplex space. The results obtained with this method were identical to MFlash without any acceleration technique, and demonstrated a good robustness. In terms of performances, 3 PB re-uses the information from previous equilibrium calculations to pass by the computation of further stability analysis. A lot of stability testing can be avoided when using this procedure. Here, the use of 3 PB improves the time of computation from 8 up to 12 times as compared with MFlash. For the whole simulation, with 3 PB the averaged computation time is 1.6 times faster than   with MFlash alone. However, using a flash based on global optimization, or with a higher number of initial guesses, or with a more accurate and computationally expensive EoS model, the differences may increase significantly.
To check the robustness of the 3 PB, a comparison was performed for the synthetic oil SAGD test case. Table 13 gives the computational results obtained with MFlash with reduction variables and 3 PB for 90 days of simulation. Taking the flash in stand-alone as the reference, the relative error on the oil production rate for the 3 PB simulation remained always under 6%, taking large steps (system ill-conditionned where convergence problems occur). By taking smaller steps, the error becomes negligible. In this case, the total computation time with 3 PB was much smaller than for MFlash: (2529s versus 3914s), decreasing the fraction of the equilibrium calculations from 45% to 10% of the total simulation. 3 PB seems quite robust and capable of handling simulations of difficult processes quite efficiently.
Note that all the tests were carried out using only one core. In each cell, the equilibrium calculation does not depend on properties from other cells (local computations), therefore equilibrium calculations can be efficiently parallelized. Using a parallel procedure, the percentage of the total computation time spent in the equilibrium calculations would become smaller, and the CPU time required for equilibrium calculations will represent a smaller proportion of the total simulation time. Finally, analyzing all the simulations (except for CO 2 injection where more initial guesses are required), one can see that the time of the equilibrium calculations with or without the use of 3 PB always remains under 50% of the total time, even for the most computationally expensive SAGD test case.
Finally, several points on the advantages and limitations of the presented approach and possible future work follow. In this paper, a general purpose algorithm to solve multiphase flash calculations has been developed. All phase equilibrium calculations are thermodynamically consistent (unlike the K-value approaches used for thermal compositional simulations in most commercial simulators). Besides, the approach is valid over a wider range than for other models such as Raoult's law or Henry constant (also used in some simulators). A comparison with commercial simulators could possibly show that the methodologies in this paper can reveal more mechanisms than the conventional ones and would be an interesting future work. In this paper, cubic equations of state have been used (as giving relative accuracy within an acceptable computational time, and also as the one implemented in standard reservoir simulators). Those kinds of equations of state do not perfectly model water solubility in hydrocarbon phases (especially the significant solubility of water in bitumen phase at high temperature) or CO 2 solubility in water. An extension of this work would be, at first, to implement a Søreide and Whitson (1992) approach. A cubic-plus-association equation of state (CPA-EoS) would greatly improve hydrocarbon-water solubility modelling (Kontogeorgis et al., 1996;Zirrahi et al., 2012). Viscosity of the bitumen cannot be entirely captured by conventional viscosity correlations; another improvement of the physical model can be obtained by using a double-log viscosity correlation for bitumen, as suggested by Zirrahi et al. (2017). Compositional thermal simulations are also highly critical for bitumen production, especially with the current interests on the new solvents with mutual solubility in aqueous and oleic phases (Zirahi et al., 2020).

Conclusions
In this work, improvements in equilibrium algorithms made possible the realization of fully compositional simulations of steam injection with heavy oils. The phase equilibrium calculations code is based on both Trust-Region and reduction methods that improve both robustness and efficiency. For the stability testing problem, in reservoir simulation purposes, the proposed reduction method has proved to be the most efficient while keeping the same convergence path as in the conventional formulation. The good reproduction of various phase diagrams (LLV-LLLV) illustrates the capability of the program to predict the true phase behavior of complex systems.
Simulation test were performed with three-and four-phase (when water is added) CO 2 injection problems, and with several steam injection cases, including a 3D steam flooding simulation for the highly heterogeneous SPE10 reservoir and a fully compositional reservoir simulation of the SAGD process for a realistic extra heavy oil mixture. Reservoir simulations cover a great variation of p-T-z conditions, including difficult situations, such as crossing phase boundaries and/or near-critical conditions. Any failure in the convergence to the correct minimum can lead to the divergence of a simulation. The capability to carry a whole simulation in all cases demonstrates the robustness of the phase equilibrium code.
The proposed thermodynamic code was also used with different acceleration procedures, namely, based on the bypass method by Rasmussen et al. (2006) and the three-phase bypass by Zaydullin et al. (2016). Without any acceleration technique, the program has taken around 45% of the total CPU time for steam injection problems. In presence of water for multiphase systems, the bypass method did not perform quite well. However, with 3 PB, the CPU time was decreased and reached around 10% of the total time, keeping the same solution as the thermodynamic program alone.