Pressure equation for single-phase flow

Quick refresher

The mass balance of a cell is an exact statement for a cell of any size and for any step. Approximations come later: when we express the flow through a face by the pressures at the centres of the neighbouring cells, and when we choose the moment whose pressures give the flows during the step.
Darcy's law and compressibility turn the balance into an equation for the pressure alone. The equation is nonlinear (its coefficients change with the unknown), because the density, the viscosity and the porosity, and sometimes the permeability, depend on the pressure we are looking for.
For a slightly compressible fluid (one whose density changes only a little with pressure) the equation becomes linear, but the heterogeneity stays: the permeability remains under the derivative. The diffusivity equation is the special case of a uniform formation.
The equation of a cell has two “parameters”: the storage capacity $S_i = V_{p,i}c_t$ , which is how much fluid the cell gives up when the pressure falls by 1 MPa, and the coefficients of the links with the neighbours, $a_{ij} = 0.0853\,T_{ij}/\mu$ .
The explicit scheme calculates the flows from the pressures at the start of the step. It is certain to be stable (errors do not grow) if the step is not larger than the characteristic time of the “fastest” cell, $\tau_i = S_i/\sum_j a_{ij}$ . With a noticeably larger step it is unstable (the errors grow from step to step). The implicit scheme calculates the flows from the pressures at the end of the step. It is stable at any step, but at every step we must solve a system of equations.
Both schemes conserve mass exactly. Stability and accuracy are separate properties: an implicit calculation with a large step does not “fall apart”, but it smooths the result and lags behind.

Key relationships

\phi\rho\,c_t\,\frac{\partial p}{\partial t} = \nabla\cdot\left[\frac{\rho}{\mu}\,\mathbf{K}\left(\nabla p - \rho g\,\nabla D\right)\right] - \rho\,\tilde q

pressure equation for single-phase flow (SI); $D$ is the depth, $\tilde q$ is the withdrawal per unit volume of the formation

S_i\,\frac{p_i^{n+1} - p_i^{n}}{\Delta t} = \sum_j a_{ij}\,(p_j - p_i) - q_i

equation of a cell: the pressures on the right are those at the start of the step (explicit scheme) or at the end of the step (implicit scheme)

S_i = V_{p,i}\,c_t,\qquad a_{ij} = \frac{0.0853\,T_{ij}}{\mu}

m³/MPa and m³/(d·MPa) for $T$ in mD·m and $\mu$ in mPa·s

\Delta t \le \tau_i = \frac{S_i}{\sum_j a_{ij}}\quad\text{for all cells}

a sufficient condition for the stability of the explicit scheme (if it holds, the scheme is stable); for a uniform one-dimensional grid $\tau = \Delta x^2/(2\eta)$

Common pitfalls

Taking a stable calculation for an accurate one: an implicit scheme with a large step smooths the change of pressure and lags behind.
Believing that a mass balance that closes confirms the pressures: the balance also holds in an unstable calculation.
Averaging the permeability of neighbouring cells arithmetically: for flow through a face we need the harmonic average.
Putting the compressibility of the oil alone into the storage capacity instead of the total compressibility $c_t$ , and mixing rates at reservoir and standard conditions.
Taking the pressure of the cell with a well for the bottomhole pressure.

After this lesson you can

derive the pressure equation from the mass balance of a cell (a control volume) and Darcy's law, and name which of its coefficients depend on pressure
obtain the linear form for a slightly compressible fluid in a heterogeneous formation (one whose properties change from place to place), and link it to the diffusivity equation
write the equation of a cell with transmissibilities, a well and boundary conditions, in metric oilfield units
tell the explicit and implicit time schemes apart (the flows are taken from the pressures at the start or at the end of the step), and estimate the largest time step of the explicit scheme from the characteristic times of the cells
solve by hand a problem on a grid of four cells with both schemes, and check the mass balance

What the simulator looks for: the pressure in every cell

In a single-phase problem, one fluid moves in the pores: undersaturated oil (oil above its bubble point pressure, with no free gas), water or gas. One field, the pressure, describes the state of the formation at every moment. Once we know it, we find the velocities and the flows from Darcy's law, the density of the fluid from the equation of state (a rule that gives the density from the pressure), and the bottomhole pressure from the well model. That is why the equation of single-phase flow is called the pressure equation. In a multiphase model, equations for the saturations are added to it (Two-phase flow equations). The pressure equation itself is built in the same way, and it is easiest to see all the numerical methods on the single-phase case.

In the lesson Diffusivity equation the same equation is built for a thin ring around a well in a uniform formation, and solved analytically (with a formula). A simulator needs a different form. The formation is not uniform, the geometry is complex, there are many wells, and the solution is found numerically (by calculation, step by step), for cells of finite size and time steps of finite length. Here we go by the “cell” route: from the mass balance of a cell to the differential equation, and back again to the equations for cells that the computer solves.

The mass balance of a cell

Take cell $i$ of a one-dimensional grid: a bar of formation of length $\Delta x_i$ , cross-section area $A$ and volume $V_i = A\,\Delta x_i$ . Fluid enters through the left face (we label it $i - 1/2$ ) and leaves through the right face ( $i + 1/2$ ), and a well withdraws some of it. The mass of fluid in the cell is the density multiplied by the pore volume: $\rho\phi V_i$ . A volume for which we count everything that enters it and leaves it is called a control volume.

Mass balance of a cell over a time step

V_i\left[(\rho\phi)_i^{n+1} - (\rho\phi)_i^{n}\right] = \Delta t\left[(\rho q)_{i-1/2} - (\rho q)_{i+1/2} - (\rho q_w)_i\right]

where:

$V_i$ volume of the cell (rock together with pores), $A\,\Delta x_i$ , m³
$\rho,\ \phi$ density of the fluid at reservoir conditions and porosity, kg/m³, —
$n,\ n + 1$ numbers of the moments of time: the start and the end of the step, —
$q_{i-1/2},\ q_{i+1/2}$ volumetric flow through the left and right faces, positive in the direction of the $x$ axis, m³/s
$q_{w,i}$ withdrawal by a well from the cell at reservoir conditions (negative for injection), m³/s
$\Delta t$ time step, s

Units: SI (any consistent system)

Holds when

The flows in the square brackets are averages over the step. Read this way, the equality is exact.
The flow is along one axis. In two and three dimensions, the flows through the other faces are added.
The density at a face is taken from the values in the neighbouring cells.

The left side is the change of the mass of fluid in the cell over the step. The right side is the mass that entered through the left face, minus the mass that left through the right face and the mass withdrawn by the well. This is a conservation law, with not a single physical assumption yet about how the fluid moves.

We divide the balance by $V_i\,\Delta t$ and let the cell size and the time step go to zero. The difference of the mass flows at the faces, divided by $\Delta x$ , becomes the derivative with respect to the coordinate. The change of mass over the step, divided by the step, becomes the derivative with respect to time. The withdrawal by the well becomes the withdrawal per unit volume of the formation, $\tilde q$ . In three dimensions the differences of the flows along the three axes add up to the divergence (the net outflow from a small volume), in the same way as in the lesson on the diffusivity equation.

Continuity equation with a source

\frac{\partial(\rho\phi)}{\partial t} = -\frac{\partial(\rho u)}{\partial x} - \rho\,\tilde q\quad\longrightarrow\quad \frac{\partial(\rho\phi)}{\partial t} + \nabla\cdot(\rho\,\mathbf{u}) + \rho\,\tilde q = 0

where:

$u,\ \mathbf{u}$ Darcy velocity $q/A$ along the $x$ axis (the flow rate per unit area), and its vector, m/s
$\tilde q$ withdrawal by wells per unit volume of the formation, $q_w/V$ ; it is non-zero only in cells with wells, 1/s
$\rho,\ \phi$ density of the fluid and porosity, kg/m³, —

Units: SI

Holds when

The limit as $\Delta x \to 0$ and $\Delta t \to 0$ . The properties change smoothly inside the formation.

In the second form, the first term is the rate at which mass accumulates in a unit volume of the formation, the second is the “divergence” of the mass flow out of this volume, and the third is the withdrawal by a well. They add up to zero. Like the balance of a cell, the equation is exact. It does not yet say by what law the fluid moves.

Darcy's law and compressibility: the pressure equation

To get to a single unknown, we need two more links. We express the Darcy velocity with Darcy's law in its general form, with the permeability tensor (a set of numbers that describes how permeability depends on direction) and gravity (Darcy's law in general form). We express the accumulation through the pressure. The density of the fluid and the porosity depend on the pressure. By the chain rule (the rule for the rate of change of a quantity that depends on another changing quantity), the rate of change of the product $\rho\phi$ equals $\rho\phi\,c_t$ multiplied by the rate of change of the pressure. Here $c_t = c + c_f$ is the total compressibility: the compressibility of the fluid plus the compressibility of the pores. Substitution gives:

The pressure equation for single-phase flow

\phi\rho\,c_t\,\frac{\partial p}{\partial t} = \nabla\cdot\left[\frac{\rho}{\mu}\,\mathbf{K}\left(\nabla p - \rho g\,\nabla D\right)\right] - \rho\,\tilde q,\qquad c_t = \frac{1}{\rho\phi}\,\frac{\dd(\rho\phi)}{\dd p}

where:

$p(x, y, z, t)$ pressure, the only unknown, Pa
$\mathbf{K}$ permeability tensor; if the model axes point along the principal directions, it is diagonal: $k_x$ , $k_y$ , $k_z$ , m²
$\mu$ viscosity of the fluid, Pa·s
$c_t$ total compressibility: compressibility of the fluid $c$ plus compressibility of the pores $c_f$ , 1/Pa
$g,\ D$ acceleration of gravity and the depth, measured downwards, m/s², m
$\tilde q$ withdrawal by wells per unit volume of the formation (negative for injection), 1/s

Units: SI

Holds when

One mobile fluid. Darcy's law (no inertial losses).
Isothermal flow (the temperature does not change).
The density, the viscosity, the porosity and the permeability may depend on the pressure. The equation does not simplify anything yet.

The left side is how much mass a unit volume of the formation stores when the pressure rises. The right side is how much mass the flows bring in and how much the well takes out. The minus sign before the gradient in Darcy's law is already taken into account: fluid flows towards lower potential, so a region of lower pressure is “filled” from its neighbours.

Density $\rho(p)$ appears both in the accumulation and in the flow. For a liquid it changes by fractions of a per cent, and for a gas it changes almost in proportion to the pressure.
Viscosity $\mu(p)$ of a gas rises noticeably with pressure. For oil above the bubble point pressure it changes little.
Porosity $\phi(p)$ and the total compressibility change through the compressibility of the pores.
Permeability $\mathbf{K}(p)$ in fractured and poorly cemented rocks falls when the pore pressure falls (Rock compressibility and the effect of stress on core properties).

The coefficients depend on the pressure we are looking for, so the equation is nonlinear, and it has no general analytical solution even in one dimension. The simulator does not simplify it “on paper”. It writes the nonlinear equations for the cells and solves them at every step by iterations (repeated rounds of calculation), updating the properties from the latest estimate of the pressure (see the advanced section and the lesson IMPES and fully implicit schemes).

To understand numerical schemes, and to test them on exact solutions, we need a linear form. For a slightly compressible fluid, such as oil above the bubble point pressure or water, the product of the compressibility of the fluid and the change of pressure is small ( $c\,\Delta p \ll 1$ ; estimates are in the lesson on diffusivity). Then we neglect the change of density in the flow and divide the equation by $\rho$ . We take the viscosity and the total compressibility as constant. We account for gravity by replacing the pressure by the potential $\Phi = p - \rho g D$ (the pressure corrected for depth).

The linear pressure equation for a heterogeneous formation

\phi\mu c_t\,\frac{\partial p}{\partial t} = \frac{\partial}{\partial x}\!\left(k_x\frac{\partial p}{\partial x}\right) + \frac{\partial}{\partial y}\!\left(k_y\frac{\partial p}{\partial y}\right) + \frac{\partial}{\partial z}\!\left(k_z\frac{\partial p}{\partial z}\right) - \mu\,\tilde q

where:

$k_x,\ k_y,\ k_z$ principal permeabilities; they may change from point to point, m²
$\phi,\ \mu,\ c_t$ porosity, viscosity, total compressibility (constants), —, Pa·s, 1/Pa
$\tilde q$ withdrawal per unit volume of the formation at reservoir conditions, 1/s
$p$ pressure; if gravity is taken into account, the potential $p - \rho g D$ with constant $\rho$ , Pa

Units: SI

Holds when

A slightly compressible fluid: $c\,\Delta p \ll 1$ . The change of density in the flow is not taken into account.
The viscosity and the total compressibility do not depend on pressure. The porosity changes little.
The coordinate axes point along the principal directions of permeability.

The equation is linear, but the permeability stays under the derivative, so it describes a heterogeneous and anisotropic formation. If the formation is uniform and isotropic, $k$ comes out from under the derivative, and after division by $k$ we get the diffusivity equation with $\eta = k/(\phi\mu c_t)$ .

The initial and boundary conditions in the model

As in analytical problems (conditions of the problem), the equation needs an initial condition and boundary conditions. In the model, each of them turns into a term in the equation of some cell.

Initial condition: The pressure in every cell at the initial moment. We calculate it from hydrostatic equilibrium (the balance of the weight of the fluid column): the pressure grows with depth by $\rho g$ for every metre, starting from the pressure at a reference depth. This is part of the initialisation of the model (What a reservoir simulation model is and why it is needed). If the initial pressure does not agree with the weight of the fluid column, the fluid in the model will start to flow between cells without any wells.
No-flow boundary: The outer faces of the model are closed by default: there is no connection through them, and the transmissibility is zero. We set sealing faults and the pinching out of a formation (the place where it thins to zero thickness) in the same way (no-flow boundary).
Constant pressure: The boundary cell is linked with a fictitious cell whose pressure does not change. The distance from the centre of the cell to the face is half a cell, so the transmissibility of such a connection, $T_b = 2k_iA/\Delta x_i$ , is twice as large as between two such cells. This is a rough way to describe a powerful aquifer (a large water-bearing rock next to the reservoir that keeps the pressure up). More accurate models are in the lesson Aquifer models: water influx.
Well with a given rate: The withdrawal $q_i$ is known and enters the equation of the cell as a known quantity.
Well with a given bottomhole pressure: The withdrawal is proportional to the difference between the cell pressure and the bottomhole pressure: $q_i = a_w\,(p_i - p_{wf})$ . We calculate the coefficient $a_w$ from the well index, which is the subject of the well model (Well model in a grid block). In the implicit scheme the unknown cell pressure enters the withdrawal, and the coefficient is added to the diagonal of the system of equations.

Cells instead of derivatives: the equation of a cell

Now we go the other way, from the differential equation to the equations that the computer solves. The formation is divided into cells, and the pressure of each cell is attached to its centre. Such a grid is called block-centred. We integrate the linear equation over the volume of the cell (that is, we add up its terms over the whole cell). The accumulation gives the storage capacity of the cell multiplied by the rate of change of its pressure. The divergence of the flow turns into the sum of the flows through the faces. We express the flow through a face with Darcy's law for the two halves of the neighbouring cells joined one after the other, through the transmissibility $T_{ij}$ introduced in the lesson on the simulation model. The derivatives disappear, and one algebraic equation (an equation with no derivatives) remains for every cell. It is the equation of a cell.

Equation of a cell (single-phase, slightly compressible fluid)

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

where:

$p_i^{n},\ p_i^{n+1}$ pressure in cell $i$ at the start and at the end of the step, MPa
$S_i$ storage capacity of the cell: the volume of fluid that it gives up when the pressure falls by 1 MPa, m³/MPa
$V_{p,i}$ pore volume of the cell, m³
$c_t$ total compressibility, 1/MPa
$a_{ij}$ coefficient of the link between cells $i$ and $j$ : the flow between them for a pressure difference of 1 MPa, m³/(d·MPa)
$T_{ij}$ transmissibility of the connection (geometric part): $A/\big(\Delta x_i/(2k_i) + \Delta x_j/(2k_j)\big)$ , mD·m
$q_i$ withdrawal by a well from cell $i$ at reservoir conditions (negative for injection), m³/d
$\Delta t$ time step, d

Units: metric oilfield units: $k$ in mD, lengths in m, $\mu$ in mPa·s, $c_t$ in 1/MPa, $p$ in MPa, $t$ in d, volumes at reservoir conditions

Holds when

A slightly compressible fluid, with constant viscosity and compressibility (a linear equation).
The sum runs over all cells with which cell $i$ has a common face or a non-neighbour connection.
The pressures on the right side are taken at the start of the step (explicit scheme) or at the end (implicit scheme).
If the cells lie at different depths, we put $(p_j - p_i) - \rho g\,(D_j - D_i)$ in place of $p_j - p_i$ , with $\rho g$ in MPa/m.

A cell behaves like an elastic tank. Its pressure changes at a rate equal to the net inflow from the neighbours and the well, divided by its capacity. The coefficient $a_{ij}$ is the same for the pair $(i, j)$ and for $(j, i)$ : whatever left one cell arrived in the other.

In other unit systems (1)

V_{p,i}\,c_t\,\frac{p_i^{n+1} - p_i^{n}}{\Delta t} = \sum_j \frac{T_{ij}}{\mu}\left(p_j - p_i\right) - q_i

$V_p$ in m³, $c_t$ in 1/Pa, $p$ in Pa, $\Delta t$ in s, $T_{ij}$ in m³ (permeability in m²), $\mu$ in Pa·s, $q$ in m³/s.

There is no constant: all quantities are in consistent units.

In two and three dimensions the equation is the same. A cell has four or six neighbours that share a face with it (the five-point and the seven-point schemes). In a corner-point geometry (cells that follow the dip of the layers) with displacements along faults, it has as many connections as it has common faces and non-neighbour connections (Grids in reservoir simulation models). The system of “tanks and pipes” repeats exactly the chain of capacitors and resistors from the lesson on diffusivity, only with a finite number of links.

The time step: explicit and implicit schemes

The main choice remains: the pressures of which moment should give the flows during the step. If we take the pressures at the start of the step, we get the explicit scheme. If we take the pressures at the end of the step, we get the implicit scheme. The lesson on the simulation model discussed this choice for one “tank” with a well. A grid adds something new: the cells are linked with each other.

Explicit scheme: the new pressure of a cell

p_i^{n+1} = p_i^{n} + \frac{\Delta t}{S_i}\left[\sum_j a_{ij}\left(p_j^{n} - p_i^{n}\right) - q_i\right]

where:

$p_i^{n},\ p_j^{n}$ pressures of the cell and its neighbours at the start of the step (known), MPa
$S_i,\ a_{ij},\ q_i,\ \Delta t$ as in the equation of a cell, m³/MPa, m³/(d·MPa), m³/d, d

Units: metric oilfield units, as in the equation of a cell

Holds when

The flows and the withdrawal over the whole step are calculated from the pressures at the start of the step.

Each cell is recalculated separately from known numbers, so no system of equations is needed. The price of this simplicity is that in one step a disturbance moves by exactly one cell, and with too large a step the calculation loses stability.

The largest step of the explicit scheme

\Delta t \le \min_i \tau_i,\qquad \tau_i = \frac{S_i}{\sum_j a_{ij}}

where:

$\tau_i$ characteristic time of the cell: the time in which the cell would equalise its pressure with its neighbours if the flows did not change, d
$S_i,\ a_{ij}$ storage capacity of the cell and the coefficients of its links, m³/MPa, m³/(d·MPa)

Units: metric oilfield units (any consistent units)

Holds when

A linear equation. The condition is sufficient: if it holds, the scheme is stable. Under it the new pressure of the cell (leaving the withdrawal aside) is a weighted average of the pressures of the cell itself and its neighbours at the start of the step, with non-negative weights. So the calculation cannot “jump over” the neighbours.
For a long uniform grid the condition is exact: a larger step is unstable. For a short or non-uniform grid the allowed step is a little larger (see the advanced section).

The step is limited by the “fastest” cell, the one with the smallest capacity and the most conductive links. For a uniform one-dimensional grid $\tau = \Delta x^2/(2\eta)$ : the largest step is proportional to the square of the cell size and inversely proportional to the hydraulic diffusivity (the number that says how fast pressure changes spread through the formation).

In the implicit scheme the flows and the withdrawal are calculated from the pressures at the end of the step, which are not yet known. Then the equation of each cell has three unknowns: its own pressure and the pressures of its neighbours. The equations of all cells have to be solved together.

Implicit scheme: the equation of a cell in the system

-a_{i,i-1}\,p_{i-1}^{n+1} + \left(\frac{S_i}{\Delta t} + \sum_j a_{ij}\right)p_i^{n+1} - a_{i,i+1}\,p_{i+1}^{n+1} = \frac{S_i}{\Delta t}\,p_i^{n} - q_i

where:

$p_{i-1}^{n+1},\ p_i^{n+1},\ p_{i+1}^{n+1}$ unknown pressures of the cell and of its neighbours at the end of the step, MPa
$a_{i,i\pm1}$ coefficients of the links with the neighbours (an end cell has no neighbour on one side, and the coefficient there is zero), m³/(d·MPa)
$S_i,\ q_i,\ \Delta t,\ p_i^{n}$ as in the equation of a cell, m³/MPa, m³/d, d, MPa

Units: metric oilfield units

Holds when

A one-dimensional grid. In two and three dimensions a row has five or seven unknowns.
A linear equation. In a nonlinear one the coefficients depend on the unknown pressures, and the system is solved by iterations.

The matrix of the system (the table of the coefficients of its equations) is tridiagonal: in each row only the diagonal and the two neighbouring elements are non-zero. The diagonal is strictly larger than the sum of the neighbouring coefficients, so the system can be solved at any step, and the calculation does not produce oscillations. The capacity of the cell enters the diagonal as $S_i/\Delta t$ : the larger the step, the weaker the “memory” of the start of the step.

A tridiagonal system is solved by the Thomas algorithm (also called the “sweep” method). In the forward pass, from top to bottom, we eliminate the previous unknown from each equation. In the backward pass, from bottom to top, we find the unknowns. The number of operations is proportional to the number of cells, so for a one-dimensional problem this is almost as fast as an explicit step. In two and three dimensions the matrix stays sparse (it has five or seven non-zero diagonals, and nearly all its other entries are zero), but the Thomas algorithm no longer fits. Systems of millions of equations are solved by iterative methods with preconditioning (they improve a guess step by step and use a cheap approximate solution to speed this up) (Solving systems of linear equations in simulators).

The Thomas algorithm for a system with a tridiagonal matrix

\begin{aligned} &w_i = b_i - l_i\,e_{i-1},\qquad e_i = \frac{u_i}{w_i},\qquad f_i = \frac{r_i + l_i\,f_{i-1}}{w_i}\\ &x_N = f_N,\qquad x_i = f_i + e_i\,x_{i+1}\quad (i = N - 1,\ \ldots,\ 1)\end{aligned}

where:

$x_i$ unknowns (pressures or their changes over the step), MPa
$b_i$ diagonal coefficient of equation $i$ , m³/(d·MPa)
$l_i,\ u_i$ coefficients of the left and right neighbours, taken with the opposite sign ( $a_{i,i-1}$ and $a_{i,i+1}$ ); $l_1 = 0$ , $u_N = 0$ , m³/(d·MPa)
$r_i$ right-hand side of equation $i$ , m³/d
$w_i,\ e_i,\ f_i$ auxiliary quantities of the forward pass; $e_0 = f_0 = 0$ , m³/(d·MPa), —, MPa

Units: any consistent units

Holds when

A system of the form $-l_i x_{i-1} + b_i x_i - u_i x_{i+1} = r_i$ with a dominant diagonal (as in the implicit scheme): there is no division by zero.

The forward pass turns each equation into a relation $x_i = f_i + e_i x_{i+1}$ with a single next unknown. The last equation gives $x_N$ at once, and then the rest are found in the reverse order.

Explicit and implicit schemes for the pressure equation

Property: Flows during the step
Explicit: from the pressures at the start of the step
Implicit: from the pressures at the end of the step

Property: What is solved
Explicit: each cell separately
Implicit: a system of equations for all cells

Property: Stability
Explicit: guaranteed for $\Delta t \le \tau_{\min}$
Implicit: at any step

Property: Error in time
Explicit: proportional to the step
Implicit: proportional to the step; a large step smooths the result

Property: Use
Explicit: simple and teaching calculations
Implicit: the pressure in almost all simulators

Property	Explicit	Implicit
Flows during the step	from the pressures at the start of the step	from the pressures at the end of the step
What is solved	each cell separately	a system of equations for all cells
Stability	guaranteed for $\Delta t \le \tau_{\min}$	at any step
Error in time	proportional to the step	proportional to the step; a large step smooths the result
Use	simple and teaching calculations	the pressure in almost all simulators

There are also intermediate schemes, for example flows from the average pressure of the step. Their accuracy and stability, and the general theory of finite differences, are in the lesson Finite differences: from equations to computation.

Example: a strip of formation with four cells

Now we follow both schemes on a grid that we can calculate by hand. A strip of formation 1000 m long, 200 m wide and 15 m thick is divided into four cells of 250 m each. Both ends of the strip are impermeable. The rock and the oil are the same as in the lesson Flow regimes. Cell 3 contains a clay-rich zone with a lower permeability. A well works in cell 4 at a rate of 100 m³/d at standard conditions.

SchematicA strip of formation with four cells (a sketch, not to scale): the permeabilities of the cells, the transmissibilities of the connections between them, the well in cell 4 and the impermeable ends. The data are those of the examples below.

Worked exampleTransmissibilities, capacities and the largest step of the explicit scheme

Given (assumed data)

Cells: 4 cells of $\Delta x$ = 250 m; width 200 m; thickness $h$ = 15 m
Permeability of cells 1–4: 80; 80; 20; 100 mD
Rock and oil (lesson Flow regimes): $\phi$ = 0.20; $\mu$ = 1.2 mPa·s; $c_t = 2.1125\cdot10^{-3}$ MPa⁻¹; $B$ = 1.3
Well in cell 4: 100 m³/d at standard conditions

Solution

Face area $A = 200\cdot15 = 3000$ m². The distance from the cell centre to the face is 125 m.
Transmissibilities: $T_{12} = \dfrac{3000}{125/80 + 125/80} = 960$ ; $T_{23} = \dfrac{3000}{125/80 + 125/20} = 384$ ; $T_{34} = \dfrac{3000}{125/20 + 125/100} = 400$ mD·m.
Link coefficients $a = 0.0853\,T/1.2$ : $a_{12} = 68.240$ , $a_{23} = 27.296$ , $a_{34} = 28.433$ m³/(d·MPa).
Pore volume of a cell $V_p = 250\cdot200\cdot15\cdot0.20 = 150{,}000$ m³. Storage capacity $S = 150{,}000\cdot2.1125\cdot10^{-3} = 316.875$ m³/MPa.
Characteristic times $\tau_i = S/\sum_j a_{ij}$ : $\tau_1 = 316.875/68.240 = 4.64$ ; $\tau_2 = 316.875/(68.240 + 27.296) = 3.32$ ; $\tau_3 = 316.875/(27.296 + 28.433) = 5.69$ ; $\tau_4 = 316.875/28.433 = 11.14$ d.
Withdrawal at reservoir conditions: $q_4 = 100\cdot1.3 = 130$ m³/d.

Answer: The explicit scheme is guaranteed to be stable with a step of not more than 3.32 d. Cell 2, which has two well-conducting links, sets this limit.

What it means. The clay-rich cell 3 lowered the transmissibility of both its connections by 2.4–2.5 times compared with the connection of two cells of 80 mD. This is how harmonic averaging works. The “fastest” cell is not the cell with the well but cell 2: it has the most transmissibility per unit of capacity. The exact stability limit for this strip is a little larger, 4.1 d. The condition $\Delta t \le \tau_{\min}$ keeps a safety margin below it. The data in this example are assumed for illustration.

Worked exampleExplicit scheme: two steps of 2 days and steps of 10 days

Given (assumed data)

Strip and coefficients: from the previous example; initial pressure 20 MPa in all cells (above the bubble point pressure)
Step: $\Delta t$ = 2 d, less than $\tau_{\min}$ = 3.3 d; then 10 d

Solution

Step 1. The pressure is the same everywhere, so there are no flows. Only cell 4 changes: $p_4 = 20 - \dfrac{2\cdot130}{316.875} = 19.18$ MPa.
Step 2. The flow from cell 3 into cell 4 is $28.433\cdot(20 - 19.18) = 23.3$ m³/d. Cell 3: $p_3 = 20 - \dfrac{2\cdot23.3}{316.875} = 19.85$ MPa.
Cell 4: $p_4 = 19.18 + \dfrac{2\,(23.3 - 130)}{316.875} = 18.51$ MPa. Cells 1 and 2 do not “see” the disturbance yet: their pressure is 20 MPa.
Comparison with the exact solution of the same cell equations (a calculation with a very small step) after 4 days: $p_3 = 19.79$ and $p_4 = 18.60$ MPa. The error of the explicit scheme is hundredths of a megapascal.
Step of 10 days (above the limit). First step: $p_4 = 20 - 10\cdot130/316.875 = 15.90$ MPa, against the exact 17.04 MPa. After that the calculation goes out of control: after 30 days the pressure in cell 3 rises to 18.7 MPa, although in cell 2 it falls, and after 40 days the pressure in cell 2 is 25.3 MPa.

Answer: With a step of 2 days the explicit scheme gives

p_4

= 18.51 MPa after 4 days (exact: 18.60 MPa). With a step of 10 days the pressure in cell 2 after 40 days is 25.3 MPa, above the initial value.

What it means. Oil is only withdrawn from the closed strip, so the pressure in no cell can rise above the initial 20 MPa. Values above 20 MPa are a sign of instability: each step “throws” fluid across a neighbour, and the oscillations grow. With a step smaller than $\tau_{\min}$ this does not happen, but even then the disturbance passes only one cell in a step. The data in this example are assumed for illustration.

Worked exampleImplicit scheme: one step of 10 days with the Thomas algorithm

Given (assumed data)

Strip and coefficients: from the first example; initial pressure 20 MPa
Step: $\Delta t$ = 10 d, three times larger than $\tau_{\min}$

Solution

The unknowns are the changes of pressure over the step, $\delta_i = p_i^{n+1} - p_i^{n}$ . At the start of the step the pressure is the same everywhere, so the right-hand side is non-zero only for cell 4: $r = (0;\ 0;\ 0;\ -130)$ .
The diagonal is $b_i = S/\Delta t + \sum_j a_{ij}$ , where $S/\Delta t = 31.6875$ : $b_1 = 99.93$ , $b_2 = 127.22$ , $b_3 = 87.42$ , $b_4 = 60.12$ m³/(d·MPa).
Forward pass: $e_1 = 68.240/99.93 = 0.6829$ ; $w_2 = 127.22 - 68.240\cdot0.6829 = 80.62$ , $e_2 = 27.296/80.62 = 0.3386$ .
$w_3 = 87.42 - 27.296\cdot0.3386 = 78.18$ , $e_3 = 28.433/78.18 = 0.3637$ ; $w_4 = 60.12 - 28.433\cdot0.3637 = 49.779$ . All $f_i$ up to the third are zero, and $f_4 = -130/49.779$ .
Backward pass: $\delta_4 = -130/49.779 = -2.612$ ; $\delta_3 = 0.3637\cdot\delta_4 = -0.950$ ; $\delta_2 = 0.3386\cdot\delta_3 = -0.322$ ; $\delta_1 = 0.6829\cdot\delta_2 = -0.220$ MPa.
Pressures after 10 days: 19.78, 19.68, 19.05 and 17.39 MPa. Exact solution of the same equations: 19.92, 19.82, 19.11 and 17.04 MPa.
Balance check: $S\sum\delta_i = 316.875\cdot(-4.10) \approx -1300$ m³, which is as much as the well withdrew in 10 days: $130\cdot10 = 1300$ m³ (without rounding the equality is exact).

Answer: After 10 days the pressure in the cell with the well is 17.39 MPa (exact: 17.04 MPa), and in the far cell 1 it is 19.78 MPa (exact: 19.92 MPa).

What it means. The implicit step, three times larger than the limit of the explicit scheme, went through without oscillations, and the mass was conserved exactly. But the large step “smeared” the disturbance. The cell with the well is depleted less than it should be (by 0.35 MPa), and the far cells are depleted more, as if the pressure had had time to level out. Stability does not replace accuracy: we still choose the step so that the result hardly changes when it is reduced. The data in this example are assumed for illustration.

CalculatedPressure in cell 2 of the strip of the examples: the exact solution of the cell equations, the explicit scheme with steps of 2 and 10 days, and the implicit scheme with a step of 10 days. The dots are the calculated steps. The explicit scheme with a step of 10 days raises the pressure above the initial value, which is a sign of instability.

Why simulators calculate the pressure implicitly

The largest step of the explicit scheme is proportional to the square of the cell size. In a formation with large cells it is measured in days. But models make the grid finer near wells and divide the formation into thin layers, and there the limit falls to minutes and seconds. The step is the same for the whole model, so a single small cell near a well would force us to calculate millions of steps. For a uniform grid it is convenient to write the limit with the hydraulic diffusivity in each direction:

The largest step of the explicit scheme for a uniform grid

\Delta t_{\max} = \frac{1}{2}\left(\frac{\eta_x}{\Delta x^2} + \frac{\eta_y}{\Delta y^2} + \frac{\eta_z}{\Delta z^2}\right)^{-1},\qquad \eta_x = \frac{0.0853\,k_x}{\phi\mu c_t}

where:

$\Delta t_{\max}$ the largest step at which the explicit scheme is stable, d
$\eta_x,\ \eta_y,\ \eta_z$ hydraulic diffusivity along the axes (with the permeability $k_x$ , $k_y$ , $k_z$ ), m²/d
$\Delta x,\ \Delta y,\ \Delta z$ dimensions of the cell, m
$k,\ \phi,\ \mu,\ c_t$ permeability, porosity, viscosity, total compressibility, mD, —, mPa·s, 1/MPa

Units: metric oilfield units: $k$ in mD, $\mu$ in mPa·s, $c_t$ in 1/MPa, dimensions in m, $\Delta t$ in d

Holds when

A uniform grid inside the formation, and a linear equation. In one dimension only one term remains, and $\Delta t_{\max} = \Delta x^2/(2\eta)$ .

This is the same condition $\Delta t \le \tau_i$ , written for an inner cell of a uniform grid. A thin layer with a low vertical permeability can limit the step more strongly than large cells in the horizontal direction.

In other unit systems (1)

Hydraulic diffusivity in metric oilfield units

\eta\,[\text{m}^2/\text{d}] = 0.0853\,\frac{k\,[\text{mD}]}{\phi\;\mu\,[\text{mPa·s}]\;c_t\,[\text{MPa}^{-1}]}

$k$ in mD, $\mu$ in mPa·s, $c_t$ in 1/MPa, $\eta$ in m²/d.

$0.0853 = 24\cdot3.553\cdot10^{-3}$ : the diffusivity constant from the lesson on diffusivity in m²/h, multiplied by 24 h.

Worked exampleThe largest step of the explicit scheme: from large cells to cells near the well

Given (assumed data)

Formation (lesson Flow regimes): $k$ = 50 mD; $\phi$ = 0.20; $\mu$ = 1.2 mPa·s; $c_t = 2.1125\cdot10^{-3}$ MPa⁻¹
Grids: square cells (seen from above) of 250, 50 and 10 m; near the well, a cell of 10 × 10 × 1 m with a vertical permeability of 5 mD

Solution

Hydraulic diffusivity: $\eta = \dfrac{0.0853\cdot50}{0.20\cdot1.2\cdot2.1125\cdot10^{-3}} = 8412$ m²/d. In the vertical direction (5 mD) it is 841.2 m²/d.
A grid with two horizontal directions: $\Delta t_{\max} = \Delta x^2/(4\eta)$ . For cells of 250 m: $250^2/(4\cdot8412) = 1.86$ d.
For cells of 50 m: 0.074 d, which is 1.8 h. For cells of 10 m: 0.0030 d, about 4.3 min.
A three-dimensional cell of 10 × 10 × 1 m: $\dfrac{8412}{100} + \dfrac{8412}{100} + \dfrac{841.2}{1} = 1009$ d⁻¹. In seconds, $\Delta t_{\max} = \dfrac{86{,}400}{2\cdot1009} = 43$ s.
A 30-year forecast ( $30\cdot365\cdot86{,}400$ s) with such a step would need about 22 million steps.

Answer: The largest step falls from 1.86 d for cells of 250 m to 43 s for the thin cell near the well.

What it means. An explicit calculation of pressure in a real model is practically impossible: its step would be set by the smallest cells. So simulators calculate pressure implicitly. They do this either for the pressure alone (the IMPES scheme: implicit pressure, explicit saturation) or for the pressure and the saturations together (the fully implicit method). The step is then chosen by the accuracy needed and by how fast the operating modes of the wells change, usually from fractions of a day to a month (IMPES and fully implicit schemes). The data in this example are assumed for illustration.

An implicit calculation is not free either. At every step we solve a system with millions of unknowns, and for a nonlinear equation we do it several times per step. But the step is now set not by stability but by the accuracy we need, and in the end we gain by orders of magnitude (by factors of ten). The step is usually made variable: small after wells are started or stopped, when the pressure changes fast, and large in quiet periods.

Terms and symbols

Pressure equation: The equation of single-phase flow for the pressure, obtained from the mass balance, Darcy's law and compressibility. In a simulator it is written for every cell.
Control volume: A volume (usually a grid cell) for which we count the mass of fluid that enters, leaves and accumulates.
Storage capacity of a cell $S_i$: $V_{p,i}c_t$ , m³/MPa: the volume of fluid that the cell gives up when the pressure falls by 1 MPa.
Link coefficient $a_{ij}$: $0.0853\,T_{ij}/\mu$ , m³/(d·MPa): the flow between neighbouring cells for a pressure difference of 1 MPa.
Characteristic time of a cell $\tau_i$: The ratio of the storage capacity of the cell to the sum of the coefficients of its links. The explicit scheme is stable for a step not larger than the smallest $\tau_i$ .
Explicit and implicit schemes: The flows during a step are calculated from the pressures at the start of the step (each cell is recalculated separately) or at the end of the step (a system of equations is solved).
Thomas algorithm (sweep method): A direct method (one that gives the answer in a fixed number of operations, without iterations) for solving a system with a tridiagonal matrix. The forward pass eliminates the previous unknown, and the backward pass finds the unknowns. The number of operations is proportional to the number of equations.

Common mistakes and misconceptions

Wrong

“An implicit scheme is stable, so we can take any step with it.”

In fact

Stable, yes; accurate, no. In the example of the lesson, a step of 10 days left the pressure in the cell with the well 0.35 MPa above the exact value, and it depleted the far cells more. We test the step for convergence: the result must not change noticeably when the step is reduced.

Wrong

“The mass balance in the model closes, so the pressures are right.”

In fact

The equations of the cells conserve mass at any step, including in an unstable calculation. In the example, the explicit scheme with a step of 10 days raised the pressure above the initial value, while the sum of $S\,\delta p$ over the cells still equalled the withdrawal. The mass balance checks the solution of the systems of equations, not the accuracy of the scheme.

Wrong

“The pressure equation of a simulator is the diffusivity equation.”

In fact

The diffusivity equation is a special case: a uniform formation, a slightly compressible fluid and constant properties. The simulator solves a nonlinear equation with variable permeability, gravity, and properties that depend on pressure.

Wrong

“The pressure in the cell with a well is the bottomhole pressure.”

In fact

The cell pressure is one number for its whole volume, and the cell is hundreds of times wider than the wellbore. It equals the pressure of radial flow at a distance of about 0.2 of the cell size from the well, and the bottomhole pressure is lower than it by the pressure drop around the well. The well model links them (Well model in a grid block).

In depth: advanced detail

Why the grid does not lose mass

We add the equations of all cells. Each flow enters the sum twice: with a plus sign in the equation of the cell the fluid came into, and with a minus sign in the equation of the cell it left. This is because the link coefficient $a_{ij} = a_{ji}$ . The flows cancel each other, and what remains is:

The sum of the cell equations is the material balance

\sum_i S_i\left(p_i^{n+1} - p_i^{n}\right) = -\Delta t\sum_i q_i

where:

$S_i$ storage capacity of the cell, m³/MPa
$p_i^{n},\ p_i^{n+1}$ cell pressures at the start and at the end of the step, MPa
$q_i$ withdrawal by wells from the cell at reservoir conditions, m³/d
$\Delta t$ time step, d

Units: metric oilfield units

Holds when

Closed outer boundaries. If there is a connection with an aquifer, the inflow through the boundary remains in the sum.
It holds at any step, for both the explicit and the implicit scheme.

This is the material balance of an accumulation as a “tank”: the sum of the pressure drops, weighted by the pore volumes of the cells, equals the volume withdrawn divided by the total compressibility. A grid of many cells contains the “tank” model as a special case and cannot contradict it.

So the balance error that simulators print says nothing about the accuracy of the scheme. It says how accurately the equations were solved at the step: how many iterations were made, and with what tolerance on the residual (the remaining mismatch of the equations). A noticeable balance error is a reason to reduce the step or the tolerances of the solver. A balance with no error says nothing about whether the grid and the step were chosen correctly.

In the strip of the lesson, after a transition period, all cells are depleted at the same rate: $130/(4\cdot316.875) \approx 0.103$ MPa/d. This is the “tank”. The pressure differences between the cells are then constant. The withdrawal of three cells, $97.5$ m³/d, goes through the face between cells 3 and 4, and $97.5/28.433 \approx 3.4$ MPa is lost across it. Across the face between cells 2 and 3 the loss is $65/27.296 \approx 2.4$ MPa, and across the face between cells 1 and 2 it is $32.5/68.24 \approx 0.5$ MPa. The two connections of the clay-rich cell take more than 90% of the pressure drop. This is the same picture as with the skin factor at a well (a number that shows how much the zone right around the well hinders flow).

Where the limit of the explicit step comes from

For one “tank” with a well, an explicit step multiplies the deviation of the pressure from the bottomhole pressure by $1 - \Delta t/\tau$ , and the calculation is stable for $\Delta t < 2\tau$ (stability of the step). On a grid the worst case is a “checkerboard” pattern, in which the pressures of neighbouring cells deviate from the average in opposite directions. Then the pressure difference between a cell and its neighbour is twice the deviation of the cell itself, and over an explicit step the deviation is multiplied by $1 - 2\Delta t/\tau$ . The absolute value of this factor is not larger than one for $\Delta t \le \tau$ . This gives the condition for a uniform grid and the factor 2 in $\Delta x^2/(2\eta)$ .

On a non-uniform or short grid, the “checkerboard” pattern cannot develop in full, and the allowed step is a little larger. In the strip of the lesson $\tau_{\min}$ = 3.3 d, and the exact limit is about 4.1 d. It is given by the largest eigenvalue of the matrix of links divided by the capacities. (An eigenvalue is a number that measures the fastest way in which the pressures in the strip can level out.) The condition $\Delta t \le \tau_i$ is simpler and safer, because we check it for each cell separately. The general theory of the stability and the error of finite-difference schemes is in the lesson Finite differences: from equations to computation.

The nonlinear equation: iterations within a step

In the full pressure equation, the density, the viscosity, the formation volume factor and the porosity depend on the pressure. So the coefficients of the cell equations, the capacity $S_i$ and the link coefficients $a_{ij}$ , also depend on the unknown pressures at the end of the step. The simplest way to solve this is iterations with lagged coefficients. We calculate the coefficients from the pressures at the start of the step and solve the linear system. Then we recalculate the coefficients from the pressures we obtained and repeat. We stop when the changes of pressure between iterations become smaller than a tolerance.

Simulators more often use Newton's method. They approximate the residuals of the equations of all cells (the amounts by which the equations are not yet satisfied) by straight lines in the changes of the unknowns. At each iteration they solve a linear system whose matrix is a table of derivatives (the Jacobian matrix). The method converges in a few iterations (it reaches the needed accuracy) even when the properties change strongly with pressure, for example for a gas whose density is almost proportional to the pressure. In analytical problems we replace the pressure of a gas by the real gas pseudo-pressure. A numerical model does not need such a replacement, because it solves the nonlinear equations directly (IMPES and fully implicit schemes).

Connections to other topics

Diffusivity equation — the same equation for a uniform formation and its analytical solutions; Darcy's law in general form.
What a reservoir simulation model is and why it is needed — the balance of a cell, the transmissibility, and the explicit and implicit step for one “tank”.
Grids in reservoir simulation models — what kinds of grids there are and how the cell size affects the result.
Well model in a grid block — how to link the pressure of the cell with a well to the bottomhole pressure and the rate.
Interblock flow and transmissibility — averaging of the permeability between cells and weighting of the mobility.
Finite differences: from equations to computation — the general theory of finite-difference schemes: error, stability, intermediate schemes.
Two-phase flow equations — the pressure and saturation equations for two phases.
IMPES and fully implicit schemes — the IMPES and fully implicit schemes, Newton iterations.

Check yourself

1
What makes the single-phase pressure equation nonlinear?
1. A.The dependence of density, viscosity and porosity on the pressure we are looking for
2. B.The withdrawal by a well
3. C.The term with gravity
4. D.The derivative of pressure with respect to time
Show the answer
Correct answer: A. The dependence of density, viscosity and porosity on the pressure we are looking for
An equation is nonlinear when its coefficients depend on the unknown function. The density and viscosity of the fluid, the porosity and sometimes the permeability change with pressure. So solutions cannot be added together, and there is no general analytical solution. The time derivative, gravity with a constant density, and a given withdrawal all enter linearly.
2
Estimate the largest step of the explicit scheme for a one-dimensional uniform grid with cells of 100 m: $k$ = 200 mD, $\phi$ = 0.25, $\mu$ = 0.8 mPa·s, $c_t = 1.5\cdot10^{-3}$ MPa⁻¹.
Show the answer
The hydraulic diffusivity is $\eta = 0.0853\cdot200/(0.25\cdot0.8\cdot1.5\cdot10^{-3}) = 17.06/3\cdot10^{-4} \approx 56{,}900$ m²/d. For one dimension, $\Delta t_{\max} = \Delta x^2/(2\eta) = 10{,}000/113{,}700 \approx 0.088$ d, about 2.1 h. For a grid with two horizontal directions and the same square cells the limit is half as large.
3
In a closed model with one producing well, the sum $\sum_i S_i\,(p_i^{n+1} - p_i^{n})$ over all cells for one step equals:
1. A.minus the withdrawal of the well only in the implicit scheme
2. B.zero: the flows between cells balance the withdrawal
3. C.minus the withdrawal of the well only for a step smaller than the limit
4. D.minus the withdrawal of the well over the step, for any scheme and any step
Show the answer
Correct answer: D. minus the withdrawal of the well over the step, for any scheme and any step
Each flow enters the equations of two cells with opposite signs and disappears in the sum, so only the withdrawal of the well remains. This holds for both the explicit and the implicit scheme and for any step, even in an unstable calculation. So the mass balance does not prove that the pressures are calculated accurately.
4
The end cell of a model (permeability 120 mD, length 200 m, face area 2000 m²) is linked with a region of constant pressure. What is the transmissibility of this connection, and why is it twice as large as the transmissibility between two such cells?
Show the answer
The distance from the centre of the cell to its face is half a cell, 100 m: $T_b = kA/(\Delta x/2) = 120\cdot2000/100 = 2400$ mD·m. Between two equal cells the flow passes two halves, $T = A/(100/120 + 100/120) = 1200$ mD·m. The constant pressure is set right at the face, so the resistance is half as large.
5
The cells of a model in plan, 200 × 200 m, were replaced by cells of 40 × 40 m. How does the largest step of the explicit scheme for the pressure change?
1. A.It becomes 25 times larger
2. B.It does not change
3. C.It becomes 25 times smaller
4. D.It becomes 5 times smaller
Show the answer
Correct answer: C. It becomes 25 times smaller
The largest step is proportional to the square of the cell size: $(200/40)^2 = 25$ . This is exactly why small cells near wells make an explicit pressure calculation impractical, and simulators calculate pressure implicitly.
6
Why did the implicit step of 10 days in the example of the lesson give a pressure above the exact value in the cell with the well, and below it in the far cell?
Show the answer
The implicit scheme calculates the flows over the whole step from the pressures at the end of the step, when the drawdown cone has already spread. At the start of the step the flows towards the well are in reality smaller, and the far cells are hardly touched yet. The scheme, as it were, spreads the withdrawal over the whole strip at once: the cell with the well is depleted less, and the far cells more. As the step gets smaller, this smoothing disappears.

Updated: 7 October 2026