Numerical dispersion

Quick refresher

The upstream scheme replaces the derivative with respect to $x$ (the slope along the flow) by a difference with the cell upstream (the cell the flow comes from). A Taylor series (a way to write a function as a sum of terms with its derivatives) shows that the difference equation is the original equation plus a diffusion term. So the scheme in fact solves an equation with numerical diffusion, $D_{num}(S) = (vf_w'\Delta x/2)(1 \mp C)$ .
Numerical diffusion is proportional to the cell size and to the speed $vf_w'$ at which a saturation moves. ( $f_w'$ is the slope of the fractional flow curve, which gives the water share of the flowing fluid.) Here $C$ is the Courant number: the fraction of a cell that the saturation travels in one step. An explicit scheme with the step at the stability limit (the largest step that keeps the errors from growing) smears less: the factor is $1 - C$ . An implicit scheme with a large step smears more: the factor is $1 + C$ .
The front of oil displaced by water sharpens itself. The saturations that numerical diffusion pushes ahead of the front hardly move, and the front catches up with them. So on a grid the front does not get wider and wider, as a tracer front does. Its width stays constant: about two cells with an explicit scheme, and more with a large implicit step.
The results: first water (the first arrival of water at the producer) comes earlier in the model, by about the time the front needs to cross one or two cells. The water cut (the share of water in the produced liquid) rises smoothly instead of jumping. The oil produced for the same injection is too low. The errors fall roughly in proportion to the cell size.
The physical smearing of the front in the reservoir, caused by capillary forces (the pull of surface tension in the narrow pores), is from tens of centimetres to metres. The numerical smearing on 50-m cells is about a hundred metres. The width of the front in a model is a property of the grid, not of the rock.
We reduce numerical dispersion by making the cells smaller in the direction of flow, by high-order schemes with limiters (TVD schemes, which do not create false oscillations), and by streamline simulation (carrying the saturation along the paths of the fluid). In coarse models we account for it with pseudo functions (curves adjusted for coarse cells).

Key relationships

\frac{\partial S}{\partial t} + v\,\frac{\partial f_w}{\partial x} = \frac{\partial}{\partial x}\!\left(D_{num}\,\frac{\partial S}{\partial x}\right)

the modified equation: what the upstream scheme really solves

D_{num}(S) = \frac{v f_w'(S)\,\Delta x}{2}\,\bigl(1 \mp C(S)\bigr),\qquad C(S) = \frac{v f_w'(S)\,\Delta t}{\Delta x}

the numerical diffusion coefficient: “−” for the explicit scheme, “+” for the implicit one

S_{i+1/2} = S_i + \tfrac12\,\Psi(r_i)\,(S_i - S_{i-1}),\qquad r_i = \frac{S_{i+1} - S_i}{S_i - S_{i-1}}

the saturation at the face in a TVD scheme; $\Psi$ is the limiter

R_\infty \approx R_{4N} + \frac{R_{4N} - R_{2N}}{2^{p} - 1}

extrapolation of the result from three grids (Grids in reservoir simulation models)

Common pitfalls

Taking a smooth rise of the water cut in the model for a sign of heterogeneity of the reservoir, without checking the grid.
Thinking that an implicit scheme is more accurate than an explicit one because it is stable at any step. A large implicit step smears the front more.
“Correcting” an early breakthrough on a coarse grid by changing the relative permeability curves or the permeability. The numerical error hides in the reservoir parameters and moves into the forecast.
Expecting that cells half as large cost only twice as much. With an explicit scheme the time step also gets smaller, so a one-dimensional calculation costs four times more, and a three-dimensional one 16 times more.
Using a high-order scheme without a limiter. False waves (oscillations) appear at the front, and some saturations fall outside the physical limits.

After this lesson you can

derive the modified equation of the upstream scheme for the saturation equation (the equation that the scheme really solves), and find the numerical diffusion coefficient, which depends on the saturation, the cell size and the time step
explain why a front of oil displaced by water is smeared over a few cells on a grid but does not keep spreading with time, unlike the front of a tracer (a marker substance carried by the flow)
estimate how the smearing and the breakthrough time depend on the cell size, the time step and the choice between explicit schemes (new values calculated directly from the old ones) and implicit schemes (new values found by solving equations together)
compare a grid calculation with the exact Buckley–Leverett solution (the time of first water, that is, when water first reaches the producer, and the oil produced up to breakthrough and after it), and find the observed order of convergence (how fast the error falls when the cells get smaller)
describe the grid orientation effect (the result changes with the direction of the flow relative to the grid lines) and the ways to reduce numerical dispersion: finer cells, high-order schemes with limiters (more accurate schemes that avoid false waves at the front), pseudo functions (relative permeability curves adjusted for coarse cells)

What smears the front in a model

The lesson on grids shows numerical dispersion on the simplest example, the transport of a tracer (a marked substance that the flow carries without changing it). A substance that enters a cell at once “dissolves” in its whole volume. At the next step part of it passes into the neighbouring cell, and a sharp front spreads out. Its smearing grows as $\sqrt{\Delta x\,x_f}$ , where $x_f$ is the distance the front has travelled. Here we study waterflooding. We see how numerical diffusion interacts with the Buckley–Leverett jump of saturation (the sharp step in water saturation at the front, Buckley–Leverett theory and the Welge method), what it depends on, and what to do about it.

The smearing of a front can also be physical. Capillary forces stretch the front over about one capillary length, the scale of the smearing that they cause (capillary term). Small-scale heterogeneity, such as laminae (very thin layers) and lenses (small lens-shaped bodies of different rock), mixes water and oil over centimetres and metres. In miscible displacement (the injected fluid mixes with the oil) there is also molecular and mechanical dispersion: mixing by the movement of molecules and by the winding paths of the flow in the pores. For a model to show the true width of the front, its cells must be smaller than the physical zone of smearing, that is, centimetres or metres. For a model of a core sample such a grid is normal. For a field, where wells are hundreds of metres apart, it is out of reach. There the numerical smearing is almost always larger than the physical one.

Physical and numerical smearing of a displacement front

Cause: Capillary forces
What it depends on: permeability, capillary pressure difference, oil viscosity, flow velocity (Darcy velocity)
Scale: in the reservoir, from tens of centimetres to metres

Cause: Small-scale heterogeneity
What it depends on: layering, lenses, their contrast and size
Scale: centimetres to metres; larger heterogeneity is described by the geological model

Cause: Numerical dispersion, waterflooding
What it depends on: cell size, time step, scheme
Scale: about two cells with an explicit scheme, more with a large implicit step

Cause: Numerical dispersion, tracer
What it depends on: cell size and distance travelled
Scale: grows as $\sqrt{\Delta x\,x_f}$

Cause	What it depends on	Scale
Capillary forces	permeability, capillary pressure difference, oil viscosity, flow velocity (Darcy velocity)	in the reservoir, from tens of centimetres to metres
Small-scale heterogeneity	layering, lenses, their contrast and size	centimetres to metres; larger heterogeneity is described by the geological model
Numerical dispersion, waterflooding	cell size, time step, scheme	about two cells with an explicit scheme, more with a large implicit step
Numerical dispersion, tracer	cell size and distance travelled	grows as $\sqrt{\Delta x\,x_f}$

The scales are approximate. The calculations for an element of a line-drive pattern (rows of injection wells and rows of production wells) are in the sections below.

The modified equation of the upstream scheme

Take the one-dimensional saturation equation without capillary pressure and gravity (water balance in a layer). Divide it by $\phi A$ (the porosity times the cross-section area): $\partial S/\partial t + v\,\partial f_w/\partial x = 0$ . Here $v = q_t/(\phi A)$ is the average velocity of the fluid in the pores, and $q_t$ is the total flow rate. In the model the reservoir is divided into cells of size $\Delta x$ . With upstream weighting, the water fraction at a face is taken from the cell the flow comes from. The explicit scheme is the saturation step of the IMPES method (implicit pressure, explicit saturation) for a strip of cells, from the lesson on the IMPES and fully implicit methods:

Explicit upstream scheme for the saturation equation

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

where:

$S_i^{n}$ water saturation of cell $i$ at time $t^n$ , —
$v$ average velocity of the fluid in the pores, $q_t/(\phi A)$ , m/d
$\Delta x,\ \Delta t$ cell size in the flow direction and time step, m, d
$f_w$ fractional flow of water (the share of water in the flowing fluid); for the first cell, instead of $f_w(S_0)$ we use the injected water, $f_w = 1$ , —

Units: any consistent units: $v$ in m/d, $\Delta x$ in m, $\Delta t$ in d

Holds when

The flow is one-dimensional, in the $+x$ direction. The phases and the rock are incompressible, and the cells are equal.
$v\Delta t/\Delta x = \Delta t\,q_t/V_p$ is the same factor $K$ as in the lesson on the IMPES method. In the implicit scheme the water fractions are taken at time $t^{n+1}$ .

A Taylor series shows how far this difference equation is from the original one (truncation error). (A Taylor series writes a smooth function near a point as a sum of terms with its derivatives. The truncation error is what we lose by cutting that sum short.) The forward difference in time and the backward difference in space are of first order: their error is proportional to the step. The main terms of the error are second derivatives. We write the second derivative in time through the one in space, by differentiating the equation itself (the advanced section shows how). The result: apart from small terms of second order, the scheme solves a different equation. It is the modified equation:

The modified equation and the numerical diffusion of the upstream scheme

\frac{\partial S}{\partial t} + v\,\frac{\partial f_w}{\partial x} = \frac{\partial}{\partial x}\!\left(D_{num}\,\frac{\partial S}{\partial x}\right),\qquad D_{num}(S) = \frac{v f_w'(S)\,\Delta x}{2}\,\bigl(1 \mp C(S)\bigr),\qquad C(S) = \frac{v f_w'(S)\,\Delta t}{\Delta x}

where:

$D_{num}(S)$ numerical diffusion coefficient for the water saturation $S$ (it works like a real diffusion coefficient), m²/d
$f_w'(S)$ slope of the fractional flow curve; $vf_w'$ is the speed at which the saturation $S$ moves (frontal advance equation), —
$C(S)$ Courant number for the saturation $S$ : the fraction of a cell that it travels in one step (Courant–Friedrichs–Lewy condition), —
$\mp$ the sign “−” is for the explicit scheme, “+” for the implicit one, —

Units: any consistent units: $v$ in m/d, $\Delta x$ in m, $\Delta t$ in d, $D_{num}$ in m²/d

Holds when

The upstream scheme on a uniform grid. The main error terms are included. The other terms are of the order of $\Delta x^2$ and $\Delta t^2$ .
The solution is smooth. At a jump of saturation the analysis gives only a qualitative picture, and the width of the smeared front is found by calculation (next section).

The scheme adds to the equation a diffusion that does not exist in the reservoir. Its coefficient depends on the saturation. It is largest where the $f_w$ curve is steep, that is, in the front zone. It is zero where only one phase flows. For a tracer ( $f_w = S$ , $f_w' = 1$ ) we get the coefficient from the lesson on grids. The explicit scheme carries a tracer exactly when $C = 1$ . But in waterflooding $C(S)$ is different for different saturations, so we cannot remove the diffusion by choosing the step. The IMPES stability limit is set by the largest slope of $f_w$ , and at the front the slope is smaller.

Worked exampleNumerical diffusion at the front: explicit and implicit schemes

Given (assumed data)

Element of a line-drive pattern (Buckley–Leverett theory and the Welge method): length 500 m, cross-section 300 × 12 m, $\phi$ = 0.20, injection 255 m³/d at reservoir conditions; oil 3.0 mPa·s
Front: $S_{wf}$ = 0.564, $f_w'(S_{wf})$ = 2.525; the largest slope of the $f_w$ curve is 5.222
Grid: 10 cells of 50 m

Solution

Fluid velocity in the pores: $v = 255/(300\cdot12\cdot0.20) = 0.3542$ m/d. The front speed is $0.3542\cdot2.525 = 0.8944$ m/d.
IMPES step limit: the pore volume of a cell is $36{,}000$ m³, so $\Delta t \le 36{,}000/(255\cdot5.222) = 27.0$ d. We take an explicit step of 13.5 d, half the limit.
Courant number at the front: $C = 0.8944\cdot13.5/50 = 0.2415$ . Then $D_{num} = (0.8944\cdot50/2)(1 - 0.2415) = 22.36\cdot0.7585 = 17.0$ m²/d.
Implicit scheme with a step of 135 d: $C = 2.4149$ and $D_{num} = 22.36\cdot(1 + 2.4149) = 76.4$ m²/d. This is 4.5 times more.
The length over which such diffusion smears a moving front is about $D_{num}/v_f$ : 19 m for the explicit scheme and 85 m for the implicit one. The capillary length of the same rock is about 0.6 m (capillary term).

Answer: At the front the scheme adds a diffusion of 17.0 m²/d with an explicit step of 13.5 d, and of 76.4 m²/d with an implicit step of 135 d. The numerical smearing is tens of times stronger than the capillary one.

What it means. The front in the model is smeared by the grid, not by the rock. On cells of 50 m the numerical “smearing length” is 30 to 140 times the capillary length. Even an explicit step right at the stability limit ( $C$ at the front is about 0.48) reduces the numerical diffusion at the front only by half. A large implicit step increases it several times. This is the price of stability at any step. The data in this example are assumed for illustration.

A self-sharpening front: the smearing does not grow

For a tracer, the numerical diffusion coefficient is the same at all concentrations, and the smeared front spreads more and more, as with real diffusion. For a displacement front it is different. Numerical diffusion pushes small water saturations forward, into the oil-saturated zone. But for them the slope of the $f_w$ curve is almost zero, so they move slowly. The jump moves at the speed of the chord, $v\,f_{wf}/(S_{wf} - S_{wi})$ ( $v$ times the slope of the straight line from the initial saturation to the front point on the $f_w$ curve), and it catches up with them. Two effects balance each other: numerical diffusion pushes saturations ahead of the front, and the front catches up with them. On the grid the front moves as a wave of constant shape, a few cells wide. Such a front is called self-sharpening, because it keeps its own sharpness.

CalculatedWater saturation in the element of the line drive after 0.10, 0.20 and 0.30 pore volumes (fractions of the volume of all the pores of the element) have been injected. The exact Buckley–Leverett solution is shown with dashed lines. The solid lines are a calculation on a grid of 20 cells of 25 m (values at the cell centres; explicit upstream scheme, step at half the limit). The smeared front is shifted forward, but its width does not grow.

Worked exampleWidth of the smeared front on two grids

Given (assumed data)

Element and oil: as in the previous example; the exact front is a jump from 0.22 to 0.564
Grids: 20 cells of 25 m and 50 cells of 10 m; explicit upstream scheme, step at half the limit
Measure of width: the distance between the points where the saturation has risen by 10% and by 90% of the jump: 0.254 and 0.530

Solution

The exact front after 0.10, 0.20 and 0.30 pore volumes: $x_f = W_{iD}\cdot2.525\cdot500$ = 126, 253 and 379 m.
Grid of 25 m: the front width is 46, 47 and 48 m, about 1.9 cells. The middle of the smeared front is ahead of the exact one by 11–17 m.
Grid of 10 m: 23, 23 and 20 m, about 2.2 cells. The middle is ahead by 5–9 m.
A tracer that has travelled the same 379 m on the 25-m grid would be smeared by $\sigma_x \approx \sqrt{25\cdot379} = 97$ m in each direction ( $\sigma_x$ is the standard deviation of the smeared profile, a measure of its spread). And this smearing would keep growing.

Answer: The displacement front is smeared over about two cells: 48 m on the 25-m grid and 20 m on the 10-m grid. Its width hardly changes as it advances.

What it means. The width of the smeared front is almost proportional to the cell size and does not depend on the distance travelled. Cells 2.5 times smaller give a front 2 to 2.4 times narrower. The convergence is close to first order, better than for a tracer (whose smearing falls only as the square root of the cell size). The forward shift of the smeared front also falls with the cell size. After 0.30 pore volumes, on cells of 50 and 100 m it is 22 and 31 m, with front widths of 110 and 139 m. The data in this example are assumed for illustration.

A grid against the exact solution: breakthrough and production

Now we compare the whole waterflood calculation for the element with the exact solution. The grids have 5, 10, 25 and 50 cells, and we use the explicit upstream scheme with a step of half the stability limit. Breakthrough in the model is not instant: the water fraction in the production grows from zero gradually. So by convention we say that the first water arrives when the water fraction in the outlet flow reaches 1%.

An element of a line-drive pattern on different grids, and the exact solution

Cells: 5
$\Delta x$ , m: 100
First water, $W_{iD}$: 0.293 (413 d)
$N_{pD}$ at the exact breakthrough: 0.364
$N_{pD}$ after 1 pore volume: 0.421

Cells: 10
$\Delta x$ , m: 50
First water, $W_{iD}$: 0.338 (477 d)
$N_{pD}$ at the exact breakthrough: 0.376
$N_{pD}$ after 1 pore volume: 0.428

Cells: 25
$\Delta x$ , m: 20
First water, $W_{iD}$: 0.370 (522 d)
$N_{pD}$ at the exact breakthrough: 0.386
$N_{pD}$ after 1 pore volume: 0.433

Cells: 50
$\Delta x$ , m: 10
First water, $W_{iD}$: 0.382 (539 d)
$N_{pD}$ at the exact breakthrough: 0.390
$N_{pD}$ after 1 pore volume: 0.435

Cells: exact solution
$\Delta x$ , m: —
First water, $W_{iD}$: 0.396 (559 d)
$N_{pD}$ at the exact breakthrough: 0.396
$N_{pD}$ after 1 pore volume: 0.438

Cells	$\Delta x$ , m	First water, $W_{iD}$	$N_{pD}$ at the exact breakthrough	$N_{pD}$ after 1 pore volume
5	100	0.293 (413 d)	0.364	0.421
10	50	0.338 (477 d)	0.376	0.428
25	20	0.370 (522 d)	0.386	0.433
50	10	0.382 (539 d)	0.390	0.435
exact solution	—	0.396 (559 d)	0.396	0.438

Explicit upstream scheme with the step at half the stability limit. First water means 1% of water in the outlet flow. $W_{iD}$ is the water injected and $N_{pD}$ is the oil produced, both in pore volumes of the element. The data are assumed for illustration.

CalculatedWater fraction in the outlet flow of the element against the water injected: the exact solution (a jump at 0.396 pore volumes) and calculations on grids of 5, 10 and 50 cells. The larger the cells, the earlier and the smoother the rise.

Time step and scheme

The modified equation shows that the time step acts on numerical diffusion through the Courant number. An explicit scheme with a large step smears less (the factor is $1 - C$ ), and an implicit scheme smears more (the factor is $1 + C$ ). At the front the Courant number is limited. The IMPES stability limit is set by the largest slope of $f_w$ (5.22 for the curves of the lesson). At the front the slope is 2.525. So even right at the limit, $C$ at the front is not more than $2.525/5.222 \approx 0.48$ . The fully implicit scheme is stable at any step (IMPES and fully implicit schemes), and a simulator willingly takes steps tens of times larger, at the price of more smearing.

An element on 20 cells of 25 m: time step and scheme

Scheme and step: explicit, step at half the limit
$\Delta t$ , d: 6.8
First water, $W_{iD}$: 0.364
$N_{pD}$ at the exact breakthrough: 0.384
Front width, cells: 1.9

Scheme and step: explicit, step at the limit
$\Delta t$ , d: 12.8
First water, $W_{iD}$: 0.366
$N_{pD}$ at the exact breakthrough: 0.387
Front width, cells: 1.7

Scheme and step: implicit, step equal to the explicit limit
$\Delta t$ , d: 13.5
First water, $W_{iD}$: 0.347
$N_{pD}$ at the exact breakthrough: 0.376
Front width, cells: 3.6

Scheme and step: implicit, step 5 times larger
$\Delta t$ , d: 67.2
First water, $W_{iD}$: 0.293
$N_{pD}$ at the exact breakthrough: 0.359
Front width, cells: 6.9

Scheme and step: implicit, step 10 times larger
$\Delta t$ , d: 134.5
First water, $W_{iD}$: 0.286
$N_{pD}$ at the exact breakthrough: 0.343
Front width, cells: 9.4

Scheme and step	$\Delta t$ , d	First water, $W_{iD}$	$N_{pD}$ at the exact breakthrough	Front width, cells
explicit, step at half the limit	6.8	0.364	0.384	1.9
explicit, step at the limit	12.8	0.366	0.387	1.7
implicit, step equal to the explicit limit	13.5	0.347	0.376	3.6
implicit, step 5 times larger	67.2	0.293	0.359	6.9
implicit, step 10 times larger	134.5	0.286	0.343	9.4

Exact solution: breakthrough after 0.396 pore volumes, and oil production up to it 0.396. The front width is measured between the 10% and 90% levels of the jump, after 0.30 pore volumes have been injected. The data are assumed for illustration.

Worked exampleWhat a large implicit step costs

Given (assumed data)

Grid: 20 cells of 25 m; the element and the oil are the same as in the examples above
Steps: explicit, at the limit (12.8 d), and implicit, 10 times the limit (134 d)

Solution

Courant number at the front: $0.1818\cdot2.525 = 0.459$ for the explicit step and $1.905\cdot2.525 = 4.81$ for the implicit one ( $K = \Delta t\,q_t/V_p$ of the cell).
Factors of numerical diffusion at the front: $1 - 0.459 = 0.541$ and $1 + 4.81 = 5.81$ . With the implicit step it is 11 times larger.
Calculation: the front width is 1.7 and 9.4 cells; first water comes after 0.366 and 0.286 pore volumes, against 0.396 for the exact breakthrough.

Answer: The implicit step ten times the limit widened the front from 1.7 to 9.4 cells and moved the first water earlier, from 0.366 to 0.286 pore volumes.

What it means. The factor $1 + C$ predicts that the smearing grows 11 times. The calculation gave 5.5 times. The modified equation is an estimate for a smooth solution, and at a jump the higher-order terms (the smaller terms that it leaves out) also matter. The conclusion is the same: a large step of the fully implicit scheme saves calculation time, but it smears the front more than a calculation with smaller cells does. So we limit the step by a target change of saturation per step (choice of the step), and we keep the step small in the front zone and near wells. The data in this example are assumed for illustration.

CalculatedWater fraction in the outlet flow of the element on a grid of 10 cells of 50 m: the explicit upstream scheme (step at half the limit), the fully implicit scheme with a step 5 times the limit, and the TVD scheme with a van Leer limiter, compared with the exact solution.

A two-dimensional grid: the orientation effect and vertical smearing

In two and three dimensions the same analysis gives numerical diffusion along each axis of the grid. A front that runs along a row of cells is smeared only along the flow. A front that runs along a diagonal is smeared across the flow too, and in a different way from the smearing along an axis. So the same model gives different results for wells that lie on a grid line and for wells that lie on a diagonal. This is the grid orientation effect (more). It is strongest at an unfavourable mobility ratio, when the displacing fluid moves much more easily than the oil: gas displacing oil, or water displacing viscous oil. Such a front is unstable: it forms “fingers” (narrow tongues of the displacing fluid), and only numerical diffusion smooths them. The scheme then directs the fingers along the grid lines. Several things help: nine-point schemes (a cell also exchanges fluid with its diagonal neighbours), unstructured grids (cells of different shapes that do not form rows), finer cells, and a simple check, which is a calculation on a grid rotated by 45°.

Numerical diffusion in the vertical direction mixes the saturations of neighbouring layers of the model. It acts as an artificial crossflow between thin layers (flow of fluid from one layer to another). It also smears a tongue of water at the bottom of the formation, or of gas at the top, when the tongue is thinner than one layer of cells. In a layered reservoir this can replace the real crossflow (crossflow and vertical equilibrium) and distort the role of gravity (Segregated flow and gravity displacement). So waterflood models need tens of layers in the vertical direction wherever the vertical cross-section of the reservoir calls for it.

How numerical dispersion is reduced

Making the cells smaller in the flow direction (grid refinement). This is reliable but expensive: with an explicit scheme, cells half as large also need half the time step. We make the cells smaller only where the fronts move and where the wells are (local grid refinement, LGR). We also check the convergence with the grid, that is, that the result stops changing as the cells get smaller (how many cells are needed).
High-order schemes (schemes whose error falls faster than that of the upstream scheme when the cells get smaller) take the saturation at the face not from one cell upstream but from the slope of the profile inside it. Simple second-order schemes give oscillations at the front and saturations outside the physical limits. A linear scheme that creates no new peaks or dips can only be of first order. So we limit the slope. This is how TVD schemes (total variation diminishing: they create no new oscillations) with limiters (functions that cut the slope near a jump) are built.
Two-point upstream weighting. The mobility at the face (how easily the phase flows through it) is extrapolated from the two cells upstream: their trend is continued to the face. Simulators use it as a cheap improvement, also with a limiter.
Transport along streamlines. The saturation is carried along the streamlines (the paths of the fluid) on a very fine one-dimensional grid (streamlines).
Pseudo functions in coarse models are chosen so that coarse cells reproduce the calculation on a fine grid. At the same time they compensate for numerical dispersion (Two-phase upscaling: pseudo-functions).

Saturation at the face in a TVD scheme with a limiter

S_{i+1/2} = S_i + \tfrac12\,\Psi(r_i)\,(S_i - S_{i-1}),\qquad r_i = \frac{S_{i+1} - S_i}{S_i - S_{i-1}},\qquad \Psi_{\text{vL}}(r) = \frac{r + |r|}{1 + |r|}

where:

$S_{i+1/2}$ saturation at the face between cells $i$ and $i + 1$ , used to find the water fraction in the flow between them (the flow goes from $i$ to $i + 1$ ), —
$r_i$ ratio of the saturation differences downstream and upstream, a measure of how smooth the profile is, —
$\Psi(r)$ slope limiter; $\Psi_{\text{vL}}$ is the van Leer limiter, —

Units: dimensionless

Holds when

The flow is in the $+x$ direction. With $\Psi = 0$ we get the upstream scheme, and with $\Psi = 1$ a linear second-order reconstruction (a straight-line profile inside each cell).
The limiter lies in the TVD region, $0 \le \Psi(r) \le \min(2r,\ 2)$ , and at extremes (peaks and dips of the profile, $r \le 0$ ) $\Psi = 0$ .
The time step has two stages (Heun's scheme: a trial step and then a correction), and the Courant number is at most 0.5.

Where the profile is smooth ( $r \approx 1$ ), the scheme takes the slope and works as a second-order scheme. At a jump and at extremes the limiter reduces the slope to zero, and the scheme turns into the upstream scheme. The total variation $\sum|S_{i+1} - S_i|$ (the sum of the saturation changes from cell to cell) does not grow, so no new oscillations appear.

Terms and symbols

Numerical dispersion (diffusion): The smearing of fronts in a model. The difference scheme causes it: the finite size of the cells, the time step and upstream weighting.
Modified equation: The equation that the difference scheme follows more closely than the original one. It is the original equation plus the main terms of the truncation error.
$D_{num}$: The numerical diffusion coefficient. For the upstream scheme it is $(vf_w'\Delta x/2)(1 \mp C)$ .
Self-sharpening front: A front toward which the characteristics converge (the characteristics are the lines along which each saturation value travels). The smearing forward is compensated, and on the grid the front keeps a width of a few cells.
TVD scheme: A scheme that does not increase the total variation of the profile. It is second order where the solution is smooth, and it becomes the upstream scheme at jumps and extremes.
Limiter: A function $\Psi(r)$ that reduces the slope of the reconstruction at jumps. Examples are van Leer, minmod and others.
Grid orientation effect: The dependence of the result on the direction of the flow relative to the grid axes. It is stronger at an unfavourable mobility ratio.

Common mistakes and misconceptions

Wrong

“An implicit scheme is stable at any step, so it is also more accurate.”

In fact

Stability is not accuracy. A large implicit step increases numerical diffusion by the factor $1 + C$ . In the example of the lesson, a step ten times larger widened the front from 1.7 to 9 cells.

Wrong

“A waterflood front on a coarse grid spreads more and more, like a tracer front.”

In fact

A displacement front sharpens itself and keeps a width of about two cells. It is the front of a tracer or of a miscible agent that spreads more and more. But two cells of 50–100 m are still 100–200 m.

Wrong

“Breakthrough in the model comes earlier than in the field, so the reservoir has a high-permeability thin layer.”

In fact

First check the grid. On cells of 50 m the first water in the example of the lesson arrives 15% earlier than in the exact solution, with no heterogeneity at all.

Wrong

“We will use a second-order scheme, and the smearing will disappear.”

In fact

Without a limiter, a second-order scheme gives oscillations and unphysical saturations at the front. We need a scheme with a limiter (TVD), and even it smears the toe of the front over a cell or two.

In depth: advanced detail

Derivation of the modified equation

We assume that the solution of the scheme is a smooth function $S(x, t)$ that satisfies the difference equation at the grid nodes. We substitute it into the explicit scheme and expand it in a Taylor series around the node $(x_i, t^n)$ : $S_i^{n+1} = S + \Delta t\,S_t + \tfrac12\Delta t^2 S_{tt} + \ldots$ and $f_{i-1} = f - \Delta x\,f_x + \tfrac12\Delta x^2 f_{xx} - \ldots$ . The subscripts $t$ and $x$ mean partial derivatives. The scheme becomes $S_t + v f_x = \tfrac12 v\Delta x\,f_{xx} - \tfrac12\Delta t\,S_{tt} + O(\Delta x^2, \Delta t^2)$ .

We express the second derivative in time through the space derivatives. We use the fact that to leading order (keeping only the biggest terms) $S_t = -v f_x = -v f'S_x$ . Then $S_{tt} = -v\,(f_x)_t = -v\,(f'S_t)_x = v^2(f'f_x)_x = v^2\left(f'^2S_x\right)_x$ . Also $f_{xx} = (f'S_x)_x$ . Putting these together we get $S_t + vf_x = \left[\tfrac12 v\Delta x\,f'\left(1 - v f'\Delta t/\Delta x\right)S_x\right]_x$ . This is the modified equation of the main section, with the sign “−”.

In the implicit scheme the flows are taken at time $t^{n+1}$ , and it is easier to expand around $(x_i, t^{n+1})$ . Then the difference in time gives $S_t - \tfrac12\Delta t\,S_{tt}$ , and the sign of the term with the step changes to “+”. This gives the factor $1 + C$ : the implicit scheme adds a “temporal” numerical diffusion to the spatial one. For linear transport, the analysis of harmonics (waves of different lengths into which an error can be split) in the lesson on the IMPES schemes gave the same result. The amplification factor $|G| < 1$ (the factor by which a harmonic is multiplied in one step) of the short harmonics is their damping, that is, diffusion.

Order of convergence and extrapolation from grids

Worked exampleFirst water on three grids, and extrapolation

Given (assumed data)

Calculations: an element of a line-drive pattern on 10, 20 and 40 cells, explicit upstream scheme, step at half the limit
First water (1% in the outlet flow): $R_{10}$ = 0.3381, $R_{20}$ = 0.3643, $R_{40}$ = 0.3785 pore volumes

Solution

Changes when the cells get smaller: $0.3643 - 0.3381 = 0.0262$ and $0.3785 - 0.3643 = 0.0142$ .
Observed order: $p = \log_2(0.0262/0.0142) = \log_2 1.845 = 0.88$ .
Extrapolation by the formula of the lesson on grids: $R_\infty \approx 0.3785 + 0.0142/(2^{0.88} - 1) = 0.3785 + 0.0169 = 0.395$ .

Answer: Extrapolation gives 0.395 pore volumes, against the exact value of 0.396. The order of convergence is about 0.88, slightly below one.

What it means. For problems with jumps, the order of convergence of a first-order scheme is usually a little less than one, and it is not constant. For the next pair of grids (40 and 80 cells) it comes out at about 0.8. Extrapolation from three grids (estimating what the result would be on an ever finer grid) gives the error of the working grid even without an exact solution. This is how we use it for two- and three-dimensional models, for which there is no exact solution. The data in this example are assumed for illustration.

Why high order needs limiters

The Lax–Wendroff scheme and other linear second-order schemes remove the main term $\tfrac12 v\Delta x\,S_{xx}$ . But the next error term contains the third derivative, and it is not diffusive but dispersive: harmonics of different lengths travel at different speeds, and “ripples” appear behind the jump. Godunov's theorem says that a linear scheme that keeps the profile monotone (with no new peaks or dips) cannot be more accurate than first order. The only way round it is a nonlinear scheme, in which the slope of the reconstruction depends on the solution itself.

A TVD (total variation diminishing) scheme must not let the total variation $\sum_i|S_{i+1} - S_i|$ grow from step to step. For this, the limiter $\Psi(r)$ must lie in the region $0 \le \Psi \le \min(2r, 2)$ . The minmod limiter ( $\Psi = \max[0, \min(1, r)]$ ) is the most cautious, superbee is the “sharpest”, and van Leer is between them. A saturation profile is smooth ahead of the front and behind it, and it jumps at the front. For such a profile, a TVD scheme gives second order almost everywhere and no oscillations at the jump.

Numerical diffusion in place of physical diffusion

Sometimes numerical diffusion is used on purpose. Suppose the physical dispersion is known, for example the coefficient of longitudinal dispersion (mixing along the flow) in miscible displacement, and the numerical dispersion is certainly larger on the working grid. Then we choose the cell size so that $v\Delta x/2$ is close to the physical coefficient. A coarse model then reproduces the width of the mixing zone. We check this against a calculation on a fine grid, where the diffusion is set explicitly. For a waterflood front this trick is not reliable. Numerical diffusion depends on the saturation differently from capillary diffusion, and it depends on the time step. So we leave its compensation to pseudo functions (Two-phase upscaling: pseudo-functions).

Connections to other topics

Grids in reservoir simulation models — numerical dispersion of a tracer, the grid orientation effect, convergence with the grid, and streamlines.
Buckley–Leverett theory and the Welge method — the exact solution that grid calculations are compared with.
IMPES and fully implicit schemes — time-step schemes, the Courant condition and the choice of the step.
Finite differences: from equations to computation — truncation error, order, and the analysis of harmonics.
Two-phase flow equations — the saturation equation and upstream weighting.
Fractional flow function and mobility ratio — the slope of the fractional flow curve and the capillary smearing of the front.
Displacement in heterogeneous reservoirs — how to tell numerical smearing from heterogeneity of the reservoir.
Two-phase upscaling: pseudo-functions — pseudo functions that compensate for the smearing in coarse models.
History matching — why numerical error must not be hidden in reservoir parameters.

Check yourself

1
How does the numerical diffusion of the upstream scheme at the front change if the cells are made half as large and the Courant number is kept the same?
1. A.It doubles
2. B.It does not change
3. C.It falls to a quarter
4. D.It halves
Show the answer
Correct answer: D. It halves
$D_{num} = (vf_w'\Delta x/2)(1 \mp C)$ . At the same Courant number (the step is reduced together with the cell) the coefficient is proportional to $\Delta x$ . The width of the smeared waterflood front, about two cells, also halves.
2
A grid has cells of 40 m, the fluid velocity in the pores is 0.5 m/d, and the slope of the fractional flow curve at the front is 2.0. Find the numerical diffusion at the front for an explicit step of 10 d and for an implicit step of 100 d.
Show the answer
The front speed is $0.5\cdot2.0 = 1.0$ m/d. Explicit step: $C = 1.0\cdot10/40 = 0.25$ , and $D = (1.0\cdot40/2)(1 - 0.25) = 20\cdot0.75 = 15$ m²/d. Implicit step: $C = 2.5$ , and $D = 20\cdot3.5 = 70$ m²/d. This is almost five times more, for a step ten times larger.
3
Why is a waterflood front smeared over a constant width on a grid, while a tracer front spreads more and more?
Show the answer
The saturations that numerical diffusion pushes ahead of a displacement front are small, and the slope of the fractional flow curve at them is almost zero. They hardly move, and the jump moves at the speed of the chord and catches up with them. Smearing and sharpening balance each other, and the profile moves as a wave of constant shape. For a tracer, all concentrations move at the same speed. Nothing brings back what was pushed forward, and the smearing grows as the square root of the distance travelled.
4
A model with 100-m cells shows breakthrough 25% earlier than the field data. What is it reasonable to do first?
1. A.Increase the time step of the implicit scheme
2. B.Repeat the calculation on a finer grid and estimate how much of the early breakthrough is numerical dispersion
3. C.Change the water relative permeability curves
4. D.Reduce the permeability of the high-permeability thin layer, to delay the water
Show the answer
Correct answer: B. Repeat the calculation on a finer grid and estimate how much of the early breakthrough is numerical dispersion
In the one-dimensional example of the lesson, a grid of 5 cells of 100 m gives the first water 26% earlier than the exact solution, with no heterogeneity at all. Until the numerical error is estimated, adjusting the permeability or the relative permeability would hide it in the reservoir parameters. A large implicit step would only increase the smearing.
5
How is a TVD scheme better than a second-order scheme without a limiter?
1. A.It is more accurate on smooth parts
2. B.It creates no oscillations or unphysical saturations at a jump, and it stays second order where the profile is smooth
3. C.It is stable at any time step
4. D.It removes numerical dispersion completely
Show the answer
Correct answer: B. It creates no oscillations or unphysical saturations at a jump, and it stays second order where the profile is smooth
A linear second-order scheme gives “ripples” and saturations outside the limits at a jump, and a monotone linear scheme can be only of first order (Godunov's theorem). A TVD scheme is nonlinear: the limiter reduces the slope at jumps and extremes, and second order stays on smooth parts. It does not become stable at any step: in the explicit version its step is limited by the Courant condition even more strictly than the step of the upstream scheme. And the toe of its front is also smeared over a cell or two.

Updated: 15 October 2026