Accepted Manuscript

Operator-based linearization approach for modeling of multiphase multi-component flow in

porous media

Denis V. Voskov

PII: S0021-9991(17)30144-4
DOI: http://dx.doi.org/10.1016/j.jcp.2017.02.041
Reference: YJCPH 7173

To appear in: Journal of Computational Physics

Received date: 23 March 2016

Revised date: 15 February 2017
Accepted date: 16 February 2017

Please cite this article in press as: D.V. Voskov, Operator-based linearization approach for modeling of multiphase multi-component flow in
porous media, J. Comput. Phys. (2017), http://dx.doi.org/10.1016/j.jcp.2017.02.041

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing
this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is
published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all
legal disclaimers that apply to the journal pertain.
 Operator-based linearization approach for modeling of
multiphase multi-component flow in porous media

Denis V. Voskov
Delft University of Technology
Department of Geoscience & Engineering
Stevinweg 1, 2628 CN Delft, NL
d.v.dvoskov@tudelft.nl

Introduction

Reservoir simulation of complex physical processes requires on the solution

of the nonlinear model equations. These include partial differential equations
describing the conservation of mass, momentum, and energy as well as different
types of local constraints that define phase behavior and/or chemical reactions
in the system. A general purpose process of robustly solving these equations
typically includes highly implicit (fully or adaptively) approximation in time
to avoid numerical limitations (such as a CFL limit for a multi-component
transport) and complex flux approximations to preserve the accuracy in spatial
discretization. Both of these approximations introduce nonlinearity into the
system of equations.
In reservoir simulation, different formulations are used for solving the system
of governing equations. These formulations usually differ in types of nonlinear
unknowns used for the solution, see [1] for an extensive overview. For the
purpose of classification, all nonlinear formulations can be separated into two
classes: mass-based and phase-based formulations. The mass-based formulation
[2, 3] uses mass-related variables (overall composition or molar mass) for solving
a nonlinear system of equations. In the phase-based (also called natural) for-
mulation [4], variables related to a phase (saturation and phase fractions) serve
as nonlinear unknowns.

Preprint submitted to the Journal of Computational Physics February 20, 2017

 Recently, several advanced nonlinear formulations were developed for com-
positional simulations [5, 6]. All of these approaches propose a better way of
dealing with phase changes. However, these formulations were never tested for
complex realistic problems. The idea described in [5] was adapted for a general
purpose simulation and tested against state of the art approaches [7]. For the
problems of practical interest, this formulation demonstrated only insignificant
improvements with respect to conventional methods.
After that equations are discretized in space and time, they require a lin-
earization step. Usually, a Newton-based method is used for the linerization
implying an assembly of a Jacobian and a residual for the combined system
of equations. In general purpose reservoir simulation this is a challenging task.
Both value and derivatives of different properties involved in the governing equa-
tions need to be computed and stored. This requires fixing nonlinear unknowns
and the formulation for the implementation of a simulator code [8, 9]. It is
still possible to derive a Jacobian for another formulation on a linear level using
a transformation matrix (see [8] for example), but this approach often lacks
robustness.
The development of an Automatic Differentiation (AD) technique helps to
improve the situation. In several research areas, the utilization of AD makes
the solution of the system of nonlinear equations much more flexible. In reser-
voir simulation, the first attempt to design an AD-based general-purpose sim-
ulator was accomplished based on the Automatic Differentiation Expression
Template Library (ADETL) developed by Younis [10]. An Automatic Differen-
tiation General Purpose Research Simulator (ADGPRS) was developed using
AD-capabilities and the special data-structures of ADETL [11]. This simulator
utilizes the idea of having different nonlinear formulations in a single framework
providing a consistent platform for the comparison of formulations [7, 12]. ADG-
PRS creates opportunities for different research directions, including advanced
discretization [13], extended physics [14, 15] and inverse capabilities [16].
Using the ADGPRS framework, the new nonlinear strategy called Compo-
sitional Space Parametrization (CSP) was developed for compositional simula-

2
tion [12, 17]. This approach is based on an analysis of compositional displace-
ment in the hyperbolic limit, which indicates that the solution is localized only
around the displacement path in tie-line space [18]. In this approach, tie-lines
are adaptively parameterized along the displacement solution and the composi-
tional problem is solved using tie-line parameters as nonlinear unknowns. The
new tie-line-based approach improves the nonlinear convergence of composi-
tional simulations and the general performance due to a significant reduction in
phase equilibrium computations. Recently, the parametrization approach based
on sparse grids was suggested for improving the performance of phase behav-
ior [19]. The authors reports significant improvements in CPU time, but their
conventional flash algorithm demonstrates a unrealistically high reference time.
While enough attention has been paid to the linearization technique and the
solution of the phase behavior problem, the application of the nonlinear solver
has also been the object of considerable interest. In the course of the simulation,
nonlinear unknowns are updated based on the solution of a linear system. It
is a common practice in reservoir simulation to modify the update. Different
versions of global or local chopping procedures exist for the nonlinear update.
One of the most commonly used in practice is the Appleyard chop [9], where
the saturation updates are corrected locally according to the end-points of a
relative-permeability function. Another approach is to control unknowns from
crossing an inflection point of the fractional flow function [20], which becomes
more complex when compositional effects [21], gravity [22], or capillarity with
upstream changes [23] are present. All of these methods are trying to resolve the
most severe nonlinearities, which are introduced by a combination of properties
in the governing equations after the conventional linearization is applied.
In reservoir simulation, we always deal with approximations of discretiza-
tions in space and time while resolving all the features of the physical descrip-
tion very accurately. In a lot of cases, this forces the nonlinear solver to struggle
and requires advanced methods. In this paper, we present a completely different
approach to the linearization of the governing equations that follows the ideas
developed by Zaydullin et al. [12]. The discretized version of the conserva-

3
tion equations is written in an operator form where each term is presented as
a product of two operators. The first type of operators depends on the phys-
ical properties of the rock and the fluid while the second type depends on the
properties altered in space. The first type of operators is parametrized in the
physical space of a simulation problem in pre-processing stage. In the course
of the simulation, multilinear interpolation is applied for physics-based opera-
tors while the second type of operators is evolved conventionally. It is shown
that using a limited number of points in the proposed physical representation,
the main nonlinearity of a simulation model can be captured quite precisely.
The numerical results of an approximation in the linearized space demonstrate
better nonlinear convergence with the error in the approximation controlled by
the accuracy of the interpolation. In addition, any expensive property eval-
uations required by the physical model are limited to a few parametrization
points, which significantly improves the performance of the simulation. Several
additional advantages of the proposed Operator-Based Linearization (OBL) ap-
proach are discussed in the last section.

Modeling approach

In this section, we describe the governing equations and nonlinear formula-

tion for a reservoir simulation problem.

Conservation equations

The transport equations for an isothermal system containing nc components

and np phases can be written as:
⎛ ⎞
np np
∂ ⎝  
φ xcj ρj sj ⎠ + div xcj ρj vj (1)
∂t j=1 j=1
np

+ xcj ρj q̃j = 0, c = 1, . . . , nc .
j=1

Here, we typically define the coefficients of the equations as functions of spatial

coordinate ξ and physical state ω:

4
 • φ(ξ, ω) – porosity,

• xcj (ω) – the mole fraction of component c in phase j,

• sj (ω) – phase saturations,

• ρj (ω) – phase molar density,

• vj (ξ, ω) – phase velocity,

• q̃j (ξ, ω) – phase rate per unit volume.

To describe the flow of each phase, we use Darcy’s law:

 k 
rj
vj = − K (∇pj − γj ∇d) , j = 1, . . . , np , (2)
μj

where

• K(ξ) – permeability tensor,

• krj (ω) – relative permeability,

• μj (ω) – phase viscosity,

• pj (ω) – vector of pressures in phase j,

• γj (ω) – gravity term,

• d(ξ) – vector of depths (positive downwards).

By applying a finite-volume discretization on a general unstructured mesh

and backward Euler approximation in time, we can transform the conservation
equations into
⎛ ⎞
 
V ⎝(φ xcj ρj sj )n+1 − (φ xcj ρj sj )n ⎠ (3)
j j
⎛ ⎞
  
−Δt ⎝ xlcj ρlj Tjl Δψ l ⎠ + Δt ρp xcj qj = 0.
l∈L j j

where V is the volume of a control volume and qj = q̃j V the source of a phase.
Here we neglected capillarity, gravity and used a Two-Point Flux Approximation

5
(TPFA) with upstream weighting introducing the summation over all interfaces
L connecting the control volume with another grid blocks. Based on these
simplifications, Δψ l becomes a simple difference in pressures between blocks a
and b, where Tjl is a phase transmissibility. These assumptions are not required
by the method, but help to simplify the further description.

Nonlinear unknowns

The system of equations (3) is the discretized form of flow and transport
equations for general multi-component fluid. Here, we used a Fully Implicit
Method (FIM) time approximation. It requires the (xlcj ρlj Tjl Δψ l ) flux term to
be defined based on nonlinear unknowns at a new timestep (n+1) and introduces
nonlinearity to the system of governing equations. Another source of nonlinear-
ities comes from the additional assumption on instantaneous thermodynamic
equilibrium, which is required to close the system.
Several different strategies exists for the nonlinear solution of the resulting
system, see [1], and [11] for an extensive description and examples. Here, we used
the overall molar formulation suggested by Collins et al. [3]. In this formulation,
thermodynamic equilibrium is assumed at every nonlinear iteration of solution
of Eq. 3. For the control volume at multiphase conditions with np -phases, the
following system of equations needs to be solved:
np

Fc = zc − νj xcj = 0, (4)
j=1
Fc+nc = fc1 (p, T, x1 ) − fcj (p, T, xj ) = 0, (5)
nc

Fj+nc ×np = (xc1 − xcj ) = 0, (6)
c=1
np

Fnp +nc ×np = νj − 1 = 0. (7)
j=1

Here zc = j xcj ρj sj / j ρj sj is overall composition and fcj (p, T, xcj ) is the
fugacity of component c in phase j. This procedure is called a multiphase flash
[24]. For a given overall composition zc , the solution of (4)-(7) provides molar
fractions for each component xcj and phase fractions νj .

6
 In the overall molar formulation, the unknowns are p and zc . They fully
define the physical state ω for a given control volume. Once multiphase flash is
solved, it can provide derivatives of all properties in (3) with respect to nonlinear
unknowns using the inverse theorem, see [11] for details.

General form of governing equations

We can re-write Eq. 3 as the component of a residual vector in an algebraic

form

rc (ξ, ω, u) = a(ξ) (αc (ω) − αc (ωn )) (8)


− βcl (ω)bl (ξ, ω) + θc (ξ, ω, u) = 0.
l

Here, we defined

αc (ω) = (1 + cr (p − pref )) xcj ρj sj , (9)
j
a(ξ) = V (ξ)φ0 (ξ), (10)
 krj
βc (ω) = xcj ρj , (11)
p
μj
b(ξ, ω) = ΔtT ab (ξ)(pb − pa ), (12)

θc (ξ, ω, u) = Δt ρj xcj qj (ξ, ω, u). (13)
j

In addition, cr is the rock compressibility, T ab is the geometric part of transmis-

sibility (which involves permeability and the geometry of the control volume),
ω and ωn are nonlinear unknowns at the current and the previous timestep
respectively, and u is a vector of well control variables.
The operator αc is dependent on the properties of rock and fluid, and in-
dependent of spatially distributed properties (initial porosity) as in the case of
the operator a. Similarly, the divergence operator is present as a fluid-related
operator βc independent of spatial distributed properties (permeability) and the
discretization-related operator b. The same approach can be applied for the well
source/sink operator θc , but for simplicity we did not apply it here.

7
Linearization methods

In this section, we describe different types of linearization using the general

algebraic form of governing equations (8).

The standard linearization approach

To solve nonlinear system (8), we need to linearize it. The conventional

approach in reservoir simulation is based on the application of the Newton-
Raphson method. In each iteration of this method, we need to solve a linear
system of equations of the following form

J(ω k )(ω k+1 − ω k ) = −r(ω k ), (14)

where J is the Jacobian defined at nonlinear iteration step k.

The typical approach requires a sequential assembly of the residual and the
Jacobian based on numerical approximations of the analytic relations in (9)-(13).
This may demand an interpolation in tables (for standard PVT correlations or
relative permeabilities), or a solution of the highly nonlinear equations (for
EoS-based properties). Each property evaluation requires a storage space for
both values of the property and its derivatives with respect to the nonlinear
unknowns.
Most reservoir simulation software performs numerical- [25]or hand-differentiation
[9] of each property with respect to nonlinear unknowns. In this work, we uti-
lized the ADGPRS framework [11] for the implementation of conventional and
newly proposed linearization procedures.

Operator-based linearization

Here, we propose a new strategy for the linearization of the reservoir simu-
lation problem described by Eq. 8. As can be seen from the structure of each
operator in (9)-(12), this system is based on a complex combination of different
nonlinear properties and relations. Since we fixed our space and time approx-
imation, the discretization error can be controlled only by the variation of the

8
timestep size Δt and the characteristic size of the mesh embedded in the T ab
term. Both of these errors are controlled by the operators ψ and θc .
The operators αc and βc represent the physics-based terms. The accuracy of
the nonlinear physics representation is controlled by these two operators (and a
part of θc ). In conventional simulation, we introduce all the nonlinear properties
into the conservation equation as is. Next, the nonlinear solver tries to resolve all
the details of the nonlinear description, struggling sometimes with unimportant
features due to the numerical nature of the property representations.
The Operator-Based Linearization (OBL) strategy, proposed in this work,
is based on the simplified representation of the nonlinear operators αc and βc
in the parameter-space of a simulation problem. For simplicity, let us assume
in the following derivations and implementation that the system of equations
(8) describes flow and transport in a binary system. In this case, we only have
two independent variables in the overall molar formulation - pressure p and
one independent overall composition z (the second is dependent based on the
constraint z1 + z2 = 1). For an isothermal reservoir simulation, the range in
the parameter space for pressure is fully defined by the injection/production
conditions on wells. The range of z is obviously restrained by [0, 1].
In the further numerical examples, we uniformly discretize the parameter
space with a fixed number of points N . The interpolation intervals are defined as
[P1 , P2 , . . . , PN ] and [Z1 , Z2 , . . . , ZN ]. Next for p ∈ [Pi , Pi+1 ] and z ∈ [Zj , Zj+1 ]
we define

p − Pi z − Zj
pi = , zj = , (15)
Pi+1 − Pi Zj+1 − Zj

and the auxiliary relation fi,j = f (Pi , Zj ). Based on pi , zj and fi,j , the inter-
polant of function Ff can be defined as

Ff (p, z) = (1 − zj )[(1 − pi )fi,j + pi fi+1,j ] + zj [(1 − pi )fi,j+1 + pi fi+1,j+1 ], (16)

9
with corresponding gradients

∂Ff (1 − zj )(fi+1,j − fi,j ) + zj (fi+1,j+1 − fi,j+1 )

= , (17)
∂p Pi+1 − Pi
∂Ff (1 − pi )(fi,j+1 − fi,j ) + pi (fi+1,j+1 − fi+1,j )
= . (18)
∂z Zj+1 − Zj

Similarly to [17], the error between an interpolant and function can be evaluated
based on the following relation

2

Vω2
∂ f

|f − Ff | ≤ sup

, (19)
4 ω ∂ω

where Vω = ||Δωm ||.

To evaluate α̂c (p, z) and β̂c (p, z) in the course of the simulation, we apply
an interpolation

α̂c (p, z) = Fαc (p, z), β̂c (p, z) = Fβc (p, z). (20)

This representation helps to provide a continuous description of the physics-

based operators in the proposed approach. In the examples below, the dis-
cretization of the entire parameter space of the problem is applied at the pre-
processing stage. In Fig. 1, an example of all operators parameterized at N = 64
is shown for a binary compositional system described below. For a general pur-
pose simulation, it is possible to apply this approach adaptively, following the
idea of a compositional space parametrization [18, 12]. The example of adaptive
interpolation applied to the Operator-Based Linearization is described by Khait
and Voskov [26].
In the proposed approach, the number of points in the interpolation con-
trols the accuracy of the approximation of the nonlinear physics, similar to the
accuracy of the approximation in space and time being controlled by the grid
size. The error described by (19) is similar to the truncation error in the dis-
cretization of equation (8). The nonlinear solver deals here with a simplified
representation directly expressed as a piece-wise linear combination of nonlinear
unknowns. Also, the relation (17)-(18) provides direct derivatives with respect
to nonlinear unknowns, which significantly simplifies the evaluation and assem-

10
 Figure 1: All operators for a binary compositional system parameterized at N = 64.

bly of Jacobians. In the next section, we provide several numerical examples to

illustrate the proposed approach.

Solution method

Once each operator in (8) is linearized, the residual vector r and the Jacobian
J can be assembled. The overall-molar formulation does not require a secondary
set of equations [3, 11] and we can apply the linear solver directly to system (14).
In this work, we employed GMRES with the two-stage Constrained Pressure
Residual (CPR) preconditioner [27] as a linear solver. For details on the linear
solution of general-purpose simulation problems, see [8]. Once the linear solution
is found, we need to update the nonlinear unknowns. Here we applied a standard
Newton-Raphson update with the maximum allowed local change in the overall
composition Δz = 0.1.
An important part of the nonlinear solver is a timestep control. This control
may include different types of heuristics to connect timestep size changes in

11
different variables from the previous timestep [28]. Here, we employ a simple
strategy:

• We start with a pre-defined minimal timestep Δtmin ;

• if the nonlinear solver converges in a given number of iterations Ni , for

the current timestep, we multiply the next timestep by a fixed ratio γ;

• if the nonlinear solver fails to converge, we divide the next timestep by

the same constant γ;

• if the maximum step Δtmax is reached, we keep it for further simulation;

• finally, if the timestep is cut to the minimal acceptable value (10−8 days),
the simulation is stopped.

In the examples below, we use a fixed set of parameters, where Δtmin = 10−3
days, Ni = 30, γ = 2 and Δtmax = 10 days. The convergence of the nonlinear
solver is based on the following criterion: maxi |Ri /(αai )| < 10−5 , where i is a
grid block number.

Numerical results

In this section, we compare the numerical results obtained by employing

OBL approach against the conventional linearization, used as a reference solu-
tion. The comparison covers different physical processes, modeled in 1D and
3D reservoirs.

One-dimensional homogeneous reservoir

Here we introduce two types of physical descriptions corresponding to dif-

ferent oil recovery processes: waterflooding (immiscible displacement) and gas
injection (miscible displacement). Physical properties for these simulations are
described in Appendix A.
A homogeneous 1D reservoir of 1000 m length with one injection and one
production well in the first and last grid blocks was modeled first. Constant
porosity φ0 = 0.2 and permeability k = 100 mD are used in this model. We

12
apply a finite-volume discretization on a standard Cartesian grid with block size
Δx = 1, Δy = 10 and Δz = 10 m. For the well discretization, the Peaceman
formula [29] is used. The injection well is controlled by a water rate qw = 20
m3 /day, and the production well is controlled by a bottom hole pressure pw =
100 bars.
We run a set of simulations and compare the reference solution, based on
a conventional linearization method, with results performed at different resolu-
tions of interpolation tables in the OBL approach. For simplicity, we use the
equal number of points for each unknown (p and z) with uniformly distributed
values in the range between pmin and pmax for pressure, and between 0 and 1
for the composition.
In the immiscible displacement problem, pmin = 50 and pmax = 650 were
used to cover the parameter space completely. The corresponding average maxi-
mum CFL number for this problem is equal to 13.0, which can be interpreted as
the ratio between the timestep size of the performed (Fully Implicit) simulation
and the timestep restricted by an explicit transport approximation (e.g., in the
Implicit Pressure Explicit Saturation (IMPES) method). The total number of
timesteps in this simulation is Nt = 63.
The results of the first simulation are displayed in Table 1. The number of
points used by each unknown for representation of operators αc and βc is in the
first column. The second column indicates the number of nonlinear iterations
used in the solution. The third and fourth columns correspond to errors in
pressure and composition solution, compared to the reference simulation, that
are introduced by the method. The last column shows the CPU time for each
simulation.
In Fig. 2, we present a spatial distribution of the water overall molar fraction
and pressure by the end of the simulation at time t = 500 days. In this figure
you can see the difference between the solution with conventional linearization of
nonlinear physics vs. the parametrized solutions with four different resolutions
of parameter space. The finest resolution shown in the figure contained only
16 points for each unknown. However, it manages to reproduce standard sim-

13
 Table 1: Results of 1D immiscible displacements simulation
Resolution Iters. Ep , % Ez , % CPU, sec
Standard 528 - - 2.0
64 × 64 521 0.02 0.01 2.0
32 × 32 468 0.06 0.04 1.8
16 × 16 422 0.27 0.19 1.6
8×8 365 1.04 0.75 1.4
4×4 316 6.40 3.66 1.3
2×2 146 51.23 22.68 0.8

ulation results quite accurately with a smaller number of nonlinear iterations.

The simulation time is not significantly different for this simulation since it is
relatively short. However, you can see the reduction in simulation time almost
proportional to the reduction in the number of nonlinear iterations.
Next, we describe the simulation results for the miscible displacement prob-
lem. Unlike in the first case, where all the fluid properties were computed based
on a simplified correlation (see Appendix A), in the compositional case, fluid-
based properties were computed based on a solution of a highly-nonlinear cubic
Equation of State (EoS). Notice that between the injection gas phase and the
initial liquid phase conditions, the solution passes through the two-phase region.
Here, the derivatives of the phase fractions xcj with respect to nonlinear un-
knowns are obtained from (4)-(7), which makes the combination of properties
even more nonlinear.
The results of the miscible displacement simulation are displayed in Table 2.
Similar to the immiscible displacement case, we used the same resolutions for
intervals in pressure and composition. However, in the case of gas injection,
the interval of pressure changes is lower and pmin = 50 and pmax = 150 bars.
Well controls were set to a gas rate q = 2000 m3 /day in the injection well and a
pressure pw = 70 bars in the production well. For this case, the corresponding
average CFLmax = 11.5 and Nt = 60 timesteps performed in the simulation.
Table 2 displays the same tendency as in the immiscible displacement case.

14
Figure 2: Comparison between overall molar fraction (upper) and pressure (lower) solutions
for different resolutions of parametrization in the immiscible displacement test case

With the increasing resolution, the error between the parametrized and the con-
ventional solutions decreases while the number of nonlinear iterations increases.
Notice that in this case, the conventional simulation requires significantly larger
time than any simulation based on OBL. It can be explained by a more expen-
sive linearization step in the conventional approach which requires multiphase
flash and solution of EoS to the evaluate properties and their derivatives in each
control volume at every nonlinear iteration. In contrast, the OBL approach re-
quires these calculations only in a limited number of interpolation points at
the pre-processing stage. This explains a significant difference in performance
between conventional and OBL approaches.
Displacement profiles for both overall molar fraction of the first component
(C1 ) and pressure are shown in Fig. 3 at time t = 500 days. Again, the finest
resolution shown in the figure with only 16 points is capable of reproducing the
conventional simulation results with the reduced number of nonlinear iterations.
In Fig. 4, the most nonlinear operator for this system β2 is shown for different

15
 Table 2: Results of 1D miscible displacement simulation
Resolution Iters. Ep , % Ez , % CPU, sec
Standard 448 - - 10.8
64 × 64 386 0.03 0.01 2.7
32 × 32 364 0.12 0.04 1.9
16 × 16 343 0.65 0.16 1.6
8×8 337 2.36 0.82 1.4
4×4 335 3.74 3.34 1.2
2×2 145 62.24 16.39 0.8

resolutions of the parameterized space N for two limiting pressures p = 60 and

p = 160. A full shape of the operator β2 can be seen in Fig. 1. It is clear that
at a low resolution of parameterized space (N = 2 and N = 4), the nonlinearity
of the operator is not resolved, which explains a large errors Ez in Table 2.
However, starting at N = 8, the nonlinearity of the operator reproduced quite
accurately, which drops the error Ez < 1%.

Three-dimensional heterogeneous reservoir

Here we present the simulation results for a three-dimensional test based on

the upper 5 layers of the SPE10 reservoir [30]. Buoyancy is neglected, and the
vertical displacement is driven by convective forces only. The injector well was
placed in the middle with four production wells in each corner of the model. All
wells were perforated into all layers and operated at the same controls as in the
previous model. Both porosity and permeability are highly heterogeneous in
this case (15 orders of magnitude difference) which explains the very scattered
results in Fig 5. The number of timesteps in this simulation Nt = 113. Again,
we compare the conventional solution with a solution based on an interpolation
of operators. It can be seen that the error introduced by all resolutions is
mostly located at the displacement interfaces and still provide a reasonable
approximation of the reference solution. With a finer resolutions, the difference
in the results is almost indistinguishable as can be seen in both the error map
of the solution with 8 × 8 resolution and cross-section distributions.

16
Figure 3: Comparison between overall molar fraction (upper) and pressure (lower) solutions
for different resolutions of parametrization in the miscible displacement test case

The performance of the full resolution study is shown in Table 3. The

improvement in the number of nonlinear iterations in a highly heterogeneous 3D
reservoir model looks similar to 1D homogeneous case. Notice that this type of
displacement is characterized by a larger average CFLmax = 223.1, which often
makes the nonlinear convergence complicated when the physical parameters are
fully refined. The improvement in CPU time is also similar and proportional to
the number of nonlinear iterations.
Similar observations can be carried out for a 3D miscible displacement case.
In this simulation, the average CFLmax = 167.6 and the number of timesteps
Nt = 214. The comparison between the reference solution and the solution,
based on interpolated operators, is shown in Fig. 6. Again, the solution based
on the coarse interpolation tables still provides a reasonable approximation for
the reference solution with an error distributed along the displacement fronts.
The full convergence study for this case is shown in Table 4. The improvement in

17
Figure 4: Operator β2 for miscible displacement case at end-point pressures and different
resolutions in composition

Table 3: Results of 3D immiscible displacement simulation

Resolution Iters. Ep , % Ez , % CPU, sec
Standard 583 - - 651.6
64 × 64 577 0.03 0.01 629.5
32 × 32 562 0.10 0.02 604.2
16 × 16 554 0.41 0.10 593.4
8×8 529 2.54 0.47 563.4
4×4 473 6.49 2.05 528.0
2×2 320 101.35 13.13 362.6

simulation time for the OBL approach is still quite significant in comparison to
the conventional approach due to the reduction in expensive EoS computations
and more efficient assembly of the Jacobian.
The approach proposed here translates all the nonlinear relations into lin-
earized versions. This approximation is directly based on nonlinear unknowns,
which significantly improves the behavior of the nonlinear solver. One can bal-
ance between the more accurate representation of physics and the performance
of the nonlinear solver in the simulation. A significant improvement in simula-
tion time can be achieved due to the reduction in expensive property evaluations
in the course of simulation.

18
Figure 5: Overall molar fractions: the reference solution based on conventional approach,
the error between reference solution and OBL solution with 8x8 interpolation table and cross-
section solution at several OBL resolutions for immiscible displacement simulations in reservoir
contains 5 upper layers of the SPE10 problem

Consistency of numerical solution

In this section, we demonstrate the consistency of the proposed linearization

method assuming that the original problem described by (1) has a numerical
solution. To simplify analysis, we assume that the model is limited by a 1D
reservoir with Cauchy boundary conditions on left and right side. This simplifies
the spacial discretization, which yields to the following equation in vector form
for the block i:

ri (ω) = (αi (ω) − αi (ωn )) (21)

− γ (βi (ω)(pi+1 − pi ) − βi−1 (ω)(pi − pi−1 )) = 0,

19
Figure 6: Overall molar fractions: the reference solution based on conventional approach, the
error between reference solution and OBL solution with 8x8 interpolation table and cross-
section solution at several OBL resolutions for miscible displacement simulations in reservoir
contains 5 upper layers of the SPE10 problem

where γ = ΔtT /(φ0 V ).

The internal Jacobian row of the equation can be written as:
⎡ ⎤T
γBi−1 (pi − pi−1 ) − γβi−1 × ep
⎢ ⎥
J =⎢ ⎥
⎣Ai − γBi (pi+1 − pi ) + γ(βi + βi−1 ) × ep ⎦ (22)
−γβi × ep

where
   
∂αi ∂βi
A= , B= , i, j = 1, . . . , nc (23)
∂ωj ∂ωj

and ep = [1 0 . . . 0]T is a nc -vector. We can present the Jacobian as a sum of

20
 Table 4: Results of the 3D miscible displacement simulation
Resolution Iters. Ep , % Ez , % CPU, sec
Standard 1025 - - 2074.4
64 × 64 951 0.06 0.01 1101.7
32 × 32 922 0.30 0.03 1079.2
16 × 16 869 1.40 0.11 1015.1
8×8 807 4.28 0.47 950.1
4×4 783 22.81 3.15 966.3
2×2 488 70.49 6.77 612.7

two matrix
⎡ ⎤T ⎡ ⎤T
−γβi−1 × ep γBi−1 (pi − pi−1 )
⎢ ⎥ ⎢ ⎥
J = J p + Jz = ⎢ ⎥ ⎢ ⎥
⎣γ(βi + βi−1 ) × ep ⎦ + ⎣Ai − γBi (pi+1 − pi )⎦ (24)
−γβi × ep 0

where Jp corresponds to the mostly elliptic part of the problem and Jz corre-
sponds to the mostly hyperbolic part. For the consistent solution of the original
problem, the matrix Jp expressed in the conventional form of βi should be
diagonally dominated. It is clear that the diagonal dominance of Jp will not
be affected by the errors in interpolation of βi since these errors will cancel
each other in diagonal and off-diagonal part. This conclusion is also valid in a
multi-dimensional case.
It is more complicated to perform the analysis of matrix Jz in general case.
In the incompressible limit, the error bound defined by (19) will guarantee that
the approximation schema is contractive for interpolation operators when the
original scheme is contractive which is sufficient for the stability. The eigenval-
ues of sub-matrices Âi and B̂i related to derivatives of composition will define
the hyperbolicity of solution. Notice, that transport in binary systems is always
hyperbolic since the only eigenvalue is always real. For systems with a larger
number of components, the simulation for a coarse resolution of parameter space
with uniform parametrization may fail to converge, as was observed by Khait
and Voskov in [26]. This can be explained by the potential loss of hyperbol-

21
icity in systems with a large number of components. As the solution for this
problem, the parametriation can be adjusted to follow tie-lines which describe
the phase equilibrium in a system. As demonstrated in [12], the convexity in
tie-line parametrization provides a consistent numerical solution for a gas in-
jection problems with multiple components (up to nc = 8) even in the case of
an extremely coarse parametrization and a significant compressibility. We leave
these aspects for future investigation.

Conclusion and Discussion

We presented a new approach for the assembly of the residual and the Ja-
cobian in reservoir simulation. In this approach, the governing equations are
transformed into an operator form where each term is presented as a product of
two operators. The first type of operators is fully defined by the physical state
of the problem and usually corresponds to the most significant nonlinearities.
The second operator depends on both spatial and state variables. Next, we pa-
rameterize the first type of operators using a uniformly distributed mesh in the
parameter space while treating the second type of operators in a conventional
manner. By example of two different binary systems, we demonstrated that a
simple linear interpolation in the space of state variables (p and z) can repro-
duce the results of the reference system quite accurately even in the case of a
coarse interpolation table. The resolution study in physical space for simplified
1D homogeneous and realistic 3D highly heterogeneous models demonstrated a
good convergence for the two types of nonlinear physics related to immiscible
and miscible displacement processes. The performance of the nonlinear solver
demonstrated an improvement due to operation in a linearized physical space.
The proposed approach has several benefits. The first of them is the sim-
plicity of the application. The fully implicit approximation is a very attrac-
tive method due to its unconditional stability. However, the complexity of the
evaluation and the storage of all the nonlinear properties and their derivatives
makes development and further extension of fully implicit code quite challeng-
ing. The proposed approach simplifies the residual and the Jacobian assembly

22
significantly and provides a unique opportunity for the extension of the phys-
ical model. The performance of the Jacobian and the residual assembly can
be improved and adjusted to the modern computational architectures by using
vectorization of the interpolation operator.
Another benefit of the proposed method is related to nonlinear solvers. Most
of them in reservoir simulation are based on a version of Newton’s method. The
interpolation of residual operators as an explicit function of nonlinear unknowns
provides an opportunity for better nonlinear strategies. Advanced nonlinear
solvers for the proposed formulation can be designed based on a trust-region
method in interpolation tables, similar to the approach proposed in [21] for
conventional molar formulation. Following this nonlinear strategy, the potential
problems introduced by rapid changes of gradients in the coarse tables can be
resolved as well. The fact that the physical kernel of a problem has a multi-linear
representation based on nonlinear unknowns may also help to achieve a better
preconditioning on a linear level similar to Algebraic Multi-Scale methods [31].
The general-purpose application of the introduced method requires dealing
with multi-linear interpolation in a highly dimensional space, especially for com-
positional problems. The adaptive parametrization approach, developed in [12]
for representation of compositional properties, can help in the extension of the
method for problems with a large number of components. The accuracy of the
approach can be controlled by the error estimator and adaptive resolution in pa-
rameter space. An example of a simple error estimator applied for a limited case
of property interpolation can be found in [17]. The extension to systems with
more than two-phases can be performed based on a tie-simplex parametrization
approach [32, 33].
In general, reservoir simulation is always based on the approximation of spa-
tial and temporal discretization. In a lot of cases, this approximation has only a
first order of accuracy. At the same time, the description of complex nonlinear
physics, embedded into a model, is always very precise. The current realization
of the approach, described in this paper, can be seen as a compromise between
the accuracy of the nonlinear representation of physics and the performance

23
of the nonlinear solver. Finding the balance among the different coarsening
strategies is the target for the future research.

Acknowledgement

We would like to acknowledge Rustem Zaydullin for the research ideas that
inspired this work and Mark Khait for the code optimization. We also thank
the Stanford University Petroleum Research Institute for Reservoir Simulation
(SUPRI-B) program for the permission to use ADGPRS in this research.

References

[1] K. Aziz, T. Wong, Considerations in the development of multipurpose reser-

voir simulation models, in: The 1st and 2nd International Forum on Reser-
voir Simulation, Alpbach, Austria, 1989, pp. 77–208.

[2] G. Acs, S. Doleschall, E. Farkas, General purpose compositional model.,

Soc. Petrol. Eng. J. 25 (4) (1985) 543–553.

[3] D. Collins, L. Nghiem, Y.-K. Li, J. Grabenstetter, Efficient approach to

adaptive-implicit compositional simulation with an equation of state, Soc.
Petrol. Eng. J. 7 (2) (1992) 259–264.

[4] K. Coats, An equation of state compositional model, Soc. Petrol. Eng. J.

20 (1980) 363–376.

[5] A. Abadpour, M. Panfilov, Asymptotic decomposed model of two-phase

compositional flow in porous media: Analytical front tracking method
for riemann problem, Transport Porous Med. 82 (3) (2009) 547–565.
doi:10.1007/s11242-009-9428-8.

[6] A. Lauser, C. Hager, R. Helmig, B. Wohlmuth, A new ap-

proach for phase transitions in miscible multi-phase flow in
porous media, Adv. in Water Resour. 34 (8) (2011) 957 – 966.
doi:http://dx.doi.org/10.1016/j.advwatres.2011.04.021.

24
 [7] D. Voskov, An extended natural variable formulation for compositional sim-
ulation based on tie-line parameterization, Transport Porous Med. 92 (3)
(2012) 541–557. doi:10.1007/s11242-011-9919-2.

[8] H. Cao, Development of Techniques for General Purpose Simulators, PhD

Thesis, Stanford University, 2002.

[9] Geoquest, Eclipse Technical Description 2005A, Schlumberger, 2005.

[10] R. Younis, Modern advances in software and solution algorithms for reser-
voir simulation, PhD Thesis, Stanford University, 2011.

[11] D. V. Voskov, H. A. Tchelepi, Comparison of nonlinear formulations for

two-phase multi-component eos based simulation, J. Petrol. Sci. Eng. 82-
83 (2012) 101–111.

[12] R. Zaydullin, D. Voskov, H. Tchelepi, Nonlinear formulation based on

an equation-of-state free method for compositional flow simulation, Soc.
Petrol. Eng. J. 18 (2) (2013) 264–273.

[13] Y. Zhou, H. Tchelepi, B. Mallison, Automatic differentiation framework for

compositional simulation on unstructured grids with multi-point discretiza-
tion schemes, in: SPE Reservoir Simulation Symposium 2011, Vol. 1, 2011,
pp. 607–624.

[14] R. Zaydullin, D. Voskov, S. James, A. Lucia, Fully compositional and

thermal reservoir simulation, Comput. Chem. Eng. 63 (2014) 51 – 65.
doi:10.1016/j.compchemeng.2013.12.008.

[15] T. Garipov, M. Karimi-Fard, H. Tchelepi, Discrete fracture model for cou-

pled flow and geomechanics, Computational Geosciences 20 (1) (2016) 149–
160. doi:10.1007/s10596-015-9554-z.

[16] O. Volkov, D. Voskov, Effect of time stepping strategy on adjoint-

based production optimization, Computat. Geosci. 20 (3) (2016) 707–722.
doi:10.1007/s10596-015-9528-1.

25
[17] R. Zaydullin, D. Voskov, H. Tchelepi, Formulation and solution of compo-
sitional displacements in tie-simplex space, in: SPE Reservoir Simulation
Symposium 2013, Vol. 2, 2013, pp. 1268–1282. doi:doi:10.2118/163668-MS.

[18] D. Voskov, H. Tchelepi, Compositional space parametrization for miscible

displacement simulation, TIPM 75 (1) (2008) 111–128. doi:10.1007/s11242-
008-9212-1.

[19] Y. Wu, C. Kowitz, S. Sun, A. Salama, Speeding up the flash

calculations in two-phase compositional flow simulations – the ap-
plication of sparse grids, J. Comput. Phys. 285 (2015) 88 – 99.
doi:http://dx.doi.org/10.1016/j.jcp.2015.01.012.

[20] P. Jenny, H. A. Tchelepi, S. H. Lee, Unconditionally conver-

gent nonlinear solver for hyperbolic conservation laws with s-shaped
flux functions, J. Comput. Phys. 228 (20) (2009) 7497 – 7512.
doi:http://dx.doi.org/10.1016/j.jcp.2009.06.032.

[21] D. V. Voskov, H. A. Tchelepi, Compositional nonlinear solver based on

trust regions of the flux function along key tie-lines, in: SPE Reservoir
Simulation Symposium, SPE 141743, 2011. doi:doi:10.2118/141743-MS.

[22] X. Wang, H. Tchelepi, Trust-region based solver for nonlinear transport

in heterogeneous porous media, J. Comput. Phys. 253 (2013) 114–137.
doi:10.1016/j.jcp.2013.06.041.

[23] B. Li, H. Tchelepi, Nonlinear analysis of multiphase transport in porous

media in the presence of viscous, buoyancy, and capillary forces, J. Comput.
Phys. 297 (2015) 104–131. doi:10.1016/j.jcp.2015.04.057.

[24] M. Michelsen, The isothermal flash problem: Part ii. phase-split calcula-
tion, Fluid Phase Equilibr. 9 (1982) 21–40.

[25] K. Pruess, S. Finsterle, G. Moridis, C. Oldenburg, Y.-S. Wu, General-

purpose reservoir simulators: The tough2 family, Bulletin. Geothermal Re-
sources Council 26 (2) (1997) 53–57.

26
[26] M. Khait, D. Voskov, Operator-based linearization for non-isothermal mul-
tiphase compositional flow in porous media, in: ECMOR XIV-15th Euro-
pean Conference on the Mathematics of Oil Recovery, 2016, pp. 1–15.

[27] J. Wallis, R. Kendall, T. Little, Constrained residual acceleration of con-

jugate residual methods, Soc Pet Eng AIME J (1985) 415–428.

[28] K. Aziz, T. Settari, Petroleum Reservoir Simulation, Applied Science Pub-

lishers, 1979.

[29] D. Peaceman, Interpretation of well-block pressures in numerical reservoir

simulation., Soc Pet Eng AIME J 18 (3) (1978) 183–194.

[30] M. Christie, M. Blunt, Tenth spe comparative solution project: A com-

parison of upscaling techniques, SPE Reservoir Eval. & Eng. 4 (4) (2001)
308–316.

[31] Y. Wang, H. Hajibeygi, H. A. Tchelepi, Algebraic multiscale solver for flow

in heterogeneous porous media, J. Comput. Phys. 259 (2014) 284 – 303.
doi:10.1016/j.jcp.2013.11.024.

[32] A. Iranshahr, D. Voskov, H. Tchelepi, Gibbs energy analysis: Com-

positional tie-simplex space, Fluid Phase Equilibr. 321 (2012) 49–58.
doi:10.1016/j.fluid.2012.02.001.

[33] A. Iranshahr, D. Voskov, H. Tchelepi, Tie-simplex based compositional

space parameterization: Continuity and generalization to multiphase sys-
tems, AIChE J 59 (5) (2013) 1684–1701. doi:10.1002/aic.13919.

[34] D. Peng, D. Robinson, A new two-constant equation of state, Ind. Eng.

Chem. Fundam. 15 (1976) 59–64.

[35] J. Lohrenz, B. Bray, C. Clark, Calculating viscosities of reservoir fluids from

their compositions, Soc. Petrol. Eng. J. (1964) 1171–1176doi:10.2118/915-
PA.

27
Appendix A. Fluid and rock properties

The parameters defined in Table A1 are common for all the models. In this
table index p means water in the immiscible displacement model and gas in
the miscible displacement model. The Corey-type relative permeabilities are

Table A1: Rock-fluid parameters

Parameter Value Description
−5
cr 10 1/bar rock compressibility
Spr 0 phase residual saturation
Sor 0 oil residual saturation
np 2 phase exponents
no 3 oil exponents

defined as

Spe = (Sp − Spr )/(1 − Spr − Sor ), (A1)

krp = (Spe )np , kro = (1 − Spe )no . (A2)

The parameters of correlations used for properties in the immiscible dis-

placement example are defined in Table A2. Here we use standard correlations

Table A2: Immiscible displacement parameters

Param. Oil Water Description
ρ0p 850 kg/m3 1000 kg/m3 surface density
cp 10−4 1/bar 10−5 1/bar compressibility
μ0p 0.8 cP 1 cP viscosity
βp 3.2 10−3 1/bar 0 viscosibility

for density and viscosity of dead-oil:

Rρ = cp (p − pref ), ρp = ρ0p (1 + Rρ + Rρ2 /2), (A3)

Rμ = β(p − pref ), μp = μ0p (1 + Rμ + Rμ2 /2). (A4)

The component mole fractions for this model are defined as xoo = xww = 1 and
xow = xow = 0. Initial saturation Sini = 0 and pressure pini = 200 bars were

28
applied for the initialization of the models.
The next table defines the main parameters for the compositional model de-
scribing the gas injection case. In this model, Peng and Robinson [34] Equation

Table A3: Miscible displacement parameters

C. Tc K Pc bar ω Mw g/mol ki,C1
C1 190.6 46.04 0.013 16.04 -
C8 575.0 28.79 0.312 107.0 0.037

of State was used to compute density ρp . For phase viscosity, the LBC correla-
tion [35] was applied. Phase fractions were defined based on a flash calculation
[24]. For the initialization of the compositional model, the initial composition
zini = {0.1, 0.9} and pressure pini = 100 bars was used.

29