IMPES and fully implicit schemes

Quick refresher

In the two-phase cell equations the flows depend on the pressures and, steeply, through the relative permeabilities, on the saturations. A time scheme decides the values of which moment to use: the end of the step (implicit) or the start (explicit).
IMPES: the mobilities come from the saturations at the start of the step, the pressures are found implicitly, by one linear system with one unknown per cell, and the saturations are found explicitly, from the new flows. The water balance holds exactly.
The explicit recalculation of saturation is stable as long as the Courant number of the cell, $C = \Delta t\,Q\,f_w'/V_p$ , is not larger than one. (The Courant number is the distance that a saturation change travels in one step, measured in cells.) The largest step is proportional to the pore volume of the cell and inversely proportional to the flow through it. It is set by the small cells near wells and by the cells in the front zone.
The fully implicit method is stable at any step. The nonlinear system for all the unknowns is solved by Newton's method: we replace the residuals (the amounts by which the equations are not yet satisfied) by their linear part and solve a linear system with a sparse Jacobian matrix (nearly all its entries are zero). Near the solution the number of correct digits doubles at every iteration.
With a large step Newton's method may fail to converge. The curve of the fractional flow of water (the share of water in the flowing fluid) is S-shaped, and the tangent (the straight-line approximation) throws the iteration beyond the root (the answer). We limit the corrections, and if the iterations do not converge, we cut the step.
Stability does not replace accuracy: a large implicit step smears the front. We choose the step from the target changes of saturation and pressure. The sequential and adaptive implicit schemes are compromises between IMPES and the fully implicit method.

Key relationships

V_{p,i}\,c_{t,i}^{n}\,\frac{p_i^{n+1} - p_i^{n}}{\Delta t} = \sum_j a_{ij}^{n}\left(p_j^{n+1} - p_i^{n+1}\right) - q_{t,i},\qquad a_{ij}^{n} = 0.0853\,T_{ij}\,\lambda_{t,ij}^{n}

IMPES: the pressure equation with the mobilities at the start of the step

C_i = \frac{\Delta t\,Q_i\,f_w'}{V_{p,i}} \le 1

Courant condition for the explicit recalculation of saturation

\mathbf{J}\,\delta\mathbf{x} = -\mathbf{R},\qquad \mathbf{x}^{(\nu+1)} = \mathbf{x}^{(\nu)} + \delta\mathbf{x}

Newton iteration: the Jacobian matrix, the residuals, the correction of the unknowns

f_w' = \frac{f_w\,(1 - f_w)}{1 - S_{wi} - S_{or}}\left(\frac{n_w}{S_n} + \frac{n_o}{1 - S_n}\right)

slope of the fractional flow curve for the Corey model

\Delta t^{n+1} = \Delta t^{n}\cdot\min\!\left(\frac{\Delta S_{\text{target}}}{\max|\Delta S|},\ \frac{\Delta p_{\text{target}}}{\max|\Delta p|}\right)

choice of the next step from the target changes per step

Common pitfalls

Taking the IMPES step limit for a property of the whole model: it is set by one cell with the largest flow per unit pore volume, and it changes together with the saturations.
Repeating an IMPES step with updated mobilities in the hope of getting around the limit: such simple iterations (recalculating again and again with the latest values) converge only if the Courant number is below one.
Loosening the convergence tolerances (how small the residual must be before we stop iterating) to speed up the run: iterations that have not converged break the mass balance.
Taking converged Newton iterations for a sign of an exact solution: the errors of the step and of the grid remain.
Taking a large implicit step during fast processes, such as the front approaching a well or wells being started and stopped, without checking the sensitivity to the step.

After this lesson you can

write the IMPES step for the two-phase cell equations (an implicit pressure equation with the mobilities at the start of the step, and an explicit recalculation of saturation), and check the water balance
derive the stability condition of the explicit recalculation of saturation (the Courant condition), use it to find the largest step for the cells of a model, and explain why this step is set by the cells near wells and in the front zone
write the residuals of the cell equations of the fully implicit method (the amounts by which the equations are not yet satisfied), explain the structure of the Jacobian matrix, and carry out Newton iterations with a convergence check
explain why Newton's method may fail to converge (to settle on an answer) with a large step and how it is protected, and tell the convergence of iterations from the accuracy of the result
choose the time step from the target changes of saturation and pressure and from the progress of the iterations
describe the sequential and adaptive implicit schemes, and choose a scheme for a task

What makes the time step difficult

The cell equations for two phases are written in the lesson on the two-phase flow equations. For each cell there is an oil balance and a water balance. The flow of a phase through a face equals the transmissibility multiplied by the mobility of the phase in the upstream cell (the cell the phase flows out of) and by the difference of the phase potentials. A cell has two unknowns, the oil pressure $p_o$ and the water saturation $S_w$ . One question has stayed open. In the single-phase problem it was settled by a choice between the explicit and the implicit scheme (time step): the values of which moment, the start or the end of the step, should we put into the flows?

Here the question has two parts. The flows depend on the pressures linearly, and on the saturations through the relative permeabilities (the share of the rock's permeability that a phase can use), and this dependence is steep. In the middle of the two-phase zone a change of water saturation by 0.05 changes the relative permeabilities of water and oil by one and a half to two times, and near its edges by several times. So different combinations are possible: the pressures at the end of the step with the mobilities at the start, everything at the end of the step, or implicit treatment in only some of the cells. The schemes that simulators use are built on these combinations.

Time-step schemes for the equations of multiphase flow

Scheme: IMPES
Pressure: end of step
Saturations in the flows: start of step; then an explicit recalculation
The step is limited by: the Courant condition

Scheme: sequential implicit
Pressure: end of step
Saturations in the flows: end of step, with the total flows known
The step is limited by: accuracy

Scheme: fully implicit
Pressure: end of step
Saturations in the flows: end of step, together with the pressure
The step is limited by: accuracy and convergence of the iterations

Scheme: adaptive implicit
Pressure: end of step
Saturations in the flows: end of step in some cells, start of step in the others
The step is limited by: the Courant condition in the explicit cells

Scheme	Pressure	Saturations in the flows	The step is limited by
IMPES	end of step	start of step; then an explicit recalculation	the Courant condition
sequential implicit	end of step	end of step, with the total flows known	accuracy
fully implicit	end of step	end of step, together with the pressure	accuracy and convergence of the iterations
adaptive implicit	end of step	end of step in some cells, start of step in the others	the Courant condition in the explicit cells

In a multiphase model the pressure is almost always calculated implicitly. For small cells an explicit pressure calculation is limited to a step of seconds or minutes (why pressure is calculated implicitly). The schemes differ in how they treat the saturations.

IMPES: implicit pressure, explicit saturation

IMPES divides the step into two parts. First, from the saturations at the start of the step, we calculate the phase mobilities at the faces, taken from the upstream cell as in the two-phase lesson, and we “freeze” them for the whole step. With the mobilities known, the pressure equation becomes linear. It is the sum of the phase balances, weighted by the formation volume factors (pressure equation). It is built in the same way as the single-phase one. In the link coefficient the total mobility (the sum of the mobilities of the phases) stands in place of $1/\mu$ .

IMPES: the pressure equation of a cell with the mobilities at the start of the step

V_{p,i}\,c_{t,i}^{n}\,\frac{p_i^{n+1} - p_i^{n}}{\Delta t} = \sum_j a_{ij}^{n}\left(p_j^{n+1} - p_i^{n+1}\right) - q_{t,i},\qquad a_{ij}^{n} = 0.0853\,T_{ij}\,\lambda_{t,ij}^{n}

where:

$p_i^{n},\ p_i^{n+1}$ oil pressure in the cell at the start and at the end of the step, MPa
$V_{p,i}$ pore volume of the cell; $V_{p,i}c_{t,i}^{n}$ is its storage capacity at the saturations at the start of the step, as in the single-phase lesson, m³
$c_{t,i}^{n}$ total compressibility $S_oc_o + S_wc_w + c_f$ at the start of the step, 1/MPa
$a_{ij}^{n}$ link coefficient of cells $i$ and $j$ with the mobilities at the start of the step, m³/(d·MPa)
$\lambda_{t,ij}^{n}$ total mobility at the face per unit permeability, $(k_{ro}/\mu_o)_\uparrow + (k_{rw}/\mu_w)_\uparrow$ : each phase is taken from its own upstream cell, 1/(mPa·s)
$T_{ij}$ transmissibility of the connection (the geometric part with the absolute permeability), mD·m
$q_{t,i}$ withdrawal of liquid by a well at reservoir conditions (negative for injection), m³/d

Units: metric oilfield units: $T$ in mD·m, $\mu$ in mPa·s, $p$ in MPa, $t$ in d, volumes at reservoir conditions; 0.0853 is the constant of Darcy's law

Holds when

The mobilities and compressibilities are taken at the values at the start of the step, so the equation is linear in the pressures at the end of the step.
The phases are slightly compressible, and the formation volume factors change little over a step. Capillary pressure (the pressure difference between the phases that comes from surface forces in the pores) and gravity are taken into account on the right-hand side, at the values at the start of the step.

This is the single-phase cell equation (equation of a cell) with coefficients recalculated from the saturations at the start of the step. It is solved in the same way, by one system for all cells with one unknown per cell. If the rock and the fluids are incompressible, the left side disappears: the pressure adjusts instantly to the withdrawals and injection.

Then, from the new pressures, we find the water flows, with the water mobility at the start of the step from the upstream cell, and we recalculate the water saturation of each cell separately, without a system of equations:

IMPES: explicit recalculation of the water saturation

V_{p,i}\,\frac{(S_w/B_w)_i^{n+1} - (S_w/B_w)_i^{n}}{\Delta t} = \sum_j 0.0853\,T_{ij}\left(\frac{k_{rw}}{\mu_w B_w}\right)_{\!\uparrow}^{\!n}\left(p_j^{n+1} - p_i^{n+1}\right) - q_{w,i}

where:

$V_{p,i}$ pore volume of the cell, m³
$S_w,\ B_w$ water saturation and formation volume factor of water; $B_w$ at the end of the step is taken at the new pressure, —, m³/m³
$\left(k_{rw}/(\mu_w B_w)\right)_\uparrow^{n}$ water mobility at the start of the step from the upstream cell, 1/(mPa·s)
$p_j^{n+1} - p_i^{n+1}$ pressure difference at the end of the step, MPa
$q_{w,i}$ withdrawal of water by a well at standard conditions; the share of water in it follows from the mobilities at the start of the step, m³/d

Units: metric oilfield units, as in the pressure equation; volumes of water at standard conditions

Holds when

Without capillary pressure and gravity. With them, we put the difference of the water potentials $\Phi_w = p_o - p_c - \gamma_w D$ in place of the pressure difference.
The pore volume is constant during the step.
The oil saturation is found from $S_o = 1 - S_w$ . With compressible phases, the linearisation of the pressure equation leaves a small error in the oil balance.

Everything on the right-hand side is known, and each cell is recalculated separately. A flow through a face enters the equations of the two neighbouring cells with opposite signs, so the water balance of the whole model holds exactly at any step, like the mass balance in the single-phase lesson.

The example of the lesson on the two-phase equations, a step of 200 days of waterflooding of a strip of four cells, is exactly one IMPES step for the incompressible case. The pressure drops are found from the total mobility at the start of the step, the water flows from the fractional flows of the upstream cells, and the saturations explicitly. Further in this lesson the same strip, with the same relative permeability curves, is used to test both schemes. The rock and the fluids are incompressible, so the total flow through each face is known in advance, 25.5 m³/d at reservoir conditions, and the water balance of a cell over a step takes a simple form:

IMPES on a one-dimensional strip with a given injection

S_i^{n+1} = S_i^{n} + K\left[f_w\!\left(S_{i-1}^{n}\right) - f_w\!\left(S_i^{n}\right)\right],\qquad K = \frac{\Delta t\,q_t}{V_p}

where:

$S_i$ water saturation of cell $i$ , —
$f_w$ fractional flow of water $\lambda_w/\lambda_t$ ; for the first cell, the injected water takes the place of $f_w(S_0)$ , $f_w = 1$ , —
$q_t$ total flow through each face at reservoir conditions, m³/d
$V_p$ pore volume of the cell, m³
$K$ the share of the pore volume of a cell that the total flow passes in a step, —

Units: $\Delta t$ in d, $q_t$ in m³/d, $V_p$ in m³

Holds when

One-dimensional horizontal flow from cell 1 to the last one. Incompressible rock and fluids. Injection at a given rate.
No capillary pressure and no gravity. The fractional flow of water at a face is taken from the upstream cell.

The new saturation is the old one plus the difference of the fractional flows of water at the inlet and the outlet, multiplied by $K$ . Pressures do not enter this form: in a one-dimensional incompressible problem they do not affect the transport of saturation. They decide the injectivity of the wells (how much water an injection well takes for a given pressure difference), which the lesson on the two-phase equations covers.

The IMPES step limit: the Courant condition

Why can the explicit recalculation of saturation not use an arbitrarily large step? We can write the difference of the fractional flows of water at the inlet and the outlet of a cell through the slope of the curve. It equals the slope $f_w'$ at some saturation $\tilde S$ between the saturations of the cell and its upstream neighbour, multiplied by the difference of these saturations. Then the new saturation of the cell is a mixture of its old saturation and the saturation of the neighbour:

The explicit recalculation as a mixture of the saturations of the cell and its neighbour

S_i^{n+1} = \left[1 - K f_w'(\tilde S)\right]S_i^{n} + K f_w'(\tilde S)\,S_{i-1}^{n}

where:

$\tilde S$ a saturation between $S_i^n$ and $S_{i-1}^n$ at which the slope of the curve equals the slope of the chord (the straight line joining the two points of the curve; this is the mean value theorem), —
$f_w'$ slope of the fractional flow curve $\dd f_w/\dd S_w$ , —
$K$ $\Delta t\,q_t/V_p$ , as in the IMPES equation for the strip, —

Units: dimensionless

Holds when

An exact regrouping of the IMPES equation for the strip. $\tilde S$ is unknown, but it lies between the saturations of the cell and its neighbour.

As long as the weight $Kf_w'$ is not larger than one, both weights are non-negative, and the new saturation lies between the old one and the saturation of the neighbour. The calculation cannot “jump over” the neighbour and go beyond the physical limits. This is the same argument as for the explicit pressure scheme in the single-phase lesson.

From this we get the stability condition of the explicit recalculation, the Courant–Friedrichs–Lewy condition (CFL, Courant condition). For a cell with several neighbours and wells, we take the whole flow into the cell in place of $q_t$ .

The Courant condition for the explicit recalculation of saturation

C_i = \frac{\Delta t\,Q_i\,f_w'}{V_{p,i}} \le 1\quad\Longrightarrow\quad \Delta t \le \min_i \frac{V_{p,i}}{Q_i\,f_w'}

where:

$C_i$ Courant number of the cell: the fraction of the cell that the saturation passes in a step, —
$Q_i$ total flow entering the cell from the neighbours and wells, at reservoir conditions, m³/d
$f_w'$ the largest slope of the fractional flow curve between the saturation of the cell and the saturations of its upstream neighbours; to be safe for any saturations, the largest slope on the whole curve, —
$V_{p,i}$ pore volume of the cell, m³
$\Delta t$ time step, d

Units: any consistent units: $Q$ in m³/d, $V_p$ in m³, $\Delta t$ in d

Holds when

A sufficient condition for upstream weighting in incompressible flow: if it holds, the calculation is stable. Capillary pressure and gravity add their own limits (see the advanced section).
In a one-dimensional strip $Q_i = q_t$ ; in two and three dimensions it is the sum of the incoming flows.
With a noticeably larger step the calculation loses monotonicity (the saturations no longer change steadily in one direction): they oscillate from cell to cell and go beyond the physical limits.

The quantity $f_w'u_t/\phi$ is the speed at which a saturation moves. The condition requires that in one step it passes not more than one cell. This is the same Courant number as in the lesson on numerical diffusion, but with the speed of the saturation. The limit is proportional to the pore volume of the cell and inversely proportional to the flow through it. It depends linearly on the cell size, not on its square, as for the explicit pressure scheme.

It is convenient to calculate the slope of the fractional flow curve for the Corey model with a ready formula. We get it by differentiating $f_w = \lambda_w/(\lambda_w + \lambda_o)$ , using the fact that $\lambda_o\lambda_w/\lambda_t^2 = f_w(1 - f_w)$ :

The slope of the fractional flow curve for the Corey model

f_w' = \frac{\dd f_w}{\dd S_w} = \frac{f_w\,(1 - f_w)}{1 - S_{wi} - S_{or}}\left(\frac{n_w}{S_n} + \frac{n_o}{1 - S_n}\right)

where:

$f_w$ fractional flow of water $\lambda_w/\lambda_t$ (without gravity and capillary pressure), —
$S_n$ normalised water saturation $(S_w - S_{wi})/(1 - S_{wi} - S_{or})$ , —
$S_{wi},\ S_{or}$ irreducible water saturation and residual oil saturation, —
$n_w,\ n_o$ Corey exponents for water and oil, —

Units: dimensionless

Holds when

The Corey model (lesson on relative permeability). The viscosities are constant.
At the ends of the mobile interval, where only one phase flows, the slope is zero.

The slope is small where almost one phase flows, and it is largest in the two-phase zone, at the inflection point of the $f_w$ curve (where the curve turns from bending up to bending down). For the curves of the lesson it reaches 5.22 at $S_w \approx 0.52$ (see the figure). Saturations near this value are carried fastest, and the cells they pass through, in the front zone, limit the step.

CalculatedThe slope of the fractional flow curve

f_w'

for the relative permeability curves of the lesson Relative permeability, for an oil viscosity of 1.2 mPa·s and a water viscosity of 0.5 mPa·s, and the saturations of the strip cells at the start of the step of 200 days. The largest slope is 5.22, at a water saturation of about 0.52. It lies between the saturation of cell 3 (0.30) and that of its upstream neighbour, cell 2 (0.55), so the step of the strip is limited by cell 3.

Worked exampleThe IMPES step limit: the strip and a cell near a well

Given (assumed data)

Strip (lessons Pressure equation for single-phase flow and Two-phase flow equations): 4 cells, the pore volume of each $V_p$ = 150,000 m³; total flow through the faces $q_t$ = 25.5 m³/d at reservoir conditions
Relative permeabilities and fluids: Corey model: $S_{wi}$ = 0.22, $S_{or}$ = 0.25, $k_{ro}'$ = 0.85, $k_{rw}'$ = 0.25, $n_o$ = 2.5, $n_w$ = 3; $\mu_o$ = 1.2 mPa·s, $\mu_w$ = 0.5 mPa·s
Water saturation of cells 1–4 at the start of the step: 0.70; 0.55; 0.30; 0.22, the start of the step of 200 days of the lesson on the two-phase equations
A cell near an injection well of a three-dimensional model: 10 × 10 m in plan, thickness 1.5 m, $\phi$ = 0.20; the same injection passes through it

Solution

The time in which the total flow passes through the pore volume of a cell: $V_p/q_t = 150{,}000/25.5 = 5882.4$ d.
The slope $f_w'$ in cell 2 ( $S_w$ = 0.55): $S_n = 0.62264$ , $f_w = 0.66077$ . The factor is $f_w(1 - f_w)/0.53 = 0.42293$ , the bracket is $3/0.62264 + 2.5/0.37736 = 11.4432$ , and $f_w' = 4.840$ . In the same way, in cell 1 ( $S_w$ = 0.70) it is 0.2902, and in cell 3 ( $S_w$ = 0.30) it is 0.1562. In cell 4 the water is immobile and the slope is zero. The largest slope on the whole curve is 5.218 at $S_w$ = 0.52.
In the stability condition the slope is taken at an unknown point between the saturation of the cell and that of its upstream neighbour. So to be safe we take the largest slope on this interval: cell 1 (from 0.70 to 0.75, the injected water) 0.2902; cell 2 (0.55–0.70) 4.840; cell 3 (0.30–0.55, the interval contains the inflection point) 5.218; cell 4 (0.22–0.30) 0.1562.
Courant numbers for a step of 200 d, $C_i = 200\,f_w'/5882.4$ : 0.0099, 0.165, 0.177 and 0.0053. All are much smaller than one.
The step limit of each cell is $5882.4/f_w'$ : about 20,300, 1215, 1127 and 37,700 d. At present the step of the strip is limited by cell 3, 1127 d. This is also the guarantee for any position of the front, $5882.4/5.218$ .
The cell near the well: pore volume $10\cdot10\cdot1.5\cdot0.20 = 30$ m³, limit $30/(25.5\cdot5.218) = 0.2255$ d, that is 5.4 h. A 30-year forecast with such a step is about 48.6 thousand steps.

Answer: A step of 200 d on the strip is stable with a large margin (the largest Courant number is 0.177), and the limit is 1127 d. For the cell of 30 m³ near the well, with the same flow, the limit is 5.4 h.

What it means. The IMPES limit is set not by the formation as a whole but by one cell, the one with the largest flow per unit pore volume and with saturations on the steep part of the $f_w$ curve. In the large cells of the strip it is measured in years, and in a small cell near a well in hours. A 30-year forecast would need tens of thousands of steps, although in the rest of the formation the step could be thousands of times larger. For pressure, the explicit limit of such cells is seconds (why pressure is calculated implicitly). The limit for saturation is less strict, because it is proportional to the cell size and not to its square. This is why IMPES calculates the pressure implicitly and the saturation explicitly. The data in this example are assumed for illustration.

Worked exampleAn IMPES step of 3000 days, above the limit

Given (assumed data)

Strip, relative permeabilities and initial state: from the previous example: $S_w$ = 0.70; 0.55; 0.30; 0.22; $V_p$ = 150,000 m³; $q_t$ = 25.5 m³/d
Step: 3000 d, 2.7 times larger than the guaranteed limit of 1127 d

Solution

$K = 3000\cdot25.5/150{,}000 = 0.51$ : in a step 51% of its pore volume passes through each cell.
Fractional flows of water at the start of the step: 0.99481, 0.66077, 0.00364 and 0, as in the lesson on the two-phase equations.
Explicit recalculation: $S_1 = 0.70 + 0.51\cdot(1 - 0.99481) = 0.70265$ ; $S_2 = 0.55 + 0.51\cdot(0.99481 - 0.66077) = 0.72036$ ; $S_3 = 0.30 + 0.51\cdot(0.66077 - 0.00364) = 0.63514$ ; $S_4 = 0.22 + 0.51\cdot0.00364 = 0.22186$ .
Courant numbers with the slopes of the previous example, $C_i = 0.51\,f_w'$ : 0.148, 2.47, 2.66 and 0.080. For cells 2 and 3 they are larger than one.
The reference is the same cell equations calculated with a step of 0.25 d: after 3000 days the saturations are 0.7025, 0.6268, 0.5654 and 0.3828.
Water balance: the sum of the saturation increments, 0.51, equals $K$ . All the water injected in the step, $25.5\cdot3000 = 76{,}500$ m³, stayed in the strip, because the share of water in the production was taken from the saturation of cell 4 at the start of the step, and it is zero.

Answer: The explicit step gave cell 2 a higher water saturation than cell 1, the cell that feeds it with water (0.72036 against 0.70265). Almost no water got into cell 4 (0.22186 against 0.3828 in the reference).

What it means. While the saturation profile is monotonic (it only falls from the injection end to the far end), and with injection at one end it is, a cell cannot get a higher water saturation than its upstream neighbour, which feeds it. The value 0.72036 is a sign that the step went beyond the stability limit: the Courant number of cell 2 is 2.47. The balance does not show the error, because it holds exactly. Only the physically impossible saturations reveal it. A second property of the explicit scheme can also be seen: in one step the water goes no farther than the neighbouring cell, however long the step is, so cell 4 stayed almost dry. If we repeat such steps, the saturations begin to oscillate from cell to cell and go beyond the physical limits (see the figure below). The data in this example are assumed for illustration.

CalculatedThe same strip on a grid of 20 cells of 50 m (the IMPES limit is 225 days): the water saturation after 5400 days of injection. The dots are the reference, the same cell equations with a step of 0.25 d. Top, IMPES: with a step of 450 d (

C \approx 2

) the profile has a “saw-tooth” shape, and with a step of 900 d (

C \approx 4

) the saturation of 0.985 exceeds the limit

1 - S_{or}

= 0.75. Bottom, the fully implicit method with the same steps: the profiles are smooth, but the front is smeared and has moved ahead. The larger the step, the farther ahead it has moved.

The fully implicit method: residuals and the Jacobian matrix

In the fully implicit method we put the values at the end of the step, both pressures and saturations, into the flows, the mobilities and the well rates. The cell equations become nonlinear in the unknowns, because the mobilities at the end of the step depend on the saturations we are looking for. We write them as residuals: the difference between the left and right sides, which must become zero in the solution:

Residual of the balance of a phase in a cell (fully implicit method)

R_{\alpha,i} = V_{p,i}\,\frac{(S_\alpha/B_\alpha)_i^{n+1} - (S_\alpha/B_\alpha)_i^{n}}{\Delta t} - \sum_j 0.0853\,T_{ij}\left(\frac{k_{r\alpha}}{\mu_\alpha B_\alpha}\right)_{\!\uparrow}^{\!n+1}\left(\Phi_{\alpha,j} - \Phi_{\alpha,i}\right)^{n+1} + q_{\alpha,i}^{n+1} = 0

where:

$R_{\alpha,i}$ residual of the balance of phase $\alpha$ in cell $i$ : how much the current approximation breaks the balance, m³/d (standard conditions)
$\alpha$ phase: $o$ is oil, $w$ is water, —
$\Phi_\alpha$ phase potential $p_\alpha - \gamma_\alpha D$ ; the water pressure is $p_w = p_o - p_c(S_w)$ , MPa
$(\cdot)_\uparrow^{n+1}$ mobility of the phase at the end of the step from the upstream cell, 1/(mPa·s)
$q_{\alpha,i}^{n+1}$ withdrawal of the phase by a well, calculated from the pressure and the mobilities at the end of the step, m³/d
$V_{p,i},\ T_{ij},\ \Delta t$ pore volume of the cell, transmissibility of the connection, step, m³, mD·m, d

Units: metric oilfield units; volumes of the phases at standard conditions

Holds when

The same cell equation as in the two-phase lesson, with all the quantities of the right-hand side taken at the end of the step.
The unknowns of a cell are $p_o$ and $S_w$ , and the other quantities are expressed through the closure relations.
In the black oil model with gas, a cell has three equations and three unknowns.

For $N$ cells there are $2N$ nonlinear equations with $2N$ unknowns. The solution of the step is a set of pressures and saturations at which all the residuals are zero to a given accuracy. Then the balance of each phase in each cell holds for the values at the end of the step, and the scheme is stable at any step.

Such a system is solved by Newton's method, also called the Newton–Raphson method (Newton's method). We approximate the residuals by their linear part, the first terms of the Taylor series around the current approximation, and find the correction that makes this linear part zero. The coefficients of the linear part, the partial derivatives of the residuals with respect to the unknowns, form the Jacobian matrix, or simply the Jacobian. It is a table that shows how each residual changes when each unknown changes.

A Newton iteration for a system of equations

\mathbf{J}\!\left(\mathbf{x}^{(\nu)}\right)\delta\mathbf{x} = -\mathbf{R}\!\left(\mathbf{x}^{(\nu)}\right),\qquad \mathbf{x}^{(\nu+1)} = \mathbf{x}^{(\nu)} + \delta\mathbf{x},\qquad J_{kl} = \frac{\partial R_k}{\partial x_l}

where:

$\mathbf{x}$ vector of the unknowns at the end of the step $(p_{o,1},\ S_{w,1},\ \ldots,\ p_{o,N},\ S_{w,N})$ , MPa, —
$\mathbf{R}$ vector of the residuals of all equations, m³/d
$\mathbf{J}$ Jacobian matrix of size $2N \times 2N$ , in units of residuals per unit of the unknowns
$\delta\mathbf{x}$ correction of the unknowns at the iteration, MPa, —
$\nu$ iteration number; the initial approximation $\mathbf{x}^{(0)}$ is the values at the start of the step or a forecast from the previous steps, —

Units: in units of the residuals and the unknowns

Holds when

The residuals are smooth functions of the unknowns. With upstream weighting the derivatives change abruptly when a flow changes direction.
At each iteration we solve a linear system. The matrix $\mathbf{J}$ is not inverted.

The method replaces a nonlinear system by a sequence of linear ones. Near the solution it converges quadratically: the error of the next iteration is proportional to the square of the error of the previous one. If it was of the order of $10^{-3}$ , it becomes of the order of $10^{-6}$ . The price is that at each iteration we must calculate all the derivatives and solve a linear system with $2N$ unknowns.

The Jacobian matrix consists almost entirely of zeros. The residuals of a cell depend only on its own unknowns and on the unknowns of the neighbours with which it shares faces: on three cells in a one-dimensional strip, and on seven in a three-dimensional Cartesian grid. A pair “equations of a cell, unknowns of a cell” corresponds to a block of 2 × 2, and in a one-dimensional strip the matrix is block-tridiagonal. Upstream weighting makes it non-symmetric. The saturation of the downstream cell does not enter the mobility at their common face, and without capillary pressure the derivatives of the flow with respect to it are zero (see the figure; more details are in the advanced section).

SchematicThe structure of the Jacobian matrix of a strip of four cells during waterflooding from cell 1 to cell 4, without capillary pressure (a sketch). The rows are the oil and water balances of the cells, and the columns are the pressure and saturation of the cells. Non-zero elements appear only for the cell itself and its neighbours. The balances of a cell do not depend on the saturation of its downstream neighbour.

In a one-dimensional incompressible problem with a given injection everything is simpler. The total flow through each face is known, the pressures drop out of the water balances, and the Jacobian matrix for the saturations becomes lower triangular (all its entries above the diagonal are zero): the balance of a cell depends only on the cell itself and on the cell upstream of it. Then we solve the equations one at a time, along the flow: first cell 1, then cell 2, and so on. It is convenient to refer the residual to the pore volume of the cell, and then it becomes dimensionless:

Implicit water balance of a strip cell

R_i(S) = \left(S - S_i^{n}\right) - K\left[f_w\!\left(S_{i-1}^{n+1}\right) - f_w(S)\right] = 0,\qquad R_i'(S) = 1 + K f_w'(S)

where:

$S$ unknown water saturation of cell $i$ at the end of the step, —
$R_i$ residual of the water balance over the step, as a fraction of the pore volume of the cell, —
$R_i'$ derivative of the residual with respect to $S$ , a Jacobian of size 1 × 1, —
$K$ $\Delta t\,q_t/V_p$ , —

Units: dimensionless

Holds when

A one-dimensional strip, incompressible rock and fluids, injection at a given rate. The fractional flow of water at a face is taken from the upstream cell.
The cells are solved along the flow: by the time we calculate cell $i$ , the saturation $S_{i-1}^{n+1}$ is already found.

The derivative $R'$ is not smaller than one: the residual grows monotonically with $S$ , and the root is unique. The Newton iteration is $S^{(\nu+1)} = S^{(\nu)} - R/R'$ .

Newton iterations and convergence criteria

Graphically, a Newton iteration is a tangent (a straight line that touches the curve). We replace the residual curve at the current point by its tangent, and we take the point where the tangent crosses zero. If the curve does not change the direction of its bend between the approximation and the root, the iterations approach the root faster and faster. But the $f_w$ curve has an inflection point (the point where its bend changes direction), and so does the residual. At low water saturations the curve is flat, and near the inflection point it is steep. A tangent drawn on the flat part throws the approximation far beyond the root, and from there the next tangent throws it back. The iterations can swing to and fro without converging.

CalculatedThe residual of the water balance of cell 3 over a step of 3000 days against its water saturation (the example below). The tangent at the starting point 0.30 throws the Newton iteration beyond the root, to 0.715, and the next tangent throws it back, to 0.275. Without a limit the iterations swing to and fro. If we do not let the correction pass the inflection point at 0.52, the following iterations converge to the root, 0.512.

Worked exampleThe same step of 3000 days, fully implicit: Newton iterations

Given (assumed data)

Strip, state and step: as in the example of the IMPES step of 3000 d; $K$ = 0.51
Residual of cell $i$: $R = (S - S_i^n) - 0.51\,[f_{\text{in}} - f_w(S)]$ , $R' = 1 + 0.51\,f_w'(S)$ , where $f_{\text{in}}$ is the fractional flow of water of the upstream cell at the end of the step
Convergence: we iterate until $|R| < 10^{-12}$ to show how the iterations proceed. The usual simulator tolerance (the residual at which it stops) is of the order of $10^{-3}$ , a thousandth of the pore volume of the cell

Solution

We solve the cells along the flow. Cell 1 ( $f_{\text{in}}$ = 1) converges in 3 iterations to 0.70232, and cell 2 converges in 5 to 0.608005. The inflow into cell 3 is $f_{\text{in}} = f_w(0.608005) = 0.881721$ .
Cell 3. The initial approximation is the saturation at the start of the step, 0.30: $f_w$ = 0.00364, $f_w'$ = 0.1562. The residual is $R = -0.51\cdot(0.881721 - 0.00364) = -0.44782$ , the derivative is $R' = 1.07966$ , and the correction $0.44782/1.07966 = 0.41478$ leads to 0.71478.
Without a limit the iteration goes beyond the root. At 0.71478: $f_w$ = 0.99802, $f_w'$ = 0.1521. The residual is already positive, $R = 0.41478 - 0.51\cdot(0.881721 - 0.99802) = +0.4741$ , $R' = 1.0776$ , and the next tangent returns the iteration to $0.71478 - 0.4741/1.0776 = 0.2748$ . Then come 0.7347, 0.2489, 0.7460 and so on. The approximations swing between ≈ 0.24 and ≈ 0.75, the residual stays near ±0.5, and the step does not converge.
With a limit, the correction does not cross the inflection point of $f_w$ in one iteration, and the approximation stops at it, at 0.52180. There $f_w$ = 0.51724 and $f_w'$ = 5.2177: $R = (0.52180 - 0.30) - 0.51\cdot(0.881721 - 0.51724) = 0.0359$ , $R' = 3.6610$ , and the new approximation is $0.52180 - 0.0359/3.6610 = 0.51199$ .
At 0.51199: $f_w$ = 0.466210, $R = (0.51199 - 0.30) - 0.51\cdot(0.881721 - 0.466210) \approx 7.9\cdot10^{-5}$ , already smaller than the usual tolerance $10^{-3}$ . One more iteration gives 0.511968 and $|R| \sim 10^{-9}$ . The exponent of the residual roughly doubles at each iteration, which is quadratic convergence.
Cell 4: the inflow is $f_w(0.511968) = 0.46610$ , and the solution 0.41097 comes in 5 iterations.
Water balance at reservoir conditions: the saturation increments add up to 0.46326, that is $150{,}000\cdot0.46326 = 69{,}489$ m³. The injected volume is 76,500 m³. The produced volume, with the fractional flow of water of cell 4 at the end of the step, is $76{,}500\cdot0.09165 = 7011$ m³, and the difference is 69,489 m³.
Comparison with the reference (the same cell equations with a step of 0.25 d), to four digits: the reference is 0.7025, 0.6268, 0.5654, 0.3828; the implicit step gives 0.7023, 0.6080, 0.5120, 0.4110; the deviations are −0.0002, −0.0188, −0.0534 and +0.0282.

Answer: The implicit step of 3000 d gives 0.70232, 0.60801, 0.51197 and 0.41097: the saturations decrease downstream, the water has reached cell 4, and the balance holds. In cell 3 Newton's method without a limit does not converge. With a limit it converges in 4 iterations to

10^{-12}

, and it reaches the usual tolerance of

10^{-3}

already in two.

What it means. The fully implicit method went through a step on which IMPES gave physically impossible saturations, but stability does not mean accuracy. Compared with the reference, the front is smeared: cells 2 and 3 hold less water than they should, and more water has come into cell 4. A large implicit step, as it were, spreads the injection over all cells at once. The error falls in proportion to the step (first order in time), and how to choose the step for accuracy is described below. The difficulty of Newton's method is not special to this example: when the Corey exponents are larger than one, the $f_w$ curve is S-shaped, and with a large step the tangent from the flat part throws the iteration beyond the root. The data in this example are assumed for illustration.

Simulators protect Newton's method from such oscillations in several ways. They limit the change of saturation per iteration (for example to 0.1–0.2), they do not let the correction cross the inflection point of the fractional flow curve, and they reduce the correction if the residual grows. If the iterations still do not converge in the given number (usually from ten to twenty), the step is discarded and repeated with a smaller one. When should we consider the iterations converged? Usually we need the residuals to be small, both in each cell and in the sum over the model:

Convergence criteria of the iterations within a step

\max_i\,\frac{\Delta t\,B_\alpha\left|R_{\alpha,i}\right|}{V_{p,i}} \le \varepsilon_R,\qquad \frac{\Delta t\,B_\alpha\left|\sum_i R_{\alpha,i}\right|}{\sum_i V_{p,i}} \le \varepsilon_{MB}

where:

$R_{\alpha,i}$ residual of the balance of phase $\alpha$ in cell $i$ after the iteration, m³/d (standard conditions)
$B_\alpha$ formation volume factor of the phase: it converts the residual into reservoir volumes, m³/m³
$\varepsilon_R$ tolerance on the largest normalised residual, a fraction of the pore volume of the cell, —
$\varepsilon_{MB}$ tolerance on the material balance error of the phase over the whole model, —

Units: dimensionless

Holds when

They are checked for each phase. A condition is often added that the corrections of the last iteration, in pressure and in saturation, are small.
Typical orders of the tolerances: $\varepsilon_R$ from $10^{-3}$ to $10^{-2}$ , $\varepsilon_{MB}$ from $10^{-7}$ to $10^{-6}$ . The specific values depend on the simulator and the task.

The first condition is local: in no cell has more than a given fraction of its pore volume been “lost” or “created” during the step. The second is the overall balance: the total error is small compared with the pore volume of the model. This is what simulators print as the material balance error.

How the time step is chosen

In IMPES the Courant condition limits the step from above. In the fully implicit method stability does not limit the step, and we choose it for two other reasons. The first is accuracy. The saturations and pressures must not change too much in a step, or a large implicit step smears the front and lags behind fast changes. The second is the convergence of the iterations: the larger the step, the farther the initial approximation is from the solution, and the harder it is for Newton's method. The simplest rule for choosing the step scales it by the changes over the previous step:

Choosing the next step from the target changes

\Delta t^{n+1} = \Delta t^{n}\cdot\min\!\left(\frac{\Delta S_{\text{target}}}{\max_i\left|\Delta S_i\right|},\ \frac{\Delta p_{\text{target}}}{\max_i\left|\Delta p_i\right|}\right),\qquad \Delta t^{n+1} \le f_{\text{growth}}\,\Delta t^{n}

where:

$\Delta S_i,\ \Delta p_i$ changes of the saturation and pressure of the cell over the previous step, —, MPa
$\Delta S_{\text{target}},\ \Delta p_{\text{target}}$ target changes per step, the accuracy settings, —, MPa
$f_{\text{growth}}$ largest allowed growth of the step, usually 2–3, —
$\Delta t^{n},\ \Delta t^{n+1}$ previous and next steps, d

Units: $\Delta t$ in d, $p$ in MPa

Holds when

The changes per step are roughly proportional to the step: a flow without abrupt events.
The step is also limited by the largest and smallest values, by report dates and by well events: after a start, a shut-in or a change in the way a well is run we begin again with a small step.

If over the previous step the saturation and pressure changed by less than the target, the next step is made larger, and if by more, it is made smaller. The target changes are set at several hundredths to tenths for saturation and about a megapascal for pressure. These are accuracy settings, and we check them with calculations at smaller values.

How a simulator goes through a step

1
Propose a step
By the step rule, with the limits on growth, the report dates and the well events.
2
Newton iterations
Residuals, the Jacobian matrix, the linear system, and the corrections with a limit. After each iteration we check the convergence criteria.
3
Not converged: cut the step
If the iterations did not converge in the given number, the step is discarded and repeated with a step 2 to 10 times smaller.
4
Converged: check the changes
If the changes over the step are much larger than the allowed ones, the step is also repeated with a smaller one.
5
Accept the step
Save the results, check the balance and propose the next step.

Worked exampleHow large to make the next step

Given (assumed data)

Step of 200 d (the lesson on the two-phase equations): the largest change of water saturation is 0.02234, in cell 3; after the step the saturations are 0.70018, 0.56136, 0.32234, 0.22012
Step control (assumed): target change of saturation per step 0.05; growth of the step by not more than a factor of two
Implicit step of 3000 d: from the example with Newton iterations

Solution

The proportional rule: $200\cdot0.05/0.02234 = 448$ d. The growth is limited to a factor of two, $2\cdot200 = 400$ d, so the next step is 400 d.
If the calculation is run with IMPES, we check the stability at the new saturations. The largest slope on the intervals is for cell 3: its interval from 0.32234 to 0.56136 contains the inflection point, and the slope is 5.218. $C = 400\cdot5.218/5882.4 = 0.355 < 1$ , so a step of 400 d is allowed.
The implicit step of 3000 d changed the saturation of cell 3 by $0.5120 - 0.30 = 0.2120$ . The rule would have proposed $3000\cdot0.05/0.2120 = 708$ d, which is 4.2 times smaller.
A check against the reference: the saturation of cell 3 after 3000 days is 0.5654. One implicit step gives 0.5120 (error −0.0534), and four steps of 750 d give 0.5496 (error −0.0158), which is 3.4 times smaller.

Answer: After a step of 200 d the next one is 400 d. A step of 3000 d is stable but too large for accuracy: the rule would have chosen about 708 d, and with steps of 750 d the error in cell 3 falls 3.4 times.

What it means. The step rule watches accuracy, not stability. It increases the step while little changes in a step, and it decreases the step when processes speed up: when the front approaches a well, and when wells are started or shut in. The target change of 0.05 is an assumed setting here. We check it by repeating the calculation with a smaller value and comparing the quantities for which the model is built: production, water cut, breakthrough time. The data in this example are assumed for illustration.

Sequential and adaptive implicit schemes

The sequential implicit scheme starts a step like IMPES: the pressure is implicit, with the mobilities at the start of the step. From the new pressures we find the total flows through the faces. Then we solve the saturation equations implicitly with these flows, by Newton's method, but only for the saturations. The Courant limit is removed, and the systems of equations are smaller than in the fully implicit method, because they are separate for the pressures and for the saturations. If the pressure and the saturations are strongly linked (compressible gas, gravity, capillary flows), the two stages are repeated within the step until they converge. Then the scheme gives the same solution as the fully implicit method.

The adaptive implicit method (AIM, adaptive implicit method) chooses the treatment cell by cell. Before a step we estimate for each cell the Courant number at the proposed step. Where it is above a threshold, usually near wells and in the front zone, we calculate the saturations implicitly. In the other cells we calculate them explicitly, as in IMPES. The common system of equations then has one unknown for an explicit cell and two (in the black oil model, three) for an implicit cell. The set of implicit cells changes from step to step together with the front.

Worked exampleWhich cells the adaptive scheme treats implicitly

Given (assumed data)

Strip of 20 cells of 50 m: pore volume of a cell 30,000 m³, flow 25.5 m³/d; the IMPES limit for any saturations is 225 d
State: after 5400 days of injection (the reference calculation, the figure above)
Proposed step: 450 d, twice the IMPES limit

Solution

$K = 450\cdot25.5/30{,}000 = 0.3825$ . The Courant number of a cell is $K$ multiplied by the largest slope $f_w'$ between its saturation and the saturation of its upstream neighbour.
It is larger than one for cells 9–12, with water saturations 0.604, 0.583, 0.523 and 0.299: $C$ = 1.08, 1.41, 2.00 and 2.00. These cells are treated implicitly. For the neighbouring cells 8 and 13, $C$ = 0.89 and 0.06, and for the others it is even smaller. They are treated explicitly.
Cell 12 stands ahead of the front: by its own saturation its Courant number is only 0.06. If we choose the implicit cells by their own saturations, it stays explicit, and after the step its water saturation is higher than that of its upstream neighbour: 0.579 against 0.566. If we choose by the interval (the largest slope between a cell and its upstream neighbour), cells 11–13 get 0.566, 0.473 and 0.328, almost as in the fully implicit method (0.567, 0.474 and 0.325).
Unknowns in the system of a step of a compressible model: 20 pressures and 4 saturations, 24 in total. In the fully implicit method there are 40, and in IMPES 20, but with a step of not more than 225 d.

Answer: 4 cells out of 20 are treated implicitly: the front zone and the cell ahead of it. The system of the step has 24 unknowns instead of 40.

What it means. In a one-dimensional strip the front zone is a noticeable share of the cells. In a three-dimensional model the fronts and the surroundings of the wells take a small part of the volume, and there are even fewer implicit cells. The example also shows the main subtlety: the cell to be made implicit is not only the one with a steep slope now but also the one that the front will reach during the step. So we choose the criterion with a margin, for example by the interval of saturations, as here, with a threshold below one, or taking the neighbours into account, and we review it at every step. The data in this example are assumed for illustration.

Comparison of the schemes

Scheme: IMPES
Unknowns per cell in the system: one (the pressure)
The step is limited by: the Courant condition
Where it is used: large cells, compositional models (which follow each component of the hydrocarbon mixture), streamlines (flow paths along which the saturation is carried)

Scheme: sequential
Unknowns per cell in the system: one, then the saturations separately
The step is limited by: accuracy; with strong coupling, outer iterations
Where it is used: the pressure and the transport are solved by different methods

Scheme: fully implicit
Unknowns per cell in the system: two, with gas three, in a compositional model more
The step is limited by: accuracy and convergence of the iterations
Where it is used: black oil, coning (water or gas drawn up to a well), gas, zones near wells

Scheme: adaptive implicit
Unknowns per cell in the system: one or all, depending on the cell
The step is limited by: the Courant condition in the explicit cells
Where it is used: compositional and large models

Scheme	Unknowns per cell in the system	The step is limited by	Where it is used
IMPES	one (the pressure)	the Courant condition	large cells, compositional models (which follow each component of the hydrocarbon mixture), streamlines (flow paths along which the saturation is carried)
sequential	one, then the saturations separately	accuracy; with strong coupling, outer iterations	the pressure and the transport are solved by different methods
fully implicit	two, with gas three, in a compositional model more	accuracy and convergence of the iterations	black oil, coning (water or gas drawn up to a well), gas, zones near wells
adaptive implicit	one or all, depending on the cell	the Courant condition in the explicit cells	compositional and large models

The uses of the schemes are a general guide. In industrial simulators of the black oil model the main choice is the fully implicit or the adaptive implicit scheme.

Which scheme pays off more depends on the task: on the number of unknowns per cell, the number of cells, how strict the Courant limit is near the wells, and how strongly the pressure is linked with the saturations. The work per cell in a fully implicit calculation grows faster than linearly with the number of unknowns in it (cost of a calculation). So in compositional models with ten or so components, schemes that treat fewer things implicitly pay off. A calculation on streamlines (streamlines) is also in essence IMPES: the pressure on the grid is implicit, and the transport along the streamlines is solved separately, which allows large steps. The smearing of fronts that large implicit steps and the grid itself add is described in the lesson Numerical dispersion.

Terms and symbols

IMPES method: The pressures of all cells are implicit, with the mobilities at the start of the step. The saturations are explicit, from the new flows. The step is limited by the Courant condition.
Fully implicit method: The flows, mobilities and withdrawals are taken at the values at the end of the step. The nonlinear system for all the unknowns is solved by Newton's method. It is stable at any step.
Courant number $C$: $\Delta t\,Q\,f_w'/V_p$ : the fraction of a cell that a saturation passes in a step. The explicit recalculation of saturation is stable for $C \le 1$ .
Residual: How much the current approximation breaks the balance of a phase in a cell. The solution of a step is the set of unknowns at which all the residuals are zero within a given tolerance.
Jacobian matrix: The matrix of the partial derivatives of the residuals with respect to the unknowns. It is sparse, made of blocks, and non-symmetric because of upstream weighting.
Newton's method: Iterations in each of which the residuals are replaced by their linear part, and a linear system with the Jacobian matrix is solved. Near the solution they converge quadratically.
Step cut: Repeating a step with a smaller $\Delta t$ if the iterations did not converge or the changes over the step are too large.
Sequential implicit scheme: The pressure is implicit, then the saturations are implicit with the total flows that were found. If the stages are repeated until convergence, it gives the solution of the fully implicit method.
Adaptive implicit scheme: Cells with a Courant number above a threshold are treated fully implicitly, and the others explicitly. The set of implicit cells changes from step to step.

Common mistakes and misconceptions

Wrong

“The fully implicit method is stable, so any step can be taken.”

In fact

Stability does not give accuracy. In the example of the lesson, one implicit step of 3000 days was wrong in the saturation of cell 3 by five hundredths, and the step rule would have taken a step four times smaller. In addition, with a large step Newton's method converges worse, and the simulator starts to cut the step.

Wrong

“If IMPES is unstable, it is enough to repeat the step several times with updated mobilities.”

In fact

Repeating the step is a simple iteration. It converges only if the Courant number is below one. With a step of 3000 days in the example of the lesson, the iterations swing between 0.24 and 0.75. To go beyond the limit we need an implicit calculation of saturation with derivatives, which is Newton's method.

Wrong

“The IMPES step limit is set by the pressure, as for the explicit pressure scheme.”

In fact

In IMPES the pressure is implicit and does not limit the step. The limit is set by the explicit transport of saturation. It is proportional to the pore volume of the cell and inversely proportional to the flow through it, so it is decided by the small cells near wells and the cells in the front zone. It is hours, not seconds as for explicit pressure.

Wrong

“The Newton iterations converged and the balance holds, so the result is accurate.”

In fact

The convergence criteria say only that the difference equations of the step are solved. In the example of the lesson the residual of cell 3 fell to $10^{-9}$ , and the saturation still differed from the reference by five hundredths. This is the error of a large step, and we check it with calculations at a smaller step (Finite differences: from equations to computation).

Wrong

“The Newton correction is found by inverting the Jacobian matrix.”

In fact

The inverse of a sparse Jacobian matrix is almost full, and for a large model it cannot be built. At each iteration we solve the linear system $\mathbf{J}\,\delta\mathbf{x} = -\mathbf{R}$ by direct or iterative methods (Solving systems of linear equations in simulators).

In depth: advanced detail

The Courant condition through harmonics, and the capillary limit

The Courant condition can also be obtained by the analysis of harmonics, like the stability conditions in the lesson on finite differences. Let the slope $f_w'$ be almost constant (a small disturbance of a uniform state). Then the explicit recalculation with upstream weighting multiplies a harmonic with phase $\xi$ per cell by a complex number (a number with a real and an imaginary part). The amplitude changes, and the harmonic is also shifted downstream:

Amplification factor of explicit upstream transport

G(\xi) = 1 - C\left(1 - \e^{-\mathrm{i}\xi}\right),\qquad |G|^2 = 1 - 2C(1 - C)(1 - \cos\xi)

where:

$C$ Courant number $\Delta t\,q_t f_w'/V_p$ , —
$\xi$ phase shift of the harmonic per cell, rad
$\mathrm{i}$ imaginary unit, —

Units: dimensionless

Holds when

Linearisation: a constant slope $f_w'$ , a uniform grid, boundaries not taken into account.

$|G| \le 1$ for all harmonics if and only if $0 \le C \le 1$ . The worst is the checkerboard harmonic ( $\xi = \pi$ ): $G = 1 - 2C$ , and for $C > 1$ it grows, changing sign at every step. A single cell whose own saturation changes alone would stay stable for $C < 2$ (the factor $1 - C$ ), but the “saw-tooth” of neighbouring cells would not. As for pressure, the limit is twice as strict as for a single “tank” (where the step limit comes from).

For a nonlinear $f_w$ curve the slope changes from cell to cell, and the guarantee is given by the largest slope on the interval of saturations. On the strip of the lesson the calculation gets visibly worse when the limit is exceeded by one and a half to two times (see the figure with the strip of 20 cells). But we cannot rely on such a margin: it depends on the relative permeability curves and on how fast the front passes the steep part.

The capillary term of the saturation equation is diffusive (saturation equation); it spreads the saturation out, like the diffusivity equation spreads pressure. If we calculate it explicitly, it has its own limit, built like the limit of the explicit pressure scheme, $\Delta t \le \tau_i$ :

The step limit for an explicit treatment of capillary pressure

\Delta t \le \frac{V_{p,i}}{\sum_j 0.0853\,T_{ij}\,\lambda_o f_w\left|\dd p_c/\dd S_w\right|}

where:

$\lambda_o f_w$ factor of the capillary flow $\lambda_o\lambda_w/\lambda_t$ per unit permeability, 1/(mPa·s)
$\left|\dd p_c/\dd S_w\right|$ slope of the capillary pressure curve, MPa
$T_{ij},\ V_{p,i}$ transmissibilities of the connections and the pore volume of the cell, mD·m, m³

Units: metric oilfield units: $\Delta t$ in d

Holds when

Explicit treatment of the capillary term. A sufficient condition of the same kind as $\Delta t \le \tau_i$ for pressure.
The coefficient is taken as the largest on the interval of saturations of the cell and its neighbours.

The pore volume is proportional to the cell size and the transmissibility is inversely proportional to it, so the limit falls as the square of the cell size. In large cells it is far away. In small-scale models, such as a core, a water cone (water drawn up towards a well in a cone shape) or the blocks of a fractured formation, it can be stricter than the Courant condition, and then the capillary term is treated implicitly. Gravity adds its own limit on transport, with the speed of the gravity separation of the phases.

Repeated IMPES within a step: simple iterations

Suppose that within an IMPES step we recalculate the mobilities from the new saturations and repeat the calculation of pressures and saturations, and so on until the changes become small. In the limit we get the solution of the equations with the mobilities at the end of the step, that is, the fully implicit ones. But does such a process converge? For the strip of the lesson the repetition reduces to simple iterations (repeated substitution of the new value into the same formula):

Simple iterations of the implicit balance of a cell

S^{(\nu+1)} = S_i^{n} + K\left[f_w\!\left(S_{i-1}^{n+1}\right) - f_w\!\left(S^{(\nu)}\right)\right],\qquad \left|S^{(\nu+1)} - S^{*}\right| \approx K f_w'(S^{*})\left|S^{(\nu)} - S^{*}\right|

where:

$S^{(\nu)}$ approximation at iteration $\nu$ , —
$S^{*}$ the solution, the root of the implicit balance of the cell, —
$K f_w'(S^*)$ the factor by which an iteration multiplies the error, —

Units: dimensionless

Holds when

A one-dimensional incompressible strip, as in the examples of the lesson.
The estimate on the right holds near the solution (linearisation).

The error is multiplied by $Kf_w'$ at each iteration, so the iterations converge only if $Kf_w' < 1$ . This is the Courant condition again. Repeating the IMPES step refines the solution when the step is allowed, but it does not allow us to exceed the limit.

Worked exampleSimple iterations and Newton's method for cell 3

Given (assumed data)

Strip and state: as in the examples of the lesson; a step of 1000 d, $K$ = 0.17
Cells 1 and 2: solved implicitly; the inflow into cell 3 is $f_{\text{in}}$ = 0.8020

Solution

The solution of the equation of cell 3 is $S^*$ = 0.41825. The slope there is $f_w'$ = 2.157, and the factor of the simple iterations is $Kf_w'$ = 0.37 < 1.
Simple iterations from 0.30: the absolute values of the residuals are 0.136, 0.0248, 0.00985, 0.00348 and so on. The sign alternates, and the absolute value falls about 2.7 times per iteration. To reach $|R| < 10^{-7}$ we need 14 iterations.
Newton's method from the same point: 0.136, 0.0197, 0.000561, $4.8\cdot10^{-7}$ , $3.5\cdot10^{-13}$ . It is below $10^{-7}$ after 4 iterations.
For a step of 3000 d the factor at the solution is 2.64 > 1: the simple iterations do not converge. They swing between 0.24 and 0.75, like Newton's method without a limit in the example of the lesson.

Answer: With an allowed step, repeating IMPES converges, but slowly: 14 iterations against 4 for Newton's method. With a step above the limit it does not converge.

What it means. Repeating the IMPES step with updated mobilities brings it closer to the fully implicit solution, but it does not remove the limit on the step: the factor $Kf_w'$ is the same Courant number. To go beyond the limit we need the derivatives, that is, Newton's method. The data in this example are assumed for illustration.

What the Jacobian matrix is made of

Take the flow of water into cell $i$ from the neighbour $j$ , without capillary pressure and gravity, in reservoir volumes: $F = 0.0853\,T_{ij}\lambda_{w\uparrow}(p_j - p_i)$ , where $\lambda_{w\uparrow} = k_{rw}/\mu_w$ of the upstream cell. Its derivatives with respect to the unknowns are the elements of the Jacobian matrix:

Derivatives of the water flow with respect to the unknowns

\frac{\partial F}{\partial p_j} = -\frac{\partial F}{\partial p_i} = 0.0853\,T_{ij}\,\lambda_{w\uparrow},\qquad \frac{\partial F}{\partial S_{w\uparrow}} = 0.0853\,T_{ij}\,\frac{\dd\lambda_w}{\dd S_w}\,(p_j - p_i),\qquad \frac{\partial F}{\partial S_{w\downarrow}} = 0

where:

$F$ flow of water from cell $j$ into cell $i$ at reservoir conditions, m³/d
$\lambda_{w\uparrow}$ water mobility $k_{rw}/\mu_w$ in the upstream cell (for $p_j > p_i$ , in cell $j$ ), 1/(mPa·s)
$S_{w\uparrow},\ S_{w\downarrow}$ water saturation of the upstream and the downstream cell, —
$\dd\lambda_w/\dd S_w$ slope of the water mobility curve, 1/(mPa·s)

Units: metric oilfield units

Holds when

Upstream weighting. The direction of the flow does not change during the iteration (otherwise the derivatives change abruptly).

The derivatives with respect to the pressures are the same link coefficient as in the single-phase cell equation, with the water mobility in place of $1/\mu$ . The derivative with respect to the saturation of the upstream cell is the slope of the mobility curve multiplied by the pressure difference. The derivative with respect to the saturation of the downstream cell is zero. This gives the non-symmetric structure of the Jacobian matrix in the figure.

The accumulation is added to the diagonal block of a cell: $V_p/\Delta t$ for the saturation and the compressibility terms for the pressure. A well with a given bottomhole pressure adds to the same block the derivatives of its own rate: with respect to pressure, the well index multiplied by the mobility of the phase; with respect to saturation, the index multiplied by the slope of the mobility curve and by the drawdown. This is how a well links the pressure and the saturation of its cell (phase rates). Recalculating all the derivatives at every iteration is expensive, so they are sometimes taken from earlier iterations or calculated approximately. The convergence is then slower than quadratic.

Connections to other topics

Two-phase flow equations — the two-phase cell equations, upstream weighting and the step of waterflooding the strip with which the examples of the lesson begin.
Pressure equation for single-phase flow — explicit and implicit schemes for pressure, the largest step and the Thomas algorithm.
Finite differences: from equations to computation — stability, the analysis of harmonics, and the difference between the convergence of a scheme and the convergence of iterations.
Relative permeability — the Corey model from which the fractional flow of water and its slope are calculated.
Solving systems of linear equations in simulators — how the sparse linear systems at each Newton iteration are solved.
Numerical dispersion — smearing of the front on the grid and with large implicit steps.
Types of reservoir simulation model — the number of unknowns per cell for different types of models, and the cost of a calculation.
Grids in reservoir simulation models — cell size, local grid refinement and calculation on streamlines.
Well model in a grid block — operating modes of wells and their link with the cell.
Buckley–Leverett theory and the Welge method — the exact solution of one-dimensional displacement for testing schemes.

Check yourself

1
What is calculated from the values at the start of the step in IMPES?
1. A.The phase mobilities in the coefficients of the equations
2. B.The water flows through the faces
3. C.Nothing: all quantities are taken at the end of the step
4. D.The pressures of all cells
Show the answer
Correct answer: A. The phase mobilities in the coefficients of the equations
IMPES “freezes” the mobilities at the saturations at the start of the step. With them the pressure equation is linear, and the pressures at the end of the step are found implicitly, by one system. The water flows are calculated from the new pressures but with the old mobilities, and the saturations are recalculated explicitly from them.
2
A cell next to a production well has a pore volume of 1200 m³, 80 m³/d of liquid at reservoir conditions pass through it, and the largest slope of the fractional flow curve is 4.0. What is the largest step that IMPES allows, and how will it change if the cell is divided in plan into four equal parts and the whole inflow goes through one of them?
Show the answer
By the Courant condition, $\Delta t \le V_p/(Q\,f_w') = 1200/(80\cdot4.0) = 3.75$ d. A quarter of the cell has a pore volume of 300 m³ at the same flow: $300/320 \approx 0.94$ d, about 22.5 h, which is four times smaller. So local grid refinement near wells hits the IMPES scheme hard and hardly bothers the fully implicit one.
3
The residuals at the iterations of a step are $3\cdot10^{-2}$ , $2\cdot10^{-3}$ , $1\cdot10^{-5}$ and $4\cdot10^{-10}$ . How do the iterations converge?
1. A.They do not converge
2. B.Linearly: the residual falls by a constant factor
3. C.This cannot be said from the residuals
4. D.Quadratically: the exponent roughly doubles
Show the answer
Correct answer: D. Quadratically: the exponent roughly doubles
The ratios of neighbouring residuals grow: 15, 200, 25,000. The exponent (the decimal logarithm of the residual with the minus sign) goes 1.5, 2.7, 5.0, 9.4, so it roughly doubles: each residual is of the order of the square of the previous one. This is how Newton's method converges near the solution. With linear convergence, as for simple iterations, the ratio of neighbouring residuals would be constant.
4
Why may Newton's method for the saturation equation “swing” and fail to converge with a large step, and what does the simulator do in this case?
Show the answer
The fractional flow curve is S-shaped: at low saturations it is flat, and near the inflection point it is steep. A tangent drawn on the flat part crosses zero far beyond the root, and from there the next tangent throws the iteration back. In the example of the lesson the approximations swing between 0.24 and 0.75. The simulator limits the corrections (the change of saturation per iteration, the crossing of the inflection point, a smaller correction when the residual grows). If the iterations still do not converge in the given number, it discards the step and repeats it with a smaller one.
5
Over a step of 20 d the largest change of water saturation was 0.12, with a target of 0.05. What will the next step be by the proportional rule?
1. A.20 d
2. B.About 8 d
3. C.40 d
4. D.48 d
Show the answer
Correct answer: B. About 8 d
The step is scaled by the ratio of the target change to the actual one: $20\cdot0.05/0.12 \approx 8.3$ d. We would get 48 d with the inverse ratio, and 40 d with a simple doubling. If the change exceeds not only the target but also the largest allowed value, the step is not accepted and is repeated with a smaller $\Delta t$ .
6
Which cells does the adaptive implicit scheme treat implicitly, and why are there usually few such cells?
Show the answer
Those where the Courant number at the proposed step is above a threshold: near wells, where a large flow passes through a small pore volume, and in the front zone, where the slope of the fractional flow curve is large, including the cells that the front will reach during the step. In the rest of the formation the saturations change slowly, the Courant number is small, and the explicit recalculation is stable and cheaper. The fronts and the surroundings of wells take a small part of the model, so the system of equations is much smaller than in the fully implicit method.

Updated: 12 October 2026