{'fraction_non_alphanumeric': 0.03184993405305066, 'fraction_numerical': 0.018172048263397002, 'mean_word_length': 4.798470471604588, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Smooth Formulation for Multi-component Compositional Simulation with Superior Nonlinear Convergence

Compositional simulation is challenging, because of highly nonlinear couplings between multi-component flow in porous media with thermodynamic phase behavior. The coupled nonlinear system is commonly solved by the fully-implicit scheme. Various compositional formulations have been proposed. However, severe convergence issues of Newton solvers can arise under the conventional formulations. Crossing phase boundaries produces kinks in discretized equations, and subsequently causing oscillations or even divergence of Newton iterations. The objective of this work is to develop a smooth formulation that removes all the property switches and discontinuities associated with phase changes. We show that it can be very difficult and costly to smoothing the conservation equations directly. Therefore, we first reformulate the coupled system, so that the discontinuities are transferred to the phase equilibrium model. In this way a single and concise non-smooth equation is achieved and then a smoothing approximation can be made. The new formulation with a smoothing parameter provides smooth transitions of variables across all the phase regimes. In addition, we employ a continuation method where the solution progressively evolves toward the target system. We evaluate the efficiency of the new smooth formulation and the continuation method using several complex problems. Comparing to the standard natural formulation, the developed formulation and method exhibit superior nonlinear convergence behaviors. The continuation method leads to smooth and stable iterative performance, with a negligible impact on solution accuracy. Moreover, it works robustly for a wide range of flow conditions without parameter tuning.


Introduction
Gas injection processes play an important role in enhanced oil recovery (EOR). Gas injection into oil reservoir involves a number of physical mechanisms that help in mobilizing and extracting oil. Depending on pressure, temperature, and fluid compositions, immiscible or miscible displacements could occur. The physical model required to describe the mass transfer between phases is the isothermal compositional model (Coats et al. 1998). An Equation of State (EoS) model is usually used to determine the component partitioning across phases (Michelsen 1982).
Compositional simulation continues to be a challenging problem. Complexities are mainly due to nonlinear couplings between multi-phase multi-component flow in porous media with thermodynamic phase behavior (Voskov 2012). Several temporal discretization schemes are available to solve the compositional model (Aziz and Settari 1979).
Explicit schemes pose severe restrictions on the timestep size and are impractical for large-scale detailed reservoir models, in which the Courant-Friedrichs-Lewy (CFL) numbers can vary by several orders of magnitude throughout the domain (Jenny et al. 2009). Therefore, the fully-implicit method (FIM) is preferred in practice, with the nonlinear system solved by a Newton method.
Various nonlinear formulations for the compositional model have been proposed (Aziz and Wong 1989). Two popular formulations employ two different variable sets: (1) natural variables (Coats 1980) and (2) molar variables (Acs et al. 1985). The primary unknowns for the natural-variables set include pressure, saturations, and component molar fractions. The conservation equations and thermodynamic constraints are assembled for cells with two phases. If one of the phases disappears during a nonlinear iteration, the corresponding saturation variable is removed, and the nonlinear system is reduced to the conservation equations. This process is referred to as 'variable substitution' and is an essential ingredient of the natural-variables formulation (Aziz and Wong 1989;Cao 2002). Under the molar formulation, one variable-set choice is the overall composition of each component. This approach has the advantage of avoiding the need for variable substitution, since the equations and unknowns are the same for any phase state. Numerical behaviors of nonlinear formulations were investigated and compared in some recent studies Alpak and Vink 2018).
For large timestep size, or poor initial guess, the standard Newton method suffers from serious convergence difficulties. Obtaining a suitable initial guess for the Newton method is referred to as globalization (Knoll and Keyes 2004). Damping, or safeguarding Newton updates is a class of globalization techniques (Deuflhard 2011). One simple heuristic strategy is to apply a local damping of variable to ensure that the value change is limited to a pre-defined range (ECLIPSE 2008). Physics-based trust-region solvers were also introduced to greatly improve the nonlinear convergence of discrete transport problems (Jenny et al. 2009;Wang and Tchelepi 2013). Voskov and Tchelepi (2011) developed a nonlinear solver specifically for molar-variables formulations. Trust regions are constructed based on the flux functions along key tie-lines in compositional space. The Newton updates are controlled from crossing inflection points and phase boundaries. Although locating trust-region boundaries is straightforward in two-phase problems with simple flow physics, it can be quite difficult for general compositional models with gravity and phase changes. Moreover, computations become expensive when the flux functions vary significantly during the iterations.
Recent studies revealed that the non-differentiability (kink) resulted from switching criteria in the numerical flux may cause frequent oscillations or divergence of Newton iterations (Wang and Tchelepi 2013;Lee et al. 2015;Hamon et al. 2016;Jiang and Younis 2017). Likewise, severe convergence problems of Newton solvers can also arise under the conventional compositional formulations. This is due to the kinks produced when crossing phase boundaries. Phase change leads to the corresponding switches in fluid properties and discretized equations.
Several formulations were proposed, attempting to improve the nonlinear performance of compositional simulations (Abadpour and Panfilov 2009;Lauser et al. 2011;Voskov 2012;Gharbia et al. 2015). Unified system of equations is achieved and thus the variablesubstitution process is avoided. Saturations (or phase fractions) can change continuously across phase boundaries under these formulations. However, the switches of fluid properties still occur, in one form or another. The kinks in discretized equations as the essential mechanism that causes convergence difficulties are not resolved by the previous works. Consequently, their results reported insignificant improvement in nonlinear performance comparing to the standard natural formulation.
Previous studies demonstrated that a smooth numerical flux can improve Newton behaviors (Wang and Tchelepi 2013;Lee et al. 2015;Hamon et al. 2016;Jiang and Younis 2017). The objective of this work is to develop a smooth formulation that removes all the property switches and discontinuities associated with phase changes. But as we will show, it can be very difficult and costly to make smoothing approximations directly for the conservation equations under the standard formulations. Therefore, we first reformulate the coupled system, so that the discontinuities are transferred to the phase equilibrium model. In this way a single and concise non-smooth equation is achieved and then a smoothing approximation can be made. The reformulation is based on a mixed complementarity problem (MCP) proposed for the phase equilibrium in the field of chemical process simulation (Bullard and Biegler 1993;Gopal and Biegler 1997;Sahlodin et al. 2016). The MCP model contains complementarity conditions that represent phase changes. The subsequently developed formulation with a smoothing parameter can lead to smooth transitions of variables across all the phase regimes.
For a smoothing parameter with sufficient values, the smooth coupled system exhibits superior global convergence behavior. However, solution accuracy may be degraded, with a fixed smoothing parameter. Its value needs to be adaptively determined to achieve an optimal balance between accuracy and nonlinear performance. For more robust and general applications, we employ a continuation method with the smooth system as a homotopy mapping. The smoothing parameter is progressively reduced after each Newton iteration, and the solution evolves toward the target system. In this work, the continuation method acts as a globalization stage to obtain better initial guesses for the Newton process (Jiang and Tchelepi 2018). As will be shown, the new smooth formulation is very effective for globalizing the compositional flow problem. It is worth noting that the smooth formulation and the continuation method can be used for different applications depending on specific accuracy and implementation considerations.
We evaluate the efficiency of the new approaches using several complex examples. We focus on the nonlinear behavior of the coupled conservation and phase equilibrium system. The smooth formulation produces valid and unique solutions in all the three phase regimes. For most of the cases, the standard natural formulation suffers from multiple timestep cuts and subsequent wasted Newton iterations. In contrast, the smooth formulation and the smoothing based continuation method (SBC) exhibit superior global convergence, requiring no timestep cut. Moreover, the SBC method can largely resolve the convergence issues due to phase changes, with a negligible impact on solution accuracy. We find that applying SBC for a few iterations in the globalization stage is sufficient. In addition, the developed method works robustly for a wide range of flow conditions without parameter tuning.

Isothermal compositional model
We consider compressible gas-oil flow in porous media without capillarity. We ignore water that does not exchange mass with the hydrocarbon phases.
The conservation equations for the isothermal compositional problem containing n c components are written as, where c ∈ {1, ..., n c }. x c and y c are molar fractions of component c in the oil and gas phases, respectively. φ is rock porosity and t is time. ρ l is phase molar density. s l is phase saturation. q c is well flow rate. Phase velocity u l is expressed as a function of phase potential gradient ∇Φ l using the extended Darcy's law, where k is rock permeability. p is pressure. Capillarity is assumed to be negligible. g is gravitational acceleration and h is height. Phase mobility is given as λ l = k rl /µ l . k rl and µ l are relative permeability and viscosity, respectively.
Phase velocity can also be expressed under the fractional-flow formulation, which involves the total velocity, where the total mobility λ T = l λ l . In order to close the nonlinear system, additional equations are needed. These include the thermodynamic equilibrium constraints, where p, T , and z c denote pressure, temperature, and overall molar fraction, respectively. f c,l is the fugacity of component c in phase l.
We now write the phase constraints, and the saturation constraint as, The above system of equations provide a complete mathematical statement for twophase multi-component flow. The local equilibrium constraints are enforced only when both phases are present.

Natural-variables formulation
An important aspect of any compositional formulation is the choice of dividing the equations and unknowns into primary and secondary sets. Two widely used formulations are (1) natural variables (Coats 1980) and (2) overall-composition variables (Collins et al. 1992;Voskov and Tchelepi 2012).
The primary unknowns include pressure, saturations, and molar fractions: The various coefficients can be obtained as functions of the base variables. For a two-phase cell, the molar phase fraction is related to saturation as follows, and overall molar fraction of component c is written as, Note that for single-phase (l) mixture, ν l = s l = 1, and x c,l ≡ z c .

Variable substitution
An essential ingredient of the natural-variables formulation is the 'variable substitution' process (Aziz and Wong 1989;Cao 2002;Voskov and Tchelepi 2012). A common strategy for variable-switching between Newton iterations during a time step is: 1. For any cell whose status in the previous iteration is single-phase, run the phase stability test (Michelsen 1982a) to check if the mixture becomes two-phase. For the mixture that splits into two phases, perform the flash to compute the phase compositions (Michelsen 1982b).
2. If a cell is already in the two-phase state, the thermodynamic constraints are included in the nonlinear system as part of the global Jacobian.
3. If a phase saturation, or phase fraction, becomes negative between two successive iterations, the phase disappears, and appropriate variable-switching is performed.
The system of conservation equations is solved for single-phase regimes, and the combination of conservation equations and thermodynamic constraints is solved for the two-phase regime.

Phase behavior
Phase behavior computation is usually a stand-alone procedure for detecting phase changes. For a mixture of n c components and two phases, the mathematical model describing the thermodynamic equilibrium is , where ν l is molar fraction of phase l. We assume that p, T , and z c are known. The objective is to find all the x c , y c and ν l . Phase behavior of a hydrocarbon mixture is commonly described using an Equation of State (EoS) model. The EoS-based computations are expensive and may consume a large portion of the total simulation time. But their cost can be significantly reduced through some advanced approaches (Voskov and Tchelepi 2009;Zaydullin et al. 2012;Yan et al. 2013).
In this work, we simplify the phase behavior computations using the K-value based method. We assume that components partition across phases with fixed ratios (Kvalues), which depend on pressure and temperature. Then the fugacity constraint (10) can be rewritten as, where K c is the equilibrium ratio.

Discontinuous issues crossing phase boundaries
Recent studies have revealed that the non-differentiability (kink) in the numerical flux due to switching criteria can be a major cause of nonlinear convergence difficulties (Wang and Tchelepi 2013;Lee et al. 2015;Hamon et al. 2016;Jiang and Younis 2017). The curvature of residual function changes abruptly at the kink, leading to oscillations between Newton iterations (flip-flopping) or convergence failure.
In a single-phase regime, the component compositions are no longer controlled by the thermodynamic equilibrium. Then Eqs. (12) and (13) cannot be satisfied at the same time. The compositions of the existing phase can only be obtained from Eq. (11), where ν l is set to 1 (or 0). For the missing phase, the compositions become undefined.
Phase change leads to the corresponding switches in fluid properties and discretized equations. To demonstrate this, we present the accumulation term in (1) for the different phase states as, where ρ o (x) indicates the density computed at a composition x, and ρ l (z) is the density computed in the single-phase regime at a composition z.
In the absence of gravity, the overall fractional flow function of component c derived from (1) becomes, By assigning unit values to ρ l and µ l , Eq. (15) for component 1 can be further simplified as, We consider a two-component system with {K 1 = 3.5, K 2 = 0.3} and quadratic relative permeability. Fig. 1 shows the plot of u 1 as a function of z 1 . Clearly, we observe two non-differentiable points associated with the phase changes. Within the two-phase regime, the flow curve has the typical S-shape. The flow curve becomes a straight line for the single-phase states. The discontinuous derivatives present in the flow function can largely degrade the convergence performance of Newton solvers. Several nonlinear formulations were previously proposed for compositional simulations (Abadpour and Panfilov 2009;Lauser et al. 2011;Voskov 2012;Gharbia et al. 2015). Through reformulation or certain consistence conditions, unified system of equations is achieved so that the variable-substitution process is avoided. Although saturations (or phase fractions) can change continuously across phase boundaries under these formulations, the switches of fluid properties still occur, in one form or another. The kinks in discretized equations as the essential mechanism that causes nonlinear convergence difficulties are not resolved by the previous works. Consequently, the results reported insignificant improvement in nonlinear performance comparing to the conventional natural formulation.

Smooth formulation
From the previous studies (Wang and Tchelepi 2013;Lee et al. 2015;Hamon et al. 2016;Jiang and Younis 2017), it is expected that a smooth system can provide much improved Newton behaviors. According to Eqs. (14) and (15) the nonlinearities (discontinuities) associated with phase changes are due to the switches of properties from the accumulation and flux terms. We can see that it is hard to develop smoothing approximations directly for those terms in the conservation equations.
Therefore, we seek to first reformulate the coupled system, transferring the discontinuities to the phase equilibrium model. In this way a single and concise non-smooth equation is achieved and then a smoothing approximation can be derived.

Non-smooth reformulation
In the field of chemical process simulation, a mixed complementarity problem (MCP) was proposed for the phase equilibrium system (Bullard and Biegler 1993;Gopal and Biegler 1997;Sahlodin et al. 2016).
The MCP model containing complementarity conditions that represent phase changes can be equivalently written in a non-smooth form as, where c ∈ {1, ..., n c }. The system R = 0 is solved for the given z c and K c . The mid function picks the median of its three arguments, for any (a, b, c) ∈ R 3 . A non-physical variable β is introduced to relax the equilibrium equation. We can see that the conditions for the three phase states are given as, The solutions obtained from the non-smooth system (17 − 20) are valid and unique in all the three phase regimes. We can see that the phase fractions (thus saturations) are bounded between 0 and 1. The variables for the non-existent phase can be viewed as pseudo molar fractions that also sum to one. For the two-phase state, the non-smooth system provides the same results as the conventional phase equilibrium model.
In this work the fugacity constraint is evaluated through K-values. The underlying assumption does not limit the applicability of the proposed formulation to a general EoS-based system.

Smoothing approximation
To present the smooth formulation for the MCP model, we first write the mid function as a mixed complementarity function, Then we construct a smoothing function for Eq. (23) as following (Chen and Harker 1993;Kanzow 1996;Li and Fukushima 2000), where ≥ 0 is a smoothing parameter, and Θ ( , a, b, c) is continuously differentiable with > 0. We can easily see that, whenever a > c.
There are several other ways to make smoothing approximations for the mid function. We choose Eq. (24) because the function is computationally simple and is well-behaved in terms of nonlinearity. The phase equilibrium model with (24) provides smooth transitions of the variables across all the phase regimes.
Again we consider the two-component system from the last section. We compare the solutions between the standard flash and smooth formulations. Fig. 2 shows the compositions y 1 and x 1 as a function of z 1 . Fig. 3 compares the results of ν g and u 1 .
Note that the negative flash concept is applied for allowing phase fractions (saturations) to exceed the bounds of 0 and 1 in the standard flash. Even though the compositions are constant along the tie-line, the discontinuities across the phase boundaries will still appear in the accumulation and flux terms. In contrast, we can see that the smooth formulation with a fixed provides continuously differentiable functions of all the variables.
We also study a ternary system with {K 1 = 3.5, K 2 = 1.5, K 3 = 0.2} and quadratic relative permeability. The composition z 2 = 0.3 has constant value. Fig. 4 shows the compositions y 1 and x 1 as a function of z 1 . We do not further plot overall fractional flow here because of its similar form as the two-component case.
In addition, we study a two-component fluid {C 1 (60%), C 4 (40%)} at a temperature of 250 K. The K-values for the two components are computed from Wilson's equation, where p c and T c are critical pressure and temperature, respectively. We plot the compositions y 1 and x 1 as a function of p in Fig. 5. The gas phase fraction versus pressure is shown in Fig. 6. Comparing to the standard flash, the smooth formulation can achieve smooth transitions with respect to pressure. This is very beneficial for the compositional scenario driven by pressure. Lauser et al. (2011) proposed a compositional formulation that is also based on complementarity conditions for handling phase changes. Comparing to the smooth formulation proposed in this work, their approach has several limitations: 1. For a nonexistent phase in single-phase regimes, the sum of the molar fractions becomes less than one. This can produce additional discontinuities in evaluating fluid properties. Moreover, there is no guarantee that the values of molar fractions remain positive under all conditions, especially if smoothing approximations are applied.
2. The phase equilibrium and the compositional systems are non-smooth in nature. Therefore a Newton solver will still suffer from the discontinuous issues due to phase changes. 3. The simulation studies demonstrated only insignificant improvement in nonlinear convergence, compared with the conventional natural formulation (Lauser et al. 2011;Gharbia et al. 2015). We give detailed descriptions and discussions on the approach of Lauser et al. (2011) in Appendix A.

Smooth formulation for coupled system
For the smooth formulation, the primary unknowns include: The coupled system contains: the [n c ] conservation equations (1), the phase (6) and saturation (7) constraints, and the [n c ] relaxed equilibrium equations,  (24), The smoothing parameter can be kept constant during simulation. It is expected that a suitable value of will be chosen for a target class of problems. Comparing to the naturalvariables set, an obvious advantage of the new smooth formulation is that the equations and unknowns are the same for any phase state. The complex variable-substitution process is thus avoided. For general multi-component models, phase boundaries are quite complex, which makes it difficult and costly to apply damping (safeguarding) based solution strategies. One particular type of discontinuity is the transition between two tie-lines (Orr 2007).
As previously demonstrated, the saturations are bounded to [0, 1] in the smooth system, so that the physical limiting is not necessary for the conservation equations. As a result, all the discontinuities associated with the phase changes transfer to the single and concise equation (28), and get removed through the smoothing approximation. The coupled system becomes smooth across the entire phase boundaries, bringing enormous benefits for Newton solvers.
Here K-values are used to describe the phase equilibrium. Nevertheless, the new formulation is generally applicable to EoS-based systems.

Nonlinear solver
The spatial and temporal discretization schemes used for the compositional flow model are summarized in Appendix B.

Newton method
At each timestep of a FIM simulation, given the unknown vector U n and a fixed timestep size ∆t, we intend to obtain the new state U n+1 .
The nonlinear system is cast in residual form and solved by the Newton method, The Newton method comprises a sequence of iterations, each involving the construction of a Jacobian matrix and solution of the resulting linear system, where and J (U η ) = ∂R ∂U U η denotes the Jacobian matrix of R with respect to U η . The iterative process is performed until the nonlinear system is converged.
The formulation and method developed in this work are integrated into the Automatic Differentiation General Purpose Research Simulator (AD-GPRS) (Younis and Aziz 2007;Zhou et al. 2011;Voskov and Tchelepi 2012). The nonlinear framework of AD-GPRS is built on top of Automatic Differentiation with Expression Templates Library (ADETL). The Jacobian is automatically assembled and analytically derived through ADETL.

Homotopy continuation method
The Newton process may often converge slowly, or even diverge, due to poor initial guess and large timestep size (Younis et al. 2010). The smooth formulation has nice global convergence property. However, solution accuracy may be degraded, with a fixed smoothing parameter . The value of needs to be adaptively determined for an optimal balance between accuracy and nonlinear performance.
The coupled system with the smoothing equation (28) can be viewed as a homotopy mapping H. The objective is to solve the original system R(U ) = 0 containing the non-smooth mid function. We can see that H(U, 0) = R(U ) with the modified residual H(U, ) and the smoothing (continuation) parameter ≥ 0. Consider that H(U, 0 ) = 0 is much easier to solve than the target problem H(U, 0) = 0. The continuation method can be developed by discretizing in to form a sequence of nonlinear systems H(U, η ) = 0. The target solution is reached by progressively decreasing from 0 towards 0 to globalize the flow equations. As a result, the algorithm will not have an impact on the accuracy of the final solutions.
The continuation method only acts as a globalization stage to obtain better initial guesses for the Newton process (Jiang and Tchelepi 2018). In our experience, applying the continuation for a few iterations in the globalization stage is already highly effective. A simple strategy is used to evolve during each timestep, with 0 as the initial value of a timestep. After each Newton iteration, is multiplied by a constant γ. η max is the number of iterations taken for globalization. The parameter values used for all the following cases are: 0 = 0.1, η max = 4 and γ = 0.5. min is specified to ensure that the system is smooth enough at the convergence limit. From the simulations we find that a value around min = 1.0e-4 can produce the solutions with satisfactory accuracy and nonlinear convergence. We note that the smooth formulation and the continuation method involve different complexities and efforts for implementation. Subsequently, their applications depend on specific accuracy and implementation considerations.

Results: 1D model
Nonlinear convergence is based on the following criterion: max i,c |R i,c | < 10 −4 , where R i,c is the residual of c-th equation in i-th cell (for conservation equations, we normalize by the total mass of c-th component).
A simple time-stepping strategy is employed: if the Newton method converges for the current timestep, the next timestep will be doubled; if the nonlinear solver fails to converge, the timestep is reduced by half and solved again. The solution from previous timestep n is taken to be the initial guess for a new timestep. During the iterative process, all fractions in the variable-set are kept within the physical interval [0, 1]. Also the damping strategy is employed to stabilize the Newton updates: the local chopping values of 0.2 are used for saturations and 0.1 for molar fractions.
We evaluate the efficiency of the new approaches using several complex problems. Five different fluids, and two different reservoir models are considered. The models include: (a) homogeneous 1D model; (b) heterogeneous 2D model taken from the bottom layer of the SPE 10 problem. In the following cases, simple relative permeabilities given by quadratic function are used, unless otherwise indicated. Note that both phase density and viscosity depend on pressure and compositions. Phase molar density ρ l is evaluated based on the compressibility (Z) factor from the PR EoS. Phase viscosity µ l is computed by the correlation of Lohrenz et al. (1964). The specification of the 1D base model is given in Table 1. We test a comprehensive suite of 1D problems. The first problem is the displacement of propane C 3 by methane C 1 in a horizontal domain. Cell size of 1m is used. The fluid and compositions are given in Table 2. Z i and Z inj are initial and injection compositions, respectively. Pressure is kept constant as 65 bars at the production end. A constant volumetric rate 5 m 3 /day is specified at the injection end. Initial pressure is 70 bars and temperature is 311 K.  Table 3. We report the total number of Newton iterations and timesteps. In parentheses, we also give the number of wasted iterations that correspond to the iterations spent on unconverged timesteps. 'SBC' denotes the smoothing based continuation method, and min is specified for Eq. (32). For each simulation, the maximum CFL number, averaged over the timesteps taken, is reported. The corresponding maxCFL of the case is 16.
As we can see from the results, there is no timestep cut during the simulations. Even for this simple scenario, the SBC method requires much fewer iterations to converge.  (0) 12 (0) Gas saturation and overall composition profiles are plotted in Fig. 7. The profiles show two shocks formed in the two-phase region. The solutions from SBC with min = 1.0e-4 fully match the standard natural formulation. For min = 1.0e-2, some differences around the shocks are observed.

Case 2
We also test a case with gravity. Temperature is set to 380 K. The maximum time step size reduces to 30 days, due to severe degradation of nonlinear convergence encountered by the standard natural formulation. The total simulation time becomes 400 days. The other parameters from the previous case remain unchanged.
The nonlinear performance of the case with gravity is summarized in Table 4. The maxCFL is 20. As can be seen, the standard formulation suffers from multiple timestep cuts and subsequent wasted iterations. In contrast, the SBC method does not require any timestep cut, resulting in much smaller number of iterations. We further study a case with variable K-values computed from Eq. (24). Production pressure changes to 70 bars. Initial pressure is 80 bars at a temperature of 370 K. The maximum timestep size becomes 50 days, with a total simulation time of 400 days.
Here we simply apply the smooth formulation without the continuation method. The case is to demonstrate the applicability of the smooth formulation using a fixed smoothing parameter . The value level of = 1.0e-2 is sufficient to provide much improved convergence performance. Solution profiles are not shown because the differences between the two formulations are small.
The nonlinear performance of Case 3 is summarized in Table 5. The maxCFL is 21. As can be seen, the pressure-dependent K-values make it very challenging for the standard formulation. By comparison, the smooth formulation is continuously differentiable with respect to pressure, leading to significant convergence speedup.

Three-component fluid
We test cases with three-component fluid systems. Cubic relative permeabilities are used. Initial pressure is 70 bars and temperature is 320 K. Pressure is kept constant at the both injection and production ends. Production pressure is 60 bars, with injection pressure as 130 bars. The maximum timestep size is 50 days, with a total simulation time of 400 days.

Case 1
The fluid and compositions for Case 1 are given in Table 6. The nonlinear performance of Case 1 is summarized in Table 7. The maxCFL is 6. In this case, the SBC method does not require any timestep cuts.  (0) The gas saturation and overall composition profiles are plotted in Fig. 8. Here the timestep size is reduced to 10 days, to ensure that the solutions from the different methods are compared under the same time-stepping schedule. As we can see, SBC produces solutions that are very close to the standard formulation.

Case 2
We also study a three-component system with {K 1 = 2.5, K 4 = 1.5, K 10 = 0.3}, and the other parameters remain unchanged.
The nonlinear performance of Case 2 is summarized in Table 8. The maxCFL is 6. The standard formulation suffers from a large number of timestep cuts and wasted iterations, despite a small CFL number for this case. The SBC method effectively stabilizes the iterative process, leading to superior convergence performance.

Case 3
We consider a gas-condensate mixture with the fluid and compositions given in Table  9. Quadratic relative permeabilities are used. Initial, production, injection pressures are  (0) 85, 80, 120 bars, respectively. Temperature is 325 K.
2.5 0.5 0.98 C 2 1.5 0.4 0.01 C 5 0.05 0.1 0.01 The nonlinear performance of the gas-condensate case is summarized in Table 10. The maxCFL is 54. The gas saturation and overall composition profiles are plotted in Fig.  9. The gas-condensate mixture forms a sharp front in the domain. For min = 1.0e-2, small errors are produced around the shock. Note that the solution difference near the production end between SBC and the standard formulation is due to different timestepping schedules.  We further consider a case with variable K-values. Initial pressure is 100 bars, at a temperature of 340 K. Production and injection pressures are 95 and 140 bars, respectively. Initial compositions are {C 1 (1%), C 2 (50%), C 3 (49%)}, and injection mixture is The value of = 1.0e-2 is used for the smoothing parameter. The nonlinear performance of Case 4 is summarized in Table 11. The maxCFL is 24. The standard formulation shows a poor nonlinear performance, while the iteration number is largely reduced under the smooth formulation.  The fluid and compositions are given in Table 12. Initial pressure is 75 bars, at a temperature of 410 K. Production pressure is 70 bars, with injection pressure as 140 bars. The total simulation time is 400 days. The nonlinear performance of Case 1 is summarized in Table 13. The maxCFL is 17. The case is to validate the applicability of the smooth formulation without the continuation method. We can see that the value of = 1.0e-2 brings a reduction in the iteration number, though both the formulations perform well for this case.
The gas saturation and overall composition profiles are plotted in Fig. 10. Note that the two formulations present the same solution profiles.

Case 2
We test a gas-condensate system with initial compositions as {80%, 1%, 14%, 5%} at an initial pressure of 100 bars and at a temperature 344 K. Production pressure is 95  bars. Injection pressure is 140 bars, and injection gas mixture is {1%, 97%, 1%, 1%}. The total simulation time is 500 days. The nonlinear performance of Case 2 is summarized in Table 13. The maxCFL is 60. A significant reduction in the iteration number is achieved by the smooth formulation. Permeability field of the bottom layer of the SPE 10 model is shown in Fig. 11. A uniform cell size 10 m is specified, and the porosity is 0.1. Positions of the producer and injector are (60, 1) and (1, 220), respectively. Relative permeabilities are quadratic for both phases. The fluid and compositions are the same as summarized in Table 6. Initial pressure is 70 bars, at a temperature of 320 K. Production and injection pressures are 65 and 100 bars, respectively. The total simulation time is 400 days.
The nonlinear performance of Case 1 is summarized in Table 15. The maxCFL of the case is 88. The heterogeneous model is very challenging because of a large variation in the CFL numbers across the domain. Also frequent phase changes make the fluid displacement process highly nonlinear. As can be seen, the standard formulation suffers from many timestep cuts and wasted iterations. By comparison, the SBC method shows a smooth behavior, resulting in much improved nonlinear convergence.  (0) 20 (0) The gas saturation and overall composition profiles are plotted in Fig. 12 and Fig.  13. The solutions from SBC closely matches the standard formulation, even for the level of min = 1.0e-2.

Case 2
We also study a three-component system with {K 1 = 2.5, K 4 = 1.5, K 10 = 0.3}, and the other parameters remain unchanged.  The nonlinear performance of Case 2 is summarized in Table 16. The maxCFL is 87. The iteration performance of the standard formulation becomes worse for the different K-values. Crossing phase boundaries along a single tie-line causes oscillations in the Newton solver. For multi-component systems, additional discontinuities can arise with the switches between key tie-lines. The SBC method provides smooth transitions across the entire phase boundaries, leading to superior global convergence performance.
The nonlinear performance of Case 3 is summarized in Table 17. The maxCFL is 520. The SBC method exhibits favorable nonlinear convergence for this challenging case.  (0) 13 (0) The overall composition profiles are plotted in Fig. 14. Here the timestep size is reduced to 10 days, ensuring that the solutions from the different methods are compared with the same time-stepping schedule.

Summary
Frequent phase changes and oscillations around phase boundaries can cause severe convergence problems during compositional simulations. The objective of this work is to develop a smooth formulation that removes all the property switches and discontinuities associated with phase changes. Here we first reformulate the coupled system, so that the discontinuities are transferred to the phase equilibrium model. A single and concise non-smooth equation is achieved and then a smoothing approximation can be made. The reformulation is based on a mixed complementarity problem (MCP) proposed for the phase equilibrium. The developed formulation with a smoothing parameter provides smooth transitions of variables across all the phase regimes. We also employ a continuation method with the smooth system as a homotopy mapping.
We evaluate the efficiency of the new smooth formulation and the continuation method using several complex problems. For most of the cases, the standard natural formulation suffers from multiple timestep cuts and subsequent wasted Newton iterations. In contrast, the developed formulation and method exhibit superior global convergence, requiring no timestep cut. The SBC method shows a negligible impact on solution accuracy, while providing smooth and stable iterative performance. Moreover, SBC works robustly for a wide range of flow conditions without parameter tuning. In this work the fugacity constraint is evaluated through K-values. The smooth formulation can be readily applied to the black-oil fluid model (with dissolved gas and vaporized oil). The assumption does not limit the applicability of the formulation to a general EoS-based system. x c,l = 0.
The minimum function represents the nonlinear complementarity function. The above formulation may produce the solutions with negative molar fractions, especially if smoothing approximations are applied. This limitation can cause severe numerical issues for compositional problems of practical interest.
In Lauser et al. (2011), a semismooth Newton method with locally superlinear convergence is employed to solve the non-smooth system. However, the semismooth method may exhibit poor global convergence: the iterations may diverge, when the starting point is not close enough to a solution.

Appendix B. Discretization methods
A standard finite-volume scheme is applied as the spatial discretization for the mass conservation equations. A two-point flux approximation (TFPA) is used to approximate the flux across a cell interface. The method of choice for the time discretization is the fully-implicit scheme. The discrete form of conservation equation is given as, where superscripts denote timesteps, and ∆t is the timestep size. V is the cell volume. All indexes related to the cell numeration are neglected. The accumulation term involves the total density, ρ T z c = x c ρ o s o + y c ρ g s g and ρ T = s o ρ o (x) + s g ρ g (y), two phase, ρ l (z), one phase.
where ρ o (x) indicates the density computed at a composition x, and ρ l (z) is the density computed in the single-phase regime at a composition z. The discrete phase flux across the interface (ij) between two cells is written as, where subscript (ij) denotes quantities defined at the cell interface. Υ ij is the interface transmissibility. ∆Φ l,ij = ∆p ij − g l,ij is the phase potential difference with the discrete weights g l,ij = ρ l,ij g∆h ij . The phase and compositional coefficients associated with the flux terms are evaluated using the Phase-Potential Upwinding (PPU) scheme.