Compositional Reservoir Simulation: A Review

Compositional Reservoir Simulation:
A Review
Larry C Young1*‍‍
1
Tilden Technologies LLC
Summary
Compositional reservoir models have many uses, including petroleum reservoir simulation, the simulation of flow and remediation of
nonaqueous petroleum liquids, and the study of various carbon-capture technologies. The present work describes these models in a
Downloaded from http://onepetro.org/SJ/article-pdf/doi/10.2118/208610-PA/2672568/spe-208610-pa.pdf/1 by Universidad Nacional Autonoma De Mexico user on 22 May 2022
general way but concentrates on isothermal compositional models based on an equation of state (EOS). These models are used to study
petroleum reservoirs containing gas condensate or volatile oil fluids and various gas injection processes like CO2 injection. The technol-
ogy is documented from the earliest rudimentary models which first appeared more than 50 years ago to more recent work. The system
of differential-algebraic equations is described along with the physical nature of the parameters. Model formulations are compared with
respect to choices of variables and other computational features. Descriptions of efficient and reliable computational procedures are given.
In addition to basic algorithms, we detail important side issues of (1) testing for phase appearance, (2) temporal stability, (3) material
balance, and (4) other efforts to speed up the calculations and/or reduce the calculation requirements. Information in the literature is
supplemented with some work which has not been presented previously. Correct rendering of the physics and chemistry is emphasized
which is often intertwined with the numerics. Although many improvements have been made, the work is not finished. We conclude with
a discussion of areas needing additional work.
Introduction
The usage of compositional models for flow in porous media has increased significantly since the first models began to emerge
50 years ago. The driving force behind this increase is: (1) faster computers with large memory, (2) improved algorithms and
implementation, and (3) greater need. In the petroleum industry, the greater need stems from deeper drilling and the greater
incidence of volatile oil and gas condensate discoveries, together with increased use of gas injection (e.g., CO2, N2, or enriched
hydrocarbon gases) for improved or enhanced recovery. The greater need also includes the desire to clean up the environment
from nonaqueous petroleum liquids contaminating our groundwater and removal of greenhouse gases by sequestration. Faster
computing hardware has received the greatest recognition; however, improved algorithms and implementation have been equally
important.
The earliest models used equilibrium K-value tables, inconsistent correlations for phase properties, time stepping with explicit
mobility terms, and non-Newton iteration methods for the resulting nonlinear problem. All of these items caused reliability prob-
lems. Later models replaced tabular K-values with a consistent cubic EOS. Full Newton-Raphson methods improved nonlinear
convergence and judicious variable selection improved on early formulations. Most of the early formulations used explicit
mobility terms, with the well-known limitation on timestep size. Several fully implicit formulations have been described, and
their use has increased along with computer speed and memory size. Adaptive implicit methods (AIM) strike a balance between
these two methods and are widely used in many commercial simulators. Temporal instability and running speed have not been
the only issues with an EOS. Some other examples are: (1) detecting the appearance of a phase, (2) issues when three hydrocar-
bon phases exist, and (3) relative permeability consistency. These and other more obscure problems are discussed toward the end
of the article.
The oil embargo of 1973 and the long lines at gas stations created an increased interest in enhanced oil recovery as part of a
drive for greater energy independence in the US. As a result, there was a need to accurately simulate enhanced recovery pro-
cesses. These methods use a displacing fluid composed of a relatively small amount of expensive solvent, e.g., CO2, enriched
hydrocarbon gases, or chemical surfactants. Accurate modeling of these problems is difficult. This was the environment when I
arrived at the Amoco Research Department in 1975. My assignment was to investigate improved discretization techniques for
simulating these problems. My colleague, Del Fussell, was developing a compositional simulator. He had given up on a simulator
developed in a joint industry-sponsored project. It was based on tabular K-values. Instead, he chose to use the Redlich-Kwong
EOS, because Lyman Yarborough had achieved excellent results correlating fluid properties with that EOS (Yarborough 1979).
The importance of compositional phenomena has only grown since then. In addition to gas injection for improved recovery, deeper
drilling has increased the frequency with which exotic fluids are encountered. Starting in the 1960s, Amoco devoted a major effort to CO2
flooding, and it was the only company with the foresight to plan oil-in-the-tank pilots rather than the observation and sampling wells used
by other companies. During my final years at Amoco, I was involved on the periphery of the Anschutz Ranch reservoir study. The fluid
was a very rich gas condensate, near its critical point, and with a substantial variation of composition with depth. Had the reservoir been
a few degrees cooler, it would have been an oil reservoir. It became the largest N2 injection project at that time, and it was simulated with
a compositional model which was substantially larger than any previous ones (Wendschlag et al. 1983).
Wong and Aziz (1989) have previously reviewed a wide range of models with compositional effects, including thermal, chem-
ical, as well as isothermal EOS models. Thermal-compositional models are important for simulating thermal processes like
combustion and steam injection for improved oil recovery or for remediation of contaminated soils. Other articles have compared
*Corresponding author; email: larry@tildentechnologies.com

Copyright © 2022 Society of Petroleum Engineers
Original SPE manuscript received for review 14 July 2021. Revised manuscript received for review 28 September 2021. Paper (SPE 208610) peer approved 29 September 2021. This
paper is published as part of the 2021 SPE Reservoir Simulation Conference Special Issue. Supplementary materials are available in support of this paper and have been published
online under Supplementary Data at https://doi.org/10.2118/208610-PA. SPE is not responsible for the content or functionality of supplementary materials supplied by the authors.
2022 SPE Journal 1

various aspects of compositional model formulations (Thele et al. 1983; Wong et al. 1990; Coats 2000; Cao 2002; Voskov and
Tchelepi 2012; Schmall et al. 2013).
The current article considers only isothermal models with an emphasis on those based on an EOS. These models are used for volatile oil and gas
condensate reservoirs and for simulating gas injection processes in oil reservoirs. Although thermal simulation is not directly addressed, much of
the methodology carries over to thermal processes. Restricting the type of models considered, allows for a more in-depth analysis than in previous
reviews.
This article is divided into several parts. First, the governing equations and their basic nature are described. Next, various algorithms used to
solve the equations are outlined. These are subdivided into early models which used methods and algorithms that are now considered obsolete and
those which overcome deficiencies of the early models. The later models invariably use a set of primary and secondary variables, but different
variable choices are in use. These choices are analyzed in some detail. Next, several other computational aspects are considered, such as material
balance, detecting phase appearance, temporal stability, alternative fluid correlations, efficient and reliable implementation of EOS calculations, and
high-performance computing. Over the years, there have been dramatic improvements in algorithms, implementation, and available computational
power; however, several problems remain. Last, some of these remaining problems are discussed, including (1) three hydrocarbon phase issues, (2)
consistency and accuracy of relative permeability, and (3) dispersive mixing effects.
Disclaimer. The author has worked in this field for over 40 years. In addition to the technology, the paper covers my own personal journey
in this field of study. The article heavily references my own work because it is most familiar to me. I have made every effort to integrate
it with the relevant work of others. Only key references are cited for many topics, which gives a starting point for future researchers. The
phrase and references cited therein is implied rather than continually repeated. I apologize in advance for any work I have inadvertently
overlooked.
Nature of the Problem

Compositional models are governed by differential material balances together with thermodynamic relationships governing fluid proper-
ties, creating a large and highly nonlinear differential-algebraic system of equations. Equilibrium phase behavior is usually assumed
although this assumption has been the subject of debate.
Differential Material Balance. Note on Notation. When describing the governing equations, we use the convention that two subscripts
indicate a value for a component within a fluid phase. A single subscript may indicate either an overall component value or the value for
a fluid phase. To distinguish between the two cases, a Greek subscript, α, β, and so on, is used to indicate a phase value, while a letter,
i, j, k, and so on, indicates a component value. An overall quantity for the multiphase-multicomponent mixture uses no subscript. When
referring to specific phases, the subscripts, ℓ, v, and a are used to indicate hydrocarbon liquid, vapor, and aqueous phases, respectively. At
times, the subscript w is used for the water component (numbered last) and h is used to indicate an overall quantity for the hydrocarbon
phases (i.e., excluding water).
With the assumption of local equilibrium, the governing equations are:
Np h
P i
@
@t
S˛ x˛i˛ + r u˛ x˛i˛ r D˛ r x˛i˛ = 0, (1)
‍ ˛=1 ‍
for i = 1,…,Nc. The equations are valid in any units for the quantity of a component (e.g., mass, moles, or standard volumes). Mole units
are used in this discussion. The phase velocities are governed by Darcy’s law extended to multiphase flow

u˛ = k˛ rp + rp˛ ˛ grH ,‍
˛
(2)
‍
where ‍p + p˛‍is the phase pressure with p a characteristic pressure function created from a combination of the individual phase pressures.
We leave the definition of p open for now, but if it is chosen to be an individual phase pressure, then the pα are capillary pressure functions.
The phase permeabilities and capillary pressures are normally treated as functions of saturation, although a dependence on interfacial
tension (IFT) or capillary number is sometimes included. The normal correlations are like:

k˛ = kkr˛ S
(3)
‍ p˛ = p˛ S , ‍
where S is the vector of phase saturations and the k with no subscript is the absolute permeability. k varies with spatial position
and should be represented as a general tensor, but in most cases, the cross terms are neglected. The relative permeability and
capillary pressure functions may also vary spatially. They also depend on the history of the displacement (Killough 1976). The
pressure dependence of these parameters is important in some situations, especially when coupled with rock mechanics models
to account for compaction and similar phenomena. Toward the end of the article, we discuss issues with the usual representations
of these parameters.
The dispersion term in Eq. 1 is composed of a small molecular diffusion term and a dispersion term to represent mechanical mixing.
We neglect molecular diffusion. The dispersion term is known to be a tensor function of velocity, with differences in the longitudinal and
transverse directions, like:
ˇ ˇ uˇ˛ uT˛
D ˛ = ˛O t ˇu˛ ˇ I + ˛O ` ˛O t ˇ
ˇu˛ ˇ , (4)
‍ ‍
where ‍u˛ uT˛ ‍is the velocity outer product (Bear 2013). With few exceptions, the dispersion terms are ignored; however, their neglect can
have significant consequences. They are included here for subsequent discussion in the section “Numerical and Physical Dispersion”
toward the end of the article.
Assuming local equilibrium, the differential equations are supplemented by phase or VLE (vapor/liquid equilibrium) relationships
(Prausnitz and Lichtenthaler 2004):
2 2022 SPE Journal

' (x , p)x˛i = 'ˇ (xˇ , p)xˇ i or
˛ ˛
'˛ ˛ˇ (5)
‍ 'ˇ x˛i = Ki x˛i = xˇ i , ‍
for α ≠ β and i = 1,…, Nc, where xα is the vector of fractions within phase α. The fugacity coefficients should be evaluated at the individual
phase pressure for tight reservoirs such as shales where the pores are extremely small (Rezaveisi et al. 2015, 2018; Sandoval et al. 2016).
The phase relationships are correlated by an EOS . Early
models
used
tabular K values. An EOS or correlation also provides values for
the phase specific volumes and viscosities, ‍ ˛ x˛ , p and ˛ x˛ , p ‍.
In most cases, these general equations can be simplified. When simulating petroleum recovery, the solubility of the hydrocarbons in
water is usually neglected. However, in water resource modeling, it could be important when the water is intended for human consump-
tion. The partitioning of water into the vapor phase can be important for hot reservoirs. The water solubility of nonhydrocarbon gases,
such as CO2, are often important and included for simulating processes like CO2 flooding. In chemical surfactant flooding processes,
partitioning between water, hydrocarbon, and microemulsion phases is highly important and must be included in a simulator. However,
these models typically use no vapor phase.
Eq. 1 is written for any number of flowing phases. In some processes, partitioning between the fluid and the rock is important (e.g., surfactant
flooding, ion exchange, adsorption). Eq. 1 can be applied
to these problems by treating the rock as an immobile phase (i.e., u = 0).
The problem can be thought of as having ‍Np Nc + 1 + 1‍ unknowns, x, S, and p. Eqs. 1 and 5 together provide NpNc equations. The
remainder are provided by constraints. There is one constraint on the saturations:
Np
P
S˛ = 1, (6)
‍ ˛=1 ‍
and Np on the fractions within each phase:
P
Nc
x˛i = 1 for ˛ = 1, : : : , Np . (7)
‍ i=1 ‍
Because there are a large variety of possibilities, for illustrative purposes, we will use the common case of two hydrocarbon phases and
water, with partitioning among the components in the hydrocarbon phases, but with an aqueous phase composed of pure water which does
not partition. In this case, there are 3Nc + 4 unknowns. Nc + 1 of the phase equilibria relationships, Eq. 5, are trivial, because by assump-
tion, the aqueous phase is pure inert water. The constraint, Eq. 7, on the composition of the aqueous phase is also trivial. Eliminating these
trivial equations leaves Nc material balance equations, Nc − 1 phase equilibria relationships and three constraints, which matches the 2Nc +
2 unknowns of p, S, and mole fractions x. The mole fractions number 2(Nc − 1) since water is excluded.
There are many ways of solving these equations, which is the subject of this article. In some problems, the volumetric properties of the
fluid are relatively simple, and consequently, the equations are relatively easy to solve. For example, for chemical surfactant processes,
one can usually either neglect the composition dependence of the fluid specific volume or assume additive volumes or ideal mixing.
However, for our example case with liquid and vapor hydrocarbon phases, the volumetric properties are more complex. Before proceed-
ing with the review, we wish to further elucidate the nature of the problem.
Acs et al. (1985) pointed out what should be obvious. Regardless of how the phase behavior and volumetric properties are correlated,
the overall specific volume of a fluid can be characterized by multiphase partial volumes:
" #1
Np
P P
Nc
S˛
= ˛ = i xi , (8)
‍ ˛=1 i=1 ‍
where υ (with no subscript) is the overall specific volume of the multiphase, multicomponent mixture and the υi are partial volumes given
by:

@ n
i = . (9)
‍ @ni ‍
The partial volumes describe the change in specific volume or density when one component is added to the mixture while pressure, tem-
perature, and the amount of all other components remain fixed. These multiphase partial volumes differ from the single-phase quantities
discussed in most thermodynamic texts.
ni and xi (with a single-letter subscript) represent the mole number and mole fraction of a component over all phases. They are of course
related:
Np
P S˛ ni
xi = ˛ x˛i = n
˛=1
Np
P
ni = n˛i (10)
˛=1
P
Nc
n = ni .
‍ i=1 ‍
Computation of the partial volumes and other derivatives of the overall υ is complex because they account for both expansion
of the individual phases and changes in the amount and composition of the phases. Advanced degrees have been granted and
peer-reviewed articles published on calculation of just the pressure derivative alone. Discussion of these derivative calculations
is reserved for later in the article. For now, we simply state that they exist and can be calculated with some effort.
The partial volumes not only describe the compositional dependence of the overall fluid volume but also the individual phase
volumes:
2022 SPE Journal 3

P
Nc
˛ = i x˛i . (11)
‍ i=1 ‍
Several other quantities can be defined in terms of the partial molar volumes. The volume fraction of a component within a phase and an
overall component volume fraction are given by:
i x˛i
S˛i = ˛ and
Np
P (12)
i xi
Si = = S˛ S˛i .
‍ ˛=1 ‍
Using these relationships and definitions, Eq. 1 can be rewritten in terms of overall mole fractions, xi, or volume fractions, Si. The former
mole fraction case is:
Np h i
P
@
@t
xi + r u˛ x˛i˛ r D ˛ r x˛i˛ = 0. (13)
‍ ˛=1 ‍
The later volume fraction form is:
P Np h i
@
@t
Sii + r u˛ S˛ii r D˛ r S˛ii = 0. (14)
‍ ˛=1 ‍
The mole fraction form is the basis for most solution algorithms, whereas the volume fraction form makes the equations analogous to
those for immiscible displacement without component partitioning.
The total fluid volume and overall specific volume can be linearized in terms of these same quantities. For an isothermal
system:
P
Nc
d(n) = i dni + n @
@p
dp or
i=1
PNc
d = i dxi + @
@p
dp (15)
i=1
PNc P Nc
= i dxi + xi di = i dxi + Sii di .
‍ i=1 i=1 ‍
The problem is simplified by differentiation of Eq. 14 with the chain rule and multiplication by υi:
Np h
P i
u˛ S˛i
@
@t
Si Si @
@t
i
+ r u˛ S˛i i ri i r D ˛ r S˛i
i = 0. (16)
‍ i
˛=1 ‍
By virtue of Eq. 15, the summation of these component balance equations gives:
P Np h
Nc P i
u˛ S˛i S˛i
C @p
@t
+r u i ri i r D ˛ r i = 0, (17)
‍ i=1 ˛=1 ‍
where the total compressibility and total velocity are:
1 @ 1 @
C= @p @p and
Np
P (18)
u= u˛ .
‍ ˛=1 ‍
The summation terms on the right of Eq. 17 are zero if υi is constant and are relatively small in other cases. Authors of streamline-
type models discard these terms without justification and without stating the true nature of the assumption made. This equation
is an overall continuity equation, but it is usually called the pressure equation. For the special case of a simple black oil fluid
treatment, Eq. 17 was described long ago, but its extension to a fully compositional fluid is relatively more recent (Martin 1959;
Young 2001).
Eq. 17 makes the nature of the equations clear. The equations are like those for flow of immiscible fluids without phase trans-
fer or miscible flow of a single phase such as the transport of tracers, salt water, or radionuclides in an aquifer or watered-out oil
reservoir. Eq. 17 provides an equation for the velocity field. It is a parabolic equation which can be almost quasistatic or essen-
tially elliptic in nature for the time scales of interest. It is primarily a function of pressure and is coupled to the component bal-
ances, Eq. 14, which are convection dominated and nearly hyperbolic.
Partial Volumes. The partial volumes defined by Eq. 9 are a key property fundamental to the problem. These quantities have
a simple physical interpretation. The value for each component determines the change in fluid volume as a small amount of the
component is added to the mixture. The overall fluid molar volume is just a weighted average as shown in Eq. 8. The partial
volumes often vary in unexpected ways and can even be negative. Although these quantities have been discussed for years
in the literature, examples illustrating their behavior have not been presented. For this reason, Supplementary file 1 contains
examples for three fluid systems: (1) a two-component black oil, (2) the ternary C1–C3–C12 system, and (3) a well-documented
N2-gas condensate characterized with 11 components. Because the effects of pressure are more widely recognized, the examples
concentrate on the effects of composition. A few results from Supplementary file 1 are reproduced here because they are central
4 2022 SPE Journal

Fig. 1—C1–C3–C12 phase diagram and lines.
to the discussion. Fig. 1 is a ternary diagram for the C 1–C3–C12 system calculated at 12.4 MPa (1,800 psi) and 121.1°C (250oF).
It includes the phase envelope and tie-lines along with three lines along which properties are examined. Fig. 2 shows the partial
volumes along Line A in Fig. 1, while Fig. 3 shows them along Lines B and C.
Line A is a tie-line far removed from the critical, where the behavior is basically like that for a black oil. Line B is a tie-line
near the critical, while Line C is in the dense phase region but runs very near the critical point. Along Line B, the discontinuities
at the phase boundaries are smaller than along Line A. The overall variations (maximum to minimum values) along all the lines
are similar, but away from the critical, the changes occur abruptly at the phase boundaries, whereas near the critical, the changes
Fig. 2—C1–C3–C12 partial volumes, Line A.
2022 SPE Journal 5

Fig. 3—C1–C3–C12 partial volumes, Lines B and C.
are greater within the single-phase and two-phase regions. Except for the small two-phase region along Line B, the differences
with Line C are minor.
Along Lines B and C, the partial volume of C12 goes to −0.5 m3/(kg-mole) when present in small amounts. With a negative partial
volume, the overall specific volume of the mixture will decrease when a small amount of that component is added. This anomalous type
of behavior has been noted previously for steam in thermal simulations (Coats 2000), but these results show it can also occur for hydro-
carbons in isothermal problems.
The main point of these calculations is that all the volume derivatives vary smoothly except near a critical point and at phase bound-
aries where the derivatives are discontinuous. The behavior of the total fluid compressibility is well known, but here we see that partial
volumes also have jump discontinuous by a factor of 5 or more. They can vary from one component to another by more than a factor of
10 with either a molar or mass basis. They can even be negative. Of course, the properties (e.g., ‍x˛ , ˛ , ˛)‍ for one phase are undefined
on one side of a phase boundary. This fact has important consequences for discretization and solution of these equations. The accuracy
and convergence rates of numerical approximations are always limited by the continuity of the coefficients and the solution. Unless the
limited continuity is somehow accounted for, one can expect no better than a first-order convergence rate regardless of a method’s order
for smooth problems.
The discontinuities also impact convergence rates of nonlinear algebraic equations. If the problem is solved by a Newton-Raphson
procedure, the elements of the Jacobian are discontinuous at a phase boundary, slowing the convergence rate. For simulations of repres-
suring by water injection, often one Newton iteration per timestep is required except when a phase transition occurs somewhere within
the grid, then invariably a second iteration is required. When some simplified Newton methods are used, the equations can fail to converge
when crossing a phase boundary. These discontinuities are not a new discovery, but they are often forgotten.
Flow Terms. To take a closer look at the flow terms, the following are defined using the same subscripting convention. The component
mobilities are defined in terms of volume fractions, S, rather than mole fractions, so they will be like saturation in a strictly immiscible
problem. First, the mobility terms are define:
kr˛
˛ = ˛
˛i = ˛ S˛i
Np
P
i = ˛i (19)
˛=1
PNp P
Nc
= ˛ = i ,
‍ ˛=1 i=1 ‍
and then the fractional flow functions:
f˛ = ˛ /
f˛i = ˛i / = f˛ S˛i . (20)
‍ fi = i / ‍
In terms of these quantities, the overall velocity is:
6 2022 SPE Journal

!
Np
P
u = k rp + f˛ rp˛ grH
N , (21)
‍ ˛=1 ‍
P
where ‍N = ˛ f˛ ˛.‍ Ideally, the characteristic pressure, p, should be defined such that the capillary pressure term is zero or small in the
overall velocity equation. For two-phase systems without component interchange, it can be made identically zero by use of Chavent’s
global pressure function (Chavent and Jaffré 2014).
Phase velocities are defined using deviations from the average velocity:

‍ u˛ = f˛ u + ıu˛ ,‍ (22)
where:
!
Np
P
ıu˛ = k rp˛ fˇ rpˇ ˛ N grH . (23)
‍ ˇ =1 ‍
These deviation velocities incorporate gravity and capillary forces but have no direct dependence on global pressure. The volumetric form
of the component balances, Eq. 14, can be expressed using these phase velocity deviations as:
P Np h i
@
@t
Sii + r u + ıu˛ f˛ii r D ˛ r S˛ii = 0. (24)
‍ ˛=1 ‍
The problem is usually formulated in terms of mobilities, λ, but formulating it with fractional flow functions has been advocated by some
(Young 1984; Chen and Ewing 1997).
As we shall see, this form shows that the fractional flow functions multiplying the velocities are key to the temporal stability of explicit discret-
ization methods.
Discretization
Discretization of the compositional equations will not be described here in detail. Models in widespread use employ either finite differ-
ence, finite volume, control volume or finite element methods. The terms block, cell, point, or node are used interchangeably in this
discussion.
In compositional models, the dispersion terms are almost universally ignored. Here, simplifications are made to describe the solution
algorithm more clearly. The dispersion terms and velocity deviations, δu, are not displayed in the description, but these terms are easy to
carry along. With these assumptions, the component balances are:

‍
@
@t
xi = r i k rp.‍ (25)
After discretization, these equations take the form:

P
Vj
t O xi = (i )u Tjk (pk pj ), (26)
‍ j k ‍

where Δ indicates the change over a timestep and ‍O = /0.‍ For a slightly compressible formation usually ‍O = 1 + Cr p p0 /0.‍
Throughout the article, we will call the left-hand side of Eq. 26 the capacity terms and the right-hand side the flow terms. The summation
is over all the connections to point j, and u is usually the upstream value, either j or k. This representation is valid for a number of differ-
ence schemes including standard finite difference grids, grids of triangles or tetrahedra (sometimes called Voronoi or PEBI grids), and
quadrilateral or hexahedral grids treated rigorously with 27-point or multipoint flux approximations. The most common discretization for
3Duses seven points or two-point flux approximations, even when applied to distorted hexahedral grids where the approximations are
‍O 1 .‍ The T are called transmissibility coefficients and V are the pore volumes associated with the nodes or blocks. Note that the original
component balances have units of moles or mass per unit volume per unit time (M/V/T), whereas the discretized equations, Eq. 26, have
units of moles or mass per unit of time (M/T). Stock tank or standard volumes could be used because they are directly related to moles
and mass.
For the time scales of interest in reservoir simulation, the pressure terms must be evaluated implicitly. If the flow coefficients, λ, are
evaluated at the old time level, it is said to have explicit flow coefficients or the abbreviation IMPES. If they are evaluated at the new time
level, it is a fully implicit method or FIM.
The introduction of the pore volume array, V, into the discrete form of the equations has caused confusion in the interpretation of var-
ious solution algorithms and how the equations and independent variables are perceived. A case in point is the saturation constraint equa-
tion, Eq. 6. Treatment of this constraint has been viewed differently by different authors. It can be expressed as either:
Np
P Np
P ˛ n˛
S˛ = O
= 1 or
˛=1 ˛=1 V
Np (27)
P P
Nc
O =
V ˛ n˛ = n = i n i .
‍ ˛=1 i=1 ‍
Some view the constraint as one on saturation, while others have called it a volume balance. These are two different ways to perceive the
same equation. We make no distinction because the difference is trivial.
Another example is in the perception of the variables actually being solved for (i.e., the independent variables). The capacity terms on
the left-hand side of Eq. 26 are:
2022 SPE Journal 7

V O xi = V
O xi = nxi = ni . (28)
‍ ‍
Some view the independent variables as ‍x or x/ ‍ , while others view them as n. Again, the difference is minor and one of perception.
Many authors make an important distinction between extensive and intensive variables, but here we see the difference is not significant.
It is natural to switch within the steps of the algorithm and some authors do precisely that.
Eq. 26 or some variation of this form of the equations are used by most authors. However, many separate the hydrocarbon balances
from the water balance. For our standard example case where the aqueous phase is pure inert water, the component balances take the form:
8 9
Vj
:O Sa ; = P (a )u Tjk (pk pj ) and
t a
P
k
(29)
Vj O 1 Sa zi =
t h j i u Tjk pk pj .
‍ k ‍
Occasionally the hydrocarbon equations are summed to create a hydrocarbon overall balance:
P
Vj O 1Sa
t h j = h u Tjk pk pj . (30)
‍ k ‍
Let water be the last component, so for the above i = 1,…, Nc − 1. zi and υh are properties of the hydrocarbon phases defined on a water-
free basis:
" #1
P S˛
NP
c 1
h = ˛ = i zi . (31)
‍ ˛=`,v i=1 ‍
Also define the hydrocarbon vapor mole fraction:
h Sv nv nv
v= = = nh . (32)
‍ v nnw ‍
The relationships for the hydrocarbon overall fractions, zi, are:
ni ni
zi = (nnw ) = nh
P S˛
= h ˛ x˛i (33)
˛=`,

‍ = 1 v xì + vxvi .‍
The equations for phase equilibria and the constraints on phase fractions, Eqs. 5 and 7, can be combined to give the Rachford-Rice equa-
tion (Rachford and Rice 1952):
NPc 1

Ki 1 zi
= 0, (34)
1+ Ki 1 v
‍ i=1 ‍
where K are the equilibrium K-values, Eq. 5.
Review of Formulations
Early Models. The first full flow models with component interchange between phases were of the black oil type with its simplified two
hydrocarbon component treatment. These models used independent variables of pressure and phase saturations. They were formulated in
terms of saturation, pressure, formation volume factors, and solution gas/oil ratio, Rs. Some call these natural variables. However, what
is natural to one person may be unnatural to another. The variables are natural in the sense these variables are used in material balance
calculations and are familiar to most petroleum engineers. Density and equilibrium K-values are more familiar (or natural) to other
engineers.
Most early black oil models used explicit flow coefficients in Eq. 26. However, even with this assumption, there remain time deriva-
tives of saturation that couple the equations. In an article describing effective variables for multiphase pressure transient analysis, (Martin
1959) used a combination of the equations which eliminates the saturation derivatives. At the SPE Annual Technical Conference in 1960
(Sheldon et al. 1960; Stone and Garder 1961) described black oil models using the same elimination procedure. Sheldon et al. state the
procedure is “implicit with regard to pressure and explicit with regard to saturation.” This designation was later shortened to the acronym
IMPES. The same IMPES procedure was used in many early simulators (Fagin and Stewart 1966; Breitenbach and Thurnau 1968;
Sheffield 1969). The pressure equation derived by the IMPES method is usually solved by an iterative method, then two saturations are
updated with the individual component balance equations, and the last saturation is calculated by difference, Eq. 6.
The weighting factors needed to create the equation for pressure are the partial volumes of Eq. 9 although they were not recognized as
such at the time. The conserved units are standard volumes rather than moles, which differ by simple conversion factors. The basic for-
mulation was refined to treat problems with repressuring. The relationships for repressuring are given in the waterflooding monograph
(Craig 1971). This is called a variable bubble point treatment (Thomas et al. 1976). When the pressure exceeds the bubble point pressure,
the vapor phase saturation becomes zero, so it is no longer a valid independent variable. In that case, the independent variables are pres-
sure, saturation, and solution gas/oil ratio, Rs, which is equivalent to solving for hydrocarbon liquid phase composition. This switch is
called variable substitution. One may consider this either natural or unnatural, depending on one’s point of view.
FIM was described by Blair and Weinaug (1969). This method was the second most common solution method. It has become more
popular along with improvements in computer speed and memory. When the flow coefficients, that is, the λi in Eq. 26 are treated explicitly,
there is the well-known CFL (Courant-Friedrichs-Lewy) limitation on timestep size for difference methods. For a simple two-phase
immiscible system with no phase transfer and negligible capillary pressure, the limitation is:
8 2022 SPE Journal

Q df
t Vjj dS j = Cj , (35)
‍ ‍
where ‍C ‍< 2 for stability and ‍C ‍< 1 insures no overshoot. f is the fractional flow function, Qj is the volumetric flow out of grid cell j, and
Vj is the volume associated with the cell. For an immiscible system, f has the usual S shape. For a strictly miscible displacement (also
called first contact miscible or FCM) with two components in a single-phase mixture, the analogy between miscible and immiscible dis-
placements gives ‍df/dS = 1.‍ Because the derivative is constant for a miscible displacement, the fundamental nature of when and where
instability can occur is different for miscible and immiscible systems. The limitation is basically a throughput condition, so it is most
severe when the flow is high, or the grid is fine. Stability issues for compositional models are discussed in the section “Temporal Stability
and Adaptive Implicit Methods.” FIM circumvents the timestep limitation by treating the flow terms implicitly. Although FIM is uncon-
ditionally stable, timesteps are still limited by nonlinear convergence issues. In the early days, owing to the greater computations, FIM
was reserved for specialized simulators used for near wellbore calculations with cylindrical coordinate grids. These were often called
coning models.
The earliest true compositional models, with more complex phase behavior, considered the calculations for only a single cell, that is,
a tank model (Jacoby and Berry 1957). These were sometimes combined to mimic displacement processes. For example, after the phase
equilibria calculations for one cell, the gas could be mixed with fresh liquid, phase split calculation performed, this new gas mixed again
with fresh oil and so forth. In some cases, relative permeability is used to move the liquid and gas forward, approximating a 1D displace-
ment at constant pressure. These are called cell-to-cell models (Attra 1961; Fussell et al. 1976).
In the late 1960s, several articles describe solution of the full set of compositional flow equations (Kniazeff and Naville 1965; Price
and Donohue 1967; Roebuck et al. 1968, 1969a, 1969b, 1969c; Wattenbarger 1970; Van-Quy et al. 1972). Often, simplifications were
made (e.g., 1D flow, three components, simplified phase behavior). For more complex multidimensional problems, the solution methods
were heavily influenced by the IMPES method for black oil simulation, so saturation was chosen as a primary variable.
In a compositional context, the term IMPES is not fully descriptive and is confusing, because there is more to the calculation than just
saturation. Some have used acronyms IMPEM and IMPEC, but that terminology gets cumbersome. These designations do not distinguish
between methods that fully couple the phase equilibria and those that do not.
With the computer speed and memory then available, compositional problems were too intractable to be tackled by FIM, so a simpli-
fied algorithm like IMPES was needed. In a compositional setting, creation of an equation with pressure as the only variable is signifi-
cantly more difficult. The simple derivation of Eq. 17 above was not apparent for some time. Also, calculation of the partial volumes is
far more complex than the simple factors for a black oil representation.
The early compositional simulators used a pressure equation, but it was not derived in a rigorous way. Some of these early articles state
that a pressure equation is used, but they do not describe how it is derived. A good description of the approach is given by (Kazemi et al.
1978). These formulations were called IMPES, but the time derivative or capacity terms, Eq. 26, were not fully linearized unlike the
rigorous linearization used for an IMPES black oil model. The underlying assumption in the Kazemi et al. formulation is that all partial
volumes are equal (Watts 1986). The partial volumes illustrated in the Supplementary file 1 show this is not a good assumption. The
assumption of equal partial volumes is in addition to the explicit treatment of flow coefficients. For this reason, the term IMPES is not a
good description for compositional models, but its usage is so thoroughly ingrained that we will continue using it here. The reader should
keep in mind that the term IMPES means the mobility terms are treated explicitly.
These formulations were relatively easy to implement. The steps in the calculation for each timestep are like:
1. Calculate properties
2. Construct approximate pressure equation
3. Solve for pressure
4. Update compositions
5. Flash the mixture to get new properties
6. Return to Step 1 if not converged
The term Flash is commonly used to denote the VLE calculations, but Step 5 is more involved than it appears. Good initial estimates are
needed. At a minimum, one would save the K-values for updating. One must determine whether multiple hydrocarbon phases exist and if
so, how the components partition.
All of the early articles used tabular equilibrium K-values to describe the phase behavior. Empirical relationships were used for den-
sity or specific volume, separate correlations for liquid and vapor. (Nolen 1973) points out many pitfalls with some of the early fluid
descriptions. For example, there must be consistent specific volume calculations for the hydrocarbon phases. A fluid may go from liquid
to vapor by passing through the dense phase region, never splitting into two phases. If different correlations are used for vapor and liq-
uid, a discontinuity occurs when switching correlations. Nolen describes the result as like an explosion occurring within the simulation
grid.
In summary, these early simulators suffered from three problems:
1. Inconsistent properties
2. Loosely coupled fluid properties
3. Explicit flow coefficients
The solution to the first problem is to use one consistent EOS for phase equilibria and the volumetric properties of both liquid and vapor.
Advances in the development of cubic equations of state for high pressure systems provided a viable alternative to the earlier fluid repre-
sentations. The Redlich-Kwong EOS (Redlich and Kwong 1949) was used more initially, but the Peng-Robinson EOS (Peng and Robinson
1976) is now the most popular. Liquid density is improved with a volume translation parameter (Péneloux et al. 1982; Jhaveri and
Youngren 1988). This article is not specialized to a specific EOS because many have been proposed. The same correlation is used for both
liquid and vapor viscosity (Lohrenz et al. 1964). Nghiem et al. (1981) describe an EOS model with consistent fluid properties, based on
a simplified pressure equation, so it overcomes the first problem.
The second problem listed leads to inefficiency, slow convergence, and in extreme cases, a failure to converge [see Mansoori’s
discussion of Nghiem et al. (1981)]. The number of iterations required for the simplified treatment is typically three to eight per
timestep (Kazemi et al. 1978; Nghiem et al. 1981). It is important to distinguish between methods with fully coupled fluid properties
or a rigorously derived pressure equation and those which use a simplified approach. (Young and Stephenson 1983) called the sim-
plified methods non-Newton-Raphson, because the pressure equation is not the result of a full linearization of the capacity terms, Eq.
26. It could also be called a simplified or approximate Newton method. Others have used the terms fully coupled and not fully
coupled.
2022 SPE Journal 9

Core Developments. Here we describe solution methods which overcome the major deficiencies of the early methods. Most simulators
in widespread use employ a similar formulation to one of these. The organization of this section is not strictly chronological, but rather a
combination of topical and chronological.
The governing equations form a differential-algebraic system. They consist of Nc differential material balances governing flow and
transport, Eq. 26, together with several nonlinear algebraic constraints, Eqs. 5–7 governing the physics. The transport equations involve
pressure and mobility terms at neighboring cells, whereas the constraints involve only variables within a single cell. Solution methods use
this problem structure in various ways.
The descriptions again use the special case of pure inert water phase together with several other components which partition between
two hydrocarbon phases. This case was considered in most articles. To make the description explicit, the three phases are designated as ℓ,
v, and a for hydrocarbon liquid, vapor, and aqueous, respectively. The water component is numbered last and designated by subscript w.
Because the total number of components is denoted by Nc, the hydrocarbon phase components are the first Nc − 1. This differs from many
authors who use a similar symbol to denote only hydrocarbon components. We assume an EOS is used where fugacities are like Eq. 5 and
the composition dependence cannot be simplified.
When I arrived at Amoco in 1975, my colleagues, Del and Lynne Fussell and others, were already working with a version of the
Redlich-Kwong EOS developed by (Yarborough 1979). This work led to the first model which circumvented the first two limitations
described above. Two years later, Keith Coats published the first EOS model based on FIM.
Fussell and Fussell (1979) describe the first EOS model and the first model with a rigorous linearization of the capacity terms, Eq. 26.
It is an IMPES method, that is, mobility terms are treated explicitly. The algorithm is called minimum variable Newton Raphson, an
extension of an earlier one for flash calculations (Fussell and Yanosik 1978). It departs from the usual procedure of using saturations as
independent variables. Instead, they used moles of each phase, phase mole fractions and pressure as independent variables (Table 1). This
choice simplifies linearization of the capacity terms in Eq. 1. When saturation is an independent variable, linearization of the densities or
specific volumes in the capacity terms create dependencies which do not occur explicitly in their formulation. This feature is discussed
further below in conjunction with Eq. 40.
Component balances, Eq. 29 ‍nw , n` , xì , i = 1, ..., Nc 2‍

Phase equilibria, Eq. 5 ‍nv , xvi , i = 1, ..., Nc 2‍
Saturation constraint, Eq. 27 ‍p‍
Table 1—Equation arrangement used by Fussell and Fussell

(1979).
Fussell and Fussell organize the equations differently depending on whether the hydrocarbon is predominantly liquid or vapor. When it is pre-
dominantly liquid, the ordering shown in Table 1 is used. When the fluid is predominantly vapor, the variables are switched around to avoid ill
conditioning. Another drawback to their method is that because the flow equations are numbered first, the flow coefficients fill in during elimination
of variables from the phase equilibria and saturation constraint equations, increasing both calculation and storage requirements.
When the equations are rigorously linearized, all of the variables are updated in the back-substitution phase of the calculation. For the earlier
simplified solution methods, the values are updated using a flash calculation, which is far less efficient because it does not use prior information.
The method of Fussell and Fussell achieved greater efficiency owing to both more accurate updating and fewer iterations. This 2D simulator was
used in several studies at Amoco, but was restricted to problems of a few hundred cells owing to available computer speed and memory.
Coats (1980) describes an implementation of FIM with a cubic EOS. The equations were rigorously linearized and solved by a Newton-
Raphson iteration. The order of equations and independent variables is summarized in Table 2. Phase saturations, phase compositions, and pressure
are used as primary variables. Because phase saturations were selected as independent variables, the composition dependence of specific volume
or density creates a dense matrix for the capacity terms at each node (see discussion below Eq. 40). However, because FIM was used, this is less of
a problem than in an IMPES formulation. Also, variable substitution is required when a vapor phase is not present at a given cell.
Phase equilibria, Eqs. 5 and 7 ‍xì , i = 1, ..., Nc 1, xv1 , xv2‍

Saturation constraint, Eq. 27 ‍Sa ‍
Component balances, Eq. 29 ‍xvi , i = 3, ..., Nc 1, S` , Sv , p‍
Table 2—Equation arrangement of Coats (1980) and Coats et al. (1998).
At the time of Coats’ paper, FIM was prohibitive for problems of significant size owing to computer memory and speed. For this reason, the
simulator was not actively marketed by Coats’ company (Intercomp). However, they later marketed two versions of a simulator based on IMPES.
Although not published in open literature, those models used the same basic variable set with explicit flow coefficients to rigorously derive a pres-
sure equation. This selection has been called natural variables; however, we will use the abbreviation SXP for this type of formulation.
The variables in the first two rows are aligned with the algebraic equations and are called secondary variables, designated by ‍.‍ Those in the
last row are aligned with the differential equations and are called primary variables, designated by ‍ ‍. There is another level of tertiary variables,
governed by simple relationships, which do not appear because they are directly substituted. The phase equilibria and saturation constraint are lin-
earized in terms of primary and secondary variables. Because the secondary equations are algebraic and do not involve interblock flow, they can be
used to eliminate secondary variables from the primary equations which are differential. In matrix notation, the secondary equations can be treated
as follows:
A + A = Ar or

‍ = A1
A Ar . ‍ (36)
10 2022 SPE Journal

Because there are Nc primary variables and Nc +2 secondary ones, ‍A‍ is a square matrix (Nc +2) × (Nc +2). As discussed below in the
section “Selection of Primary Variables,” further manipulation is required to produce a pressure equation for the IMPES case (Coats
2000). After the secondary variables are eliminated with Eq. 36, the remaining system of equations is solved simultaneously. The number
of simultaneous equations is NcNg for FIM and Ng for IMPES.
A procedure like Eq. 36 is used to eliminate secondary variables in all the procedures. This explanation is conceptual. First, one would
never want to actually invert the matrix, because it is more efficient to achieve the same result by Gaussian elimination or equivalently by
3 ‍A‍into upper and lower components (LU factorization). These calculations are important because they are usually the only ones
factoring
‍ Nc ,‍ at least for the IMPES case.
O
The calculations needed to achieve the reduction of Eq. 36 depends on the structure of these matrices. In turn, their structure depends
on the selection of primary and secondary variables. These choices and consequently the matrix structure vary widely for the different
formulations. In some formulations, the matrices are sparse. There may be submatrices which are diagonal or symmetric. These special
features can be exploited in the implementation. For this reason, the implementation details are discussed below in the sections “Selection
of Primary Variables” and “Choice of Secondary Variables” when discussing the pros and cons of the various variable choices. For now,
the discussion is in conceptual terms.
In the mid-1990s, Coats et al. (1998) developed another compositional model for EOS compositional and black oil fluid representa-
tions. In addition, they could use pressure dependent K-values in place of the EOS. Both FIM and IMPES are implemented. The ordering
of the equations and unknowns are identical to Coats (1980) formulation, see Table 2. They describe some innovative procedures such as
correcting flow reversals during the nonlinear iterations, even with IMPES. Some improved methods for treating wells, such as wellbore
crossflow, are also described. In this article, the term relaxed volume is used to describe Coats’ method for making the equations conser-
vative. This method is discussed below in the section “Material Balance.”
Harper et al. (1985) and Quandalle and Savary (1989) describe formulations using variables like those in Table 2. However, rather than
treating the flow terms by FIM, only relative permeability is treated implicitly (often called IMPSAT). The compositional dependence of
the component mobility terms, Eq. 19, is neglected (see Eq. 42 and related discussion). Similar variables have been selected by others.
Cao and Aziz (2002) considered IMPSAT and an adaptive implicit IMPSAT procedures. Branco and Rodriguez (1996) describe an
IMPSAT-like procedure where the composition dependence of the mobilities is lagged during the iterations. Fernandes et al. (2017) and
several others describe FIM procedures. The efficacy of the IMPSAT approach is discussed in the section “Temporal Stability and
Adaptive Implicit Methods.”
Chien et al. (1985) describe an alternate FIM for a general compositional formulation. The variables selected and corresponding equa-
tions are listed in Table 3. They use the equilibrium K-values as secondary variables. This is a better choice than previous work because
K-values are well defined at phase boundaries and depend only on pressure for systems away from a critical point. The saturation con-
straint is used to eliminate the dependence on one component composition. With a strategy similar to pivoting in linear algebra (Forsythe
and Moler 1967), they select as a key component, the one with the greatest influence on the saturation constraint equation. This key
component is designated as the first one in Table 3, but it could be any of the components. The component balances, including the flow
coefficients, λ, are linearized in terms of the remaining composition variables and pressure. The secondary variables in the first two rows
are eliminated by a method like Eq. 36. Then the remaining system of NcNg equations is solved simultaneously. This procedure is arguably
a better FIM procedure due both to conditioning and the choice of K-values as variables.
Phase equilibria, Eq. 5 ‍Ki , i = 1, ..., Nc 1‍

Saturation constraint, Eq. 27 ‍x1 / ‍
Component balances, Eq. 26 ‍xi /, i = 2, ..., Nc , p‍
Table 3—Equation arrangement of Chien et al. (1985).
About 1980, I was assigned the task of developing a simulator specifically for CO2 flooding. At the time, most simulators for miscible gas
flooding were modified black oil models, often called four-component models (Lantz 1970). These models assume FCM between the CO2 and oil,
so they ignore all the complexities of phase equilibria and the mechanism by which miscibility is generated or lost. My assignment was to improve
on this approach but with a simulator which would permit grids of at least 10,000 nodes. Somewhat larger black oil model grids were in use in sector
models for simulating waterfloods of the Amoco CO2 pilots (Ader and Stein 1984). These simulations were large at the time. Typical EOS models
were smaller by at least an order of magnitude. Some companies had purchased Cray supercomputers, but at Amoco, we were stuck with IBM
mainframes which were less powerful than a typical early 1990s workstation.
In the early 1980s, the viability of an EOS-based simulator for CO2 flooding was uncertain for two reasons: (1) adequacy of the EOS
for simulating the displacement process and (2) available computer speed and memory. For these reasons, an algorithm that would allow
a flexible representation of the fluids was needed.
Generality was a key feature of the method developed by Young and Stephenson (1983). They describe a formulation which was designed to
permit interchangeable fluid property correlations within one simulator. Three correlation methods were demonstrated: (1) K-values dependent
only on pressure, (2) K-values dependent on pressure and overall composition, z, using Hand’s rule, and (3) cubic EOS. IMPES was used with
implicit production rates (Spivak and Coats 1970). Like the variables of (Fussell and Fussell 1979), saturation was not chosen as an independent
variable. For the EOS case, the order of the selected independent variables and equations are summarized in Table 4. This selection produces a
sparse matrix at each node, so elimination is simple. The formulation separates the fluid properties portion of the calculations so they can be tailored
to the correlation used. This was an object-oriented idea before the term had been invented.
Phase equilibria, Eq. 5 ‍v, xvi , i = 1, ..., Nc 2‍

Component balances, Eq. 29 ‍zi , nh , nw , i = 1, ..., Nc 2‍
Table 4—Equation arrangement of Young and Stephenson (1983).
2022 SPE Journal 11

At the 1982 SPE Symposium on Reservoir Simulation, the papers by Acs et al. (1985) (henceforth ADF) and Young and Stephenson
(henceforth YS) were presented back-to-back in the same session. Acs, Doleschall, and Farkas are from Hungary, and several fields there
were under active CO2 flood. The motivation and even the titles of the two papers are similar—“General Purpose Compositional …” and
“Generalized Compositional ….” The benefits of a generalized approach are echoed in numerous later publications.
In the algorithm of YS, the pressure equation results from the elimination process in a Newton-Raphson iteration. Owing to the sparse
nature of the capacity terms, the elimination is much simpler and more efficient than for Coats’ formulation. It is ultimately a weighted
sum of the component balances, like the construction of Eq. 17 from Eq. 14. ADF introduced the concept of multiphase partial volumes
and described a procedure like the construction of Eq. 17 from Eq. 14. If a rigorous linearization is used, the weighting factors in the
summation are unique (within a multiplicative constant), so the procedures must be related. This fact was apparent to me even in 1982,
and it was simple enough to show. The relationship was presented in later publications, which were largely ignored (Young 1987, 1991)
The relationship was not included in the final 1983 publication because it was considered unethical to change the paper based on what
was learned at the conference.
The YS simulator was originally conceived for simulating West Texas CO2 projects. These were invariably simulated with multiple
noncommunicating layers (i.e., 2½D). However, when Amoco discovered the large Anschutz Ranch East rich gas condensate reservoir, it
prompted a quick extension to full 3D capability. This discovery was an early one with a significant compositional gradient (Metcalfe
et al. 1988). Pressure maintenance with N2 injection was contemplated. At the time, the displacement mechanism was thought to be an
immiscible gas cycling process. However, the YS paper and later work (Renner et al. 1989) show the process is clearly multicontact-
miscible (MCM). The simulator was used for a full field 3D simulation with nine hydrocarbon phase components and 8,000 blocks
(Wendschlag et al. 1983). Because commercially available supercomputers did not have sufficient memory, the calculations were carried
out with overnight runs on an IBM mainframe. It was by far the largest compositional simulation up until that date and signaled that these
simulators were no longer restricted to tiny problems in one or two dimensions.
The relationship between the YS and ADF formulations was described in later articles (Young 1987, 1991). Nevertheless, the industry
continued to view them as different methods (Wong et al. 1987, 1990). Because there has been so much confusion over the years, we will
go through this comparison in some detail.
The phase equilibria and constraint equations apply at each node and involve no neighboring nodes, so for convenience no subscript
is used for the node number. YS work with variables ‍F = nh /V and W = nw /V ‍, the rest are defined above. They did not consider rock
compressibility,
˚ so ‍O = 1‍. To describe
the procedure, first break
˚ the variables into primary and secondary ones. The primary variables are
‍ = z z
1, 2, , : : : , z Nc 2 , F, W, p ‍and the secondary ones, ‍ = v, xv1, xv2, , : : : , xv,Nc 2 .‍
The saturation constraint or volume balance equation, Eq. 27, is a function of all the variables. It is linearized and then the secondary
variables are eliminated with Eq. 36 to give an equation of the form [see Fig. 4 of Young and Stephenson (1983)]:
NP
c 2
Ci zi + Cf n nw
V + Cw V + Cp p = Cr .
h
(37)
‍ i=1 ‍
Because the component balances involve numerous flow terms, it is desirable to minimize the operations on the flow terms, so a natural
P
solution procedure is as follows. First, modify the coefficients C O O O O O
‍ i = Ci V/nh , Cf = Cf Ci zi , Cp = Cp V/t, and Cr = Cr V/t,‍ so
(neglecting a second order z
‍ i h‍
n term) the equation becomes:
NP
c 2
C O f nh + Cw nw + C
O i ni + C O p p = C
O r. (38)
t t t
‍ i=1 ‍
Then the pressure equation is created by substitution of Eqs. 29 and 30, which amounts to just a weighted sum. YS describe the operations
to produce Eq. 38 in the context of Gaussian elimination.
So, how are the coefficients of Eq. 38 related to the partial volumes? The original equation, Eq. 37, is the linearization of the saturation
constraint, so it is dimensionless. However, because of the modifications, the final Eq. 38 has units of volume per unit of time. It is the
same equation but can now be viewed as a volume balance, like the variations in Eq. 27. The relationship of these coefficients can be
found by comparison with Eq. 15. Accounting for the fact that the overall was used in place of the last component, Eq. 38 is:
c 2
NP V
ni
+ Nc 1 n nw O
t + w t C t p = Cr . (39)
h
i Nc 1 t
‍ i=1 ‍
YS were well aware the pressure coefficient gives the total fluid compressibility, see Eq. 18, but they did not interpret the other terms as
partial volumes until the method was compared with that of ADF soon after the 1982 symposium. This was a case of not seeing the forest
for all the trees. The realization does not change the algorithm but provides a clear physical interpretation of the method. Obviously, many
readers were also confused because some reported on a comparison of these “two fundamentally different approaches” (Wong et al. 1987).
This relationship was explained long ago (Young 1987, 1991), but perceived differences continued to linger (Wong et al. 1990; Mifflin
et al. 1991; Farkas 1998; Coats 2000).
There are some minor differences between YS and ADF, but the method for deriving the pressure equation is identical. One difference
is that ADF do not iterate and claim it is not necessary with their formulation. Although the description above outlines only the first iter-
ation, YS allow for additional iterations. Most timesteps converge in one iteration. When more than one iteration is required, it is normally
caused by a phase transition somewhere within the grid, because the coefficients of the linearization are discontinuous at a phase bound-
ary. When a phase transition causes excessive error, it makes sense to take a second iteration rather than cut timesteps. The number of
iterations is greatly reduced from three to eight for methods using an approximate pressure equation (Kazemi et al. 1978; Nghiem et al.
1981).
Independently, Watts (1986) developed a solution procedure based on partial volumes and gives a clear physical interpretation of their
meaning. He contrasts the method to an earlier one using an approximate pressure equation (Kazemi et al. 1978) and explains that the
underlying assumption in that method is that all partial volumes are equal. The calculations presented in the Supplementary file 1 and
Figs. 2 and 3 show this is not a good assumption. He describes an optional sequential IMPSAT extension to the basic method. As
explained below in the section “Temporal Stability and Adaptive Implicit Methods,” this approach will not improve temporal stability for
most compositional gas injection problems.
12 2022 SPE Journal

Watts describes how to calculate the partial volumes only for the black oil case. He states that the method was in use at Exxon with an
EOS, but that formulation was not reported until later (Subramanian et al. 1987). Subramanian et al. selected the liquid phase mole num-
bers as secondary variables. Table 5 shows the organization used in the Exxon simulator. After linearization, the secondary variables can
be eliminated, like Eq. 36, to produce partial volumes and then a pressure equation.
Phase equilibria, Eq. 5 ‍nì , i = 1, ..., Nc 1‍

Component balances, Eq. 26 ‍ni , i = 1, ..., Nc ‍
Volume balance, Eq. 27 ‍p‍
Table 5—Equation arrangement of Subramanian et al. (1987).
Others describe procedures which use the same variable set as in Table 5 (Wong et al. 1990; Collins et al. 1992; Dogru et al. 2009);
however, some use vapor phase mole numbers as secondary variables—not a significant difference.
Farkas (1996) sees differences in the use of extensive vs. intensive variables. In the methods ˚ of YS, ADF, and Watts, what are the
independent variables? The YS algorithm is usually portrayed as having primary variables of ‍ z1, z2, , : : : , zNc 2 , F, W, p ‍as listed above.
These are intensive variables, because they do not depend on amount. The procedures of ADF and Watts are described as having extensive
primary variables of n and p. However, the modifications which ˚produce Eq. 38 or Eq. 39 from Eq. 37 can be viewed as a change from
the original intensive variables to the equivalent extensive ones, ‍ n1 , n2, , : : : , nNc 2 , nh , nw , p ‍, and from a saturation constraint to a vol-
ume balance, see discussion of Eqs. 27 and 28. In all three articles, the variables are transitioned in a similar way. In their article, ADF
start with variables r which is related to our notation by, ‍ri = xi / ‍(Acs et al. 1985). Later, they transition to the variables n after the equa-
tions are discretized. Watts does the same thing—his zm is equivalent to x in Eq. 10 (Watts 1986). For some reason, the YS method got
labeled with the intensive form of the variables, whereas the other two were labeled with the corresponding extensive form. The only real
difference is that YS use the overall hydrocarbon moles, nh, in place of one component mole number. This difference is minor. The point
here is that one’s view of the independent variables and type of constraint is largely a matter of perception.
Wong et al. (1990) and Mifflin et al. (1991) describe differences between a Newton-Raphson and a volume balance method. However,
all methods employ a volume balance or equivalent saturation constraint, see Eq. 27. It would be a misrepresentation to imply one method
uses a volume balance and the other does not. Furthermore, as we have seen, the equation linearization is identical and leads to an identical
pressure equation. Whether one always stops after one iteration or allows for multiple iterations is a minor difference. Given these iden-
tical representations, it seems silly to label one method for the constraint and another for the iterative procedure. These are entirely differ-
ent aspects of the same basic method.
The idea of partial volumes is clearly and simply presented by ADF and Watts, but neither describes how to calculate them analytically
with an EOS. ADF describe a clever numerical differentiation procedure, which is far less efficient than an analytical procedure like that
of YS or later by others (Subramanian et al. 1987; Young 1987; Wong et al. 1990; Collins et al. 1992; Zick 1992; Wang et al. 1997). Also,
an analytical procedure linearizes all the phase equilibria and constraint equations which are needed for efficient updating. Because ADF
and Watts did not describe these details, they were able to produce a clearer presentation.
The difference between the method of ADF and Watts vs. that of YS is one of philosophy or view of the approach used. The view of
the former is—We need partial volumes to create a pressure equation, whereas the view of the later is—We’ve got this messy set of equa-
tions, let’s linearize them and reduce them to a pressure equation. The two approaches lead to the same result, and the partial volumes are
clearly identified in the YS algorithm by comparing Eqs. 38 and 39. The former view is elegant, cleaner, and more direct. Of the three
papers though, YS is the only one to describe a way to calculate the partial volumes with an EOS even though they did not interpret their
coefficients as such at the time.
Young (1987, 1991) made several refinements to the basic methods of YS, ADF, and Watts. The YS formulation supported only an inert
water phase, and there is no significant advantage to use an overall hydrocarbon balance, so the primary variables were generalized to be
like those of ADF and Watts. The simulator used templates to describe the phases present and which phases a component could reside in,
so the framework was general.
An improved set of secondary variables was selected for the EOS. The variables and equations are listed in Table 6. Like (Chien et al.
1985), Table 3, K-values are selected as secondary variables, but instead their logarithms are used. Away from critical points, the K-values
are inversely proportional to pressure and independent of composition, so their logarithms are linear in pressure.
YS achieved efficiency improvements by selective Jacobian updating. The equilibria relations, Eq. 5, should linearize
nicely with
‍ln K ‍and rarely need updating. However, the Rachford-Rice equation, Eq. 34, contains singularities at ‍v = 1/ 1 Ki ‍. The singularities
can cause problems. The problem usually occurs when the amount of a heavy component (Ki → 0) nearly disappears, so that term in the
summation approaches zero divided by zero (Zick 1992; Wang et al. 1997). For that reason, the problem was formulated with a separate
Rachford-Rice equation, which is always updated together with phase equilibria relationships which are selectively updated. Although
this idea seems reasonable, it did not always work well in practice. Zick found the equations can get out of sync, so they should be updated
together (Zick 1992). Further discussion of this aspect of the calculation and the use of ‍ln K ‍as a secondary variable is provided later in
the article in the section “Choice of Secondary Variables.”
In Young (1987), the number of nonlinear iterations was reduced from that of YS. During a timestep, YS update the fluid properties,
and if convergence is not achieved, new coefficients are calculated (at nonconverged nodes only) and then another iteration is performed.
On the other hand, if needed, Young performs inner iterations on the fluid properties until all fluid property calculations are within
Phase equilibria, Eq. 5

‍ln Ki , i = 1, ..., Nc 1‍
Rachford-Rice, Eq. 34 v
Component balances, Eq. 26 ‍ni , i = 1, ..., Nc ‍
Table 6—Equation Arrangement of (Young 1987, 1991).
2022 SPE Journal 13

tolerance. Another outer iteration is performed only if the resulting saturation constraint or volume balance is not within tolerance. YS
mention an average of 2.0 iterations per timestep with an EOS, while Young reports 1.0 iteration per timestep for his tests. Inner iterations
seem to be worthwhile because they are performed only at troublesome nodes like those near a critical point, whereas the outer iterations
involve all the nodes. This idea is used by others (Collins et al. 1992).
The algorithmic improvements described in these articles were overshadowed by the performance of the code (Killough and Kossack
1987). The first article (Young 1987) describes the first benchmarks of the simulator on a Cray supercomputer. It had previously been run
only on an Intel 80286-based PC. The results were rushed, and there were a few bottlenecks that were corrected by the time of the second
article a year later (Young 1991). In the second article, the same benchmarks, based on the third and fifth SPE comparative solution prob-
lems (Kenyon and Behie 1987; Killough and Kossack 1987), were executed in less than half the time.
These articles are known primarily for efficiency improvements using vector calculations. On the average, the time for computing the
fifth comparative problem was roughly 100 times less than others reported, all on a Cray supercomputer. For this tiny problem, only a
factor of 4 can be attributed to the vector processing hardware, so the ratio was about 25 if strictly scalar calculations were used. Other
comparisons have sometimes revealed similar large differences between simulators. CPU times varied by a factor of about 20 for the black
oil models in the first SPE comparative solution problem (Odeh 1981).
The later article shows results for a nine-component, 10,000 gridblock problem solved in as little as 75 CPU seconds (Young 1991).
These results should not have been surprising given the Amoco work (Young and Stephenson 1983; Wendschlag et al. 1983), it was nev-
ertheless a shock to many, especially potential competitors. There was disparaging talk about cheating and skipping important calcula-
tions. In reality, the discrepancy was because others had not put as much effort into compositional simulation. It prompted some to work
harder. When added to the Amoco work, these results helped to make compositional modeling a practical reservoir engineering tool.
An efficient black oil model was later developed within the same simulator framework (Young and Hemanth-Kumar 1992). It was
implemented so that the type of model created was determined by which version of fluid routines was specified during the link step of the
compilation. The bulk of the simulator was completely independent of the fluid property calculations. It was demonstrated on a problem
with 1,000,000 blocks, at a time when 50,000 blocks was considered large.
Wang et al. (1997) describe a simulator also based on the structure in Table 6. They discuss many implementation details (e.g.,
improved methods for testing phase stability and fully implicit option). Many of these issues are discussed below.
Implementation Details
By the mid-1990s, the formulations of EOS compositional models fell primarily into two groups. There were those based on primary
variables of saturation, individual phase composition and pressure, and those based on mole, mass, or overall fractions and pressure. The
approach of Coats and others typifies the first choice (Coats 1980; Harper et al. 1985; Quandalle and Savary 1989; Branco and Rodriguez
1996; Coats et al. 1998; Cao 2002; Fernandes et al. 2017), while the formulations of YS, ADF, Watts, Chien, et al., and others are of the
second type (Young and Stephenson 1983; Acs et al. 1985; Chien et al. 1985; Watts 1986; Subramanian et al. 1987; Young 1991; Collins
et al. 1992; Litvak 1994; Wang et al. 1997; Dogru et al. 2009). Coats calls the first choice natural variables, but we will use the abbrevi-
ation SXP. The second choice will henceforth be MP. Here we discuss these choices and related issues. Although these two choices are
the principal ones in use, others have been discussed. Haukas et al. (2007) investigated the use of isochoric variables, while Gharbia et al.
(2015) evaluate fugacity as an independent variable.
Selection of Primary Variables. Chien et al. (1985) discuss the need for diagonal dominance of the Jacobian, while Acs et al. (1985)
discuss asymmetry of early formulations, although they do not define the term. They are both alluding to the condition of algebraic
equations (Forsythe and Moler 1967). Most likely in their use of the term asymmetric, ADF are referring to the condition of the
procedures published before 1982 (Fussell and Fussell 1979; Coats 1980). Both methods required some variable switching. When one
uses a heterogenous set of variables and equations and some variables become zero or undefined at phase boundaries, it is difficult to
maintain a well-conditioned system of equations. MP formulations are naturally well conditioned because the capacity terms in Eq. 26
are represented by the identity matrix.
Another consideration is the generality of the formulation. For example, in many cases, injected gases can cause three hydrocarbon
phases. The most common are CO2 —oil systems at relatively low temperature, <120°F (49°C). For an MP approach, a general structure
is simpler because the calculations are naturally segregated into the algebraic part (fluid property, physics) and differential calculations
which involve interblock flow. The simulator described in (Young 1987) used a general template structure to define the number of phases
and the partitioning of the components. A general structure is also possible with SXP variables using a similar template arrangement [see
Table 1 in Cao et al. (2009)]. However, with SXP variables, a generalized implementation is more cumbersome because the primary
variables are a function of the problem type. Also, with SXP variables, the generalization affects both the flow terms and capacity terms,
whereas with an MP formulation the capacity terms are always the same.
Many state, without justification, that with MP variables the saturation or volume constraint, Eq. 27, is “highly nonlinear,” whereas
with SXP variables it is linear. However, Figs. 2 and 3 and those in the Supplementary file 1 show that the partial volumes and fluid
compressibility vary slowly with composition and pressure, except at phase boundaries and near a critical point. ADF claim that iteration
is not necessary to converge the constraint (Acs et al. 1985). Experience has shown that it usually requires one iteration except when a
phase transition occurs somewhere in the grid. Be assured the discontinuities at phase boundaries will be felt by either choice of
variables.
There is an obvious difference in linearization of the capacity terms for hydrocarbon components in the transport equations, Eq. 26.
For the MP approach, these terms are especially simple because mass is the quantity conserved, see Eq. 28. The SXP approach requires
the following expansion for the hydrocarbon components:
8 9 h 8 9i
V : O xi ; = V
O :S` xì` + Sv xviv ;
t t
(40)
O O dO
= V
t ` xì S` + S` xì S` xì
`
`
+ v xvi Sv + Sv xvi Sv xvi
v
v
+ dp p S` xì` + Sv xviv .
‍ ‍
‍ ` and v must

‍ be expanded further in terms of the primary and secondary variables. The secondary variables can then be eliminated
using Eq. 36. When all the component balances are included and residual Ar terms are negligible, the resulting capacity terms look like:
14 2022 SPE Journal

" #
8 9 P
Nc NP
c +2
V : O xi ; = V
D + Dik k
t t ik k
k=1 k=1
(41)
V P
Nc
= t Dik k.
‍ k=1 ‍
Using the Schur complement of Aχ, ‍D = D D A1

A ‍. D is a dense matrix at each node. For the IMPES case, (Coats 2000) describes
how a pressure equation can be constructed. The procedure he describes requires factoring Aχ, the calculation of ‍A1 A ,‍ a matrix-matrix
multiply to calculate D, then D is factored and the equivalent of a solve produces the desired parameters. These parameters
are proportionl
1 3
to the partial volumes, but Coats calls them IMPES reduction variables. These calculations require a total of
‍ O 5 3 c ‍ floating point
N
operations. In the next section, we show that with MP variables the same result can be achieved with O 13 N3c operations. Also, after the
‍ ‍
pressure solution with IMPES, the SXP approach requires the solution of a matrix problem of the form ‍D = Dr ‍at each node to update
the primary variables.
Some argue that SXP variables make for simpler linearization of the mobility terms for FIM (Cao 2002). Similar statements are echoed
in several articles; they go something like—“SXP variables simplify use of the chain rule to linearize mobility terms.” The argument has
merit if only a subset of the variables are implicit, for example, Np − 1 saturations or IMPSAT (Harper et al. 1985; Quandalle and Savary
1989), but as explained below in the section “Temporal Stability and Adaptive Implicit Methods,” the IMPSAT procedure does not
improve stability for many important problems. Superficially, the linearization appears to be simpler, but let us look at the issue in more
detail.
For the general case, there are a total of Np(Nc +1) mobility variables. For the example case with inert water, there are 2Nc variables in
the mobility terms if one saturation is eliminated with the saturation constraint. Excluding pressure, the number of primary composition
variables is only Nc − 1. Therefore, flow terms involving the secondary variables are unavoidable. The secondary variables must be elim-
inated using Eq. 36 to reduce the dependence to primary variables only. These eliminations are straight forward, but the number of calcu-
lations is significant. No interblock flow is involved, so eliminations can be executed on a node-by-node basis.
The mobility terms depend on saturation and all the individual phase compositions, Sα, xαi. Capillary pressure terms are also functions
of phase saturation, so the linearization issues are similar and are not discussed. If the dependence is on a primary variable, the computa-
tions are simple. It is common to neglect the dependence on viscosity and density, so the following is required:
NP
p 1
kr˛ x˛i x˛i @kr˛ kr˛
˛ ˛ ˛ ˛ @Sˇ
Sˇ + ˛ ˛ x˛i . (42)
‍ ˇ =1 ‍
The first term is simple because saturations are independent variables, see Table 2. The second term is neglected for the so-called IMPSAT
procedure. For FIM, the second term is simple for the phase fractions that are primary variables (i.e., α = v and i = 3,…, Nc − 1); however,
for the others, we have from Eq. 36:
h i
x˛i = A1 A Ar , (43)
‍ k ‍
where again ‍ ‍ are the primary variables listed in Table 2, and k = 1,…,Nc +1 corresponds to the secondary variables
‍
= fxì , i = 1, ...,
Nc3 1, xv1 , xv2 g‍in that order. Given the LU factors of Aχ, the calculations are equivalent to a matrix-matrix
multiply, requiring O
‍ 2Nc ‍floating point calculations.
Watts (1986) claims the simplicity of the elimination is an advantage for the MP approach. Using MP variables, the phase equilibria
relationships must be factored to calculate the partial volumes. The fluid volume is linearized in terms of primary and secondary variables
and then the phase equilibria relationships are used to eliminate the dependence on secondary variables, like Eq. 36. For implicit calcula-
tions, the phase saturations and individual phase compositions linearize just like the fluid volume, and the secondary variables can be
eliminated the same way. The coefficients and computer code needed for these other eliminations can be reused (Subramanian et al. 1987;
Young and Russell 1993).
As discussed in the next section, factorization of the phase equilibria relationships, Aχ, for MP requires O 13 N3c floating point calcu-
‍ ‍
lations. Two phase saturations and all 2(Nc − 1) phase composition variables are reduced to a dependence on primary variables with the
calculations outlined in the Appendix. The result is:
NP
c 1
S˛ = a˛j nj + ap˛ p
j=1
NPc 1
(44)
x˛i = b˛ij nj + bp˛i p.
‍ j=1 ‍
3

‍ 2Nc ‍floating point operations are required to calculate these coefficients at points with two hydrocarbon phases. This is the only cal-
O
culation of that order. Substitution into Eq. 42 produces the desired result. As shown in the Appendix, there are many common subexpres-
sions in these calculations which lead to improved efficiency. Note, as described in the Appendix, one could easily relax the assumption
of constant phase molar volume in Eq. 42 by calculating ‍ x˛i ‍in Eq. 44. This calculation would require
‍ x˛i /˛ ,‍ and using it in place of
little additional effort.
From
3
this simple analysis, it appears that the calculations needed to linearize flow terms with both SXP and MP variables for FIM are
‍ 2Nc ‍, not counting the factorization step. This simple analysis is based on the dominant calculations and the mobility terms only, but
O
it suggests there is no advantage with SXP variables even for FIM. This result is not apparent with a superficial comparison.
Some have claimed superiority of SXP variables with regard to solution of the nonlinear algebraic problem (e.g., Voskov et al. 2009;
Zaydullin et al. 2012). This issue is most important for FIM. In the section “Temporal Stability and Adaptive Implicit Methods,” we dis-
cuss tradeoffs between FIM and other methods and show that for FIM to be competitive with other methods, the timesteps must be large.
The claimed superiority for SXP variables centers around the ability to dampen Newton-Raphson iterations using procedures like the
Appleyard chop (Appleyard et al. 1981). Advantages have been demonstrated for immiscible displacements. For immiscible fluids (e.g.,
2022 SPE Journal 15

water/oil), there is little difference between the variable choices,
because when the phases are pure and inert, the component and phase are
synonymous, so combining Eqs. 10 and 27, S˛ /˛ = n˛ / V O . The two variable choices differ only by terms dependent on phase and
‍ ‍
rock volume which are generally weak functions of pressure. Also, (Alpak and Vink 2018) show that Eq. 44 can be regarded as a trans-
form between the two variable sets and how this information can be used to achieve equivalent damping in either formulation. These same
articles have found MP variables to be superior for MCM displacements (see Cao et al. 2017). An MCM displacement involves near
critical calculations, see Fig. 5 below. Near a critical point, small composition changes produce large saturation changes, see Fig. 3, caus-
ing ill conditioning for an SXP formulation. In summary, for the nonlinear problem there appears to be no advantage to the use of SXP
variables.
In this section, we have compared many computational features of SXP and MP approaches. Most of the advantages claimed for SXP
formulations do not stand up to scrutiny while MP variables appear to offer advantages in several areas. What about direct comparisons?
There have been a few comparisons (Thele et al. 1983; Killough and Kossack 1987; Voskov and Tchelepi 2012). We will let the reader
evaluate these results, but we offer a few words of caution. Comparisons like these are highly dependent on test problem selection and
implementation details. The next section shows how the choice of secondary variables can make a difference and how some procedures
can produce ill-conditioned algebraic systems near phase boundaries, causing the potential for convergence problems. In the past, large
differences have been found for simulators using the same basic formulation. The section “Efficient Implementation” discusses imple-
mentation efficiency and gives examples showing large differences unrelated to the basic formulation.
Choice of Secondary Variables. Many different secondary variables have been used with both SXP and MP approaches. Some of the
choices used are listed in Table 7.
(Coats 1980; Coats et al. 1998; Fernandes et al. 2017) ‍xì , i = 1, ..Nc 1, xv1 , xv2 , Sa ‍
(Young and Stephenson 1983) ‍v, xvi , i = 1, ..Nc 2‍
(Subramanian et al. 1987; Wong et al. 1990; Collins et al. 1992; Dogru et al. 2009) ‍nvi , i = 1, ..Nc 1‍
(Chien et al. 1985; Litvak 1994) ‍Ki , i = 1, ..Nc 1‍
(Young 1987; Wang et al. 1997)
‍v, ln Ki , i = 1, ..Nc 1‍
Table 7—Secondary variables.
Some may use slightly different but similar variables. For the case on Line 3, some use vapor mole numbers while others use liquid
mole numbers—not a significant difference. Line 2 lists extensive variables, while Line 3 lists the corresponding intensive ones. As dis-
cussed above, this difference is not significant. The last two cases are also similar. The vapor fraction, v, can be thought of as either a
secondary or tertiary variable depending on one’s view. As shown below and in the Appendix, it is relatively simple to eliminate v by
direct substitution.
When making this selection, proper scaling of the equations is of paramount importance to produce a well-conditioned system of
equations. The first three cases listed use variables which become zero or undefined at phase boundaries. These choices can cause a loss
of accuracy near a phase boundary (Zick 1992; Michelsen and Mollerup 2004) and lead to convergence difficulties.
The last two cases use either equilibrium K-values or their logarithms, which are well defined at all phase
boundaries. Rather than
describe the pros and cons of each choice, we will make the case for ‍ln K ‍and let the reader decide. ‍ln K ‍is arguably better because it

is more symmetrical, that is, it does not depend on how the phases are ordered, because ‍ln K = ln 1/K ‍. Fugacity coefficients are more
directly expressed in terms of their logarithms (Prausnitz and Lichtenthaler 2004). Away from the critical point, the K-value logarithms
are linear in pressure and independent of composition, so their linearization should often be accurate and rarely need updating.
Treating K-values as unknowns rather than as data does require reorienting one’s thinking. When the ‍ln K ‍are variables, the equations
for the phase equilibria, Eq. 5, look like:

‍ !i = ln xvi /xì = ln 'ì x` , p ln 'vi xv , p ,‍ (45)
where ω are the logarithms of the K-values. These equations are linearized with respect to the primary and secondary variables. The
Appendix describes the linearization of the phase equilibria and Rachford-Rice equations and the overall hydrocarbon-specific volume.
Respectively, they can be represented in matrix notation as:
Aw w + Av v + An + Ap p = Ar
eTw Bv v + BTn = 0
nh
‍ n = CTw w + Cv v + CTn + Cp p,‍ (46)

where for convenience, the intermediate variables ‍wi = xvi xì /zi !i and i = ni /n‍ are defined, along with the unit vector ei = 1. All
T
vectors and matrices are in bold face. The matrix Aw is symmetric, so it can be factored into the form LDL
, where
D is a diagonal matrix
1 3
and L is unit lower triangular. By taking advantage
of the symmetry, the matrix can be factored in O 3 Nc ‍ floating point calculations,
‍
whereas a nonsymmetric matrix requires O 23 N3c operations. So, by exploiting the symmetry, both storage and calculations associated
‍ ‍
with Aw are cut nearly in half. The original work (Young 1991) did not exploit this symmetry but later work did (Zick 1992).
These calculations are often described using the generic term flash calculations. However, they are more general. For a true flash cal-
culation, the pressure and feed are specified, whereas here they are variables determined by the transport equations. The equations for a
true flash calculation are produced when Δp and η are set to zero.
Michelsen (1982b) describes a Newton-Raphson flash calculation using ‍ln K ‍as variables and requiring only symmetric matrix cal-
culations. Several authors claim they use this procedure, but we discovered it becomes singular at phase boundaries (Zick 1992). Near
16 2022 SPE Journal

phase boundaries, the singularity can cause ill conditioning, a loss of accuracy, and convergence issues for implementations using the
procedure. Owing to the singularity, we formulate the problem in a different way from Michelsen but the ideas are similar.
A problem with exploiting the symmetry is that the phase equilibria relationships are constrained by the second equation, the Rachford-
Rice equation. If it is used to eliminate the vapor fraction, the phase equilibria relationships are:
A Q n + Ap p = Ar ,
Q ww + A (47)
‍ ‍
where:
Av eT
‍ AQ w = Aw + Bv .‍ (48)
And the overall molar volume linearization is reduced to:

nh Q Tw w + CQ Tn + Cp p.
‍ n =C ‍ (49)
The matrix ‍AQ w ‍is not symmetric. Faddeeva (1959) treats border matrices like this or equivalently the Sherman-Morrison theorem can be
used (Press and Vetterling 1999). The inverse is:
1
1 A1 Av eT T
AQ w = I + w T 1 A1
w = I + Ue Aw . (50)
‍ Bv +e Aw A v ‍
2 1 T
The formula requires an O ‍ 2Nc ‍ solve for ‍Aw Av ‍, calculated using the LDL factors. With this approach, the savings of a symmetric
factorization are realized. In addition to reducing calculations, Eq. 50 offers some flexibility in the calculation of the vector U to create a
properly scaled, well-conditioned system of equations. Alternate ways to calculate U are discussed briefly in the Appendix.
Once Aw is factored, the partial molar volumes and total hydrocarbon fluid compressibility, Eq. 15, can be calculated:

‍ nh = T n nh Ch p + r ,‍ (51)
T 1
where after calculating Ri = CQ w AQ w with Eq. 50:
‍ i‍
T
T = C
Q RT AQ n
T n
Ch = R Ap Cp /h (52)
‍ r = RT Ar .‍
2
The O 13 N3c factorization of Aw is the only operation of that order, but there are a few O ‍ Nc ‍operations (e.g., in the calculation of U
‍ ‍
and R). This implementation using Eq. 50 is a major improvement over the original one used in Young (1991) which did not exploit sym-
metry and performed eliminations on the entire An array resulting in O 2 23 N3c calculations. In addition to the reduced computation, the
‍ ‍
equations are scaled better to reduce rounding errors. Details are given in the Appendix.
This procedure, developed by (Zick 1992), is far and away the most efficient one for calculating the partial volumes and total fluid
3 on An. In that
compressibility. If the linearizations, Eq. 44, are also required for FIM, it is then most efficient to perform the eliminations
case, minor improvements to the partial volume calculation may be possible. However, for IMPES, it saves the O ‍ 2Nc ‍calculations to
perform eliminations on An. The savings are possible because the MP approach yields a diagonal matrix for the capacity terms, so no
eliminations are needed for those terms.
Because the implementation of this method involved modifications to the Jacobian, the procedure for periodic Jacobian updates was
also changed. In the original implementation (Young 1991), the Rachford-Rice equation was treated separately. It was always updated,
whereas the rest of the Jacobian was updated only periodically. For reasons discussed above, this procedure was replaced by one where
the entire Jacobian is updated together, producing about a 14% increase in the speed of the calculations. The other changes leading to Eq.
51 netted an improvement of only 2%. The improvement was small because on the average only about 2 to 3% of the two-phase points
require a full relinearization. This is a testament to the effectiveness of using the ‍ln K ‍as secondary variables.
Material Balance. An advantage of the MP approach is that it is naturally conservative. After solving for pressure, it is substituted into
the component balances to update the primary variables. The differential material balance equations are exact within machine precision.
When this was first advanced (Young and Stephenson 1983), many thought there must be some cheating involved. The equations are
nonlinear, so there are errors in the saturation constraint or volume balance and in the phase equilibria calculations. The primary stopping
criteria is based on volume errors. A maximum error of 0.005 in the saturation constraint is usually sufficient. It is natural to carry
this volume or saturation constraint error forward from timestep to timestep. For each timestep, it is iterated (if necessary) down to an
acceptable tolerance. Coats et al. (1998) and Coats (2009) calls it relaxed volume and writes as if it is somehow contrived, but it is the
natural and correct treatment.
The SXP approach will usually produce material balance errors. First, a material balance error will occur if the saturation constraint,
Eq. 27, is used to update one saturation by difference. Black oil models normally use this procedure. Errors also occur owing to expansion
of the capacity terms. For example, with saturations as primary variables, consider the following term which is embedded in Eq. 40:

‍ S` xì = xì S` + S` xì + S` xì .‍ (53)
The last term is second order and normally neglected when updating. Nevertheless, it will create a small error in material balance. Error
terms of this type abound throughout the calculations with the SXP approach, see Eq. 40. Material balance is not often discussed by pro-
ponents of the SXP formulation.
2022 SPE Journal 17

YS update the ‍ zi nh ‍ term so that mass is conserved (Young and Stephenson 1983). One can make any method conservative by
careful updating. Even with implicit calculations, a conservative method can be achieved by ensuring the material balance equations are
satisfied exactly (Young and Russell 1993). With the SXP approach, material balance errors are unavoidable unless the moles or mass are
updated first and then the other variables are calculated from them. Coats et al. (1998) do precisely that in the so-called relaxed volume
method. This procedure is not truly SXP but rather a hybrid, because effectively n is used as an intermediate variable. SXP or natural
variables lead to this decidedly unnatural method of updating, whereas an MP formulation is naturally conservative. Coats describes some
other advantages of a conservative updating scheme (Coats et al. 1998; Coats 2009). Mass conservation in a sequential SXP formulation
is discussed in (Møyner and Tchelepi 2018).
Detecting Phase Appearance. An important issue when using an EOS for compositional modeling is determining the number of phases
that exist at a given point. The easiest case is when a node is two-phase. In that case, it is simply updated, and if it has become single
phase, the resulting vapor fraction will usually be outside the interval [0,1]. Alternatively, the update may converge to the trivial condition
where all K → 1, requiring further examination. You would think nodes with a single hydrocarbon phase would be easy, but that is not the
case. They must be monitored to determine when a transition to two-phase occurs. Coats (1980) and Coats et al. (1998) use a saturation
pressure calculation, so if the node is above the saturation pressure, it is single phase. Coats describes corrective actions he must take
when the procedure fails. Nghiem et al. (1981) describes a method based on a flash calculation. (Young and Stephenson 1983) describe
their experience with these and other methods and did not find any to be satisfactory in all cases.
Young and Stephenson saved computations by testing only well blocks and blocks adjacent to two-phase blocks. This method, called
the neighborhood method, was adopted by (Wang et al. 1999). An alternative procedure is based on a shadow region (Rasmussen et al.
2006; Belkadi et al. 2013). The objective of these methods is to reduce the number of blocks where testing is needed. One is based on the
proximity in physical space, while the other is based on the proximity in composition and pressure space. The earlier method is not men-
tioned in later work, so the most effective method is yet unknown. Of course, with the neighborhood method, one must iterate when a new
two-phase condition is found, because all nodes can go to two phases on a single timestep.
The phase stability test, based on Gibbs free energy analysis (Michelsen 1982a), has become the standard method to test for phase
appearance. Its first use in reservoir simulation was reported by Young (1987). It is an expensive calculation. For the benchmarks in Young
(1991), it consumed over 40% of the total calculation time, when only about 10% of the blocks were tested. However, it is very reliable
when initiated with K-values from the Wilson equation (Whitson and Brulé 2000). With millions of these calculations over the course of
a simulation, reliability is extremely important. A side benefit of the method is that it provides estimates for the K-values which can be
used to start a Newton-Raphson flash calculation, and it too will converge except when in extremely close proximity to a critical point.
Several others perform slower successive substitution iterations before switching to a faster Newton-Raphson iteration, but both stability
and flash calculations can be made very reliable without the slower iterations (Zick 1992). Their elimination makes for a much cleaner
more efficient code.
In the original work by Young (1987), the stability test calculation was not implemented efficiently for many of the reasons discussed
in Belkadi et al. (2013). Most importantly, the implementation in Young (1987) did not exploit the symmetry of the matrix problem. In
Young (1989), a procedure like that of Belkadi et al. was implemented, and the improvements were consistent with their results. Our
implementation differs in that only Newton-Raphson iterations are performed. By exploiting the symmetry, the matrix problem can be
both created and solved in almost half the time (Faddeeva 1959). Because this makes up only part of the overall calculation, the calcula-
tion time for stability testing was reduced about 25% and overall simulation time by 11% (Young 1989). This reduction was relatively
easy to achieve and is significant given the simulator was already more efficient than most. However, testing of single-phase points still
consumed 35% of the overall time, about the same as the calculations for a much larger number of two-phase blocks.
With the phase stability test by Gibbs free energy analysis, one must test from both dewpoint and bubble point sides (i.e., two separate
tests). In addition, when it fails to find a two-phase condition, it often converges to a trivial condition where all K → 1. Consequently, it
does not provide any information that can be used in the next time step. It must be restarted from scratch every time. In contrast, the lin-
earization at two-phase points gives excellent starting estimates. With ‍ln K ‍as secondary variables, the linearization rarely needs updat-
ing. It seems logical to attempt replacing the stability tests with flash calculations.
Whitson and Michelsen (1989) show that so-called negative flash calculations provide a definitive test for phase stability. A negative
flash is one where the overall vapor or liquid mole fraction, v or 1 – v is negative. It is much like a normal flash calculation, except v is
allowed to stray outside the interval 0 to 1. However, there are three issues with the method. First, sometimes there is only a trivial solution
and not a valid negative flash solution. Generally, a valid negative flash solution requires the overall composition to lie on a tie-line exten-
sion. For the ternary system displayed in Fig. 1, no valid negative flash solutions exist to the right of the critical tie-line. Second, there is
the question of how to start the flash. A normal flash is initiated using the K-values determined during a stability test. For a negative flash,
this option is often not available. Third, it is more difficult to maintain a well-conditioned algebraic system of equations. For example,
rounding errors can become prevalent when the overall composition of the heaviest component approaches zero. The problem is related
to the one which can occur with the Rachford-Rice equation (Zick 1992; Wang et al. 1997), because that equation is embedded in those
for a Newton-Raphson flash.
Zick (1992), Wang et al. (1997), and Heidari and Maini (2014) discuss the use of negative flash calculations to test for phase appear-
ance. Wang et al. do not discuss the limitations of the method or how it was implemented, nor do they present any computational experi-
ences, such as the gains in efficiency. Later articles from the University of Texas state that stability tests are used in this code, so it is
unclear how thoroughly negative flashes were implemented by Wang et al. (Heidari and Maini 2014) use negative flash calculations in a
thermal simulator but do not mention the method will not always work.
Zick (1992) analyzed the simulator in Young (1991) after the more efficient symmetric stability test calculations were implemented.
He found that the computation time to perform a stability test was roughly 8 to 10 times the calculations for a typical two-phase point.
Because stability tests consumed about 35% of the computation time, a potential saving of perhaps 25 to 30% appeared possible if nega-
tive flash calculations replace most of the stability tests. All the stability tests cannot be replaced because the negative flash method works
only for compositions that lie on a tie-line extension.
Zick implemented several options for initiating negative flash calculations. One option is to use K-values from stability tests which
produce a single-phase result, but do not converge to a trivial solution. Another, more aggressive, option is to use K-values from neigh-
boring two-phase cells to start a negative flash. The method was tested on several problems, including the third and fifth comparative
solution problems and CO2 floods with 16,000 to 256,000 blocks. The results were both problem dependent and computer dependent with
reduced CPU times ranging from 10 to 36%. The largest reductions were for the CO2 floods which spent 47% of the time doing stability
18 2022 SPE Journal

tests with the original code. On the average, about 80% of the stability tests could be converted to negative flash calculations. After ana-
lyzing the results, negative flash calculations appear to run 2.7 times faster on the average. This differs from the original estimate of 8 to
10, suggesting that negative flash calculations are slower than the calculations for a typical two-phase point. The savings are large never-
theless and yield a substantial reduction in CPU time.
Temporal Stability and Adaptive Implicit Methods. FIM is more reliable than IMPES because it is not subject to a timestep limitation
like Eq. 35. However, it requires more calculation. The tradeoffs between IMPES and FIM are different for black oil and compositional
simulations. The differences in the computations per timestep are much larger for a compositional simulator because of the larger number
of coupled equations. The computations per timestep for FIM vs. IMPES is about 3 to 6 for a black oil simulator and perhaps 10 to 30
for a compositional model.
The temporal stability characteristics are also different. Consider a two-component system, either completely miscible (FCM) or
immiscible. Eq. 35 determines the stability for IMPES. The fractional flow functions for the two cases are shown in Fig. 4. The miscible
function is linear, while the immiscible curve has the usual S shape. For immiscible flow, the curve has a zero derivative at both extremes.
The derivatives are large only at intermediate saturation. These saturations are within the displacement front, so each point quickly passes
through these conditions. For a typical water/oil displacement, stability issues are isolated in both space and time (Young and Russell
1993). They occur only at producing wells for intermediate water cuts. Injection wells quickly reach a high saturation for which ‍f 0 0‍,
so stability is not an issue there.
Fig. 4—Fractional flow functions.
In contrast, for an FCM displacement, the fractional flow curve is linear, so ‍f 0 = 1.‍ Consequently, when throughput (i.e., Q/V) is high,
stability is important for all times. The third comparative solution problem is a case in point. Like most compositional problems, it has a
gas injection well . Stability limits timesteps to relatively small values for the entire injection period (Young and Russell 1993). These
basic differences in fractional flow characteristics suggest that compositional simulation will benefit more from the use of FIM than black
oil simulators, but the computational expense is considerably greater.
Many seem to be confused about the differences between variable changes and stability. They are related only in the sense that large
fluctuations result after temporal instability has occurred. However, changes are not predictive of temporal instability. My first simulator
was for FCM problems (Young 1981). I quickly discovered that instability could crop up at a production well long before the flood front
reached the well.
The stability restrictions with IMPES can be reduced with implicit well calculations (Spivak and Coats 1970). These calculations can
be implemented at production wells without altering the matrix structure. The implementation involves modifications to the partial vol-
ume terms used to create the pressure equation, so its cost is insignificant. A theoretical fourfold increase in stable timestep size is possible
for 2D areal flow. Smaller but significant improvements are observed in practice (Young 1991). IMPES should never be implemented
without this feature.
In addition to IMPES and FIM, the other options are sequential implicit calculations and adaptive implicit calculations. Watts (1986)
describes a sequential IMPSAT procedure, that is, the second term in Eq. 42 is neglected. Others have applied the same IMPSAT approach
but use a simultaneous rather than sequential solution (Harper et al. 1985; Quandalle and Savary 1989; Cao and Aziz 2002). The IMPSAT
approach will improve stability for immiscible displacements but will do nothing to improve the stability of compositional gas injection
problems with their undying restrictions. For this reason, IMPSAT is not attractive for the majority of compositional simulations which
involve gas injection.
2022 SPE Journal 19

Another approach is the adaptive implicit method (AIM). Unfortunately, this name has been used for two different methods. The first
uses simplified iterative procedures to solve the equations for FIM. These include the original method of Thomas and Thurnau (1983) and
related methods (Bertiger and Kelsey 1985; Tan 1988; Watts and Rame 1999). We will call these quasi-Newton methods. The other
method is true AIM, where individual connections in the simulation grid are treated implicitly while others are treated explicitly. AIM
does not converge to the FIM solution but rather something between FIM and IMPES. AIM procedures are described for black oil
(Forsyth and Sammon 1986; Fung et al. 1989), compositional (Collins et al. 1992; Farkas 1993; Young and Russell 1993; Cao and Aziz
2002; Moortgat 2017), and thermal (Tan 1988; Moncorgé and Tchelepi 2009) simulation.
After making plausible assumptions about typical distributions of Q/V and the computational demands of explicit vs. implicit calcula-
tions, (Young and Russell 1993) (henceforth YR) developed a heuristic model to examine the likely tradeoffs between IMPES, FIM, and
AIM. Calculations are made using the parameters listed in Table 8. Most are self-explanatory, but are described more fully in YR.
Work Ratios Factor

FIM:IMPES black oil 4
FIM:IMPES compositional 15
AIM:IMPES (no implicit) 1.1
AIM:FIM (100% implicit) 2
Table 8—Parameters of heuristic

model.
As one increases timestep size with any method, the computational work per timestep increases owing to a greater number of iterations
0.3
and so forth. For the heuristic model, the work escalates in proportion to ‍ t ,‍ so we assume the work per timestep doubles with a
tenfold increase in timestep size. Throughout this discussion, we ignore truncation error, although it certainly plays a role (Lantz 1971;
Peaceman 2000; Loubens et al. 2009). The heuristic model also assumes the nonlinear problem will converge for any timestep size,
whereas timesteps for AIM and FIM are limited by nonlinear convergence issues.
Figs. 5 and 6 show results for black oil and compositional cases. The IMPES calculations are assumed to run at the stability limit, but
relatively larger timesteps are considered for FIM and AIM. The results are in terms of calculation efficiency relative to IMPES and are
plotted against relative timestep size. The AIM results depend on a heterogeneity parameter (see YR) which characterizes the distribution
of Q/V. The figures list values of the parameter from 5 to 320 together with the percentage of implicit blocks at the optimum. The model
suggests that the computations for FIM and IMPES are comparable when the timesteps with FIM are larger by a factor of 7 and 48 for
black oil and compositional cases, respectively. Timesteps of this size are often quite reasonable for the black oil case but could be exces-
sive for a compositional simulation. One would likely encounter convergence difficulties before timesteps of that magnitude are reached
with a compositional simulator.
Fig. 5—IMPES, FIM, and AIM tradeoffs black oil.
20 2022 SPE Journal

Fig. 6—IMPES, FIM, and AIM tradeoffs compositional.
As one would expect, the potential efficiency of AIM increases with the heterogeneity of the system. Figs. 5 and 6 show relative effi-
ciency values greater than 10 are predicted for cases with high heterogeneity. The optimum efficiency predicted by the model occurs with
less than 7% implicit blocks for black oil and less than 1.7% for compositional simulation. Most benchmark problems are relatively
homogenous. Relative efficiency factors of only 1.5 to 2.4 are reported for the third comparative solution problem in (Collins et al. 1992;
Young and Russell 1993). However, YR state that an efficiency of 4.5, with 1% implicit blocks, was achieved for a full field simulation.
The qualitative behavior of the heuristic model has been born out with calculations, and the values in Table 8 are based on experience
gained from calculations. One may argue with the values selected, but the overall trends are undeniable.
The model clearly shows that compositional simulation can benefit from AIM more than black oil. Not because of the efficiency
improvement relative to IMPES, but because FIM is so computationally intensive. For example, when the timestep size is 10 times the
IMPES stability limit, the model suggests that relative to IMPES, FIM is slower by a factor of 3 for compositional problems but is 25%
faster for black oil simulations.
There are two important difficulties in the use of AIM. The linear algebraic problem is dynamic and unstructured, which puts a consid-
erable burden on the linear solver. Also, one must have accurate criteria for selecting the connections to be treated implicitly. This is called
the switching criteria. The first AIM implementations based the switching criteria on threshold changes in the solution variables. This
approach is not very good because it can only detect instability after it has already occurred and then there is no basis for switching back
to explicit calculations. Later implementations have used numerical stability analysis in the switching criteria.
For AIM to be efficient, the switching criteria must be accurate, and it obviously must consider temporal stability. Several studies
investigate the temporal stability of the IMPES approach (Todd et al. 1972; Peaceman 1977; Fung et al. 1989; Russell 1989; Grabenstetter
et al. 1991; Watts and Rame 1999; Cao and Aziz 2002; Coats 2003b; Wan et al. 2005; Agarwal et al. 2006; Lu et al. 2009; Moncorgé and
Tchelepi 2009; Maes et al. 2013). One study performs computations comparing various criteria (Franc et al. 2016). Most of these articles
consider immiscible flow, often with only two components (e.g., water/oil systems). Some analyze simplified versions of the governing
equations, sometimes just a single dimension. The methods of analysis include von Neumann stability analysis and amplification matrix
analysis. The analyses are biased toward immiscible displacement. When compositional effects are considered, they are treated in a sim-
plistic way (e.g., constant K-values and volumetric considerations are mostly ignored).
This stability analysis problem is more complicated than the usual stability analysis for a partial differential equation because the
analysis must account for the implicit pressure. The pressure equation determines the flow field, which can be considered fixed when
analyzing for the instability caused by explicit mobilities. That the flow field is relatively static is the entire basis for streamline
simulators.
When phase equilibria is considered, the problem is usually described in terms of the mass or molar form of the equations, but that is
not the proper analog. Of course, the equations need to be simplified down to their essence. Rather than working with the molar or mass
form of the equations, Eq. 13, it is better to work with the volumetric form, Eq. 14. The latter case makes the equations like those for
strictly immiscible displacements, but they include all the effects of phase equilibria. We start with the volumetric form of the equations,
Eq. 24, and assume the partial molar volumes are constant.
P Np
@
@t
Si + r u + ıu˛ f˛i r D ˛ r S˛i = 0. (54)
‍ ˛=1 ‍
The pressure equation or continuity equation for this simplified model is found by summation:
‍ r u = 0.‍ (55)
Of course, as shown in Supplementary file 1 and Figs. 2 and 3, the partial volumes are not constant, but they are relatively slow changing
functions except near a critical point and at phase boundaries. They vary to a lesser degree than the specific volume of the individual
2022 SPE Journal 21

phases. The variation of partial volumes is universally ignored in streamline models. If this is not a good assumption, then those models
must be invalid. We believe this assumption is justified for the purpose of stability analysis. Because pressure is implicit, the analysis is
performed assuming a fixed flow field.
Most of the previous analyses simplify volumetric and phase equilibria considerations. Because they work with the molar form of the
equation, individual phase specific volumes are assumed constant. Phase equilibria is either ignored or K-values are assumed constant.
These assumptions are far worse. With the formulation here, all the phase equilibria constraints are rigorously accounted for in the linear-
ization of the individual phase saturation and composition, Eq. 44.
Eq. 55 is multiplied by the fractional flows, fi, and subtracted from the component balances, Eq. 54, to produce the nondivergent form
of the component balances:
P Np
@
@t
Si + u + ıu˛ rf˛i + f˛i r ıu˛ r D ˛ r S˛i = 0. (56)
‍ ˛=1 ‍
The discrete form of this equation is the one that should be analyzed for stability. The second term in the summation is usually small for
forced convection problems with significant pressure gradients.
If dispersion is neglected, then after discretization, Eq. 54 is:
PNp P
Vj
t Sij = f˛i u Tjk pk pj + p˛k p˛j ˛ Hjk , (57)
‍ ˛=1 k ‍

where k sums over all connections to point j, u is the upstream point, ρα is an average between j and k, and ‍Hjk = g Hk Hj ‍represents
the head term.
The pressure equation is found by summation over i:

P P
u Tjk pk pj + f˛ u p˛k p˛j H N jk = 0. (58)
‍ k ˛ ‍
If the pressure equation, Eq. 58, is multiplied by (fi)j and subtracted from the component balance, Eq. 57, the discrete nondivergent form
of the component balance is produced:
Np P
P
Vj
t Sij = ()k Tjk [(f˛i )k (f˛i )j ][pk pj + p˛k p˛j ˛ Hjk ]
"
˛=1 k=in # (59)
PNp P P
+ u f˛i j Tjk N ˛ Hjk + p˛k p˛j fˇ u pˇ k pˇ j ,
‍ ˛=1 k ˇ ‍
where k = in indicates a summation over all neighbors with flows into block j, so the pressure gradient term in brackets is positive. The
last term should be relatively small and is neglected for stability analysis. Note, owing to gravity, the inflow for one phase may be different
from another—countercurrent flow. Also note that no assumptions have been made regarding the structure of the grid, because that affects
only the number of neighbors and the values for pore volume, V, and transmissibility coefficients, T.
The approach advocated here is consistent with other analyses. For example, if the capillary pressure terms are neglected, and terms
are combined for simplicity, the equation which should be analyzed for stability is:
Np P h
P i
Vj
t Sij = f˛i k f˛i j Q˛jk + hij , (60)
‍ ˛=1 k=in ‍

where the flow terms, ‍Q˛jk = Tjk k pj pk + ˛ Hjk ‍, here are negative. These equations can be linearized in terms of Nc − 1 compo-
sition variables. If r denotes timestep number, then the new component volume fraction over a step is:
c 1
Np P NP
P
Sr+1 r
ij = Sij
t
Vj Srlk Sr`j @f˛i
| Q
@S` k ˛jk
+ Vj hij .
t
(61)
‍ ˛=1 k=in `=1 ‍
Clearly stability depends on the matrix of fractional flow derivatives. This adds a complication which does not appear in simple two-phase
immiscible flow problems.
Rather than using inflows, Russell (1989) suggests the problem can be analyzed in terms of outflows, which is more convenient. Eq.
55 indicates they should be roughly equivalent. Russell (1989) suggests an L∞ norm or row sum criteria can be used for the stability con-
dition for each component balance, so for a particular node stability requires:
ˇ !
c 1 ˇ ˇ
@f˛i ˇˇ
Np P
P NP
Cj = Vj max Q˛jk ˇ
ˇ @S` j ˇ < 1,
t
(62)
i ˛=1 k=out
‍ `=1 ‍
where for outflows, the values of ‍Q˛jk ‍are positive. The row sum gives an estimate of the eigenvalue of the matrix. However, that appears
to be too conservative. YR used the maximum diagonal in place of the row sum. Later, after the publication, other pathological cases were
encountered causing adoption of a method using the maximum eigenvalue of the matrix. With this approach, the condition is:
Cj = t
Vj ƒj < 1, (63)
‍ ‍
where Λj is the maximum eigenvalue of the matrix:
22 2022 SPE Journal

P Np P
fi`0 = Q˛jk @f ˛i
@S` j
|, (64)
‍ j ˛=1 k=out ‍
for i = 1,…, Nc − 1. One must include all the outflows in Eq. 64, especially when the flow is countercurrent. Dramatic changes in stability
occur when the flow changes from countercurrent to cocurrent (Young and Russell 1993).
The results with this approach are the same as reported by others using von Neumann stability analysis for the simple case of three
phases without component partitioning (Coats 2003a; Wan et al. 2005). For that special case, the component fractional flow and volume
fractions are identical to the phase fractional flow and saturations, so the derivatives in Eq. 64 are the usual ones for an immiscible dis-
placement. The eigenvalues are :
q

0 f 0 2 + 4f 0 f 0
ƒ = 12 max f110 +f0 ˙
22 f11 22 21 12 . (65)
‍ ‍
For this case of immiscible flow and neglected capillary pressure, the result is the same as Eq. 3 of Coats (2003a) and Eq. 21 of Wan et al.
(2005). However, it is more general because it can be applied to problems with more than three components and phase partitioning with-
out simplifying assumptions like constant K-values. For multiple components, the maximum eigenvalue can be estimated with a few
iterations of the power method (Press and Vetterling 1999). The approach was implemented in the code of YR (Young and Russell 1993)
and proved effective for problems without capillary pressure. The failure to consider capillary pressure was an important oversight
because it has a significant effect on stability in some problems. The approach can and should be extended to include capillary pressure
effects. Additional theoretical justification is also warranted.
YR also found that stable results were often achieved with a CFL limit of 2 instead of 1. Russell (1989) notes that eigenvalue analysis
of the equations can lead to a limit of 2, which he states is erroneous. Coats (2003a) shines further light on the subject and finds that a limit
of 2 can lead to solutions which oscillate, but the oscillations die out. On that basis, the solutions are stable. This type of behavior is fre-
quently observed for timestepping methods, such as the trapezoidal rule which are classified as A-stable (absolutely stable) even though
they may oscillate.
The equations used here to analyze stability are based on volume fractions, whereas the equations are normally solved with mole or
mass fractions using the coefficients of Eq. 44. Continuing with the assumption of constant partial volumes and given the coefficients of
the linearizations, Eq. 44, the matrix of derivatives in Eq. 64 can be approximated with the aid of Eq. 12:

@f˛i i @nj i @ f˛ x˛i
= @
@nj h˛
f˛ x˛i = V
@Sj j @nj
P 0 i
@Sj ˛
(66)
V x˛i 0
i
f˛ b˛ij + f˛ ˇ ˇ j ,
‍ ˛ j ˛j ‍
0 P @kr˛
where ˛j = 1 a and a and b are the coefficients of the linearizations in Eq. 44. The variation of specific volume and viscosity
˛ @Sˇ ˇ j
‍ ˇ ‍
are neglected as in Eq. 42. The assumption of constant phase-specific volume or density in Eq. 66 can easily be rectified as described in
the Appendix.
Note how the partial volumes enter the calculation of the matrix, Eq. 66. Supplementary file 1 and Figs. 2 and 3 show that the partial
volumes often differ from one component to another by a factor of 10 or more, and they can even be negative. This difference should be
accounted for when evaluating stability.
The problem is formulated with a pressure equation constructed using a sum of component balances weighted by partial volumes, just
like IMPES, but it has dependence on the component mole numbers for the implicit blocks. Component balances are used for only the
implicit nodes, so the matrix problem looks like:
Dp p + Dn n = Dr and
(67)
‍ Ep p + En n = Er , ‍
where Dp represents the normal pressure equation with a few pressure terms from the linearizations, Eq. 44. Dn embodies the influence of
composition-dependent flow terms on the flow field. This dependence should be weak. From Eqs. 17 and 21, it is apparent that the strength
of the coupling depends on how the characteristic pressure, p, is defined. The global pressure ideas of Chavent and Jaffré (2014) can be
used to reduce the influence of capillary pressure. Simple truncation error analysis shows the strength of the coupling also depends on the
way upwinding or upstream weighting (USW) is applied (Young 1984, 2001). Other coupling issues are discussed by Moncorgé et al.
(2019).
Ep and En represent the component balance equations for only the implicit unknowns, ‍ n‍. For each implicit block, Nc − 1 component
balances must be solved, so for a typical case with 1% of the blocks implicit and 8–10 components, the dimension of ‍ n‍is less than 10%
of Δp in Eq. 67.
Young and Russell (1993) did not solve Eq. 67 efficiently, because they used an off-the-shelf solver which took no advantage of the
matrix structure. The solution scheme should account for the knowledge that the crux of the problem is in Dp because that is the elliptic-
like part from IMPES. En comes from the nearly hyperbolic transport equations, so it is nonsymmetric. Also, the implicit points frequently
consist of disjoint regions around wells, so En is often reducible.
One could consider an iterative sequential scheme, where the first equation is solved for pressure assuming the compositions are
known, then those pressures could be used in the second equation to update compositions. This could be stopped after one iteration like a
classic sequential implicit procedure (Spillette et al. 1973), or it could be repeated until convergence. Recent work on sequential methods
is described by Moncorgé (2018), Moncorgé et al. (2019), Møyner and Tchelepi (2018), and Sheth and Younis (2017).
AIM, quasi-Newton, and sequential methods are related in the sense that they all try to exploit the nature of the equation coupling and
the local need for implicitness. Relative to a FIM solved by a Newton-Raphson iteration, AIM neglects Jacobian and residual terms at
explicit points. Quasi-Newton methods neglect some Jacobian elements but include all residual terms. An iterative sequential method can
also be viewed as a type of quasi-Newton method; all residuals are included but the coupling terms, Dn, are neglected in the Jacobian.
Several of the best iterative methods for FIM take advantage of the special structure of Eq. 67. These include the combinative method
(Behie and Vinsome 1982), constrained pressure residual (Wallis et al. 1985), and multiscale methods (Lacroix et al. 2003; Stueben et al.
2007; Al-Shaalan et al. 2009). Similar strategies are described for adaptive implicit methods (Watts and Rame 1999; Clees and Ganzer
2022 SPE Journal 23

2007). These methods use various procedures (usually algebraic) to extract the elliptic-like part of the problem. Direct computation of Dp
and Dn using partial volumes is an inexpensive alternative.
A robust and efficient implementation of AIM is more complicated than meets the eye. Originally, YR envisioned that implicit points
and flow directions would be specified at the beginning of the timestep and would not change until the next timestep. Unfortunately,
testing showed a more complex approach is required to achieve the stability of FIM. After extensive testing, several of their conclusions
were obvious and remain unchallenged:
1. For compositional simulations, optimum results usually occur when a small number (. ‍ ‍ 2%) of blocks are implicit. Taking large
timesteps, comparable to FIM, is less accurate and less efficient.
2. Both threshold changes and stability must be used to monitor connections for implicit treatment.
3. Only heterogeneous simulations with large variation in throughputs, Q/V, show large efficiency improvements over IMPES.
4. Iteration is required for the method to be effective. During the iterations, the correction of flow reversals and the ability to add new
implicit connections are essential.
For problems amenable to IMPES, YR show how stability calculations can be used to select timesteps, thus improving the reliability of
IMPES. Coats (2003a) points out a discrepancy between his results and those of YR. The IMPES version of the YR model treated the
production well terms implicitly like the method originally proposed by Spivak and Coats (1970). Although it requires more effort,
implicit injection terms are likely worthwhile. Coats did not state how his implicit wells were implemented; if he also treated the injection
well terms implicitly, it would explain the factor of 2 discrepancy he describes.
Differences in IMPES well terms could also explain some of the differences between the results of YR and those of Collins et al.
(1992). For the third comparative solution problem, Collins et al. report a maximum efficiency relative to IMPES of 2.4, whereas YR
reported 1.5. Collins et al. do not state whether implicit rates were implemented in their IMPES calculations. If the YR results are com-
pared to earlier work with explicit rates (Young 1987), the optimum efficiency would be 2.2 for that problem, which is closer to the
optimum reported by Collins et al.
Young and Russell (1993) show an example in their Fig. 8 where only connections around a well require implicit treatment. One could
envision treating all well terms and perhaps adjacent connections implicitly as a matter of course. Alternatively, some other method could
be devised to determine subdomains of the reservoir to be treated implicitly. Using implicit subdomains which are static or only occasion-
ally updated would be simpler to code than AIM and reduce the computational burden of FIM. If this approach were combined with
timestep selection using stability analysis, the method would likely work quite well. The author implemented an extended implicit region
around wells in an earlier, simpler model (Young 1981). It is much easier to implement FIM (the bigger hammer approach), even though
it may be slower.
Alternative Fluid Property Representations. The discussions so far have focused on improved efficiency and reliability with
algorithmic improvements using conventional formulations of an EOS. There are also potential improvements by modifying the fluid
property representation. One type of modification is the so-called reduced variable methods, where the EOS is formulated in terms of a
smaller number of parameters. Another modification is essentially tabular, called the tie-line-simplex parameterization.
Reduced variable methods originated with the work of Michelsen (1986). He shows that for problems without binary interaction
parameters, the VLE calculations can be performed with only three variables, regardless of the number of components considered. The
method has been generalized and extended in several ways and the variants are reviewed and compared in Wang and Barker (1995),
Haugen and Beckner (2013a), Michelsen et al. (2013), Gorucu and Johns (2014), and Petitfrere and Nichita (2015). One problem with
these studies is they typically consider flash and stability calculations which use no prior information, so they do not accurately reflect the
computations in a compositional simulation. Although there is conflicting evidence, the consensus seems to be that only modest gains can
be achieved unless the number of hydrocarbon components is fairly large, Nc > 10–15. The explanation was first given by Wang and
Barker (1995). They point out that true flash calculations are not used frequently in some simulators. One often speaks of the fluid prop-
erties being dominated by flash calculations, used as a generic term. To achieve significant gains, one must perform all the calculations in
terms of the reduced variables, including tests for phase stability and calculation of partial molar volumes and additional linearizations for
FIM or AIM. Stability testing shows the greatest potential for improvement (Petitfrere and Nichita 2015). However, there will always be
some overhead because the results must ultimately be expressed in terms of pressure and compositions (e.g., Eqs. 44 and 51).
For the IMPES case and an MP approach, the equations are given by Eq. 46. Implicit calculations also require the linearizations, Eq.
44. For a true flash calculation, only the first two of Eq. 46 are needed and the pressure and feed are specified, so ‍ = p = 0.‍ In a simu-
lator, a true flash calculation is needed only for new two-phase cells. For an entire simulation, a cell will typically transition to two phases
only once, for example, owing to drawdown below the bubble point. For cells that were two-phase in previous time steps, the full set of
Eq. 46 is solved together with the differential material balances. For negative flash cells, the same equations are used, but the volume
linearization is not needed. For each time step, the values are updated with back substitution, new residuals are calculated, and then each
cell is checked for convergence. Typically, more than 97% of the two-phase cells will pass the convergence test without iteration. For the
remainder, a full relinearization or Jacobian update is calculated, and an inner iteration is performed. Because the Jacobian is updated
infrequently, a reduced variable method is not likely to provide a significant advantage. In addition to the dissimilarity of tests to actual
simulator calculations, efficient implementation is undoubtedly important for a valid evaluation (Haugen and Beckner 2013a). One study
suggests a good simple approach is to exploit the sparse nature of binary interaction coefficients in the calculation of the EOS “A” param-
eter (Michelsen et al. 2013). That procedure was implemented in our simulator in 1989 (Young 1989) and gave roughly a 3% reduction
in computation time. The conclusions reached here are for the case of two hydrocarbon phases. Reduced methods could very well be more
beneficial with the three-phase systems discussed in the section “Three Hydrocarbon Phases.”
Another method for reducing the burden of fluid property calculations uses a tabular approach. The use of tabular correlations has a
long history in compositional modeling. Black oil models invariably use tables. The early compositional models used tabular K-values.
When many components are used, it is difficult to formulate that approach and achieve consistency at a critical point (Nolen 1973; Rowe
1978). An early paper (Van-Quy et al. 1972) achieved consistency by restricting the treatment to ternary systems using Hand’s Rule
(Treybal 1968, p. 29–30).
As mentioned above, my assignment at Amoco in 1980 was to develop a model for CO2 flooding. With the limited computing
resources and the questionable accuracy of an EOS, some sort of tabular approach seemed most promising. In addition to the work of
(Van-Quy et al. 1972), ternary treatments were used for method of characteristic and other studies of gas displacements at Shell (Helfferich
1981; Gardner et al. 1981). To test feasibility, the method was implemented in the generalized compositional model at Amoco (Young and
Stephenson 1983). A ternary approach was used to approximate two different 17-component Redlich-Kwong descriptions, a CO2-oil
system, and an N2-gas condensate system. A 1D constant pressure displacement was used to develop tabular data for the ternary
24 2022 SPE Journal

representation. In both cases, computations with the ternary systems agreed well with the 17-component calculations. The computation
time was reduced to that of a black-oil model. This approach seemed like an excellent middle ground between the modified black-oil FCM
approach (often called a four-component model) which neglects phase behavior (Lantz 1970) and an EOS-based model. Unfortunately,
the idea received little support at Amoco, so was not developed further. Similar ideas were used to develop a limited compositional model
at Arco (Tang and Zick 1993).
The extension of the tabular approach beyond four components is nontrivial. It was successfully accomplished years later and called
the tie-line-simplex method (Voskov and Tchelepi 2009a, 2009b). It has been applied primarily as a method to speed up calculations,
while retaining the same number of components. The tabular data are sometimes created on-the-fly from computations with a conven-
tional EOS. The method has found application in thermal simulation (Iranshahr et al. 2010b, 2013). The method has also been used for
problems with three hydrocarbon phases (Petitfrere et al. 2020), see discussion in the section “Three Hydrocarbon Phases.” One can
achieve only so much efficiency when using many components. The flexibility of a tabular approach should permit a significant reduction
in the number of components. Some work in this direction is reported in (Ganapathy and Voskov 2018). There are potentially some advan-
tages to the method, but more work is needed.
Efficient Implementation. An efficient implementation of a method often has a greater effect than the algorithm used. As an example,
CPU times for the first comparative solution problem (Odeh 1981) were given in an oral presentation. The times varied by about a factor
of 20 even for this simple black oil problem. The fastest and slowest simulators were both based on the IMPES method. In between were
simulators using other solution algorithms (e.g., FIM, sequential implicit). The least amount of time was required by the Amoco black
oil model (Sheffield 1969). It was improved significantly in unpublished work by Arden McCracken. The original code was improved in
several ways, but the greatest improvements were achieved by more efficient table lookup calculations. For example, rather than call the
lookup function once for each cell, it is more efficient to move the loop inside the procedure so that all lookups are performed with one
call. The improvements had nothing to do with the basic solution algorithm. Producing an efficient implementation is the hard work which
does not lead to publications, fame, and glory.
The details of efficient implementation are not normally discussed for reservoir simulation. Obviously, every effort should be made to
avoid unnecessary computations. The calculation rate is also important, and it depends, to some extent, on the specific computer architec-
ture. The table lookup issue above is a simple example where function call overhead is reduced by ensuring that the function does signif-
icant work. In our original work (Young 1991), a key issue was to exploit the vector processing hardware of supercomputers. Current
microprocessors have vector processing capabilities and other forms of low-level parallelism. Some of the programming issues for newer
microprocessors are discussed in Haugen and Beckner (2013b) and Gandham et al. (2016).
For an efficient EOS compositional model, one of the most important requirements is that the fundamental properties, for example,
molar volume, fugacity, and their derivatives, be calculated efficiently. Efficient methods for these calculations are well documented
(Michelsen and Mollerup 1986, 2004). It is likely some simulators do these calculations inefficiently. There were numerous inefficiencies
in the original implementation of our simulator (Young 1987) even though it was significantly more efficient than others (Killough and
Kossack 1987). That experience does not appear to be unique. Working with a different simulator, Belkadi et al. (2013) report significant
gains owing to efficient implementation of basic calculations.
Table 9 gives results of some of the modifications to the implementation reported in Young (1987). The first two improvements were
reported in Young (1991), but the rest have not been presented in the open literature. These improvements occurred over a period of a few
years, and owing to computer changes, it was not possible to compare original vs. final results on the same computer. The improvements
were both problem- and computer-dependent, so a typical average is listed. Also, computers have evolved considerably since then. The
results are presented to give a general indication of the magnitude of the improvements one might expect, because this type of information
is sorely missing in the literature.
The first item in Table 9 lists the improvements owing to the implementation of implicit production rate terms. This procedure is sim-
ple to implement, and its cost is insignificant. For the third and fifth comparative problems, the number of timesteps were reduced by an
average factor of about 2. The original code had some bottlenecks in the method used to track neighboring cells for the neighborhood
method and in the separator calculations. The average speedup reported in Young (1991) was about 2.6. As discussed above, exploitation
of sparse binary interaction parameters and stability test improvements yielded an overall speedup of almost 3 relative to the original code.
Improvements to the Jacobian calculations and negative flash calculations netted additional gains to produce an overall improvement
greater than 4. Keep in mind that these improvements were for a simulator with demonstrated high performance (Killough and Kossack
1987). For a simulator which is initially less efficient, much larger gains can be achieved. Differences of this magnitude are one reason
one must be careful comparing algorithms with direct calculations — implementation is key.
One difficulty with implementation of EOS fluid property calculations is that, depending on the conditions at a particular cell, there are
a variety of different calculations that are required. In Young (1991), the term fluid type was used to designate the required calculations.
In that implementation, five fluid types were defined: (1) two-phase points not requiring a full Jacobian update, (2) two-phase points
Incremental Overall
Normalized Time Speedup Speedup
Original (Young 1987) 1,000
Implicit rates (Young 1991) 500 2.00 2.00
Separators, neighbor tracking (Young 1991) 385 1.30 2.60
Sparse binary interaction parameters 375 1.03 2.67
Symmetric stability test 347 1.11 2.96
Jacobian updating 297 1.14 3.37
Symmetric Jacobian, Eqs. 46–51 291 1.02 3.44
Negative flash 236 1.23 4.23
Table 9—Implementation improvements.
2022 SPE Journal 25

requiring a Jacobian update, (3) single-phase points tested for phase appearance, (4) single-phase points not tested, and (5) aquifer points
with no hydrocarbon components. Later, the negative flash implementation required an additional type. The variety of calculations can
lead to inefficient code with many if tests to branch around unneeded calculations. An approach called stripmining was used by Killough
and Levesque (1982). With that approach, one scans the grid for contiguous cells requiring a particular calculation, then the calculations
are performed on each group with vectors of those contiguous cells. Others use gather and scatter operations to reorder the cells based on
fluid type (Chien et al. 1987; Subramanian et al. 1987; Malachowski et al. 1990).
The implementation in Young (1991) maintained two orderings of the cells—the normal grid-based order in the main simulator and a
fluid type-based order for fluid property calculations. The data supplied to calculate the fluid properties (i.e., the overall compositions,
pressure, etc.) and the data calculated (i.e., saturations, phase compositions, partial volumes, etc.) are much smaller than the data required
to perform the calculations (i.e., the Jacobian coefficients, Eq. 46, etc.). Typically, for each node, the number of inputs is 15—20, the
number of outputs is 30—40, while internally 200—250 quantities are required. The number of outputs is larger if linearizations, Eq. 44,
are required for FIM. With this approach, gather/scatter operations are required for only the inputs and the outputs, but not the internal
data. When a cell changes from one type to another, simple swaps are used to maintain the ordering by fluid type. The CPU time spent
moving data was only 2 to 4% of the computation time (Young 1991), so the procedure is effective even without vectorization. Except for
the swaps, all calculations are performed on consecutive memory locations, that is, with stride one, which is important for both scalar and
vector calculations.
A primary motivation for the dual ordering was to achieve efficient vector calculations. Vector calculations were essential for super-
computers of the 1980s and 1990s but also are important for newer microprocessors (Haugen and Beckner 2013b).
The ability to perform calculations in single precision (32-bit arithmetic) is also important. The solution of many problems would not
have been possible in double precision owing to memory limitations (Wendschlag et al. 1983; Hsu et al. 1997). Memory is no longer the
restriction it once was, but single precision is still important for computational speed. The transfer of data is always slower than calcula-
tion speed, whether it is from a processor’s central memory, from another processor in a distributed memory computer, or between a CPU
and GPU (graphical processing unit). Single precision vector calculations on CPUs are faster by a factor of 2. The difference on GPUs is
even greater (Gandham et al. 2016). If properly implemented, fluid property calculation can be performed in single precision. Modern
languages contain facilities to make precision changes easy. The simulator should be coded using this feature to allow easy precision
changes.
All of the calculations in Young (1991) are performed with vector lengths equal to the number of cells of a particular type. This is not
the most efficient procedure for most modern computers. Most of the fluid property calculations are nested. If performed with looping
operations, the innermost calculations are with vectors of grid cells, but they are nested inside loops over components. All modern com-
puters have cache memory. When data are available in the cache, it loads much faster than from the main memory. For example, to cal-
culate the EOS “A” parameter, requires a calculation like:
NP
c 1 NP
c 1
AL = aL ij xi xj . (68)
‍ i=1 j=1 ‍
If these parameters are to be calculated at a number of points, it is best if the number of points is not too large. There is a lot of data reuse
in the calculation. The calculations will execute faster when all the data fit within the cache. Haugen and Beckner (2013b) describe other
strategies for cache optimization. GPUs have large register files making locality of data references even more important (Gandham et al.
2016). Locality of data references is clearly important for an efficient implementation.
High Performance Computing—Parallel Processing. Because of the intensity of EOS computations, they have long been a target for
high performance computing. The computing techniques used have varied along with available computing hardware. In 1980, a typical
computation rate for an IBM mainframe was about 0.002 GFLOPS (billions of floating point operations per second), more than 1,000
times slower than a typical 2020 laptop running in serial mode. Some companies purchased Cray supercomputers in the early 1980s
(Killough and Levesque 1982). Even before that, there were attached array processors, which were separate hardware devices which
could be attached to a mainframe (Johnson 1970; Woo 1979). These devices were difficult to use, required nonportable code, and suffered
from a relatively slow transfer rate from the mainframe. In many ways, current GPUs are reminiscent of early array processors.
Supercomputers were a significant breakthrough, partly because they could perform ordinary scalar calculations three to four times
faster than other computers. If the vector processing hardware was brought to bear, the calculation rate could be increased another factor
of 5 for the early Cray—1 to 20 or more for later models from Cray and other vendors. However, most simulation software was not written
to take advantage of the faster vector processing hardware.
The literature on parallel and high-performance computing for compositional simulation is vast. We will touch on the high spots of the
overall development history but not try to cite every article. Many of the early articles concentrated on iterative methods for solving the
linear algebraic problems (Nolen et al. 1981; Scott et al. 1987; Killough and Wheeler 1987; Wheeler and Smith 1990). The literature on
linear solvers is not reviewed here. The consensus is that the best solvers for distributed memory parallel processing are based on domain
decomposition, where the problem domain is divided into several subdomains, one for each processor.
After extensive coding for vector processing, Nolen and Stanat (1981) achieved a speedup of 3.7 on a CDC Cyber 203 for a black oil
simulator using FIM. Early work on vector processing in compositional models is discussed briefly in the previous section. With the
simulator design described above, a speedup of 16 and a calculation rate of more than 0.15 GFLOPS was achieved on a single processor
Fujitsu machine (Young 1991). An expanded nine hydrocarbon components version of the SPE third comparison problem with almost
10,000 cells was solved in just 75 CPU seconds. Later calculations with a black oil formulation achieved a speedup of 10.5 and 0.12
GFLOPS on a single processor of a Cray Y-MP (Young and Hemanth-Kumar 1992). The black oil simulator was demonstrated by solving
for the first time a representative waterflood with one million cells. The full 20-year history required only 26 CPU minutes.
Once these calculation rates were achieved, memory became the dominant limitation. In the early 1980s, a Cray with 4 Mwords (32
Mbytes) of memory was rare. The compositional simulations for the Anschutz Ranch reservoir study with over 8,000 nodes required 3.75
Mwords or 15 Mbytes single precision (Wendschlag et al. 1983). It was run overnight on an IBM mainframe with 16 Mbytes of virtual
memory because there was not a commercially available Cray with sufficient memory.
Later supercomputers had more memory which was shared by multiple processors. Problems could be run in parallel on these proces-
sors. If all the available memory was needed to run a problem, it makes sense to use all the processors too. This was relatively easy to do
using compiler directives. If a code was written for vector processing, the loops could be broken up so that each processor calculates part
26 2022 SPE Journal

of a loop or nest of loops. Parallel/vector calculations are discussed by Chien et al. (1987). This works reasonably well provided the code
is highly vectorized, the problem is large, and there are not too many processors.
To understand these and other calculations, we need to introduce some terminology. First, Amdahl’s law applies for calculations where
one portion of the calculation is in parallel while the remainder are serial. Amdahl’s law is:
t1 N
= 1 = proc , (69)
tN
‍ fs + 1fs /Nproc 1+fs Nproc 1 ‍
where t1/tN, is the speedup, Nproc is the number of processors, fs is the serial fraction. Amdahl’s law can be applied to vector processing if
Nproc is replaced with the speedup for 100% vector calculations. Also, many use the term scalability. This term is thrown around loosely,
and usually no concise definition is given. We define a calculation to be scalable if:
lim t1 = E1 Nproc , (70)

‍ Nproc !1 tN ‍
where E∞ is a constant representing the efficiency with a large number of processors. The efficiency from Amdahl’s law is:
t1 1
EN = tN Nproc = . (71)
‍ 1+fs Nproc 1 ‍

Using the equation for Amdahl’s law, scalability requires fs Nproc = 1 E1 /E1. So, according to Amdahl’s law, the serial
lim
‍Nproc !1 ‍
fraction must be inversely proportional to Nproc, so it must approach zero in order for a calculation to be scalable. Using Eq. 71, the per-
missible serial fraction to produce a given efficiency is plotted in Fig. 7. For a serial content of 1%, an efficiency of 80% is produced with
27 processors but only 40% with 150 processors. To achieve 80% efficiency with 250 processors, the serial content must be 0.1%. This
level is difficult to attain for a reservoir simulator.
Fig. 7—Efficiency and serial fraction, Amdahl’s law.
Amdahl’s law assumes inefficiency is due exclusively to serial calculations. Based on Eq. 70, an algorithm can be scalable and not be
100% efficient. For example, suppose blocks are divided between the processors for fluid property calculations. Some processors may
have mostly aquifer blocks, while others may have near-critical fluids. The aquifer cells require less computation, so those processors will
remain idle while waiting for the others to finish, causing a loss of efficiency. However, in the limit, the speed will still double when the
number of processors is doubled. This example illustrates inefficiency owing to load imbalance between the processors, discussed in a
number of articles (Killough et al. 1996).
These concepts of serial content and scalability are key to the evaluation of articles on parallel processing and high-performance com-
puting. For example, Hemanth-Kumar and Young (1996) describe black oil and EOS simulations on a Cray C90 with 16 processors and
256 Mwords of memory. They reported complete simulations with their million cell black oil problem and an EOS problem with 11
hydrocarbon components and up to 580,000 cells. The largest problems were solved in only 2.5 and 9 minutes for black oil and compo-
sitional problems, respectively. With the compositional model, speedups of 18 and 81 were achieved for vector (1 processor) and vector/
parallel (16 processors) computations, respectively. An overall computation rate of 1.8 GFLOPS was reported for the largest EOS prob-
lem. Superficially, these results seem quite good; however, let us look in more detail.
Vector processing on the C90 used two vector pipelines, which could produce up to four floating point calculations per clock cycle. By
examining calculations which are 100% vector, one concludes that the vector to scalar ratio is about 22. Using Amdahl’s law with the
observed speedup of 18.3, the calculated serial content is 0.95%. For the black oil calculations, the observed speedup was 17.4, so the
serial content is 1.26%. The serial parts of that simulator primarily involve well calculations and input-output, which are difficult to vec-
torize. To be realistic, the largest test problems had 500 to 800 wells. Test problems in some articles have only two or a small number of
wells, so they are not representative.
2022 SPE Journal 27

The 16 processor Cray C90 is equivalent to 352 (16 × 22) scalar processors. With all 16 processors, the observed speedup of 81 relative
to scalar calculations is exactly what one would predict using Amdahl’s law and a serial content of 0.95%. The overall efficiency of only
23% is to be expected from Fig. 7. These results are not scalable according to Eq. 70, because the small serial content is constant and
becomes limiting with a large number of processors. Many authors claim their results are scalable when their results show a similar con-
stant serial content.
Such a low efficiency on an expensive supercomputer is not workable. By the early 1990s, supercomputers were being supplanted by
high performance workstations for day-to-day simulations (Young et al. 1990). The Cray C90 is an example of a computer with a rela-
tively small number of vector processors and shared memory. Workstations with multiple processors and shared memory were also intro-
duced. Now we have microcomputers with multiple cores and vector processing capability, some with attached GPUs.
The shared memory concept is not viable for a large number of processors owing to the difficulty of building a robust memory subsys-
tem which can be accessed by all the processors without slowdowns. To build computers with a truly large number of processors, memory
is attached locally to each processor. The transmission of data between processors is slow compared to the computation rates when data
are in the processor’s local memory. The same distributed memory approach can be used with a network of workstations or PCs. These
distributed memory computers are simpler/cheaper to build than a supercomputer. However, they are more difficult to program and com-
munication between processors must be minimized.
The vectorization/parallelization strategy, discussed above, for shared memory machines is a form of SIMD (Single Instruction
Multiple Data) parallelism. For distributed memory machines, a more general MIMD (Multiple Instruction Multiple Data) strategy is
used. With SIMD, each processor performs the same calculation, whereas with MIMD, each processor can do different calculations.
SIMD is better suited for low-level, fine-grained, or loop-level parallelism, whereas MIMD is better suited for a high level, coarse-grained
approach. EOS calculations are naturally parallel with an MIMD approach using domain decomposition. For example, with MIMD, one
processor might be performing an EOS stability test while another is calculating relative permeability. A SIMD approach is more difficult
to implement for EOS calculation owing to the fluid types (see discussion above), that is, the calculations are not the same for all points.
Some reorganization is needed to perform the calculations efficiently and different strategies are described (Chien et al. 1987; Subramanian
et al. 1987; Young 1987; Malachowski et al. 1990).
True scalability cannot be achieved in reservoir simulation when a problem of a given size is run on multiple processors. A form of
scalability, called weak scalability, can be achieved with domain decomposition methods if the problem size grows in proportion to the
number of processors and subdomains are constructed so the number of shared interface cells does not grow faster than the number of
processors. Even for this weak form of scalability, these restrictions are severe. The requirement of problem growth is not unrealistic,
because often a larger problem is dictated by a desire to simulate a larger reservoir area or to achieve greater detail. If the problem is large
enough, true scalability is not necessary.
Some articles describing parallel calculations for full simulators and not just linear solvers are as follows:
• Black Oil — (T. van Daalen et al. 1989; Meijerink et al. 1991; Chien et al. 1997; Shiralkar et al. 1998; Collins et al. 2003; Fjerstad
et al. 2005; Dogru et al. 2009; Beckner et al. 2015; Liu et al. 2015b; Guan et al. 2015; Vidyasagar et al. 2019)
• Compositional — (Killough and Bhogeswara 1991; Wang et al. 1999)
• Both — (Killough et al. 1996; Hemanth-Kumar and Young 1996; DeBaun et al. 2005; Bogachev et al. 2018)
These articles span several years and consider a variety of computers. Some plot their results in deceptive ways and all fail to provide
sufficient details concerning their test problems. Even though true scalability is not possible, they almost universally claim their results
are scalable. Detailed results with a large number of processors is rarely given. Some seem to be overly promotional even when from
universities and major oil company research labs.
Table 10 gives the serial content from several publications to evaluate how truly scalable these simulators are. The values in the table
are calculated with Amdahl’s law directly from results presented in the papers. The results are not intended to disparage any work, but
rather to convey a clearer picture of the scalability achieved. The easiest way to evaluate results is to plot tNNproc vs. the number of pro-
cessors. If the results are linear, then they obey Amdahl’s law, and the scalar content can be determined from a linear fit. If the results are
scalable, the slope should approach zero for a large number of processors. The results were invariably linear within the accuracy of the
data, so a least squares fit was used to determine the values for serial content listed in Table 10. Also, for domain decomposition methods,
the solver characteristics change with the number of processors, making the results difficult to interpret. In most cases, results were pre-
sented for a fixed problem size, but in a few cases, the problem size was proportional to the number of processors. Some articles do not
present their results in a form amenable to this type of analysis.
Serial
Content Calculation Computer Reference
1.46% Solver Meiko (Meijerink et al. 1991, their Table 3)
1.26% Black oil Cray C90 (Hemanth-Kumar and Young 1996, their Fig. 3)
0.95% EOS Cray C90 (Hemanth-Kumar and Young 1996, their Fig. 6)
1.54% Black oil IBM SP2 (Killough et al. 1996, their Fig. 14)
0.59% Black oil IBM SP2 (Chien et al. 1997, their Table 1)
1.47% Solver Cray T3D (Shiralkar et al. 1998, their Table 1)
2.27% Black oil IBM SP2 (Shiralkar et al. 1998, their Fig. 4)
0.38% Black oil Cray T3D (Shiralkar et al. 1998, their Fig. 5)
1.15% EOS IBM SP2H2 (Wang et al. 1999, their Fig. 2)
4.08% Black oil IBM NH2 (Dogru et al. 2002, their Fig. 2)
14.19% Black oil IBM P690 (Collins et al. 2003, their Table 3)
7.46% Black oil IBM P690 (Collins et al. 2003, their Table 4)
Table 10—Serial content for various parallel simulators.
28 2022 SPE Journal

Serial
Content Calculation Computer Reference
1.27% Solver Cluster (DeBaun et al. 2005, their Fig. 12)
0.97% Black oil Cluster (DeBaun et al. 2005, their Fig. 12)
1.58% EOS Cluster (DeBaun et al. 2005, their Fig. 17)
7.17% Black oil Cluster (Fjerstad et al. 2005, their Fig. 14)
1.27% Black oil Cluster (Fjerstad et al. 2005, their Fig. 18)
0.11% Black oil Discovery (Beckner et al. 2015, their Fig. 2) (128 node)
0.12% Black oil Tianhe-2 (Guan et al. 2015, their Table 1)
Table 10 (continued)—Serial content for various parallel simulators.
Although the analysis is admittedly flawed, it is still worthwhile, because it gives a rough idea of the scalability achieved, and the
comparison is on a more uniform basis. Only large problems are considered, because smaller ones are not representative. The serial con-
tent for most of these scalable formulations are greater than 1%, which is not scalable. Even with a 0.5% serial content, the calculations
achieve only 50% efficiency with 200 processors. Although the earlier simulators were not scalable, several later ones approach scalability
when the problem is large (Dogru et al. 2009; Beckner et al. 2015; Liu et al. 2015b; Guan et al. 2015). Although scalability is necessary
for efficient calculations, it says nothing about the efficiency of the calculations on individual processors, which is the subject of the pre-
vious section.
The most efficient way to perform compositional simulator computations has always been a moving target. With the advent of super-
computers, low-level parallelism was achieved with vector processing or SIMD. With the continued development of microprocessors,
distributed memory computers were used with an MIMD strategy using domain decomposition, but the calculations on the processors
were scalar. Owing to physical limits, microprocessor speed improvements began to plateau in the early 2000s. For this reason, virtually
all microprocessors now have multiple cores with built-in SIMD and/or vector instructions, so each one is equivalent to a share memory
multiprocessor supercomputer of the 1990s. Distributed memory calculations are needed only for the largest problems. Even when dis-
tributed memory parallel computations are used, they must exploit these low-level forms of parallelism to achieve overall efficiency. To
further confuse the picture, GPUs are now available. They can perform certain types of calculations very fast. One can perform all calcu-
lations on the main CPU cores (Beckner et al. 2015), some calculations on the CPUs, and others on GPUs (Fung et al. 2014; Bogachev
et al. 2018), or all the calculations on the GPUs (Gandham et al. 2016; Vidyasagar et al. 2019).
Programming for distributed memory systems is considerably more complex than for a normal computer. The complexity is com-
pounded by demands for unstructured grids, local grid refinement, faults, fractures, dual porosity, advanced wells, surface network inte-
gration and the incorporation of advanced physics for thermal processes, chemical reactions, and coupled rock mechanics. To incorporate
all of these features within one simulator designed for efficient parallel processing is a daunting task. Fortunately, with an MIMD approach,
many calculations are inherently parallel (e.g., rock and fluid properties and Jacobian calculation). Implementation of the domain decom-
position solver is the complicated part.
Many articles promote the concept of a generalized simulator framework which was the primary motivation of work in the early 1980s
(Young and Stephenson 1983; Acs et al. 1985; Watts 1986). The main idea is to separate the physics, that is, rock and fluid property cal-
culations, from the flow calculations. Modern languages support this object-oriented approach. Many articles describe the use of a frame-
work for the calculations (Parashar et al. 1997; Parashar and Yotov 1998; Fjerstad et al. 2005; Maliassov et al. 2013; Beckner et al. 2015;
Liu et al. 2015a; Guan et al. 2015; Qiao et al. 2017). Some of these incorporate more general public domain frameworks for the domain
decomposition linear solver (Amestoy et al. 2001; Li 2005; Heroux et al. 2005; Falgout et al. 2006; Balay et al. 2019; Notay 2019). With
a well-formulated framework, the complexity is isolated so that computation of the physics, like EOS calculations, can be implemented
without regard for the domain decomposition.
Remaining Work
Although compositional simulation technology has improved greatly over the past 50 years, there are a few areas where additional work
is needed. Here we discuss three remaining problem areas.
Three Hydrocarbon Phases. Simulations with an EOS encounter difficulties when the simulator performs VLE calculations for only
two hydrocarbon phases when three actually exist. Experiments describing multiphase phenomena were first reported in Shelton and
Yarborough (1977). My first exposure to the problem was in a group meeting at Amoco in the mid- to late 1970s. Del Fussell talked
about problems with his cell-to-cell model (Fussell et al. 1976). Most of us present did not understand what he was talking about, but
he kept saying “the EOS is trying to predict a third phase.” Later, Lynne Fussell developed a three-phase flash code to prove the point
(Fussell 1979). The associated difficulty posed for reservoir simulation was apparently first mentioned in the open literature by Young
and Stephenson (1983).
Okuno et al. (2010) give a description of the problem. When, according to the EOS, three phases exist, there may be two solutions to
a two-phase VLE calculation. If the path of the simulation, in composition-pressure space, traverses the three-phase region, at some point,
the VLE will switch from one two-phase solution to the other. For example, Lins et al. (2011) describe calculations where the three-phase
region is traversed by incrementing pressure, both increasing and decreasing. Starting at high pressure in the L1–L2 region, the calcula-
tions continue to produce liquid-liquid equilibrium until it abruptly changes to a condition where L2 disappears and a V-L1 solution is
found. For increasing pressure, the opposite trends occur. For CO2 systems, the actual progression is that the CO2 rich phase, V or L2, is
being vaporized or condensed while in the presence of the L1 phase. The specific volumes are different for the two solutions. In a simulator
with two-phase VLE, when the solution jumps from one two-phase solution to the other, the discontinuity causes an explosion like those
described by Nolen (1973). Owing to the jump discontinuity in fluid volume, cutting timesteps and other corrective actions often do not
work. On successive iterations, the solution may flip back and forth between the two solutions.
There are potentially two solutions to this problem, and both have received some attention. One could perform three-phase VLE cal-
culations or one could modify the EOS so it does not predict a third phase. The second approach is probably not reasonable for fluids like
the Kuparek fluid studied in Godbole et al. (1995), but it might be valid for the more prevalent CO2 systems where the three-phase region
2022 SPE Journal 29

is small and often well below the average reservoir pressure. For CO2 systems, the three-phase region is not usually present for tempera-
tures above 49°C (120oF). When present, it usually exists in a narrow pressure range about 8 MPa (1,160 psi) (Negahban and Kremesec
1992). Although there are undoubtedly differences in detail (Okuno and Xu 2014), laboratory data show no observable change in displace-
ment mechanism when the temperature is dropped to, say 40°C (104oF).
Reported attempts to recharacterize the fluid so the EOS does not predict three phases, while retaining accuracy for other parameters,
have not been entirely satisfactory (Whitson 1991; Fong et al. 1992; Pan et al. 2015). Because the fugacity coefficients are integrals of the
volumetric data (Prausnitz and Lichtenthaler 2004), it makes sense that it would be difficult to get the correct specific volumes or densities
while modifying the VLE. I have had experience with CO2-oil characterizations developed by different companies for various low tem-
perature West Texas oils (Hill et al. 1994; Hsu et al. 1997). Some seem to have no difficulty traversing a three-phase region, so there may
be techniques not described in the open literature.
The approach of implementing VLE calculations with three hydrocarbon phases has been considered by a few (Nghiem and Li 1986;
Perschke et al. 1989; Khan et al. 1992; Godbole et al. 1995; Guler et al. 2001; Okuno et al. 2010; Petitfrere et al. 2020). The implemen-
tations to date are research codes, sometimes restricted to 1D. The majority concern ongoing work at the University of Texas. Most arti-
cles concentrate on the flash and stability calculations and do not describe details of the formulations. Most use the IMPES approach, and
all use an EOS. Problems with simulator reliability have been reported, but no details have been described. Efficient implementation will
require that phase stability testing be restricted using some sort of shadow region. To date the reported computation times are
considerable.
In summary, a good solution to this problem is yet to be found. However, reduced variable methods show potential (Okuno et al. 2010)
as does a tie-simplex tabular approach (Iranshahr et al. 2010a, 2013; Petitfrere et al. 2020).
Relative Permeability. As discussed above, the early compositional models suffered from inconsistent properties. The use of an EOS
alleviated most problems. However, in many simulators, there remains an inconsistency in the treatment of relative permeability, Eq. 3.
For example, consider water flowing together with a hydrocarbon fluid having composition along Line C in Fig. 1. The usual convention
is to define absolute permeability to be the permeability at connate or irreducible water saturation, so ‍kr` = krv = 1‍ at irreducible water
saturation, although this condition may be inconsistent with experimental measurements. At higher water saturations, the relative
permeability of liquid and vapor are normally different. The hydrocarbon is obviously liquid when there is no C1 and vapor when there is
no C12, so oil/water and gas/water relative permeability would apply at the two extremes. The relative permeability should be continuous.
If the relative permeability is changed abruptly at some point along the line, it would not only violate physics but also would cause
difficulties for a simulator owing to the discontinuous flow coefficients.
The state of a single-phase fluid is defined by the relationship of the fluid critical temperature relative to the reservoir temperature: a
liquid Tc > T or a vapor Tc < T. This definition is strictly a convention and says nothing about the properties of the fluid and the physics of
flow. If the fluid composition lies on a tie-line extension, a negative flash calculation could be used to define the fluid state. However, for
the compositions along Line C in Fig. 1, there is no unambiguous way to define the fluid state.
This problem is sometimes called a phase labeling problem, but it is much worse. If the hydrocarbon starts out as a liquid, say the
mixture of C3 and C12 at the end of Line C, but then becomes more vapor-like owing to composition or pressure changes, relabeling it as
vapor would cause a discontinuity. However, if it is not relabeled and then later enters the two-phase region along the dewpoint line, there
is a discontinuity.
The inconsistency is also not a low IFT or high capillary number problem. It is well known that two-phase relative permeability curves
trend continuously to miscible displacements as the capillary number becomes large (low IFT). For surfactant flooding, IFT effects are
fundamental to the displacement mechanism. However, for MCM gas displacements, both numerical simulations and method of charac-
teristics theory shows it is not necessary to achieve high recovery. In the scenario described above, the fluid could enter the two-phase
region at a dewpoint far removed from the critical and with high IFT. Also, the IFT is not low between water and the fluid along Line C
in Fig. 1.
What about discontinuous flow coefficients? When simulator difficulties are encountered, corrective actions such as cutting timesteps
are usually employed. Such actions normally allow the simulator to continue. If the simulator manages to struggle through, the true source
of the difficulty is seldom determined. I encountered difficulties when first running the SPE fifth comparison problem (Killough and
Kossack 1987) with my new simulator (Young 1987). This problem uses a light fluid which is unlike typical reservoir fluids. Alternating
gas/water injection is used, so the pressure fluctuates causing many phase changes. I encountered problems at a node as the two-phase
region was entered. Backtracking revealed that at that point, the two-phase region was exited at a bubble point, but then, later it reentered
at a dew point. It was labeled as liquid, but when it reentered as a gas, a discontinuity in the flow coefficients resulted. These inconsisten-
cies are not as severe as the explosions caused by volumetric inconsistencies, but they can still cause problems. It is likely the impact will
be greatest for FIM owing to the larger timesteps and the use of derivatives of discontinuous relative permeability functions.
This inconsistency is just an obvious symptom of a larger problem, that is, that the simulator uses the wrong physics. The relative
permeability correlations in use at the time did not account for compositional effects (Stone 1970, 1973; Baker 1988). Clearly, relative
permeability must depend on composition explicitly, not just saturation.
Several empirical correlations incorporate a composition dependence. For consistency, these methods typically blend measured rela-
tive permeability curves for oil and gas using a single composition-dependent property. For the composition dependence, the correlations
have used: density (Blunt 2000; Fayers et al. 2000), parachor weighted density (Jerauld 1997), Gibbs free energy (Yuan and Pope 2012;
Beygi et al. 2014; Neshat and Pope 2018), and critical temperature (Martin et al. 2019). The correlations typically define limiting values
of the composition parameter and then interpolate curves between the limiting values. However, all have drawbacks. Density is often not
a good indicator of fluid state. For typical conditions of West Texas CO2 injection projects, the oil specific gravity is usually similar to that
of CO2, often higher near injection wells and lower near producing wells, so it is not a valid indicator. Gibbs free energy is also not a
reliable indicator (Churchwell et al. 2020), because fluids having the same Gibbs free energy can have different flow behavior. The defi-
nition of phase state based on critical temperature is strictly a convention. Like the other compositional parameters, critical temperature
does not adequately correlate the physics of multiphase flow.
Most of these empirical correlations consider not only compositional consistency, but also hysteresis, low IFT (high capillary number),
and the related correlation of capillary pressure. They appear to be continuous or consistent, but most still rely on phase labeling to some
extent. Some use rigid Corey-type functions, so they only approximate measured data. Most rely on numerous parameters which are
difficult to determine from measurements.
The compositional dependence of relative permeability is an area where additional work is needed. After comparing all the models
with then available data, Baker (1988) called for a more fundamental approach to relative permeability correlations. Those sentiments are
30 2022 SPE Journal

still valid. Rather than selecting a compositional parameter which has little relationship to flow behavior, the models need to incorporate
quantifiable physical parameters which govern the physics of multiphase flow, e.g. wettability, surface roughness, pore structure and
topology. Some progress has been made in this direction using high resolution x-ray imaging technology and pore network models (Blunt
et al. 2013; Blunt 2017). Work to date has shown promise that this approach may provide correlations founded on the physics of the pro-
cess (Khorsandi et al. 2017, 2018; Purswani et al. 2020).
Numerical and Physical Dispersion. The physical dispersion terms in Eq. 1 are usually neglected in compositional models. However,
here we summarize work showing that dispersion can significantly influence unit displacement and volumetric sweep efficiency. When
dispersion is neglected, one is left with numerical dispersion. Numerical dispersion depends on variables such as the grid configuration
and difference approximation, which are nonphysical but can still affect predictions dramatically.
The dispersivities in Eq. 4 depend on several parameters. There is obviously a dependence on viscosity or mobility ratio of the fluids
(Brigham et al. 1961; Young 1990) and heterogeneity of the rock (Lake and Hirasaki 1981; Kelkar and Gupta 1991). They also vary with
phase saturation; some write the equations with dispersive flux proportional to phase saturation, but the dependence is more complex
(Delshad et al. 1985). These relationships are not well understood. Owing to their dependence on heterogeneity, the magnitude of disper-
sion at field scale is problematic and should be viewed as an upscaling parameter. Because physical dispersion effects are usually neglected,
let us look at the consequences of neglecting physical dispersion.
When USW or upwinding is used, numerical dispersion is introduced into the solution (Lantz 1971; Peaceman 2000). Most take the
view that numerical dispersion is bad and must be minimized or eliminated. However, complete elimination of all dispersion is not desir-
able. It has long been known that dispersion or mechanical mixing can tear up small slugs of miscible fluids (Perkins and Johnston 1963).
This effect cannot be modeled unless it is correctly and accurately implemented in the reservoir simulator. There are other issues which
have been known for years but are ignored by most. For MCM displacements, dispersion reduces the displacement efficiency (Gardner
et al. 1981; Young et al. 1990). In most cases, displacements with injected gas involve adverse mobility ratios. Few realize that for these
displacements in multiple dimensions, the influence of dispersion on sweep and recovery is just as strong as the mobility ratio (Young
1981, 1984; Ewing et al. 1989).
In one dimension, there is an approximate relationship between physical dispersion and the numerical dispersion produced by USW.
Most implementations weight the mobilities, the λ terms in Eq. 25, upstream. Here we weight the fractional flows upstream, which pro-
duces a similar result. For illustration, consider Eq. 24 in one dimension, neglect gravity effects and capillary pressure, assume constant
positive velocity and partial volumes, and represent dispersion by Eq. 4:
PNp h i
@S
@t
i
+ u @f@x˛i ˛O l u @x
@
f˛ @S@x˛i = 0. (72)
‍ ˛=1 ‍
The equation could be further simplified by carrying out the summation on the convection term, but no simplification can be achieved for
the general case. When the convection term is approximated with upstream weighting or backward difference, the finite difference repre-
sentation for a uniform grid is:
PNp n h i h i h i o

Sij = t
x2
ux f˛i j1 f˛i j + ˛O ` u S˛i j+1 S˛i j f˛i j+ 1 + ˛O ` u S˛i j1 S˛i j f˛i j 1 + O x . (73)
‍ ˛=1 2 2
‍
The right-hand-side terms are evaluated at old or new time levels for explicit and implicit methods, respectively. The temporal approxi-
mation affects numerical dispersion too, but that effect is not considered in this discussion (Lantz 1971; Peaceman 2000). The truncation
error associated with the backward first derivative approximation takes the form of a dispersion term. Eq. 73 is reproduced by applying
second-order central differences to the following:
Np h
P i
Sij @2 f˛i
t = u @f@x˛i + u x
2 @x 2 + O
˛ ` u @
@x
f˛
@S˛i
@x
+ O x2
j
˛=1
Np h h ii . (74)
P
= u @f@x˛i + u x @ @f˛
2 @x S˛i @x +
@
@x
uf˛ x
2 + ˛O ` @S˛i
@x
+ O x2
‍ ˛=1 j ‍
In the second equation above, the fractional flow is expanded according to Eq. 20. This equation shows that upstream weighting intro-
duces a numerical dispersion term, proportional to uΔx/2, composed of two parts—an immiscible part (the second term) and a miscible
part (added to the dispersion term). It is common to ignore the immiscible part and claim that the approximation with no dispersion term
represents ‍˛O ` = x/2‍. This is approximately correct in 1D, but completely wrong for standard difference methods in multiple dimensions
(Fanchi 1983; Russell and Wheeler 1983).
MultiContact-Miscible Displacements. Simulations of MCM displacements show the displacement efficiency is strongly influenced
by the dispersion level (Gardner et al. 1981). With the correct dispersion level, laboratory experiments can be simulated with reasonable
accuracy (Young et al. 1990; Negahban and Kremesec 1992), including important length-scaling effects.
A 1D example of N2 displacing a rich gas condensate is taken from Young (1984), but the same behavoir occurs for CO2 and other
miscible gas displacements (Young 1990). A 17-component EOS description and upstream weighting was used with no physical disper-
sion, ‍˛O ` = 0‍. Fig. 8 shows a pseudo ternary diagram constructed from the compositions present after 0.6 hydrocarbon pore volumes
injected with 50 grid intervals. The ternary diagram was constructed using mass fractions and by lumping the original 16 hydrocarbon
components into two pseudo components, C2− and C3+. Fig. 9 shows the recovery curves for various grid sizes. Because physical disper-
sion is neglected, the results are strongly dependent on the grid size.
When dispersion is neglected, the equations are hyperbolic and can be solved by the method of characteristics (Helfferich 1981; Orr
2007). The method of characteristics solution shows that with no dispersion, the composition path intersects the phase diagram at the
tangent point and then follows the dewpoint boundary to 100% N2. For MCM displacements such as this, the solution is that of a piston
displacement with 100% recovery at breakthrough. The recovery curves in Fig. 9 show the computed results trend in that direction.
Fig. 10 plots the amount unrecovered vs. grid size because this is effectively the error for the case of no dispersion. Sure enough, the
solutions are converging to the analytical solution, that is, a piston displacement with 100% recovery or a zero final oil saturation (Sof)
after miscible flooding.
2022 SPE Journal 31

Fig. 8—Compositions for N2 displacing gas condensate: o—in place fluid, _ _ _ _ tie-lines, ···· path, from Young (1984).
Fig. 9—N2 displacing gas condensate, recovery for various grid sizes, from Young (1984).
32 2022 SPE Journal

Fig. 10—N2 displacing gas condensate, convergence with grid refinement, from Young (1984).
It should go without saying, the reason for using a compositional model is to predict performance. Laboratory and field tests show that a residual
saturation is left by the displacement. Some use four-component models with a specified Sof e.g. (Chopra et al. 1990), but these models cannot
predict loss of miscibility due to pressure or chase gas contamination nor can they predict the effect of continuing vaporization behind the displace-
ment front. The only way to predict these effects is to explicitly and accurately model dispersion in the simulator, yet this is virtually never done.
One is left with some grid dependent numerical dispersion with a nebulous relationship to true physical dispersion. Continued grid refinement gives
a solution with ‍Sof = 0‍which is unrealistic, i.e. wrong.
Adverse Mobility Ratio Displacements. For most gas injection processes, the mobility of the injected gas is greater than the mobility
of the in-place fluid. These displacements lead to the phenomenon of viscous fingering. Viscous fingering has been studied in various
laboratory experiments (Hall and Geffen 1957; Lacey et al. 1958; Blackwell et al. 1959; Habermann 1960) and by fine scale simulations
(Christie 1989; Ewing et al. 1989; Scovazzi et al. 2017). The last article referenced is a review, providing numerous references to articles
describing fine scale calculations. These problems are usually studied with 2D cross sections or five-spot patterns using fluids which are
FCM. When dispersion is modeled, viscous instability should occur in a simulator only for heterogenous systems. Many of the early
simulations and some of the later ones are suspect because the dispersion model, Eq. 4, was not properly implemented. Some have used
numerical dispersion in an effort to mimic physical dispersion, while others have used specified constant dispersion coefficients in x
and y in place of Eq. 4. These approaches are not adequate. Few have used convergence studies to demonstrate valid results in multiple
dimensions, even for homogenous systems. If a method is not accurate for a homogenous system, it is ridiculous to claim the results are
valid for a heterogenous one.
The basic features of these FCM displacements are that displacement efficiency remains high but sweep efficiency is low.
Explicit simulation of viscous fingering is not feasible on a field scale, so some sort of averaging technique is required. Two
basic macroscopic methods have been proposed involving: (1) alteration of the convection terms while neglecting dispersion
(Koval 1963; Todd and Longstaff 1972; Fayers 1988) and (2) use of the standard equations. In the first approach a mixing param-
eter is used to alter the convection term. In the later approach, dispersion is either not explicitly included (Yanosik and McCracken
1979) or effective dispersion coefficients are used (Ewing et al. 1989; Young 1990). Of course, there is always numerical dis-
persion if the effect is not explicitly modeled.
The mixing parameter approach has been extended to compositional flow (Blunt et al. 1994; Barker and Fayers 1994). 1D simulations
of immiscible gas drives match heterogenous 2D cross-section simulations. Upscaling of the approach has also been considered (Iranshahr
et al. 2014; Li and Durlofsky 2016). Areal sweep effects in pattern floods and the more common MCM displacements were not investi-
gated in these studies.
When areal sweep in FCM pattern floods is considered, the mixing parameter approach produces solutions with good sweep efficiency
and poor displacement efficiency, so it does not reflect the true physics of the process. Furthermore, in many cases, it does a poor job of
matching available recovery data (Ewing et al. 1989). Unlike the dispersion model, adequate simulation of laboratory MCM displace-
ments has not been demonstrated. For these reasons, the mixing parameter approach is not considered further here.
When a dispersion model is used, the dispersion coefficients are viewed as upscaling parameters. Most studies on the effects of hetero-
geneity and upscaling of dispersion have considered only unit mobility ratio displacement (e.g., Lake and Hirasaki 1981; Arya et al. 1988;
John et al. 2008) rather than the more common adverse mobility ratio displacements. Only a few have addressed adverse mobility ratio
problems (Kelkar and Gupta 1991; Garmeh and Johns 2010) and additional work is needed.
In addition to matching laboratory MCM displacements, a dispersion model also matches areal sweep and recovery in available exper-
imental data and fine scale simulations (Ewing et al. 1989). It is widely used in field scale simulation of MCM CO2 floods (Chopra et al.
1990; Hill et al. 1994; Hsu et al. 1997), although none of these field studies explicitly modeled dispersion. These applications were for
alternate gas/water injection which mitigates adverse mobility ratio issues.
Simulating these displacements with standard USW finite difference methods leads to the well-known grid orientation effect. For 2D
areal simulations, the effect is tested by using the two different elements of symmetry shown in Fig. 11. For a five-point difference
method with USW, Fig. 12 shows the front locations for a FCM displacement with a mobility ratio M = 10. Numerous articles have
discussed this problem with many cures claimed. The grid orientation sensitivity is circumvented by use of a nine-point finite difference
method. Many nine-point methods have been proposed (Yanosik and McCracken 1979; Potempa 1982; Coats 1983). Nine-point
2022 SPE Journal 33

Fig. 11—Two grids for five-spot pattern, from Yanosik and McCracken (1979).
methods in 2D and 11-point approximations in 3D are available in commercial simulators and are the most widely used methods for
simulating these processes. In these models, dispersion is not explicitly modeled, but a grid size depend numerical dispersion is
present.
Fig. 12—Front locations with USW five-point for M = 10, from Yanosik and McCracken (1979).
34 2022 SPE Journal

The explanation originally given for the orientation sensitivity is that the problem occurs because a five-point difference star does not
allow for flow diagonal to the grid lines. Many still believe this explanation is correct. However, long ago, it was shown to be wrong
because a five-point with physical dispersion implemented according to Eq. 4 does not suffer from the problem when the grid is fine
enough to resolve the front (Watts and Silliman 1980; Young 1981). Young gave the now accepted explanation that a five-point with USW
introduces numerical dispersion which is not rotationally invariant, leading to the problem.
The problem has been studied extensively in the petroleum industry (e.g. Russell and Wheeler 1983; Young 1984; Shubin and Bell
1984). Russell and Wheeler showed that the USW nine-point finite difference method introduces numerical dispersion which is almost
rotationally invariant with ‍˛O ` /˛O t 3.‍ The sensitivity of the solution to the dispersion level suggests that the USW nine-point method will
be sensitive to grid size. This was demonstrated by comparing the recovery of a nine-point method with various grid sizes to a dispersion
model with various dispersion levels, see Fig. 13. The dispersion level is characterized by the Peclet number, ‍Pe = L/˛O `‍ and the ratio
‍˛O ` /˛O t ‍. The results in Fig. 13 are for a dispersion model with ‍˛O ` /˛O t = 1‍, whereas nine-point results are closer to those for ‍˛O ` /˛O t = 3‍. These
results indicate the nine-point method may not exhibit a sensitivity to the orientation of the grid, but the recovery is a strong function of
grid size. The sweep efficiency is continually reduced as the grid is refined, so the results do not converge. Many alternatives to USW have
been proposed (upcoming paper titled “50 Years of Reservoir Simulation” by Lie et al. in this special issue), but virtually all introduce
some sort of numerical dispersion, so they will suffer from the same convergence problem. Without dispersion, these displacements are
inherently unstable, so the problem is not well posed. The physical dispersion term is needed to stabilize the displacement. When one
neglects this essential physical phenomenon, the predictions are influenced strongly by numerical truncation error. This result should be
no surprise. The problem is not a grid orientation problem but rather a convergence problem which is not solved by an USW nine-point
or other methods that neglect physical dispersion. The direct simulation of physical dispersion effects is the most reasonable solution to
this problem.
Fig. 13—Recovery vs. grid size o – 9-pt diagonal, ● – 9-pt parallel; and vs. 2/Pe x - dispersion model from Young (1984).
2022 SPE Journal 35

Many articles have been written describing methods for simulating dispersion for FCM displacement problems, see (Scovazzi et al.
2017) but few have addressed the problem for general compositional models which include phase equilibria, phase boundaries, and
aspects of both miscible and immiscible flow. Some throw a difference formula for dispersion into the equations but few have demon-
strated that meaningful, convergent results are achieved, an essential requirement given the number of methods shown not to converge.
(Young 1984) considered a hierarchy of problems: miscible, immiscible, a combination of the two, with and without phase changes, and
the full compositional model equations for an MCM displacement. Physical dispersion was implemented along with a rotationally invari-
ant artificial dispersion to stabilize the immiscible aspects of the flow. For each case, the convergence was tested. A few other studies have
considered direct implementation of physical dispersion in general compositional models. For these more complex problems, phase
boundaries are an issue, because the dispersion term in Eq. 1 involves gradients of a composition which is undefined on one side of the
transition. Various ad hoc procedures have been used to deal with this difficulty (Young 1990; Shiralkar and Stephenson 1991; Shrivastava
et al. 2005a, 2005b), but explicit treatment of the phase boundary has not been considered.
In summary, for MCM displacements, both unit displacement efficiency and sweep efficiency have a strong dependence on dispersion.
If physical dispersion is modeled, the results depend on its magnitude. If physical dispersion is not modeled, the results have an undesir-
able dependence on grid size, difference scheme, and other numerical truncation errors. For these reasons, the correct, reliable, and accu-
rate simulation of physical dispersion in compositional models is another area where additional work is needed. In addition to the
numerical treatment, an improved understanding is needed of dispersion or mechanical mixing at field scale (i.e., upscaling).
Conclusions
Some major conclusions documented in this work are:
1. Multiphase partial volumes are key physical parameters which are smoothly varying functions except for jump discontinuities at phase
boundaries and rapid changes near critical points.
2. Early models showed the need for property consistency, fully linearized and coupled fluid properties and attention to temporal stability.
3. The so-called volume balance (Acs et al. 1985) and Newton-Raphson methods (Young and Stephenson 1983) are essentially the same.
4. So-called natural variables (SXP) are inferior to mass/molar/pressure (MP) variables for IMPES and AIM and appear to offer no ad-
vantage even for FIM. Solutions with MP variables are naturally conservative and well conditioned.
5. Phase stability analysis reliably detects the appearance of a phase. Negative flash calculations offer a more efficient alternative for most
conditions. However, it is applicable only for compositions lying on a tie-line extension, so phase stability calculations are needed as
a backup.
6. AIM is attractive for compositional simulation, because it overcomes the stability restrictions of IMPES while mitigating the compu-
tational burden of FIM.
7. The simulation of huge problems is feasible with distributed parallel computing. Efficient implementation and exploitation of fine-
grained parallelism is still needed and often can obviate the need for distributed memory computations.
8. Although many improvements have been made to the speed and reliability of compositional models, several areas need additional
work. These include
1. calculations for problems with three hydrocarbon phases;
2. consistent and physically correct treatment of compositionally dependent relative permeability;
3. modeling and upscaling of physical dispersion (mechanical mixing) for realistic compositional problems with phase boundaries,
because current models have an ill-defined dependence on grid size.
Nomenclature
a = coefficients in saturation linearization, Eq. 44
‍aL ‍ = parameters to calculate for EOS “A” parameter
‍AL ‍ = EOS “A” parameter
A = matrix representation of secondary/phase equilibria equations, Eqs. 36 and 46
b = coefficients in mole fraction linearization, Eq. 44
B = linearization coefficients of Rachford-Rice Eqs. 34 and 46
ci = the combination of terms‍(xvi −xì )/zi ‍
‍C ‍ = Courant-Friedrichs-Lewy number
C = coefficients of linearized volume balance or saturation constraint Eq. 37
Cr = rock compressibility
‍ ‍ = total compressibility
C
di = the combination of terms‍(xvi −xì )/zi ‍
‍ ‍ = dispersion coefficient
D
D = coefficients in primary/flow equations
e = unit vector
E = coefficients in component balance equations, Eq. 67
EN = overall efficiency with N processors
f = fractional flow
fs = fraction of serial calculations
F = derivatives of fugacity coefficients
g = gravitational constant
h = constant term in Eq. 60
H = depth
I = identity matrix
k = permeability
kr = relative permeability
K = equilibrium K-value
L = length of side of 5-spot quadrant
M = mobility ratio, (displacing fluid)/(displaced fluid)
n = mole numbers for component, overall or components within a phase
36 2022 SPE Journal

Nc = total number of components
Ng = number of grid points
Nh = number of hydrocarbon components
Np = total number of phases
Nproc = number of computational processors
p= characteristic pressure
pα = difference between phase pressure and characteristic pressure
Pc = critical pressure
Pe = Peclet number, L/âl
q= constant terms, residual and source/sink terms in component balances
Q= volumetric flow out of a grid cell
R= vector defined in conjunction with Eq. 51
S= saturation or volume fraction
t= time
tN = total calculation time with N processors
T= transmissibility coefficient
Tc = critical temperature
u= velocity
υ= specific volume, molar volume or partial molar volume
U= vector defined by Eq. 50
v= vapor mole fraction, water free basis
V= pore volume
w= intermediate variables (xvixli/zi)Δωi, defined in conjunction with Eq. 46
x= mole or mass fraction
‍‍ = secondary variables
z= overall hydrocarbon mole fraction, water free basis
Zc = critical compressibility factor
Λ= maximum eigenvalue of Eq. 64
ϕ= porosity
‍O ‍ =
normalized porosity, ϕ/ϕο
ϕο = reference porosity
‍(xvi −xì )/zi ‍ = dispersivity
δu = deviation from average velocity
η= intermediate variables Δni/n, defined in conjunction with Eq. 46
λ= mobility of phase, component or overall
μ= phase viscosity
‍N ‍ = flow weight average density, Σαfαρα
ρ= mass density
φ= fugacity coefficient
ψ= primary variables
ω= logarithm of equilibrium K-values
Subscripts and superscripts

a = aqueous phase quantity
h = overall value for hydrocarbon phases
α, β = quantity for a phase
i, j, k = quantity for a component or generic index
‍`‍ = hydrocarbon liquid phase quantity or longitudinal term
t = transverse term
ν = vapor phase quantity
T = transpose
w = water component property
Acknowledgments
I thank Amoco Production company for providing a rich research environment during my formative years. Specifically, I thank Del
Fussell for the initial contact, Howard Hall for supporting my research efforts, and John Yanosik for keeping me grounded. I also thank
Rami Younis for his encouragement and council while writing this article.
References
Acs, G., Doleschall, S., and Farkas, E. 1985. General Purpose Compositional Model. SPE J. 25 (4): 543–553. SPE-10515-PA. https://doi.org/10.2118/10515-
PA.
Ader, J. C. and Stein, M. H. 1984. Slaughter Estate Unit Tertiary Miscible Gas Pilot Reservoir Description. J Pet Technol 36 (5): 837–845. SPE-10727-PA.
https://doi.org/10.2118/10727-PA.
Agarwal, A., Moncorge, A., and Tchelepi, H. A. 2006. Stability Criteria for the Thermal Adaptive Implicit Method. Paper presented at the SPE Annual
Technical Conference and Exhibition, San Antonio, Texas, USA, 24–27 September. SPE-103060-MS. https://doi.org/10.2118/103060-MS.
Al-Shaalan, T. M., Klie, H. M., Dogru, A. H. et al. 2009. Studies of Robust Two Stage Preconditioners for the Solution of Fully Implicit Multiphase Flow
Problems. Paper presented at the SPE Reservoir Simulation Symposium. The Woodlands, The Woodlands, Texas, USA, 2–4 February. SPE-118722-
MS. https://doi.org/10.2118/118722-MS.
2022 SPE Journal 37

Alpak, F. O. and Vink, J. C. 2018. A Variable-Switching Method for Mass-Variable-Based Reservoir Simulators. SPE J 23 (5): 1469–1495. SPE-182606-
PA. https://doi.org/10.2118/182606-PA.
Amestoy, P. R., Duff, I. S., L’Excellent, J.-Y. et al. 2001. MUMPS: A General Purpose Distributed Memory Sparse Solver. eds. Tor Sørevik, Fredrik Manne,
and Assefaw Hadish Gebremedhin, 121–130. Berlin, Heidelberg: Springer.
Appleyard, J. R., Chesire, I. M., and Pollard, R. K. 1981. Special Techniques for Fully-Implicit Simulators. Oxfordshire, UK: Harwell.
Arya, A., Hewett, T. A., Larson, R. G. et al. 1988. Dispersion and Reservoir Heterogeneity. SPE Res Eng 3 (1): 139–148. SPE-14364-PA. https://doi.org/
10.2118/14364-PA.
Attra, H. D. 1961. Nonequilibrium Gas Displacement Calculations. SPE J. 1 (3): 130–136. SPE-1522-G. https://doi.org/10.2118/1522-G.
Baker, L. E. 1988. Three-Phase Relative Permeability Correlations. Paper presented at the SPE Enhanced Oil Recovery Symposium, Tulsa, Oklahoma,
16–21 April. SPE-17369-MS. https://doi.org/10.2118/17369-MS.
Balay, S., Abhyankar, S., Adams, M. et al. 2019. PETSc Users Manual, Revision 3.8. Argonne, Illinois, USA: Argonne National Laboratory.
Barker, J. W. and Fayers, F. J. 1994. Transport Coefficients for Compositional Simulation With Coarse Grids in Heterogeneous Media. SPE Advanced
Technology Series 2 (2): 103–112. SPE-22591-PA. https://doi.org/10.2118/22591-PA.
Bear, J. 2013. Dynamics of Fluids in Porous Media. Mineola, New York, USA: Dover Publications.
Beckner, B. L., Haugen, K. B., Maliassov, S. et al. 2015. General Parallel Reservoir Simulation. Paper presented at the Abu Dhabi International Petroleum
Exhibition and Conference, Abu Dhabi, UAE, 9–12 November. SPE-177532-MS. https://doi.org/10.2118/177532-MS.
Behie, A. and Vinsome, P. K. W. 1982. Block Iterative Methods for Fully Implicit Reservoir Simulation. SPE J. 22 (5): 658–668. SPE-9303-PA. https://
doi.org/10.2118/9303-PA.
Belkadi, A., Yan, W., Moggia, E. et al. 2013. Speeding Up Compositional Reservoir Simulation through an Efficient Implementation of Phase Equilibrium
Calculation. Paper presented at the SPE Reservoir Simulation Symposium, The Woodlands, Texas, USA, 18–20 February. SPE-163598-MS. https://
doi.org/10.2118/163598-MS.
Bertiger, W. I. and Kelsey, F. J. 1985. Inexact Adaptive Newton Methods. Paper presented at the SPE Reservoir Simulation Symposium, Dallas, Texas,
USA, 10–13 February. SPE-13501-MS. https://doi.org/10.2118/13501-MS.
Beygi, M. R., Delshad, M., Pudugramam, V. S. et al. 2014. Novel Three-Phase Compositional Relative Permeability and Three-Phase Hysteresis Models.
SPE J. 20 (1): 21–34. SPE-165324-PA. https://doi.org/10.2118/165324-PA.
Blackwell, R. J., Rayne, J. R., and Terry, W. M. 1959. Factors Influencing the Efficiency of Miscible Displacement. Trans 217 (1): 1–8. SPE-1131-G.
https://doi.org/10.2118/1131-G.
Blair, P. M. and Weinaug, C. F. 1969. Solution of Two-Phase Flow Problems Using Implicit Difference Equations. SPE J. 9 (4): 417–424. SPE-2185-PA.
https://doi.org/10.2118/2185-PA.
Blunt, M. J. 2000. An Empirical Model for Three-Phase Relative Permeability. SPE J. 5 (4): 435–445. SPE-67950-PA. https://doi.org/10.2118/67950-PA.
Blunt, M. J. 2017. Multiphase Flow in Permeable Media: A Pore-Scale Perspective. Cambridge, UK: Cambridge University Press.
Blunt, M. J., Barker, J. W., Rubin, B. et al. 1994. Predictive Theory for Viscous Fingering in Compositional Displacement. SPE Res Eng 9 (1): 73–80.
SPE-24129-PA. https://doi.org/10.2118/24129-PA.
Blunt, M. J., Bijeljic, B., Dong, H. et al. 2013. Pore-Scale Imaging and Modelling. Adv Water Resour 51: 197–216. https://doi.org/10.1016/j.
advwatres.2012.03.003.
Bogachev, K., Milyutin, S., Telishev, A. et al. 2018. High-Performance Reservoir Simulations on Modern CPU-GPU Computational Platforms. Paper
presented at the SPE Russian Petroleum Technology Conference, Moscow, Russia, 16–18 October. SPE-187797-MS. https://doi.org/10.2118/187797-
MS.
Branco, C. M. and Rodriguez, F. 1996. A Semi-Implicit Formulation for Compositional Reservoir Simulation. SPE Advanced Technology Series 4 (1):
171–177. SPE-27053-PA. https://doi.org/10.2118/27053-PA.
Breitenbach, E. A. and Thurnau, D. H. 1968. The Fluid Flow Simulation Equations. Paper presented at the Numerical Simulation Symposium, Dallas,
Texas, USA, 22–23 April. SPE-2020-MS. https://doi.org/10.2118/2020-MS.
Brigham, W. E., Reed, P. W., and Dew, J. N. 1961. Experiments on Mixing During Miscible Displacement in Porous Media. SPE J. 1 (1): 1–8. SPE-1430-G.
https://doi.org/10.2118/1430-G.
Cao, H. 2002. Development of Techniques for General Purpose Simulators. Stanford, California, USA: Stanford University.
Cao, H. and Aziz, K. 2002. Performance of IMPSAT and IMPSAT-AIM Models in Compositional Simulation. Paper presented at the SPE Annual Technical
Conference and Exhibition, San Antonio, Texas, USA, 29 September–2 October. SPE-77720-MS. https://doi.org/10.2118/77720-MS.
Cao, H., Crumpton, P. I., and Schrader, M. L. 2009. Efficient General Formulation Approach For Modeling Complex Physics. Paper presented at the SPE
Reservoir Simulation Symposium, The Woodlands, Texas, USA, 2–4 February. SPE-119165-MS. https://doi.org/10.2118/119165-MS.
Cao, H., Zaydullin, R., and Obi, E. 2017. Nonlinear Convergence for Near-Miscible Problem: A Mystery Unveiled for Natural Variable Simulator. Paper
presented at the SPE Reservoir Simulation Conference, Montgomery, Texas, USA, 20–22 February. SPE-182633-MS. https://doi.org/10.2118/182633-
MS.
Chavent, G. and Jaffré, J. 2014. Mathematical Models and Finite Elements for Reservoir Simulation: Single Phase, Multiphase and Multicomponent Flows
through Porous Media. Burlington: Elsevier Science.
Chen, Z. and Ewing, R. E. 1997. Comparison of Various Formulations of Three-Phase Flow in Porous Media. J Comput Phys 132 (2): 362–373. https://
doi.org/10.1006/jcph.1996.5641.
Chien, M. C. H., Lee, S. T., and Chen, W. H. 1985. A New Fully Implicit Compositional Simulator. Paper presented at the SPE Reservoir Simulation
Symposium, Dallas, Texas, USA, 10–13 February. SPE-13385-MS. https://doi.org/10.2118/13385-MS.
Chien, M. C. H., Tchelepi, H. A., Yardumian, H. E. et al. 1997. A Scalable Parallel Multi-Purpose Reservoir Simulator. Paper presented at the SPE
Reservoir Simulation Symposium, Dallas, Texas, USA, 8–11 June. SPE-37976-MS. https://doi.org/10.2118/37976-MS.
Chien, M. C. H., Wasserman, M. L., Yardumian, H. E. et al. 1987. The Use of Vectorization and Parallel Processing for Reservoir Simulation. Paper
presented at the SPE Symposium on Reservoir Simulation, San Antonio, Texas, USA, 1–4 February. SPE-16025-MS. https://doi.org/10.2118/16025-
MS.
Chopra, A. K., Stein, M. H., and Dismuke, C. T. 1990. Prediction of Performance of Miscible-Gas Pilots. J Pet Technol 42 (12): 1564–1572. SPE-18078-
Christie, M. A. 1989. High-Resolution Simulation of Unstable Flows in Porous Media. SPE Res Eng 4 (3): 297–303. SPE-16005-PA. https://doi.org/
10.2118/16005-PA.
Churchwell, L., Radhakrishnan, A., and DiCarlo, D. 2020. Measurements of Three-Phase Relative Permeability as a Function of Fluid Composition.
Paper presented at the SPE Improved Oil Recovery Conference, Virtual, 31 August–4 September. SPE-200438-MS. https://doi.org/10.2118/200438-
MS.
38 2022 SPE Journal

Clees, T. and Ganzer, L. 2007. An Efficient Algebraic Multigrid Solver Strategy for Adaptive Implicit Methods in Oil Reservoir Simulation. Paper presented
at the SPE Reservoir Simulation Symposium, Houston, Texas, USA, 26–28 February. SPE-105789-MS. https://doi.org/10.2118/105789-MS.
Coats, K. H. 1980. An Equation of State Compositional Model. SPE J. 20 (5): 363–376. SPE-8284-PA. https://doi.org/10.2118/8284-PA.
Coats, K. H. 1983. A Consistent Method for Calculating Transmissibilities in Nine-Point Difference Equations. Paper presented at the SPE Reservoir
Simulation Symposium, San Francisco, California, USA, 15–18 November. SPE-12248-MS. https://doi.org/10.2118/12248-MS.
Coats, K. H. 2000. A Note on IMPES and Some IMPES-Based Simulation Models. SPE J. 5 (3): 245–251. SPE-65092-PA. https://doi.org/10.2118/65092-
PA.
Coats, K. H. 2003a. IMPES Stability: Selection of Stable Timesteps. SPE J 8 (2): 181–187. SPE-84924-PA. https://doi.org/10.2118/84924-PA.
Coats, K. H. 2003b. IMPES Stability: The CFL Limit. SPE J. 8 (3): 291–297. SPE-84924-PA. https://doi.org/10.2118/85956-PA.
Coats, K. H. 2009. Reservoir Simulation. In Petroleum Engineering Handbook, eds. H. B. Bradley. Richardson, Texas, USA: Society of Petroleum
Engineers.
Coats, K. H., Thomas, L. K., and Pierson, R. G. 1998. Compositional and Black Oil Reservoir Simulation. SPE Res Eval & Eng 1 (4): 372–379. SPE-
50990-PA. https://doi.org/10.2118/50990-PA.
Collins, D. A., Grabenstetter, J. E., and Sammon, P. H. 2003. A Shared-Memory Parallel Black-Oil Simulator with a Parallel ILU Linear Solver. Paper
presented at the SPE Reservoir Simulation Symposium, Houston, Texas, 3–5 February. SPE-79713-MS. https://doi.org/10.2118/79713-MS.
Collins, D. A., Nghiem, L. X., Li, Y.-K. et al. 1992. An Efficient Approach to Adaptive- Implicit Compositional Simulation With an Equation of State. SPE
Res Eng 7 (2): 259–264. SPE-15133-PA. https://doi.org/10.2118/15133-PA.
Craig, F. F. 1971. The Reservoir Engineering Aspects of Waterflooding. Richardson, Texas: Monograph 3, Society Petroleum Engrs.
DeBaun, D., Byer, T., Childs, P. et al. 2005. An Extensible Architecture for Next Generation Scalable Parallel Reservoir Simulation. Paper presented at
the SPE Reservoir Simulation Symposium, The Woodlands, Texas, USA, 31 January–2 February. SPE-93274-MS. https://doi.org/10.2118/93274-
MS.
Delshad, M., MacAllister, D. J., Pope, G. A. et al. 1985. Multiphase Dispersion and Relative Permeability Experiments. SPE J. 25 (4): 524–534. SPE-
Dogru, A. H., Fung, L. S. K., Middya, U. et al. 2009. A Next-Generation Parallel Reservoir Simulator for Giant Reservoirs. Paper presented at the SPE
Dogru, A. H., Sunaidi, H. A., Fung, L. S. et al. 2002. A Parallel Reservoir Simulator for Large-Scale Reservoir Simulation. SPE Res Eval & Eng 5 (1):
Ewing, R. E., Russell, T. F., and Young, L. C. 1989. An Anistropic Coarse-Grid Dispersion Model of Heterogeneity and Viscous Fingering in Five-Spot
Miscible Displacement That Matches Experiments and Fine-Grid Simulations. Paper presented at the SPE Symposium on Reservoir Simulation,
Houston, Texas, USA, 6–8 February. https://doi.org/10.2118/18441-MS.
Faddeeva, V. N. 1959. Computational Methods of Linear Algebra. New York: Dover.
Fagin, R. G. and Stewart, C. H. 1966. A New Approach to the Two-Dimensional Multiphase Reservoir Simulator. SPE J. 6 (02): 175–182. SPE-1188-PA.
Falgout, R. D., Jones, J. E., and Yang, U. M. 2006. The Design and Implementation of Hypre, a Library of Parallel High Performance Preconditioners
BT. In Numerical Solution of Partial Differential Equations on Parallel Computers, eds. A. M. Bruaset and A. Tveito, 267–294. Berlin, Heidelberg:
Springer. https://doi.org/10.1007/3-540-31619-1.
Fanchi, J. R. 1983. Multidimensional Numerical Dispersion. SPE J. 23 (1): 143–151. SPE-9018-PA. https://doi.org/10.2118/9018-PA.
Farkas, E. 1993. SPE Symposium on Reservoir Simulation. Paper presented at the Adaptive Implicit Volume Balance Techniques, New Orleans, Louisiana.
SPE-25273-MS. https://doi.org/10.2118/25273-MS.
Farkas, E. 1996. Linearization Techniques of Reservoir Simulation Equations - IMPES/IMPEM Cases. Paper presented at the 5th European Conference on
the Mathematics of Dil Recovery, Leoben, Austria, 3–6 September. https://doi.org/10.3997/2214-4609.201406879.
Farkas, E. 1998. Linearization Techniques of Reservoir-Simulation Equations: Fully Implicit Cases. SPE J. 3 (4): 316–323. SPE-52051-PA. https://doi.
org/10.2118/52051-PA.
Fayers, F. J. 1988. An Approximate Model With Physically Interpretable Parameters for Representing Miscible Viscous Fingering. SPE Res Eng 3 (2):
Fayers, F. J., Foakes, A. P., Lin, C. Y. et al. 2000. An Improved Three Phase Flow Model Incorporating Compositional Variance. Paper presented at the SPE/
DOE Improved Oil Recovery Symposium, Tulsa, Oklahoma, 3–5 April. SPE-59313-MS. https://doi.org/10.2118/59313-MS.
Fernandes, B. R. B., Sepehrnoori, K., Marcondes, F. et al. 2017. A Natural Variable Fully Implicit Compositional Reservoir Simulation. In XXXVIII
Iberian-Latin American Congress on Computational Methods in Engineering, eds. M. Noronh, P. O. Faria, and R. H. Lopez. SC, Brazil: Florianópolis.
https://doi.org/10.20906/CPS/CILAMCE2017-0447.
Fjerstad, P. A., Dasie, W. J., Sikandar, A. S. et al. 2005. Next Generation Parallel Computing for Large-Scale Reservoir Simulation. Paper presented at
the SPE International Improved Oil Recovery Conference in Asia Pacific, Kuala Lumpur, Malaysia, 5–6 December. SPE-97358-MS. https://doi.org/
10.2118/97358-MS.
Forsythe, G. and Moler, C. B. 1967. Computer Solution of Linear Algebraic Systems. Englewood Cliffs, NJ: Prentice-Hall.
Franc, J., Horgue, P., Guibert, R. et al. 2016. Benchmark of Different CFL Conditions for IMPES. Comptes Rendus Mécanique 344 (10): 715–724. https://
doi.org/10.1016/j.crme.2016.08.003.
Fung, L. S. K., Collins, D. A., and Nghiem, L. X. 1989. An Adaptive-Implicit Switching Criterion Based on Numerical Stability Analysis (Includes
Associated Paper 18727). SPE Res Eng 4 (01): 45–51. SPE-16003-PA. https://doi.org/10.2118/16003-PA.
Fung, Larry S.K., Sindi, M. O., and Dogru, A. H. 2014. Multiparadigm Parallel Acceleration for Reservoir Simulation. SPE J. 19 (04): 716–725. SPE-
Fong, W. S., Sheffield, J. M., Ehrlich, R. et al. 1992. Phase Behavior Modeling Techniques for Low-Temperature CO2 Applied to McElroy and North Ward
Estes Projects. Paper presented at the SPE/DOE Enhanced Oil Recovery Symposium, Tulsa, Oklahoma, 22–24 April. SPE-24184-MS. https://doi.org/
10.2118/24184-MS.
Forsyth, P. A. and Sammon, P. H. 1986. Practical Considerations for Adaptive Implicit Methods in Reservoir Simulation. J Comput Phys 62 (2): 265–281.
https://doi.org/10.1016/0021-9991(86)90127-0.
Fussell, D. D., Shelton, J. L., and Griffith, J. D. 1976. Effect of “Rich” Gas Composition on Multiple-Contact Miscible Displacement A Cell-to-Cell Flash
Model Study. SPE J. 16 (06): 311–316. SPE-5051-PA. https://doi.org/10.2118/5051-PA.
Fussell, D. D. and Yanosik, J. L. 1978. An Iterative Sequence for Phase-Equilibria Calculations Incorporating the Redlich-Kwong Equation of State. SPE
J 18 (3): 173–182. SPE-6050-PA. https://doi.org/10.2118/6050-PA.
2022 SPE Journal 39

Fussell, L. T. 1979. A Technique for Calculating Multiphase Equilibria (Includes Associated Papers 8734 and 8746). SPE J 19 (4): 203–210. SPE-6722-PA.
Fussell, L. T. and Fussell, D. D. 1979. An Iterative Technique for Compositional Reservoir Models. SPE J 19 (4): 211–220. SPE-6891-PA. https://doi.org/
10.2118/6891-PA.
Ganapathy, C. and Voskov, D. V. 2018. Multiscale Reconstruction in Physics for Compositional Simulation. J Comput Phys 375: 747–762. https://doi.org/
10.1016/j.jcp.2018.08.046.
Gandham, R., Esler, K., Mukundakrishnan, K. et al. 2016. GPU Acceleration of Equation of State Calculations in Compositional Reservoir Simulation.
Paper presented at the ECMOR XV - 15th European Conference on the Mathematics of Oil Recovery, Amsterdam, Netherlands. https://doi.org/
10.3997/2214-4609.201601747.
Gardner, J. W., Orr, F. M., and Patel, P. D. 1981. The Effect of Phase Behavior on CO2-Flood Displacement Efficiency. J Pet Technol33 (11): 2067–2081.
Garmeh, G. and Johns, R. T. 2010. Upscaling of Miscible Floods in Heterogeneous Reservoirs Considering Reservoir Mixing. SPE Res Eval & Eng13 (5):
Gharbia, I. B., Flauraud, E., and Michel, A. 2015. Study of Compositional Multi-Phase Flow Formulations with Cubic EOS. Paper presented at the SPE
Reservoir Simulation Symposium, Houston, Texas, USA, 23–25 February. SPE-173249-MS. https://doi.org/10.2118/173249-MS.
Godbole, S. P., Thele, K. J., and Reinbold, E. W. 1995. EOS Modeling and Experimental Observations of Three-Hydrocarbon-Phase Equilibria. SPE Res
Eng 10 (02): 101–108. SPE-24936-PA. https://doi.org/10.2118/24936-PA.
Gorucu, S. E. and Johns, R. T. 2014. Comparison of Reduced and Conventional Two-Phase Flash Calculations. SPE J 20 (2): 294–305. SPE-163577-PA.
https://doi.org/10.2118/163577-PA.
Grabenstetter, J., Li, Y.-K., Collins, D. A. et al. 1991. Stability-Based Switching Criterion for Adaptive-Implicit Compositional Reservoir Simulation.
Paper presented at the SPE Symposium on Reservoir Simulation, Anaheim, California. SPE-21225-MS. https://doi.org/10.2118/21225-MS.
Guan, W., Qiao, C., Zhang, H. et al. 2015. On Robust and Efficient Parallel Reservoir Simulation on Tianhe-2. Paper presented at the SPE Reservoir
Characterisation and Simulation Conference and Exhibition, Abu Dhabi, UAE, 14–16 September. SPE-175602-MS. https://doi.org/10.2118/175602-
MS.
Guler, B., Wang, P., Delshad, M. et al. 2001. Three- and Four-Phase Flow Compositional Simulations of CO2/NGL EOR. Paper presented at the SPE
Annual Technical Conference and Exhibition, New Orleans, Louisiana, 30 September–3 October. SPE-71485-MS. https://doi.org/10.2118/71485-MS.
Habermann, B. 1960. The Efficiency of Miscible Displacement as a Function of Mobility Ratio. Trans 219 (1): 264–272. SPE-1540-G. https://doi.org/
10.2118/1540-G.
Hall, H. N. and Geffen, T. M. 1957. A Laboratory Study of Solvent Flooding. Trans 210 (1): 48–57. SPE-711-G. https://doi.org/10.2118/711-G.
Harper, J. L., Heinemann, R. F., Ray, M. B. et al. 1985. A Compositional Simulator for Performing Large Field Studies in a Vector Computing
Environment. Paper presented at the Middle East Oil Technical Conference and Exhibition, Bahrain. SPE-13714-MS. https://doi.org/10.2118/13714-
MS.
Haugen, K. B. and Beckner, B. L. 2013a. A Critical Comparison of Reduced and Conventional EOS Algorithms. SPE J 18 (2): 378–388. SPE-141399-PA.
https://doi.org/10.2118/141399-PA.
Haugen, K. B. and Beckner, B. L. 2013b. Highly Optimized Phase Equilibrium Calculations. Paper presented at the SPE Reservoir Simulation Symposium,
The Woodlands, Texas, USA. SPE-163583-MS. https://doi.org/10.2118/163583-MS.
Haukas, J., Aavatsmark, I., Espedal, M. et al. 2007. A New IMPSAT Formulation for Compositional Simulation. Paper presented at the SPE Reservoir
Simulation Symposium, Houston, Texas, USA, 26–28 February. SPE-106235-MS. https://doi.org/10.2118/106235-MS.
Heidari, M. and Maini, B. B. 2014. Equation of State Based Simulation of Hybrid and Thermal Processes. Paper presented at the SPE Heavy Oil
Conference-Canada, Calgary, Alberta, Canada, 10–12 June. SPE-170088-MS. https://doi.org/10.2118/170088-MS.
Helfferich, F. G. 1981. Theory of Multicomponent, Multiphase Displacement in Porous Media. SPE J. 21 (01): 51–62. SPE-8372-PA. https://doi.org/10.
2118/8372-PA.
Hemanth-Kumar, K. and Young, L. C. 1996. Parallel Reservoir Simulator Computations. SPE Comp App 8 (02): 50–53. SPE-29104-PA. https://doi.org/
10.2118/29104-PA.
Heroux, M. A., Bartlett, R. A., Howle, V. E. et al. 2005. An Overview of the Trilinos Project. ACM Trans Math Softw 31 (3): 397–423. https://doi.org/
10.1145/1089014.1089021.
Hill, W. J., Tinney, T. J., Young, L. C. et al. 1994. CO2 Operating Plan, South Welch Unit, Dawson County, Texas. Paper presented at the Permian Basin
Oil and Gas Recovery Conference, Midland, Texas, USA. SPE-27676-MS. https://doi.org/10.2118/27676-MS.
Hsu, C.-F., Morell, J. I., and Falls, A. H. 1997. Field-Scale CO2 Flood Simulations and Their Impact on the Performance of the Wasson Denver Unit. SPE
Iranshahr, A., Chen, Y., and Voskov, D. V. 2014. A Coarse-Scale Compositional Model. Comput Geosci 18 (5): 797–815. S10596-014-9427-X. https://doi.
org/10.1007/s10596-014-9427-x.
Iranshahr, A., Voskov, D. V., and Tchelepi, H. A. 2010a. Generalized Negative-Flash Method for Multiphase Multicomponent Systems. Fluid Phase
Equilibria 299 (2): 272–284. https://doi.org/10.1016/j.fluid.2010.09.022.
Iranshahr, A., Voskov, D. V., and Tchelepi, H. A. 2010b. Tie-Simplex Parameterization for EOS-Based Thermal Compositional Simulation. SPE 15 (2):
Iranshahr, A., Voskov, D. V., and Tchelepi, H. A. 2013. A Negative-Flash Tie-Simplex Approach for Multiphase Reservoir Simulation. SPE J 18 (6):
Jacoby, R. H. and Berry, V. J. 1957. A Method for Predicting Depletion Performance of a Reservoir Producing Volatile Crude Oil. Trans 210 (1): 27–33.
SPE-688-G. https://doi.org/10.2118/688-G.
Jerauld, G. R. 1997. General Three-Phase Relative Permeability Model for Prudhoe Bay. SPE Reservoir Eng 12 (4): 255–263. SPE-36178-PA. https://doi.
org/10.2118/36178-PA.
Jhaveri, B. S. and Youngren, G. K. 1988. Three-Parameter Modification of the Peng-Robinson Equation of State To Improve Volumetric Predictions. SPE
John, A. K., Lake, L. W., Bryant, S. L. et al. 2008. Investigation of Field Scale Dispersion. Paper presented at the SPE Symposium on Improved Oil
Recovery, Tulsa, Oklahoma, USA, 20–23 April. SPE-113429-MS. https://doi.org/10.2118/113429-MS.
Johnson, O. G. 1970. Reservoir Simulation Using a General Processor and an Array Processor in Parallel. Paper presented at the SPE Symposium on
Numerical Simulation, Dallas, Texas, USA, 5–6 February. SPE-2813-MS. https://doi.org/10.2118/2813-MS.
Kazemi, H., Vestal, C. R., and Shank, D. G. 1978. An Efficient Multicomponent Numerical Simulator. SPE J. 18 (5): 355–368. SPE-6890-PA. https://doi.
org/10.2118/6890-PA.
40 2022 SPE Journal

Kelkar, M. and Gupta, S. P. 1991. A Numerical Study of Viscous Instabilities: Effect of Controlling Parameters and Scaling Considerations. SPE Res Eng
6 (01): 121–128. SPE-18094-PA. https://doi.org/10.2118/18094-PA.
Kenyon, D. E. and Behie, A. 1987. Third SPE Comparative Solution Project: Gas Cycling of Retrograde Condensate Reservoirs. J Pet Technol 39 (08):
Khan, S. A., Pope, G. A., and Sepehrnoori, K. 1992. Fluid Characterization of Three-Phase CO2/Oil Mixtures. Paper presented at the SPE/DOE Enhanced
Oil Recovery Symposium, Tulsa, Oklahoma. SPE-24130-MS. https://doi.org/10.2118/24130-MS.
Khorsandi, S., Li, L., and Johns, R. T. 2017. Equation-of-State Approach to Model Relative Permeability Including Hysteresis and Wettability Alteration.
Paper presented at the SPE Reservoir Simulation Conference, Montgomery, Texas, USA, 20–22 February. SPE-182655-MS. https://doi.org/
10.2118/182655-MS.
Khorsandi, S., Li, L., and Johns, R. T. 2018. A New Way of Compositional Simulation Without Phase Labeling. Paper presented at the SPE Improved Oil
Recovery Conference, Tulsa, Oklahoma, USA, 14–18 April. SPE-190269-MS. https://doi.org/10.2118/190269-MS.
Killough, J. E. 1976. Reservoir Simulation With History-Dependent Saturation Functions. SPE J. 16 (01): 37–48. SPE-5106-PA. https://doi.org/10.2118/
5106-PA.
Killough, J. E. and Bhogeswara, R. 1991. Simulation of Compositional Reservoir Phenomena on a Distributed-Memory Parallel Computer. J Pet Technol
Killough, J. E., Foster, J. A., Nolen, J. et al. 1996. Parallelization of a General-Purpose Reservoir Simulator. Paper presented at the ECMOR V - 5th
European Conference on the Mathematics of Oil Recovery, Houten, The Netherlands‎. cp-101-00003. https://doi.org/10.3997/2214-4609.201406864.
Killough, J. E. and Kossack, C. A. 1987. Fifth Comparative Solution Project: Evaluation of Miscible Flood Simulators. Paper presented at the SPE
Symposium on Reservoir Simulation, San Antonio, Texas, USA, 1–4 February. SPE-16000-MS. https://doi.org/10.2118/16000-MS.
Killough, J. E. and Levesque, J. M. 1982. Reservoir Simulation and the In-House Vector Processor: Experience for the First Year. Paper presented at the
SPE Reservoir Simulation Symposium, New Orleans, Louisiana, USA, 31 January–3 February. SPE-10521-MS. https://doi.org/10.2118/10521-MS.
Killough, J. E. and Wheeler, M. F. 1987. Parallel Iterative Linear Equation Solvers: An Investigation of Domain Decomposition Algorithms for Reservoir
Simulation. Paper presented at the SPE Symposium on Reservoir Simulation, San Antonio, Texas, USA, 1–4 February. SPE-16021-MS. https://doi.
org/10.2118/16021-MS.
Kniazeff, V. J. and Naville, S. A. 1965. Two-Phase Flow of Volatile Hydrocarbons. SPE J. 5 (01): 37–44. SPE-962-PA. https://doi.org/10.2118/962-PA.
Koval, E. J. 1963. A Method for Predicting the Performance of Unstable Miscible Displacement in Heterogeneous Media. SPE J. 3 (02): 145–154. SPE-
Lacey, J. W., Draper, A. L., and Binder, G. G. 1958. Miscible Fluid Displacement in Porous Media. Trans 213 (1): 76–79. SPE-903-G. https://doi.org/
10.2118/903-G.
Lacroix, S., Vassilevski, Y., Wheeler, J. et al. 2003. Iterative Solution Methods for Modeling Multiphase Flow in Porous Media Fully Implicitly. SIAM J
Sci Comput 25 (3): 905–926. https://doi.org/10.1137/S106482750240443X.
Lake, L. W. and Hirasaki, G. J. 1981. Taylor’s Dispersion in Stratified Porous Media. SPE J. 21 (04): 459–468. SPE-8436-PA. https://doi.org/10.2118/
8436-PA.
Lantz, R. B. 1970. Rigorous Calculation of Miscible Displacement Using Immiscible Reservoir Simulators. SPE J 10 (2): 192–202. SPE-2594-PA. https://
doi.org/10.2118/2594-PA.
Lantz, R. B. 1971. Quantitative Evaluation of Numerical Diffusion (Truncation Error). SPE J 11 (3): 315–320. SPE-2811-PA. https://doi.org/10.2118/2811-
PA.
Li, H. and Durlofsky, L. J. 2016. Upscaling for Compositional Reservoir Simulation. SPE J. 21 (03): 0873–0887. SPE-173212-PA. https://doi.org/10.2118/
173212-PA.
Li, X. S. 2005. An Overview of SuperLU: Algorithms, Implementation, and User Interface. ACM Trans Math Softw 31 (3): 302–325. https://doi.org/
10.1145/1089014.1089017.
Lins, A. G., Nghiem, L. X., and Harding, T. G. 2011. Three-Phase Hydrocarbon Thermodynamic Liquid-Liquid-Vapour Equilibrium in CO2 Process. Paper
presented at the SPE Reservoir Characterisation and Simulation Conference and Exhibition, Abu Dhabi, UAE, 9–11 October. SPE-148040-MS. https://
doi.org/10.2118/148040-MS.
Litvak, M. L. 1994. New Procedure for the Phase-Equilibrium Computations in the Compositional Reservoir Simulator. SPE Advanced Technology Series
Liu, H., Wang, K., Chen, Z. et al. 2015a. A Parallel Framework for Reservoir Simulators on Distributed-Memory Supercomputers. Paper presented at the
SPE/IATMI Asia Pacific Oil & Gas Conference and Exhibition, Nusa Dua, Bali, Indonesia. SPE-176045-MS. https://doi.org/10.2118/176045-MS.
Liu, H., Wang, K., Chen, Z. et al. 2015b. Development of Parallel Reservoir Simulators on Distributed-Memory Supercomputers. Paper presented at the
SPE Reservoir Characterisation and Simulation Conference and Exhibition, Abu Dhabi, UAE. SPE-175573-MS. https://doi.org/10.2118/175573-MS.
Lohrenz, J., Bray, B. G., and Clark, C. R. 1964. Calculating Viscosities of Reservoir Fluids From Their Compositions. J Pet Technol 16 (10): 1171–1176.
de Loubens, R., Riaz, A., and Tchelepi, H. A. 2009. Error Analysis of an Adaptive Implicit Scheme for Hyperbolic Conservation Laws. SIAM J Sci Comput
31 (4): 2890–2914. https://doi.org/10.1137/080724502.
Lu, P., Shaw, J. S., Eccles, T. K. et al. 2009. Experience With Numerical Stability, Formulation, and Parallel Efficiency of Adaptive Implicit Methods. Paper
presented at the SPE Reservoir Simulation Symposium, The Woodlands, Texas, USA, 2–4 February. SPE-118977-MS. https://doi.org/10.2118/118977-
MS.
Maes, J., Moncorgé, A., and Tchelepi, H. A. 2013. Thermal Adaptive Implicit Method: Time Step Selection. J Pet Sci Eng 106: 34–45. https://doi.org/
10.1016/j.petrol.2013.03.019.
Malachowski, M. A., Stephenson, R. E., Oppenheim, J. P. et al. 1990. Vectorization of Fluid Property Calculations within a Petroleum Reservoir Simulator.
Comput Chem Eng 14 (6): 655–663. https://doi.org/10.1016/0098-1354(90)87034-M.
Maliassov, S., Beckner, B. L., and Dyadechko, V. 2013. Parallel Reservoir Simulation Using a Specific Software Framework. Paper presented at the SPE
Martin, J. C. 1959. Simplified Equations of Flow in Gas Drive Reservoirs and the Theoretical Foundation of Multiphase Pressure Buildup Analyses. Trans
216 (1): 321–323. SPE-1235-G. https://doi.org/10.2118/1235-G.
Martin, P., Romain, D. L., and Leonardo, P. 2019. Continuous Relative Permeability Model for Compositional Reservoir Simulation, Using the True
Critical Point and Accounting for Miscibility. Paper presented at the SPE Reservoir Simulation Conference, Galveston, Texas, USA, 10–11 April. SPE-
193826-MS. https://doi.org/10.2118/193826-MS.
Meijerink, J. A., Van Daalen, D. T., Hoogerbrugge, P. J. et al. 1991. Towards a More Effective Parallel Reservoir Simulator. Paper presented at the SPE
Symposium on Reservoir Simulation, Anaheim, California, USA. SPE-21212-MS. https://doi.org/10.2118/21212-MS.
2022 SPE Journal 41

Metcalfe, R. S., Vogel, J. L., and Morris, R. W. 1988. Compositional Gradients in the Anschutz Ranch East Field. SPE Res Eng 3 (3): 1025–1032. SPE-
Michelsen, M. L. 1982a. The Isothermal Flash Problem. Part I. Stability. Fluid Phase Equilib 9 (1): 1–19. https://doi.org/10.1016/0378-3812(82)85001-2.
Michelsen, M. L. 1982b. The Isothermal Flash Problem. Part II. Phase-Split Calculation. Fluid Phase Equilib 9 (1): 21–40. https://doi.org/10.1016/0378-
3812(82)85002-4.
Michelsen, M. L. 1986. Simplified Flash Calculations for Cubic Equations of State. Ind Eng Chem Proc Des Dev 25 (1): 184–188. https://doi.org/10.1021/
i200032a029.
Michelsen, M. L. and Mollerup, J. 1986. Partial Derivatives of Thermodynamic Properties. AIChE J 32 (8): 1389–1392. https://doi.org/10.1002/
aic.690320818.
Michelsen, M. L. and Mollerup, J. 2004. Thermodynamic Modelling: Fundamentals and Computational Aspects. Holte, Denmark: Tie-Line Publications.
Michelsen, M. L., Yan, W., and Stenby, E. H. 2013. A Comparative Study of Reduced-Variables-Based Flash and Conventional Flash. SPE J. 18 (05):
Mifflin, R. T., Watts, J. W., and Weiser, A. 1991. A Fully Coupled, Fully Implicit Reservoir Simulator for Thermal and Other Complex Reservoir Processes.
Paper presented at the SPE Symposium on Reservoir Simulation, Anaheim, California, USA. SPE-21252-MS. https://doi.org/10.2118/21252-MS.
Moncorgé, A. 2018. Sequential Fully Implicit Methods for Multiscale Modeling of Compositional Flows. ETH Zurich. https://doi.org/org/10.3929/
ethz-b-000314439.
Moncorgé, A., Møyner, O., Tchelepi, H. A. et al. 2019. Consistent Upwinding for Sequential Fully Implicit Multiscale Compositional Simulation. Comput
Geosci 24 (2): 533–550. S10596-019-09835-6. https://doi.org/10.1007/s10596-019-09835-6.
Moncorgé, A. and Tchelepi, H. A. 2009. Stability Criteria for Thermal Adaptive Implicit Compositional Flows. SPE J. 14 (2): 311–322. SPE-111610-PA.
https://doi.org/10.2118/111610-PA.
Moortgat, J. 2017. Adaptive Implicit Finite Element Methods for Multicomponent Compressible Flow in Heterogeneous and Fractured Porous Media.
Water Resour Res 53 (1): 73–92. https://doi.org/10.1002/2016WR019644.
Møyner, O. and Tchelepi, H. A. 2018. A Mass-Conservative Sequential Implicit Multiscale Method for Isothermal Equation-of-State Compositional
Problems. SPE J. 23 (06): 2376–2393. SPE-182679-PA. https://doi.org/10.2118/182679-PA.
Negahban, S. and Kremesec, V. J. 1992. Development and Validation of Equation-of-State Fluid Descriptions for CO2/Reservoir-Oil Systems. SPE Res
Eng 7 (3): 363–368. SPE-19637-PA. https://doi.org/10.2118/19637-PA.
Neshat, S. S. and Pope, G. A. 2018. Three-Phase Relative Permeability and Capillary Pressure Models With Hysteresis and Compositional Consistency.
Nghiem, L. X., Fong, D. K., and Aziz, K. 1981. Compositional Modeling With an Equation of State (Includes Associated Papers 10894 and 10903). SPE
J. 21 (06): 687–698. SPE-9306-PA. https://doi.org/10.2118/9306-PA.
Nghiem, L. X. and Li, Y.-K. 1986. Effect of Phase Behavior on CO2 Displacement Efficiency at Low Temperatures: Model Studies With an Equation of
State. SPE Res Eng 1 (04): 414–422. SPE-13116-PA. https://doi.org/10.2118/13116-PA.
Nolen, J. S. 1973. Numerical Simulation Of Compositional Phenomena In Petroleum Reservoirs. Paper presented at the SPE Symposium on Numerical
Simulation of Reservoir Performance, Houston, Texas, USA, 11–12 January. SPE-4274-MS. https://doi.org/10.2118/4274-MS.
Nolen, J. S., Kuba, D. W., and Kascic, M. J. 1981. Application of Vector Processors To Solve Finite Difference Equations. SPE J 21 (4): 447–453. SPE-
Nolen, J. S. and Stanat, P. L. 1981. Reservoir Simulation On Vector Processing Computers. Paper presented at the Middle East Technical Conference and
Exhibition, Bahrain. SPE-9644-MS. https://doi.org/10.2118/9644-MS.
Notay, Y. 2019. “User’s Guide to AGMG, Revision 3.3.5.” Brussels, Belgium.
Odeh, A. S. 1981. Comparison of Solutions to a Three-Dimensional Black-Oil Reservoir Simulation Problem (Includes Associated Paper 9741). J Pet
Technol 33 (1): 13–25. SPE-9723-PA. https://doi.org/10.2118/9723-PA.
Okuno, R., Johns, R. T., and Sepehrnoori, K. 2010. Three-Phase Flash in Compositional Simulation Using a Reduced Method. SPE J 15 (3): 689–703.
Okuno, R. and Xu, Z. 2014. Mass Transfer on Multiphase Transitions in Low-Temperature Carbon-Dioxide Floods. SPE J. 19 (06): 1005–1023. SPE-
Orr, F. M. 2007. Theory of Gas Injection Processes. Copenhagen, Denmark: Tie-Line Publications.
Pan, H., Chen, Y., Sheffield, J. et al. 2015. Phase-Behavior Modeling and Flow Simulation for Low-Temperature CO2 Injection. SPE Res Eval & Eng 18
(2): 250–263. SPE-170903-PA. https://doi.org/10.2118/170903-PA.
Parashar, M., Wheeler, J. A., Pope, G. et al. 1997. A New Generation EOS Compositional Reservoir Simulator: Part II - Framework and Multiprocessing.
Paper presented at the SPE Reservoir Simulation Symposium, Dallas, Texas, USA, 8–11 June. SPE-37977-MS. https://doi.org/10.2118/37977-MS.
Parashar, M. and Yotov, I. 1998. An Environment for Parallel Multi-Block, Multi-Resolution Reservoir Simulations. Paper presented at the Proceedings of
the 11th International Conference on Parallel and Distributed Computing and Systems, 230–235, Chicago.
Peaceman, D. W. 1977. A Nonlinear Stability Analysis for Difference Equations Using Semi-Implicit Mobility. SPE J. 17 (01): 79–91. SPE-5735-PA.
Peaceman, D. W. 2000. Fundamentals of Numerical Reservoir Simulation. United States: Elsevier Science.
Péneloux, A., Rauzy, E., and Fréze, R. 1982. A Consistent Correction for Redlich-Kwong-Soave Volumes. Fluid Phase Equilibria 8 (1): 7–23. https://doi.
org/10.1016/0378-3812(82)80002-2.
Peng, D.-Y. and Robinson, D. B. 1976. A New Two-Constant Equation of State. Ind Eng Chem Fund 15 (1): 59–64. https://doi.org/10.1021/i160057a011.
Perkins, T. K. and Johnston, O. C. 1963. A Review of Diffusion and Dispersion in Porous Media. SPE J. 3 (01): 70–84. SPE-480-PA. https://doi.org/10.
2118/480-PA.
Perschke, D. R., Chang, Y., Pope, G. A. et al. 1989. Comparison Of Phase Behavior Algorithms For An Equation-Of-State Compositional Simulator.
Moscow, Kuala Lumpur: Society of Petroleum Engineers.
Petitfrere, M. and Nichita, D. V. 2015. A Comparison of Conventional and Reduction Approaches for Phase Equilibrium Calculations. Fluid Phase
Equilibria 386: 30–46. https://doi.org/10.1016/j.fluid.2014.11.017.
Petitfrere, M., Nichita, D. V., Voskov, D. V. et al. 2020. Full-EoS Based Thermal Multiphase Compositional Simulation of CO2 and Steam Injection
Processes. J Pet Sci Eng 192: 107241. https://doi.org/10.1016/j.petrol.2020.107241.
Potempa, T. 1982. Finite Element Methods for Convection Dominated Transport Problems. Houston, Texas, USA: Rice University.
Prausnitz, J. M. and Lichtenthaler, R. N. 2004. Molecular Thermodynamics of Fluid-Phase Equilibria. Taiwan: Pearson Education Taiwan.
Press, W. H. and Vetterling, W. T. 1999. Numerical Recipes in C: The Art of Scientific Computing. Cambridge: Cambridge Univ. Press.
Price, H. W. and Donohue, D. A. T. 1967. Isothermal Displacement Processes With Interphase Mass Transfer. SPE J. 7 (02): 205–220. SPE-1533-PA.
42 2022 SPE Journal

Purswani, P., Johns, R. T., Karpyn, Z. T. et al. 2020. Predictive Modeling of Relative Permeability Using a Generalized Equation of State. SPE J. 26 (01):
Qiao, C., Khorsandi, S., and Johns, R. T. 2017. A General Purpose Reservoir Simulation Framework for Multiphase Multicomponent Reactive Fluids. Paper
presented at the SPE Reservoir Simulation Conference, Montgomery, Texas, USA, 20–22 February. SPE-182715-MS. https://doi.org/10.2118/182715-
MS.
Quandalle, P. and Savary, D. 1989. An Implicit in Pressure and Saturations Approach to Fully Compositional Simulation. Paper presented at the SPE
Symposium on Reservoir Simulation, Houston, Texas, USA. https://doi.org/10.2118/18423-MS.
Rachford, H. H. and Rice, J. D. 1952. Procedure for Use of Electronic Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium. J
Pet Technol 4 (10): 19–3. SPE-952327-G. https://doi.org/10.2118/952327-G.
Rasmussen, C. P., Krejbjerg, K., Michelsen, M. L. et al. 2006. Creasing the Computational Speed of Flash Calculations With Applications for Compositional,
Transient Simulations. SPE Res Eval & Eng 9 (01): 32–38. SPE-84181-PA. https://doi.org/10.2118/84181-PA.
Redlich, O. and Kwong, J. N. S. 1949. On the Thermodynamics of Solutions; an Equation of State; Fugacities of Gaseous Solutions. Chem Rev 44 (1):
233–244, . https://doi.org/10.1021/cr60137a013.
Renner, T. A., Metcalfe, R. S., Yellig, W. F. et al. 1989. Displacement of a Rich Gas Condensate by Nitrogen: Laboratory Corefloods and Numerical
Simulations. SPE Res Eng 4 (01): 52–58. SPE-16714-PA. https://doi.org/10.2118/16714-PA.
Rezaveisi, M., Sepehrnoori, K., Pope, G. A. et al. 2015. Compositional Simulation Including Effect of Capillary Pressure on Phase Behavior. Paper
presented at the SPE Annual Technical Conference and Exhibition, Houston, Texas, USA, 28–30 September. SPE-175135-MS. https://doi.org/
10.2118/175135-MS.
Rezaveisi, M., Sepehrnoori, K., Pope, G. A. et al. 2018. Thermodynamic Analysis of Phase Behavior at High Capillary Pressure. SPE J 23 (6): 1977–1990.
Roebuck, I. F., Ford, W. T., Henderson, G. E. et al. 1968. The Compositional Reservoir Simulator; Case IV The Two-Dimensional Model. Paper presented
at the Fall Meeting of the Society of Petroleum Engineers of AIME, Houston, Texas, USA, 29 September–2 October. SPE-2235-MS. https://doi.org/
10.2118/2235-MS.
Roebuck, I. F., Henderson, G. E., Douglas, J. et al. 1969a. The Compositional Reservoir Simulator: Case II - The Linear Model. SPE J. 9: 115–130. SPE-
Roebuck, I. F., Henderson, G. E., Douglas, J. et al. 1969b. The Compositional Reservoir Simulator: Case III - The Linear Model. SPE J. 9: 115–130. SPE-
Roebuck, I. F., Henderson, G. E., Douglas, J. et al. 1969c. The Compositional Reservoir Simulator: Case I - The Linear Model. SPE J 9 (1): 115–130.
Rowe, A. M. 1978. Internally Consistent Correlations For Predicting Phase Compositions For Use In Reservoir Composition Simulators. Paper presented
at the SPE Annual Fall Technical Conference and Exhibition, Houston, Texas, USA, 1–3 October. SPE-7475-MS. https://doi.org/10.2118/7475-MS.
Russell, T. F. 1989. Stability Analysis and Switching Criteria for Adaptive Implicit Methods Based on the CFL Condition. Paper presented at the SPE
Symposium on Reservoir Simulation, Houston, Texas, USA. SPE-18416-MS. https://doi.org/10.2118/18416-MS.
Russell, T. F. and Wheeler, M. F. 1983. Finite Element and Finite Difference Methods for Continuous Flows in Porous Media. In The Mathematics of
Reservoir Simulation, 35–106. Philadelphia, Pennsylvania, USA: Frontiers in Applied Mathematics. Society for Industrial and Applied Mathematics.
https://doi.org/10.1137/1.9781611971071.
Sandoval, D. R., Yan, W., Michelsen, M. L. et al. 2016. The Phase Envelope of Multicomponent Mixtures in the Presence of a Capillary Pressure
Difference. Ind Eng Chem Res 55 (22): 6530–6538. https://doi.org/10.1021/acs.iecr.6b00972.
Schmall, L., Varavei, A., and Sepehrnoori, K. 2013. A Comparison of Various Formulations for Compositional Reservoir Simulation. Paper presented
at the SPE Reservoir Simulation Symposium, The Woodlands, Texas, USA, 18–20 February. SPE-163630-MS. https://doi.org/10.2118/163630-MS.
Scott, S. L., Wainwright, R. L., Raghavan, R. et al. 1987. Application of Parallel (MIMD) Computers to Reservoir Simulation. Paper presented at the SPE
Symposium on Reservoir Simulation, San Antonio, Texas, USA, 1–4 February. SPE-16020-MS. https://doi.org/10.2118/16020-MS.
Scovazzi, G., Wheeler, M. F., Mikelić, A. et al. 2017. Analytical and Variational Numerical Methods for Unstable Miscible Displacement Flows in Porous
Media. J Comput Phys 335: 444–496. https://doi.org/10.1016/j.jcp.2017.01.021.
Sheffield, M. 1969. Three-Phase Fluid Flow Including Gravitational, Viscous and Capillary Forces. SPE J. 9 (02): 255–269. SPE-2012-PA. https://doi.org/
10.2118/2012-PA.
Sheldon, J. W., Harris, C. D., and Bavly, D. 1960. A Method for General Reservoir Behavior Simulation on Digital Computers. Paper presented at the Fall
Meeting of the Society of Petroleum Engineers of AIME, Denver, Colorado, USA, 2–5 October. SPE-1521-G. https://doi.org/10.2118/1521-G.
Shelton, J. L. and Yarborough, L. 1977. Multiple Phase Behavior in Porous Media During CO2 or Rich-Gas Flooding. J Pet Technol 29 (9): 1171–1178.
Sheth, S. M. and Younis, R. M. 2017. Localized Linear Systems in Sequential Implicit Simulation of Two-Phase Flow and Transport. SPE J. 22 (05):
Shiralkar, G. S. and Stephenson, R. E. 1991. A General Formulation for Simulating Physical Dispersion and a New Nine-Point Scheme. SPE Res Eng 6
(01): 115–120. SPE-16975-PA. https://doi.org/10.2118/16975-PA.
Shiralkar, G. S., Stephenson, R. E., Joubert, W. et al. 1998. Falcon: A Production Quality Distributed Memory Reservoir Simulator. SPE Res Eval & Eng
Shrivastava, V. K., Nghiem, L. X., Moore, R. G. et al. 2005a. Modelling Physical Dispersion in Miscible Displacement-Part 1: Theory and the Proposed
Numerical Scheme. J Can Pet Technol 44 (5): 9. PETSOC-05-05-01. https://doi.org/10.2118/05-05-01.
Shrivastava, V. K., Nghiem, L. X., Moore, R. G. et al. 2005b. Modelling Physical Dispersion in Miscible Displacement-Part 2: Validation, Numerical Tests,
and Applications. J Can Pet Technol 44 (5). PETSOC-05-05-02. https://doi.org/10.2118/05-05-02.
Shubin, G. R. and Bell, J. B. 1984. An Analysis of the Grid Orientation Effect in Numerical Simulation of Miscible Displacement. Comput Methods Appl
Mech Eng 47 (1–2): 47–71. https://doi.org/10.1016/0045-7825(84)90047-1.
Spillette, A. G., Hillestad, J. G., and Stone, H. L. 1973. A High-Stability Sequential Solution Approach to Reservoir Simulation. Paper presented at the Fall
Meeting of the Society of Petroleum Engineers of AIME, Las Vegas, Nevada, 30 September–3 October. SPE-4542-MS. https://doi.org/10.2118/4542-
MS.
Spivak, A. and Coats, K. H. 1970. Numerical Simulation of Coning Using Implicit Production Terms. SPE J. 10 (03): 257–267. SPE-2595-PA. https://doi.
org/10.2118/2595-PA.
Stone, H. L. 1970. Probability Model for Estimating Three-Phase Relative Permeability. J Pet Technol 22 (2): 214–218. SPE-2116-PA. https://doi.org/
10.2118/2116-PA.
2022 SPE Journal 43

Stone, H. L. 1973. Estimation of Three-Phase Relative Permeability And Residual Oil Data. J Can Pet Technol 12 (4). PETSOC-73-04-06. https://doi.org/
10.2118/73-04-06.
Stone, H. L. and Garder, A. O. 1961. Analysis of Gas-Cap or Dissolved-Gas Drive Reservoirs. SPE J 1 (2): 92–104. SPE-1518-G. https://doi.org/
10.2118/1518-G.
Stueben, K., Clees, T., Klie, H. et al. 2007. Algebraic Multigrid Methods (AMG) for the Efficient Solution of Fully Implicit Formulations in Reservoir
Simulation. Paper presented at the SPE Reservoir Simulation Symposium, Houston, Texas, USA, 26–28 February. SPE-105832-MS. https://doi.org/
10.2118/105832-MS.
Subramanian, G., Trangenstein, J. A., Mochizuki, S. et al. 1987. Efficient Fluid Behavior Computations in a Sequential Compositional Reservoir
Simulator. Paper presented at the SPE Symposium on Reservoir Simulation, San Antonio, Texas, USA, 1–4 February. SPE-16024-MS. https://doi.org/
10.2118/16024-MS.
van Daalen, D. T., J. Hoogerbrugge, P., A. Meijerink, J. et al. 1989. The Parallelisation of Bosim, Shell’s Black/Volatile Oil Reservoir Simulator. In. Paper
presented at the ECMOR I - 1st European Conference on the Mathematics of Oil Recovery, Houten, The Netherlands‎. cp-234-00019. https://doi.org/
10.3997/2214-4609.201411321.
Tan, T. B. 1988. Implementation of an Improved Adaptive-Implicit Method In a Thermal Compositional Simulator. SPE Res Eng 3 (04): 1123–1128. SPE-
Tang, D. E. and Zick, A. A. 1993. A New Limited Compositional Reservoir Simulator. Paper presented at the SPE Symposium on Reservoir Simulation,
New Orleans, Louisiana, USA. SPE-25255-MS. https://doi.org/10.2118/25255-MS.
Thele, K. J., Lake, L. W., and Sepehrnoori, K. 1983. A Comparison of Three Equation-of-State Compositional Simulators. Paper presented at the SPE
Reservoir Simulation Symposium, San Francisco, California, USA, 15–18 November. SPE-12245-MS. https://doi.org/10.2118/12245-MS.
Thomas, G. W. and Thurnau, D. H. 1983. Reservoir Simulation Using an Adaptive Implicit Method. SPE J 23 (5): 759–768. SPE-10120-PA. https://doi.
org/10.2118/10120-PA.
Thomas, L. K., Lumpkin, W. B., and Reheis, G. M. 1976. Reservoir Simulation of Variable Bubble-Point Problems. SPE J. 16 (01): 10–16. SPE-5107-PA.
Todd, M. R. and Longstaff, W. J. 1972. The Development, Testing, and Application Of a Numerical Simulator for Predicting Miscible Flood Performance.
J. Pet. Technol 24 (07): 874–882. SPE-3484-PA. https://doi.org/10.2118/3484-PA.
Todd, M. R., O'Dell, P. M., and Hirasaki, G. J. 1972. Methods for Increased Accuracy in Numerical Reservoir Simulators. SPE J. 12 (06): 515–530. SPE-
Treybal, R. E. 1968. Mass-Transfer Operations. New York: McGraw-Hill.
Van-Quy, N., Simandoux, P., and Corteville, J. 1972. A Numerical Study of Diphasic Multicomponent Flow. SPE J. 12 (02): 171–184. SPE-3006-PA.
Vidyasagar, A., Patacchini, L., Panfili, P. et al. 2019. Full-GPU Reservoir Simulation Delivers on Its Promise for Giant Carbonate Fields. In Third EAGE
WIPIC Workshop: ReservoirManagement in Carbonates, Nov 2019, Volume 2019. Doha, Qatar: European Association of Geoscientists & Engineers.
https://doi.org/10.3997/2214-4609.201903118.
Voskov, D. V. and Tchelepi, H. A. 2009a. Compositional Space Parameterization: Multicontact Miscible Displacements and Extension to Multiple Phases.
Voskov, D. V. and Tchelepi, H. A. 2009b. Compositional Space Parameterization: Theory and Application for Immiscible Displacements. SPE J. 14 (3):
Voskov, D. V. and Tchelepi, H. A. 2012. Comparison of Nonlinear Formulations for Two-Phase Multi-Component EoS Based Simulation. J Pet Sci Eng
82–83: 101–111. https://doi.org/10.1016/j.petrol.2011.10.012.
Voskov, D. V., Tchelepi, H. A., and Younis, R. M. 2009. General Nonlinear Solution Strategies for Multiphase Multicomponent EoS Based Simulation.
Paper presented at the SPE Reservoir Simulation Symposium, The Woodlands, Texas, USA, 2–4 February. SPE-118996-MS. https://doi.org/
10.2118/118996-MS.
Wallis, J. R., Kendall, R. P., and Little, T. E. 1985. Constrained Residual Acceleration of Conjugate Residual Methods. Paper presented at the SPE
Reservoir Simulation Symposium, Dallas, Texas, USA. SPE-13536-MS. https://doi.org/10.2118/13536-MS.
Wan, J., Sarma, P., Usadi, A. K. et al. 2005. General Stability Criteria for Compositional and Black-Oil Models. Paper presented at the SPE Reservoir
Simulation Symposium, The Woodlands, Texas, USA, 31 January–2 February. SPE-93096-MS. https://doi.org/10.2118/93096-MS.
Wang, P., Balay, S., Sepehrnoori, K. et al. 1999. A Fully Implicit Parallel EOS Compositional Simulator for Large Scale Reservoir Simulation. Paper
presented at the SPE Reservoir Simulation Symposium, Houston, Texas, USA, 14–17 February. SPE-51885-MS. https://doi.org/10.2118/51885-MS.
Wang, P. and Barker, J. W. 1995. Comparison of Flash Calculations in Compositional Reservoir Simulation. Paper presented at the SPE Annual Technical
Conference and Exhibition, Dallas, Texas, USA, 22–25 October. SPE-30787-MS. https://doi.org/10.2118/30787-MS.
Wang, P., Yotov, I., Wheeler, M. et al. 1997. A New Generation EOS Compositional Reservoir Simulator: Part I - Formulation and Discretization. Paper
presented at the SPE Reservoir Simulation Symposium, Dallas, Texas, USA, 8–11 June. SPE-37979-MS. https://doi.org/10.2118/37979-MS.
Wattenbarger, R. A. 1970. Practical Aspects of Compositional Simulation. Paper presented at the SPE Symposium on Numerical Simulation, Dallas, Texas,
Watts, J. W. 1986. A Compositional Formulation of the Pressure and Saturation Equations. SPE Res Eng 1 (03): 243–252. SPE-12244-PA. https://doi.org/
10.2118/12244-PA.
Watts, J. W. and Rame, M. 1999. An Algebraic Approach to the Adaptive Implicit Method. Paper presented at the SPE Reservoir Simulation Symposium,
Houston, Texas, USA. SPE-51900-MS. https://doi.org/10.2118/51900-MS.
Watts, J. W. and Silliman, W. J. 1980. Numerical Dispersion and the Origins of the Grid Orientation Sensitivity. Paper presented at the 73rd Meeting of
AICHE.
Wendschlag, D. D., Stephenson, R. E., and Clark, T. J. 1983. SPE Reservoir Simulation Symposium. Paper presented at the Fieldwide Simulation of the
Anschutz Ranch East Nitrogen Injection Project With a Generalized Compositional Model, San Francisco, California, USA. SPE-12257-MS. https://
doi.org/10.2118/12257-MS.
Wheeler, J. A. and Smith, R. A. 1990. Reservoir Simulation on a Hypercube. SPE Res Eng 5 (4): 544–548. SPE-19804-PA. https://doi.org/10.2118/19804-
PA.
Whitson, C. H. 1991. Improved Simulation of Complex Reservoir Systems, Report No: 5. Reservoir Simulation Research Corp, Tulsa, Oklahoma.
Whitson, C. H. and Brulé, M. R. 2000. Phase Behavior. Richardson, Texas: Monograph 20, Society Petroleum Engrs.
Whitson, C. H. and Michelsen, M. L. 1989. The Negative Flash. Fluid Phase Equilibria 53: 51–71. https://doi.org/10.1016/0378-3812(89)80072-X.
Wong, T. W. and Aziz, K. 1989. Considerations in the Development of Multipurpose Reservoir Simulation Models. In 1st and 2nd International Forums
on Reservoir Simulation. Alpbach, Austria: Paul Steiner Publications.
44 2022 SPE Journal

Wong, T. W., Firoozabadi, A., and Aziz, K. 1990. Relationship of the Volume-Balance Method of Compositional Simulation to the Newton-Raphson
Method. SPE Res Eng 5 (03): 415–422. SPE-18424-PA. https://doi.org/10.2118/18424-PA.
Wong, T. W., Firoozabadi, A., Nutakki, R. et al. 1987. A Comparison of Two Approaches to Compositional and Black Oil Simulation. Paper presented at
the SPE Symposium on Reservoir Simulation, San Antonio, Texas, USA, 1–4 February. SPE-15999-MS. https://doi.org/10.2118/15999-MS.
Woo, P. T. 1979. Application Of Array Processor To Sparse Elimination. Paper presented at the SPE Reservoir Simulation Symposium, Denver, Colorado,
USA, 31 January–2 February. SPE-7674-MS. https://doi.org/10.2118/7674-MS.
Yanosik, J. L. and McCracken, T. A. 1979. A Nine-Point, Finite-Difference Reservoir Simulator for Realistic Prediction of Adverse Mobility Ratio
Displacements. SPE J. 19 (4): 253–262. SPE-5734-PA. https://doi.org/10.2118/5734-PA.
Yarborough, L. 1979. Application of a Generalized Equation of State to Petroleum Reservoir Fluids. In Equations of State in Engineering and Research,
Vol. 182, 21–385. N.p.: Advances in Chemistry. American Chemical Society. https://doi.org/10.1021/ba-1979-0182.
Young, L. C. 1981. A Finite-Element Method for Reservoir Simulation. SPE J. 21 (01): 115–128. SPE-7413-PA. https://doi.org/10.2118/7413-PA.
Young, L. C. 1984. A Study of Spatial Approximations for Simulating Fluid Displacements in Petroleum Reservoirs. Comput Methods Appl Mech Eng 47
(1–2): 3–46. https://doi.org/10.1016/0045-7825(84)90046-X.
Young, L. C. 1987. Equation of State Compositional Modeling on Vector Processors. Paper presented at the SPE Symposium on Reservoir Simulation, San
Antonio, Texas, USA, 1–4 February. SPE-16023-MS. https://doi.org/10.2118/16023-MS.
Young, L. C. 1989. Improved Simulation of Complex Reservoir Systems, Report 1. Reservoir Simulation Research Corp, Tulsa, Oklahoma.
Young, L. C. 1990. Use of Dispersion Relationships To Model Adverse-Mobility-Ratio Miscible Displacements. SPE Res Eng 5 (3): 309–316. SPE-14899-
Young, L. C. 1991. Full-Field Compositional Modeling on Vector Processors. SPE Res Eng 6 (1): 107–114. SPE-17803-PA. https://doi.org/10.2118/17803-
PA.
Young, L. C. 2001. Continuous Compositional Volume-Balance Equations. Paper presented at the SPE Reservoir Simulation Symposium, Houston, Texas,
Young, L. C. and Hemanth-Kumar, K. 1992. High-Performance Black Oil Computations. SPE Comp App 4 (4): 10–15. SPE-21215-PA. https://doi.org/
10.2118/21215-PA.
Young, L. C. and Russell, T. F. 1993. Implementation of an Adaptive Implicit Method. Paper presented at the SPE Reservoir Simulation Symposium, New
Orleans, Louisiana, USA, 28 February–3 March. 28SPE-25245-MS. https://doi.org/10.2118/25245-MS.
Young, L. C. and Stephenson, R. E. 1983. A Generalized Compositional Approach for Reservoir Simulation. SPE J. 23 (5): 727–742. SPE-10516-PA.
https://doi.org/10.2118/10516-PA.
Young, L. C., Hemanth-Kumar, K., and Bratvold, R. B. 1990. Reservoir Simulation on Low-Cost, High-Performance Workstations. SPE Comp App 2 (5):
Yuan, C. and Pope, G. A. 2012. A New Method To Model Relative Permeability in Compositional Simulators To Avoid Discontinuous Changes Caused by
Phase-Identification Problems. SPE J. 17 (4): 1221–1230. SPE-142093-PA. https://doi.org/10.2118/142093-PA.
Zaydullin, R., Voskov, D. V., and Tchelepi, H. A. 2012. Nonlinear Formulation Based on an Equation-of-State Free Method for Compositional Flow
Simulation. SPE J. 18 (2): 264–273. SPE-146989-PA. https://doi.org/10.2118/146989-PA.
Zick, A. A. 1992. Improved Simulation of Complex Reservoir Systems, Report No 6. Reservoir Simulation Research Corp, Tulsa, Oklahoma.
Appendix—Linearizations
The coefficients needed for linearizing the equations are described here. The equations requiring linearization are the phase equilibria,
overall molar volume, the material balance constraint, and for AIM or FIM, the saturation and composition of the phases. For the first
three, the relationships are:

!i + ln 'vi ln 'ì = 0

h = 1 v ` x` , p + vv xv , p
c 1
(A-1)
NP
xvi xì = 0,
‍ i=1 ‍

where ‍!i = ln Ki = ln xvi /xì .‍ Assume all the properties in these equations have been calculated, and the following derivatives are
available from the EOS:
ln '˛i @ ln '˛i
F˛ij = n˛ @ @n and Fp˛i = @p
@˛
˛j
@˛ (A-2)
˛i = n˛ @n and p˛ = .
‍ ˛i @p ‍

The overall composition is ‍zi = 1 v xì + vxvi ‍from which we define the following:

xì
nì = 1 v ni and nvi = xzvii vni
zi
xì 1 1 (A-3)
‍ zi = 1 + e!i 1 v and xzvii = 1 + e!i 1 1 v .‍

Basic Linearizations. For convenience, define ‍ci = xvi xì /zi , di = xì xvi /zi , wi = di !i and i = ni /n.‍ The basic quantities
linearize to:

zi xzìi = vwi xì ci v
(A-4)
zi xzvii = 1 v wi xvi ci v.
‍ ‍
The individual phase mole numbers then linearize to:
2022 SPE Journal 45

nì
n = v 1 v wi di v + 1 v xzìi i
(A-5)
nvi
n = v 1 v wi + di v + v xzvii i .
‍ ‍
The overall composition is easily expressed in terms of the overall mole numbers:
P
zi = i zi j .
‍ j ‍ (A-6)
The symmetry evident in these expressions is a direct result of choosing the K-value logarithms as variables.
Rachford-Rice Constraint. The constraint equation is basically the linearization of the Rachford-Rice, Eq. 34. Defining the unit vector
‍ei = 1,‍ the linearization is given in Eq. 46, which is repeated here:
‍ eT w Bv v + BTn = Br ,‍ (A-7)
where Br is the equation residual which is negligible after solution and:

Bni = ci Br
NPc 1 (A-8)
Bv = ci xvi xì .
‍ i=1 ‍
Phase Equilibria. The linearized phase equilibria relations from Eq. 46 are:
‍ Aw w + Av v + An + Ap p = Ar ,‍ (A-9)
where Ar is the residual and the other terms are:

ıij
Awij = + 1 v Fvij + vFìj
di
x x
Anij = zvjj Fvij z`jj Fìj
P N (A-10)
Avi = Anij xvj x`j
j=1

@ ln 'vi @ ln 'ì
Api = .
‍ @p @p ‍
The vapor fraction is readily eliminated using the constraint equation, Eq. A-7:

T A BT
Aw + ABv ev w + An + Bv v n + Ap p = Ar
(A-11)
A Q n + Ap p = Ar .
Q ww + A
‍ ‍
The decomposition of this system is calculated with the aid of the Sherman-Morrison theorem, Eq. 50, which is repeated here:

1 A1 Av eT
AQ w = I + w T 1 A1
w = I + UeT A1
w . (A-12)
‍ B v +e Aw A v ‍
An obvious way to improve the scaling and reduce potential rounding errors is to modify the way the denominator term is calculated.
Using the definition of Bv the denominator term is:

‍ Bv + eT A1 T 1
w Av = e Aw Aw + Av ,‍ (A-13)
2
where ‍ i = zi c2i = xvi xì /zi ‍. Using ‍F˛ x˛ = 0,‍ the term in brackets expands to:
zi c2i
N
P
Aw + Av i = di + Fvij x`j + Fìj xvj . (A-14)
‍ j=1 ‍
This arrangement is less prone to rounding errors than the original. It would be difficult to alter the original equations to achieve this
improvement. Zick discusses this and other equation arrangements to minimize rounding errors in the calculation of U (Zick 1992).
Molar Volume. The linearization of the volume constraint is also given in Eq. 46 and repeated here:

nh
‍ n = CTw w + Cv v + CTn + Cp p,‍ (A-15)
where:
46 2022 SPE Journal

c 1
NP
Cv = vi ì di
i=1

Cwi = v 1 v vi ì

Cni = 1 v xzìi ì + v xzvii vi

‍ Cp = 1 v p` + vpv . ‍
Eliminating the vapor fraction gives the following:

nh Cv eT Cv BTn
n = CTw + Bv w + CTn + Bv + Cp p
Q T Q Tn + Cp p.
‍ = Cw w + C ‍ (A-16)
Once Aw is factored, the partial molar volumes and total hydrocarbon fluid compressibility, Eq. 15, can be calculated:

‍ nh = T n nh Ch p + r ,‍ (A-17)
T 1
where after calculating Ri = CQ w AQ w :
‍ i‍
T = CQ T RT A
Qn
nT
Ch = R Ap Cp /h
‍ r = RT Ar . ‍
1 T
Using the fact that ‍A1
w = Aw ‍, the calculation of R is achieved with Eq. 50 as follows:
T 1 T
R = CQ w AQ w = CQ w I + UeT A1
w
T T
‍ = A1 Q Q
w Cw + Cw UAw e.
1
‍
T
‍Q w U.‍
It requires the two solves indicated and the dot productC
Linearization for AIM and FIM. For implicit calculations, the phase saturation and individual phase fractions also require linearization.
Once the molar volume has been linearized, the additional linearizations for implicit calculations are simple. The individual phase
fractions linearize as:

x˛i = zi xz˛ii + xz˛ii zi . (A-18)
‍ ‍
The linearizations for these terms are given in Eqs. A-4 and A-6. Substitution produces a dependence on v,
‍ w, and ‍. The saturation
of the vapor phase is:
nh vv
Sv = O
, (A-19)
‍ V ‍
so its linearization is given by:

1
Sv = nh v v + nh vv + vv nh Sv O .
O ‍ (A-20)
‍ V
The vapor phase molar volume, ‍ v ‍, can be linearized like total hydrocarbons molar volume,
‍ h‍, Eq. A-15. Treatment of the other terms
is straightforward.
Note, the assumption of constant υα in Eqs. 42 and 66 is easily remedied, because we could just as easily linearize:
h i
x˛i˛ = 1˛ zi xz˛ii + xz˛ii zi x˛i˛ ˛ . (A-21)
‍ ‍
Because the phase molar volume terms are already linearized for the volume constraint and saturation, Eqs. A-15 and A-2, the linear-
ization of this concentration would require little additional effort and could be represented like Eq. 44. xα does not appear alone in Eq. 1,
so its linearization is usually not needed.
In summary, the quantities needed for implicit calculations can be linearized using the same procedure as used for the total hydrocar-
bon molar volume. The dependence on ‍ v and w ‍can be eliminated the same way. The vapor fraction is eliminated as before using the
linearization of the Rachford-Rice equation, Eq. 7. The K-value terms, w, could be eliminated like Eq. A-17. However, when the elimina-
tions are to be performed for all individual phase composition variables, it is more efficient to first perform a complete elimination on Eq.
A-11, so the equation is equivalent to:
1
w + AQ w A Q 1
Q n + A Q 1
w Ap p = Aw Ar .‍ (A-22)
‍

Because ‍AQ w ‍has been factored for the partial volume calculation, these calculations require ‍O 2N3c ‍operations. The additional calcu-
lations
3 to produce Eq. 44 are then straightforward. Owing to the simplicity of Eqs. A-4, A-6, A-18, and A-2, no other calculations of
‍O Nc ‍are required.
2022 SPE Journal 47

Compositional Reservoir Simulation: A Review

Uploaded by

Copyright:

Available Formats

Compositional Reservoir Simulation: A Review

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Compositional Reservoir Simulation: A Review

Uploaded by

Copyright:

Available Formats

Compositional Reservoir Simulation:

*Corresponding author; email: larry@tildentechnologies.com

2022 SPE Journal 1

Nature of the Problem

2 2022 SPE Journal

and Np on the fractions within each phase:

2022 SPE Journal 3

4 2022 SPE Journal

Fig. 2—C1–C3–C12 partial volumes, Line A.

2022 SPE Journal 5

6 2022 SPE Journal

After discretization, these equations take the form:

2022 SPE Journal 7

8 2022 SPE Journal

2022 SPE Journal 9

Component balances, Eq. 29  ‍nw , n` , x`i , i = 1, ..., Nc  2‍

Table 1—Equation arrangement used by Fussell and Fussell

Phase equilibria, Eqs. 5 and 7  ‍x`i , i = 1, ..., Nc  1, xv1 , xv2‍

Table 2—Equation arrangement of Coats (1980) and Coats et al. (1998).

10 2022 SPE Journal

Phase equilibria, Eq. 5  ‍Ki , i = 1, ..., Nc  1‍

Table 3—Equation arrangement of Chien et al. (1985).

Phase equilibria, Eq. 5 ‍v, xvi , i = 1, ..., Nc  2‍

Table 4—Equation arrangement of Young and Stephenson (1983).

2022 SPE Journal 11

12 2022 SPE Journal

Phase equilibria, Eq. 5  ‍n`i , i = 1, ..., Nc  1‍

Table 5—Equation arrangement of Subramanian et al. (1987).

Phase equilibria, Eq. 5  

Table 6—Equation Arrangement of (Young 1987, 1991).

2022 SPE Journal 13

‍ ` and v must

14 2022 SPE Journal

Using the Schur complement of Aχ, ‍D = D  D A1

2022 SPE Journal 15

Table 7—Secondary variables.

16 2022 SPE Journal

And the overall molar volume linearization is reduced to:

2022 SPE Journal 17

18 2022 SPE Journal

Fig. 4—Fractional flow functions.

2022 SPE Journal 19

Work Ratios Factor

Table 8—Parameters of heuristic

Fig. 5—IMPES, FIM, and AIM tradeoffs black oil.

20 2022 SPE Journal

2022 SPE Journal 21

where Λj is the maximum eigenvalue of the matrix:

22 2022 SPE Journal

2022 SPE Journal 23

24 2022 SPE Journal

Table 9—Implementation improvements.

2022 SPE Journal 25

26 2022 SPE Journal

lim t1 = E1 Nproc , (70)

Fig. 7—Efficiency and serial fraction, Amdahl’s law.

2022 SPE Journal 27

Table 10—Serial content for various parallel simulators.

Component balances, Eq. 29 ‍nw , n` , x`i , i = 1, ..., Nc 2‍

Phase equilibria, Eqs. 5 and 7 ‍x`i , i = 1, ..., Nc 1, xv1 , xv2‍

Phase equilibria, Eq. 5 ‍Ki , i = 1, ..., Nc 1‍

Phase equilibria, Eq. 5 ‍v, xvi , i = 1, ..., Nc 2‍

Phase equilibria, Eq. 5 ‍n`i , i = 1, ..., Nc 1‍

Phase equilibria, Eq. 5

‍ ` and v must

Using the Schur complement of Aχ, ‍D = D D A1

PNp n h i h i h i o

‍ Aw w + Av v + An + Ap p = Ar ,‍ (A-9)