Part III: Reservoir & Production Engineering

Chapter 11

Introduction to Reservoir Simulation with Python

schedule15 min readfitness_center10 exercises

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

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

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} = \frac{k_{i+1/2}}{\mu}$ is the transmissibility between cells $i$ and $i+1$ .

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

main.py

import numpy as np
import matplotlib.pyplot as plt

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

    Solves the pressure diffusivity equation using
    finite differences with implicit time-stepping.
    """

    def __init__(self, nx, length_ft, k_md, phi, ct_psi,
                 mu_cp, Pi_psi, Boi=1.0):
        """
        Parameters:
        -----------
        nx        : int    -- number of grid cells
        length_ft : float  -- total reservoir length (ft)
        k_md      : float or array -- permeability (md)
        phi       : float  -- porosity (fraction)
        ct_psi    : float  -- total compressibility (1/psi)
        mu_cp     : float  -- viscosity (cp)
        Pi_psi    : float  -- initial pressure (psi)
        Boi       : float  -- initial formation volume factor
        """
        self.nx = nx
        self.dx = length_ft / nx                # ft
        self.k = np.full(nx, k_md) if np.isscalar(k_md) else np.array(k_md)
        self.phi = phi
        self.ct = ct_psi
        self.mu = mu_cp
        self.Bo = Boi

        # Initialize pressure uniformly
        self.P = np.full(nx, Pi_psi, dtype=float)
        self.Pi = Pi_psi

        # Conversion: md·ft / cp → RB/(d·psi) requires 6.328e-3
        self.alpha = 6.328e-3  # Darcy unit conversion

        # Storage for results
        self.time_days = [0.0]
        self.pressure_history = [self.P.copy()]

    def _build_coefficient_matrix(self, dt):
        """
        Build the tridiagonal coefficient matrix for implicit solve.
        A·P^(n+1) = b
        """
        n = self.nx
        A = np.zeros((n, n))
        b = np.zeros(n)

        for i in range(n):
            # Transmissibilities to neighbors
            if i > 0:
                k_left = 2 * self.k[i] * self.k[i-1] / (self.k[i] + self.k[i-1])
                T_left = self.alpha * k_left / (self.mu * self.Bo * self.dx**2)
            else:
                T_left = 0

            if i < n - 1:
                k_right = 2 * self.k[i] * self.k[i+1] / (self.k[i] + self.k[i+1])
                T_right = self.alpha * k_right / (self.mu * self.Bo * self.dx**2)
            else:
                T_right = 0

            # Accumulation term
            acc = self.phi * self.ct / (dt)  # 1/day

            # Fill matrix row
            A[i, i] = acc + T_left + T_right
            if i > 0:
                A[i, i-1] = -T_left
            if i < n - 1:
                A[i, i+1] = -T_right

            # Right-hand side: accumulation × P^n + source terms
            b[i] = acc * self.P[i]

        return A, b

    def add_well(self, cell_index, rate_stbd=None, bhp_psi=None, wi=None):
        """
        Add a well to a specific cell.
        Specify either rate (STB/d, negative for producer) or BHP.
        wi: well index (productivity index), STB/d/psi
        """
        if not hasattr(self, 'wells'):
            self.wells = []

        self.wells.append({
            'cell': cell_index,
            'rate': rate_stbd,
            'bhp': bhp_psi,
            'wi': wi if wi else 10.0  # Default PI
        })

    def step(self, dt_days):
        """Advance the simulation by dt_days."""
        A, b = self._build_coefficient_matrix(dt_days)

        # Apply well conditions
        if hasattr(self, 'wells'):
            for well in self.wells:
                i = well['cell']
                if well['rate'] is not None:
                    # Specified rate: add source term to RHS
                    # Negative rate = production
                    vol_cell = self.phi * self.dx * 1.0 * 1.0  # ft³ (unit cross-section)
                    b[i] += well['rate'] / (vol_cell * self.Bo)
                elif well['bhp'] is not None:
                    # Specified BHP: add well equation
                    wi = well['wi']
                    A[i, i] += wi
                    b[i] += wi * well['bhp']

        # Solve the linear system
        self.P = np.linalg.solve(A, b)
        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 the simulation for total_days with time step dt_days."""
        n_steps = int(total_days / dt_days)
        for _ in range(n_steps):
            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:
            # Plot 6 evenly spaced snapshots
            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]
            label = f"t = {t:.0f} days" if t >= 1 else "t = 0"
            ax.plot(x, self.pressure_history[idx], linewidth=1.5, label=label)

        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.py

# Create a 1D reservoir
reservoir = Reservoir1D(
    nx=50,              # 50 grid cells
    length_ft=5000,     # 5,000 ft long
    k_md=100,           # 100 md permeability
    phi=0.20,           # 20% porosity
    ct_psi=1.5e-5,      # Total compressibility
    mu_cp=2.0,          # Oil viscosity
    Pi_psi=4000         # Initial pressure
)

# Place a producer at the left end (cell 0)
reservoir.add_well(cell_index=0, bhp_psi=2000, wi=5.0)

# 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.py

# 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.py

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, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Fractional flow curve
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")
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)

# Saturation profile at a given time
# x_D = (fw'(Sw) * qi * t) / (A * phi * L)
# Simplified: just show the characteristic profile shape
x_profile = np.linspace(0, 1, 200)
Sw_profile = np.full_like(x_profile, Swi)

# Behind the front: Sw increases from front saturation to 1-Sor
front_pos = 0.4  # Fraction of reservoir length
Sw_profile[x_profile < front_pos] = 1 - Sor  # Swept zone (simplified)
# At the front: sharp jump
# Ahead: connate water

ax2.plot(x_profile, Sw_profile, "b-", linewidth=2)
ax2.axvline(x=front_pos, color="r", linestyle="--", alpha=0.5, label="Water front")
ax2.fill_between(x_profile, Swi, Sw_profile, alpha=0.1, color="blue")
ax2.set_xlabel("Dimensionless Distance (x/L)")
ax2.set_ylabel("Water Saturation, Sw")
ax2.set_title("Saturation Profile")
ax2.set_ylim(0, 1)
ax2.legend(fontsize=8)
ax2.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.")

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.

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. Understanding the fundamentals makes you a better user of those tools.

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