Part III: Reservoir & Production Engineering

Chapter 11

Introduction to Reservoir Simulation with Python

schedule19 min readfitness_center10 exercises

auto_awesomeAI Key Takeaways

Get the TL;DR and the key concepts before you dive in — or as a quick review after.

Why This Chapter Exists

Material balance tells you how much oil is in the reservoir and what is driving it out. But it treats the entire reservoir as a single tank; it cannot tell you where the oil is, how the fluids move through the rock, or what happens when you change well locations, injection rates, or completion designs.

Reservoir simulation solves this problem. It divides the reservoir into thousands or millions of grid cells, assigns rock and fluid properties to each cell, and solves the equations of fluid flow at every cell and every time step. The result is a dynamic, spatial model of the reservoir that can predict how fluids move, where pressure fronts develop, and how different development strategies compare.

Commercial reservoir simulators (Eclipse, CMG, tNavigator) cost tens or hundreds of thousands of dollars per license. Understanding what they do under the hood (the physics, the math, and the numerical methods) makes you a better user of those tools and a better engineer. This chapter builds a reservoir simulator from scratch. It will be simple compared to commercial software, but the principles are identical.

infoWhat You'll Learn

The diffusivity equation governing fluid flow in porous media
Finite difference discretization of partial differential equations
Building a 1D single-phase reservoir simulator from scratch
Implicit vs. explicit time-stepping: stability, accuracy, and cost
Buckley-Leverett two-phase displacement theory
Working with output from commercial simulators
Validating simulation results against analytical solutions

The Physics of Flow in Porous Media

Every cell in the grid you are about to build needs an equation that says how fast fluid crosses its faces and how much its pressure changes as a result. Three pieces of physics supply that equation. Fluid flow through reservoir rock is governed by three fundamental principles: conservation of mass, Darcy's Law, and an equation of state.

Darcy's Law

In 1856, Henry Darcy showed experimentally that the flow rate of a fluid through a porous medium is proportional to the pressure gradient and inversely proportional to the fluid viscosity:

u = -\frac{k}{\mu} \frac{\partial P}{\partial x}

where $u$ is the superficial velocity (flow rate per unit area), $k$ is the rock permeability, $\mu$ is the fluid viscosity, and $\partial P / \partial x$ is the pressure gradient.

The physical interpretation: fluid flows from high pressure to low pressure. The rate of flow depends on how permeable the rock is (how connected the pore spaces are) and how viscous the fluid is (how much it resists flowing). Darcy's Law is the petroleum engineering equivalent of Ohm's Law in electrical circuits: pressure difference drives flow, permeability is conductance, viscosity is resistance.

Conservation of Mass

The mass of fluid entering a volume element minus the mass leaving must equal the change in mass stored within the element:

\frac{\partial}{\partial x}\left(\frac{k}{\mu}\frac{\partial P}{\partial x}\right) = \phi c_t \frac{\partial P}{\partial t}

This is the diffusivity equation for single-phase, slightly compressible flow. It combines Darcy's Law with mass conservation and a compressibility equation of state.

$\phi$ is porosity (the fraction of rock volume that is pore space)
$c_t$ is total compressibility (how much the fluid and rock compress per unit pressure drop)
The left side represents spatial variation in flow
The right side represents temporal change in pressure

This equation is a partial differential equation (PDE). It cannot be solved analytically for most realistic reservoir geometries and boundary conditions. Numerical methods are required.

Finite Difference Discretization

The diffusivity equation has no closed-form solution for a real reservoir, so a computer cannot integrate it directly; it can only add and multiply numbers in arrays. Discretization is what turns the continuous equation into exactly that. To solve the diffusivity equation numerically, we discretize it: convert the continuous PDE into a system of algebraic equations that a computer can solve.

Spatial Discretization

Divide the reservoir into $n$ grid cells of width $\Delta x$ . Replace the spatial derivative with a central difference approximation:

\frac{\partial}{\partial x}\left(\frac{k}{\mu}\frac{\partial P}{\partial x}\right) \approx \frac{T_{i+1/2}(P_{i+1} - P_i) - T_{i-1/2}(P_i - P_{i-1})}{\Delta x^2}

where $T_{i+1/2} = \dfrac{\beta\, k_{i+1/2}\, A}{\mu\, \Delta x}$ is the transmissibility between cells $i$ and $i+1$ (RB/d/psi). $A$ is the cross-sectional area, $k_{i+1/2}$ is the harmonic mean of the two cells' permeabilities (the correct average for series flow), and $\beta = 1.127\times10^{-3}$ is the Darcy field-unit constant. Flux, accumulation, and the well source are all carried in reservoir barrels, so the formation volume factor $B_o$ enters only on the well term (converting STB/d to RB/d), not here. Multiplying the equation through by the cell bulk volume gives the finite-volume form the code solves.

Temporal Discretization

Replace the time derivative with a forward difference:

\frac{\partial P}{\partial t} \approx \frac{P_i^{n+1} - P_i^n}{\Delta t}

How you combine spatial and temporal discretization determines whether the scheme is explicit or implicit.

Building a 1D Single-Phase Simulator

The class below assembles exactly the finite-volume system from the last section. Each grid cell contributes one row to a tridiagonal matrix: the diagonal holds its accumulation term plus the transmissibilities to its neighbours, the off-diagonals hold those neighbour transmissibilities, and the right-hand side carries the old pressure and any well source. spsolve then returns the new pressure field in one step.

main.pystarting Python…

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve


class Reservoir1D:
    """1D single-phase, slightly-compressible reservoir simulator.

    Implicit finite-volume scheme: transmissibilities use a harmonic-mean
    permeability and the real cross-sectional area, and the tridiagonal system
    is solved with a sparse solver. Wells are rate-controlled (a volumetric
    source) or BHP-controlled through a Peaceman well model.
    """

    def __init__(self, nx, length_ft, width_ft, height_ft, k_md, phi, ct_psi,
                 mu_cp, Pi_psi, Bo=1.0):
        self.nx = nx
        self.dx = length_ft / nx
        self.width = width_ft
        self.height = height_ft
        self.area = width_ft * height_ft          # cross-section, ft^2
        self.vb = self.area * self.dx             # bulk volume per cell, ft^3
        self.k = np.full(nx, k_md, float) if np.isscalar(k_md) else np.asarray(k_md, float)
        self.phi = phi
        self.ct = ct_psi
        self.mu = mu_cp
        self.Bo = Bo
        self.P = np.full(nx, Pi_psi, dtype=float)
        self.Pi = Pi_psi
        self.beta = 1.127e-3                       # Darcy constant -> T in RB/(d.psi)
        self.wells = []
        self.time_days = [0.0]
        self.pressure_history = [self.P.copy()]
        self.well_rates = []                      # net voidage (RB/d) per step

    def _transmissibility(self, i, j):
        """Inter-cell transmissibility (RB/d/psi), harmonic-mean permeability.

        Reservoir-volume basis: flux, accumulation, and the rate source are all
        in reservoir bbl, so Bo appears only on the source (q_stb -> q_res), not
        here. (Putting Bo here too would double-count it.)
        """
        kh = 2 * self.k[i] * self.k[j] / (self.k[i] + self.k[j])
        return self.beta * kh * self.area / (self.mu * self.dx)

    def _well_index(self, well):
        """Peaceman productivity index (RB/d/psi) for a BHP-controlled well."""
        re = 0.2 * self.dx                        # equivalent drainage radius (1D)
        return (2 * np.pi * self.beta * self.k[well['cell']] * self.height
                / (self.mu * (np.log(re / well['rw']) + well['skin'])))

    def add_well(self, cell_index, rate_stbd=None, bhp_psi=None, rw_ft=0.354, skin=0.0):
        """Add a rate- or BHP-controlled well (rate_stbd < 0 is a producer, > 0 is an injector (STB/d))."""
        self.wells.append({'cell': cell_index, 'rate': rate_stbd, 'bhp': bhp_psi,
                           'rw': rw_ft, 'skin': skin})

    def step(self, dt_days):
        """Advance one implicit time step of dt_days."""
        n = self.nx
        # Accumulation coefficient (RB/d/psi); 5.615 converts ft^3 -> RB.
        C = self.vb * self.phi * self.ct / (5.615 * dt_days)
        lower = np.zeros(n); diag_ = np.zeros(n); upper = np.zeros(n); b = np.zeros(n)
        for i in range(n):
            t_left = self._transmissibility(i, i - 1) if i > 0 else 0.0
            t_right = self._transmissibility(i, i + 1) if i < n - 1 else 0.0
            diag_[i] = C + t_left + t_right
            if i > 0:
                lower[i] = -t_left
            if i < n - 1:
                upper[i] = -t_right
            b[i] = C * self.P[i]

        for w in self.wells:
            i = w['cell']
            if w['rate'] is not None:
                # Volumetric source: STB/d -> RB/d (sign per add_well docstring).
                b[i] += w['rate'] * self.Bo
            elif w['bhp'] is not None:
                wi = self._well_index(w)
                diag_[i] += wi
                b[i] += wi * w['bhp']

        matrix = diags([lower[1:], diag_, upper[:-1]], [-1, 0, 1], format="csc")
        self.P = spsolve(matrix, b)

        voidage = 0.0
        for w in self.wells:
            if w['rate'] is not None:
                voidage += w['rate'] * self.Bo
            elif w['bhp'] is not None:
                voidage += self._well_index(w) * (w['bhp'] - self.P[w['cell']])
        self.well_rates.append(voidage)
        self.time_days.append(self.time_days[-1] + dt_days)
        self.pressure_history.append(self.P.copy())

    def run(self, total_days, dt_days):
        """Run for total_days with time step dt_days."""
        for _ in range(int(total_days / dt_days)):
            self.step(dt_days)

    def plot_pressure(self, times_to_plot=None):
        """Plot pressure profiles at selected times."""
        x = np.linspace(self.dx / 2, self.nx * self.dx - self.dx / 2, self.nx)
        fig, ax = plt.subplots(figsize=(10, 5))
        if times_to_plot is None:
            indices = np.linspace(0, len(self.pressure_history) - 1, 6, dtype=int)
        else:
            indices = [np.argmin(np.abs(np.array(self.time_days) - t)) for t in times_to_plot]
        for idx in indices:
            t = self.time_days[idx]
            ax.plot(x, self.pressure_history[idx], linewidth=1.5, label=f"t = {t:.0f} days")
        ax.set_xlabel("Distance (ft)")
        ax.set_ylabel("Pressure (psi)")
        ax.set_title("Reservoir Pressure Profile")
        ax.legend(fontsize=8)
        ax.grid(True, alpha=0.3)
        plt.tight_layout()
        plt.show()

Running the Simulator

main.pystarting Python…

# An OML 58 sand drained toward producer OD-001 at the left edge
reservoir = Reservoir1D(
    nx=50,              # 50 grid cells
    length_ft=5000,     # 5,000 ft long
    width_ft=1000,      # 1,000 ft wide
    height_ft=100,      # 100 ft net pay
    k_md=100,           # 100 md permeability
    phi=0.20,           # 20% porosity
    ct_psi=1.5e-5,      # total compressibility (1/psi)
    mu_cp=2.0,          # oil viscosity (cp)
    Pi_psi=4000,        # initial pressure (psi)
    Bo=1.2,             # oil FVF (RB/STB)
)

# Producer OD-001 at the left end (cell 0), flowing at 2,500 psi BHP
# (the well index is the physical Peaceman PI, not a numerical penalty)
reservoir.add_well(cell_index=0, bhp_psi=2500)

# Run for 365 days with 1-day time steps
reservoir.run(total_days=365, dt_days=1)

# Plot pressure at selected times
reservoir.plot_pressure(times_to_plot=[0, 30, 90, 180, 365])

print(f"Initial pressure:  {reservoir.Pi} psi (uniform)")
print(f"Final pressure at well (cell 0): {reservoir.P[0]:.0f} psi")
print(f"Final pressure at boundary (cell {reservoir.nx-1}): {reservoir.P[-1]:.0f} psi")
print(f"Average reservoir pressure: {np.mean(reservoir.P):.0f} psi")

The pressure profiles show the characteristic shape of a drawdown: pressure drops sharply near the well and the disturbance propagates outward over time. Early time steps show the pressure front still traveling through the reservoir (transient flow). At late times, the pressure front has reached the boundary and the entire reservoir is depleting (boundary-dominated flow).

This transition from transient to boundary-dominated flow is fundamental in reservoir engineering. It affects decline curve behavior, well test interpretation, and spacing decisions.

Implicit vs. Explicit Methods

The simulator above uses implicit time-stepping, meaning it solves a linear system at each time step. The alternative is explicit time-stepping, where the new pressure is calculated directly from the old pressure without solving a system.

Explicit Method

P_i^{n+1} = P_i^n + \frac{\Delta t}{\phi c_t \Delta x^2} \left[ T_{i+1/2}(P_{i+1}^n - P_i^n) - T_{i-1/2}(P_i^n - P_{i-1}^n) \right]

The explicit method is simpler to implement but has a critical limitation: it is only stable if the time step satisfies the CFL condition:

\Delta t < \frac{\phi \mu c_t \Delta x^2}{2k \cdot 6.328 \times 10^{-3}}

If you exceed this limit, the solution oscillates and diverges. In practical terms, explicit methods require very small time steps for high-permeability reservoirs or fine grids, making them prohibitively slow.

Implicit Method

The implicit method has no stability restriction: you can take arbitrarily large time steps without divergence. The cost is that you must solve a linear system ( $A \cdot P^{n+1} = b$ ) at each step. For a tridiagonal system (1D problem), this is efficient. For 3D problems with millions of cells, it requires iterative solvers and sophisticated linear algebra.

main.pystarting Python…

# Demonstrate the stability difference
# Explicit method with a time step that violates CFL

def explicit_step(P, k, phi, ct, mu, dx, dt):
    """One explicit time step."""
    n = len(P)
    P_new = P.copy()
    alpha = 6.328e-3
    coeff = alpha * dt / (phi * ct * mu * dx**2)

    for i in range(1, n-1):
        k_right = 2*k[i]*k[i+1]/(k[i]+k[i+1])
        k_left = 2*k[i]*k[i-1]/(k[i]+k[i-1])
        P_new[i] = P[i] + coeff * (k_right*(P[i+1]-P[i]) - k_left*(P[i]-P[i-1]))

    return P_new

# Calculate the CFL limit
k_val = 100  # md
phi_val = 0.20
ct_val = 1.5e-5
mu_val = 2.0
dx_val = 100  # ft (5000/50)
dt_cfl = phi_val * mu_val * ct_val * dx_val**2 / (2 * k_val * 6.328e-3)

print(f"CFL stability limit: {dt_cfl:.2f} days")
print(f"Implicit method:     No stability limit (any dt works)")
print(f"Explicit method:     dt must be < {dt_cfl:.2f} days")
print()
print("This is why commercial simulators use implicit methods.")
print(f"For this reservoir, explicit requires dt < {dt_cfl:.2f} days,")
print(f"meaning {365/dt_cfl:.0f} time steps per year vs. 365 for implicit (dt=1 day).")

Buckley-Leverett Two-Phase Displacement

When water is injected into an oil reservoir, two immiscible fluids flow simultaneously through the same rock. The displacement process is governed by the Buckley-Leverett equation, which predicts the position of the water front as a function of time.

The key concept is fractional flow: the fraction of total flow that is water at any point in the reservoir.

f_w = \frac{1}{1 + \frac{k_{ro} \mu_w}{k_{rw} \mu_o}}

where $k_{ro}$ and $k_{rw}$ are the relative permeabilities of oil and water, which depend on water saturation.

main.pystarting Python…

def corey_kr(Sw, Swi, Sor, kro_max=0.8, krw_max=0.3, no=2.0, nw=3.0):
    """Corey relative permeability model."""
    Sw_norm = (Sw - Swi) / (1 - Swi - Sor)
    Sw_norm = np.clip(Sw_norm, 0, 1)

    kro = kro_max * (1 - Sw_norm)**no
    krw = krw_max * Sw_norm**nw

    return kro, krw

def fractional_flow(Sw, Swi, Sor, mu_o, mu_w, **kr_params):
    """Calculate water fractional flow."""
    kro, krw = corey_kr(Sw, Swi, Sor, **kr_params)
    fw = 1 / (1 + kro * mu_w / (krw * mu_o + 1e-20))
    return fw

# Parameters
Swi = 0.20   # Irreducible water saturation
Sor = 0.20   # Residual oil saturation
mu_o = 3.0   # Oil viscosity (cp)
mu_w = 0.5   # Water viscosity (cp)

# Generate fractional flow curve
Sw_range = np.linspace(Swi, 1 - Sor, 200)
fw_range = fractional_flow(Sw_range, Swi, Sor, mu_o, mu_w)

# Buckley-Leverett: find the shock front using the Welge tangent
# The tangent line from (Swi, 0) to the fw curve gives the front saturation
dfw_dSw = np.gradient(fw_range, Sw_range)

# Welge tangent: maximize (fw - 0) / (Sw - Swi)
tangent_slope = fw_range / (Sw_range - Swi + 1e-10)
idx_front = np.argmax(tangent_slope)
Sw_front = Sw_range[idx_front]
fw_front = fw_range[idx_front]

fig, ax1 = plt.subplots(figsize=(7, 5))

# Fractional flow curve with the Welge tangent
ax1.plot(Sw_range, fw_range, "k-", linewidth=2)
ax1.plot([Swi, Sw_front], [0, fw_front], "r--", linewidth=1.5, label="Welge tangent")
ax1.axvline(x=Sw_front, color="r", linestyle=":", alpha=0.5)
ax1.scatter([Sw_front], [fw_front], c="red", s=50, zorder=5)
ax1.set_xlabel("Water Saturation, Sw")
ax1.set_ylabel("Fractional Flow, fw")
ax1.set_title("Fractional Flow Curve (Welge Tangent)")
ax1.annotate(f"Sw_front = {Sw_front:.3f}", (Sw_front, fw_front),
             xytext=(15, -20), textcoords="offset points", fontsize=9)
ax1.legend(fontsize=8)
ax1.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"Shock front saturation: {Sw_front:.3f}")
print(f"Fractional flow at front: {fw_front:.3f}")
print(f"Oil displaced by water injection depends on this front behavior.")

Buckley-Leverett water saturation profile showing the displacement front.

The Buckley-Leverett solution reveals a sharp saturation discontinuity, the shock front, where water abruptly replaces oil. Behind the front, the rock is at near-residual oil saturation. Ahead of the front, only connate water is present. The front velocity depends on the derivative of the fractional flow curve, which is why the viscosity ratio between oil and water is so important for displacement efficiency.

Coupling Pressure and Saturation: an IMPES Waterflood

Fractional-flow theory is analytical: it tells you the shape of the answer. A real simulator couples it to the pressure solution. The standard scheme is IMPES (IMplicit Pressure, Explicit Saturation): each step solves the pressure equation with the current total mobility, computes the inter-cell flows, then advances water saturation explicitly with upstream weighting (the saturation of a face is taken from the cell the fluid is flowing from, central differencing of a hyperbolic transport term produces non-physical oscillations). This is the actual engine behind a waterflood forecast, and it reproduces the Buckley-Leverett front rather than a hand-drawn sketch.

main.pystarting Python…

from scipy.sparse import diags
from scipy.sparse.linalg import spsolve

# reuse corey_kr defined in the Buckley-Leverett cell above


def impes_waterflood(nx=100, length_ft=1000.0, width_ft=500.0, height_ft=50.0,
                     k_md=300.0, phi=0.25, mu_o=3.0, mu_w=0.5,
                     Swi=0.2, Sor=0.2, q_inj_stbd=400.0, Bw=1.0,
                     Pprod=3000.0, days=600.0, dt=0.5,
                     snapshots=(150, 300, 450, 600)):
    """1D incompressible waterflood by IMPES. Returns (x, {time: Sw profile})."""
    dx = length_ft / nx
    area = width_ft * height_ft
    Vp = area * dx * phi / 5.615             # pore volume per cell, RB
    beta = 1.127e-3
    Sw = np.full(nx, Swi)
    q_inj = q_inj_stbd * Bw                  # injected water, RB/d
    profiles = {}
    snap_steps = {int(s / dt) for s in snapshots}
    for step in range(1, int(days / dt) + 1):
        kro, krw = corey_kr(Sw, Swi, Sor)
        lam_t = kro / mu_o + krw / mu_w                 # total mobility per cell
        T = beta * k_md * area / dx * lam_t[:-1]        # upstream face transmissibility
        lower = np.zeros(nx); diag_ = np.zeros(nx); upper = np.zeros(nx); b = np.zeros(nx)
        for i in range(nx):
            tl = T[i - 1] if i > 0 else 0.0
            tr = T[i] if i < nx - 1 else 0.0
            diag_[i] = tl + tr
            if i > 0: lower[i] = -tl
            if i < nx - 1: upper[i] = -tr
        b[0] += q_inj                                   # water injector at cell 0
        wi = 1e7; diag_[-1] += wi; b[-1] += wi * Pprod  # producer holds BHP = Pprod
        P = spsolve(diags([lower[1:], diag_, upper[:-1]], [-1, 0, 1], format="csc"), b)
        qf = T * (P[:-1] - P[1:])                       # face flux (RB/d, rightward +)
        fw = (krw / mu_w) / (krw / mu_w + kro / mu_o + 1e-30)
        fw_face = np.where(qf >= 0, fw[:-1], fw[1:])    # UPSTREAM fractional flow
        qw = fw_face * qf                               # water flux per face
        Sw = Sw.copy()
        Sw[0] += dt / Vp * (q_inj - qw[0])
        Sw[1:-1] += dt / Vp * (qw[:-1] - qw[1:])
        Sw[-1] += dt / Vp * (qw[-1] - fw[-1] * qf[-1])
        Sw = np.clip(Sw, Swi, 1 - Sor)
        if step in snap_steps:
            profiles[round(step * dt)] = Sw.copy()
    x = (np.arange(nx) + 0.5) * dx
    return x, profiles


x_imp, profiles = impes_waterflood()

fig, ax = plt.subplots(figsize=(9, 5))
for t, prof in sorted(profiles.items()):
    ax.plot(x_imp, prof, linewidth=1.8, label=f"t = {t} days")
ax.set_xlabel("Distance from injector (ft)")
ax.set_ylabel("Water Saturation, Sw")
ax.set_title("IMPES Waterflood - Simulated Saturation Fronts")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Validate the simulated front against the Buckley-Leverett prediction
v_front = (400 * 5.615) / (500 * 50 * 0.25) * (fw_front / (Sw_front - Swi))
final = profiles[max(profiles)]
sim_front = x_imp[np.where(final > Swi + 0.02 * (1 - Swi - Sor))[0][-1]]
print(f"Buckley-Leverett front velocity:  {v_front:.2f} ft/day")
print(f"Predicted front at 600 days:      {v_front * 600:.0f} ft")
print(f"Simulated front at 600 days:      {sim_front:.0f} ft "
      f"(agrees with B-L to within grid dispersion)")

IMPES waterflood: the simulated saturation front advancing toward the producer.

The simulated front lands within a few percent of the analytical Buckley-Leverett position; the small lag is numerical dispersion, the smearing a first-order upwind scheme adds to a true shock. Refining the grid (larger nx) sharpens the front toward the analytical discontinuity, the displacement-engineering version of the grid-refinement study in Exercise 11.1.

Summary

Reservoir simulation solves the equations of fluid flow at every point in space and time, providing spatial resolution that material balance cannot.
Darcy's Law relates flow rate to pressure gradient, permeability, and viscosity. It is the fundamental flow equation in porous media.
The diffusivity equation combines Darcy's Law with mass conservation. It is a PDE that must be solved numerically for realistic reservoirs.
Finite differences convert the continuous PDE into discrete algebraic equations by approximating derivatives on a grid.
Implicit methods have no stability restriction and allow large time steps. Explicit methods are simpler but require time steps below the CFL limit.
Buckley-Leverett theory predicts the displacement of oil by water. The shock front saturation and velocity are determined by the fractional flow curve.
Commercial simulators use the same principles but extend them to 3D, multiphase flow, complex geology, and sophisticated well models. The leap to a field model is one of scale and bookkeeping, not of new physics; the mass-balance and Darcy core you built here is unchanged.

Exercises

fitness_center

Exercise 11.1Practice

: Grid Refinement Study

Run the 1D simulator from this chapter with 10, 25, 50, 100, and 200 grid cells. Plot the pressure profile at t=180 days for each case. At what grid r...

arrow_forward