Permeability Heterogeneity - A Low-k Streak Throttles the Flow

Level 3

Chapter 11: Reservoir Simulation

descriptionProblem

Modify the simulator to use a non-uniform permeability field: cells 1–25 have k=200 md, cells 26–50 have k=20 md. Run the simulation and compare the pressure profile to the uniform case. Where does pressure drop most steeply and why?

---

Real reservoirs are never a single uniform sand. The OML 58 flow unit you are modelling is a 5,000-ft-long, 200 md sand, except for a tight, low-permeability streak (a cemented or shaly interval) of just 10 md sitting in the middle of the grid. Flow has to squeeze through that streak, and Darcy's law guarantees the consequence: where transmissibility collapses, the pressure gradient must spike to keep the same flux moving. That localized pressure drop is exactly what a simulator captures and a hand-drawn straight-line gradient never would.

The setup is a steady-throughflow line drive: a constant-rate water injector at cell 0 pushes fluid toward a constant-BHP producer at the last cell, so a stable pressure gradient develops across the whole grid. You will build the non-uniform permeability field, run the simulator, and prove the streak throttles the flow.

The verified Reservoir1D simulator (harmonic-mean transmissibility, real cross-section, Peaceman wells) and all reservoir constants are embedded. Do not re-derive or edit them.

Your tasks:

Write build_k() returning a length-NX float array: every cell is K_SAND

(200 md) except the streak cells STREAK_START:STREAK_END (cells 20–24, five cells), which are K_STREAK (10 md). This is the heterogeneous permeability field.

Write run_hetero(k_array):

construct `Reservoir1D(NX, LENGTH_FT, WIDTH_FT, HEIGHT_FT, k_array, PHI, CT_PSI,

MU_CP, PI_PSI, Bo=BO)`,

add a constant-rate injector at cell 0: add_well(0, rate_stbd=INJ_RATE_STBD),
add a constant-BHP producer at the last cell:

add_well(NX - 1, bhp_psi=PROD_BHP_PSI),

run(DAYS, DT_DAYS) and return the final pressure array res.P (psi).

Write dp_across(p_psi, i0, i1) returning the absolute pressure drop

abs(p_psi[i1] - p_psi[i0]) between two cell indices.

Write harmonic_T(k_array, i, j) returning the harmonic-mean transmissibility

at the boundary between adjacent cells i and j, recomputed from scratch:

T = \frac{\beta\,k_h\,A}{\mu\,B_o\,\Delta x}$$ with `BETA = 1.127e-3`, area `A = WIDTH_FT * HEIGHT_FT`, and `dx = LENGTH_FT / NX`. (This is the same `_T` the simulator uses internally; you are confirming you understand where the throttling comes from.) 5. Run it: build `k_field = build_k()`, `p_final = run_hetero(k_field)`, then store the drop across the streak as `dp_streak_psi = dp_across(p_final, STREAK_START, STREAK_END)` and the drop across an **equal five-cell span of clean sand** as `dp_sand_psi = dp_across(p_final, SAND_LO, SAND_HI)`. > **Think about it:** the streak's boundary transmissibility is the *harmonic* mean > of 200 md and 10 md (about **19 md**, barely above the low value, not the average > of ~105 md. A series of resistances is dominated by the worst one. That is why a > thin tight streak can drop more pressure than the entire rest of the reservoir > combined, and why upscaling permeability with an arithmetic average badly > over-predicts deliverability through layered or baffled rock.

lightbulbHints (0/3)

Stuck? Reveal hints one at a time — they progress from nudge to near-solution.

codeYour solution

main.py

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve

# ── Verified 1D single-phase implicit simulator (do not edit) ────────────
class Reservoir1D:
    """1D single-phase implicit simulator (harmonic-mean T, real area, Peaceman wells)."""
    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.area = width_ft * height_ft; self.vb = self.area * self.dx; self.height = height_ft
        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, float); self.Pi = Pi_psi; self.beta = 1.127e-3; self.wells = []
    def _T(self, i, j):
        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 _wi(self, w):
        re = 0.2*self.dx
        return 2*np.pi*self.beta*self.k[w["cell"]]*self.height/(self.mu*(np.log(re/w["rw"])+w["skin"]))
    def add_well(self, cell, rate_stbd=None, bhp_psi=None, rw_ft=0.354, skin=0.0):
        self.wells.append(dict(cell=cell, rate=rate_stbd, bhp=bhp_psi, rw=rw_ft, skin=skin))
    def step(self, dt):
        n = self.nx; C = self.vb*self.phi*self.ct/(5.615*dt)
        lo = np.zeros(n); di = np.zeros(n); up = np.zeros(n); b = np.zeros(n)
        for i in range(n):
            tl = self._T(i, i-1) if i > 0 else 0.0; tr = self._T(i, i+1) if i < n-1 else 0.0
            di[i] = C + tl + tr
            if i > 0: lo[i] = -tl
            if i < n-1: up[i] = -tr
            b[i] = C*self.P[i]
        for w in self.wells:
            i = w["cell"]
            if w["rate"] is not None: b[i] += w["rate"]*self.Bo
            elif w["bhp"] is not None: wi = self._wi(w); di[i] += wi; b[i] += wi*w["bhp"]
        self.P = spsolve(diags([lo[1:], di, up[:-1]], [-1, 0, 1], format="csc"), b)
    def run(self, days, dt):
        for _ in range(int(days/dt)): self.step(dt)

# ── OML 58 flow unit: a 200 md sand with a tight 10 md streak in the middle ─
NX = 50                  # grid cells along the 5,000 ft length
LENGTH_FT = 5000.0
WIDTH_FT = 500.0
HEIGHT_FT = 50.0
PHI = 0.20
CT_PSI = 1.0e-5
MU_CP = 1.2
PI_PSI = 3500.0          # initial reservoir pressure
BO = 1.2

K_SAND = 200.0           # clean sand permeability (md)
K_STREAK = 10.0          # tight low-perm streak permeability (md)
STREAK_START = 20        # first streak cell (inclusive)
STREAK_END = 25          # one past the last streak cell -> cells 20..24 (5 cells)

SAND_LO = 5              # an equal five-cell span of CLEAN sand, away from the streak
SAND_HI = 10

INJ_RATE_STBD = 600.0    # constant-rate water injector at cell 0
PROD_BHP_PSI = 2500.0    # constant-BHP producer at the last cell
DAYS = 1095.0            # 3 years of throughflow
DT_DAYS = 15.0

BETA = 1.127e-3          # transmissibility constant (field units)

def build_k():
    """Length-NX permeability field: K_SAND everywhere, K_STREAK in the streak cells."""
    # TODO: k = np.full(NX, K_SAND, float)
    # TODO: k[STREAK_START:STREAK_END] = K_STREAK
    # TODO: return k
    pass

def run_hetero(k_array):
    """Run the line-drive (rate injector at cell 0, BHP producer at NX-1) and
    return the final pressure profile (psi)."""
    # TODO: res = Reservoir1D(NX, LENGTH_FT, WIDTH_FT, HEIGHT_FT, k_array,
    #                         PHI, CT_PSI, MU_CP, PI_PSI, Bo=BO)
    # TODO: res.add_well(0, rate_stbd=INJ_RATE_STBD)
    # TODO: res.add_well(NX - 1, bhp_psi=PROD_BHP_PSI)
    # TODO: res.run(DAYS, DT_DAYS)
    # TODO: return res.P
    pass

def dp_across(p_psi, i0, i1):
    """Absolute pressure drop (psi) between cell indices i0 and i1."""
    # TODO: return abs(p_psi[i1] - p_psi[i0])
    pass

def harmonic_T(k_array, i, j):
    """Harmonic-mean transmissibility at the boundary between cells i and j."""
    # TODO: kh = 2 * k_array[i] * k_array[j] / (k_array[i] + k_array[j])
    # TODO: area = WIDTH_FT * HEIGHT_FT;  dx = LENGTH_FT / NX
    # TODO: return BETA * kh * area / (MU_CP * dx)
    pass

k_field = None       # TODO: build_k()
p_final = None       # TODO: run_hetero(k_field)
dp_streak_psi = None # TODO: dp_across(p_final, STREAK_START, STREAK_END)
dp_sand_psi = None   # TODO: dp_across(p_final, SAND_LO, SAND_HI)

print("pressure span (psi):", None if p_final is None else round(float(p_final.max() - p_final.min()), 1))
print("dP across streak (psi):", dp_streak_psi)
print("dP across equal sand span (psi):", dp_sand_psi)

visibilityReveal reference solutionexpand_more

Try solving it yourself first — the hints walk you through it. The solution below is one valid approach; yours may differ and still be correct.

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve


# ── Verified 1D single-phase implicit simulator (do not edit) ────────────
class Reservoir1D:
    """1D single-phase implicit simulator (harmonic-mean T, real area, Peaceman wells)."""
    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.area = width_ft * height_ft; self.vb = self.area * self.dx; self.height = height_ft
        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, float); self.Pi = Pi_psi; self.beta = 1.127e-3; self.wells = []
    def _T(self, i, j):
        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 _wi(self, w):
        re = 0.2*self.dx
        return 2*np.pi*self.beta*self.k[w["cell"]]*self.height/(self.mu*(np.log(re/w["rw"])+w["skin"]))
    def add_well(self, cell, rate_stbd=None, bhp_psi=None, rw_ft=0.354, skin=0.0):
        self.wells.append(dict(cell=cell, rate=rate_stbd, bhp=bhp_psi, rw=rw_ft, skin=skin))
    def step(self, dt):
        n = self.nx; C = self.vb*self.phi*self.ct/(5.615*dt)
        lo = np.zeros(n); di = np.zeros(n); up = np.zeros(n); b = np.zeros(n)
        for i in range(n):
            tl = self._T(i, i-1) if i > 0 else 0.0; tr = self._T(i, i+1) if i < n-1 else 0.0
            di[i] = C + tl + tr
            if i > 0: lo[i] = -tl
            if i < n-1: up[i] = -tr
            b[i] = C*self.P[i]
        for w in self.wells:
            i = w["cell"]
            if w["rate"] is not None: b[i] += w["rate"]*self.Bo
            elif w["bhp"] is not None: wi = self._wi(w); di[i] += wi; b[i] += wi*w["bhp"]
        self.P = spsolve(diags([lo[1:], di, up[:-1]], [-1, 0, 1], format="csc"), b)
    def run(self, days, dt):
        for _ in range(int(days/dt)): self.step(dt)


# ── OML 58 flow unit: a 200 md sand with a tight 10 md streak in the middle ─
NX = 50                  # grid cells along the 5,000 ft length
LENGTH_FT = 5000.0
WIDTH_FT = 500.0
HEIGHT_FT = 50.0
PHI = 0.20
CT_PSI = 1.0e-5
MU_CP = 1.2
PI_PSI = 3500.0          # initial reservoir pressure
BO = 1.2

K_SAND = 200.0           # clean sand permeability (md)
K_STREAK = 10.0          # tight low-perm streak permeability (md)
STREAK_START = 20        # first streak cell (inclusive)
STREAK_END = 25          # one past the last streak cell -> cells 20..24 (5 cells)

SAND_LO = 5              # an equal five-cell span of CLEAN sand, away from the streak
SAND_HI = 10

INJ_RATE_STBD = 600.0    # constant-rate water injector at cell 0
PROD_BHP_PSI = 2500.0    # constant-BHP producer at the last cell
DAYS = 1095.0            # 3 years of throughflow
DT_DAYS = 15.0

BETA = 1.127e-3          # transmissibility constant (field units)


def build_k():
    """Length-NX permeability field: K_SAND everywhere, K_STREAK in the streak cells."""
    k = np.full(NX, K_SAND, float)
    k[STREAK_START:STREAK_END] = K_STREAK
    return k


def run_hetero(k_array):
    """Run the line-drive (rate injector at cell 0, BHP producer at NX-1) and
    return the final pressure profile (psi)."""
    res = Reservoir1D(NX, LENGTH_FT, WIDTH_FT, HEIGHT_FT, k_array,
                      PHI, CT_PSI, MU_CP, PI_PSI, Bo=BO)
    res.add_well(0, rate_stbd=INJ_RATE_STBD)
    res.add_well(NX - 1, bhp_psi=PROD_BHP_PSI)
    res.run(DAYS, DT_DAYS)
    return res.P


def dp_across(p_psi, i0, i1):
    """Absolute pressure drop (psi) between cell indices i0 and i1."""
    return abs(p_psi[i1] - p_psi[i0])


def harmonic_T(k_array, i, j):
    """Harmonic-mean transmissibility at the boundary between cells i and j."""
    kh = 2 * k_array[i] * k_array[j] / (k_array[i] + k_array[j])
    area = WIDTH_FT * HEIGHT_FT
    dx = LENGTH_FT / NX
    return BETA * kh * area / (MU_CP * dx)


k_field = build_k()
p_final = run_hetero(k_field)
dp_streak_psi = dp_across(p_final, STREAK_START, STREAK_END)
dp_sand_psi = dp_across(p_final, SAND_LO, SAND_HI)

print("pressure span (psi):", round(float(p_final.max() - p_final.min()), 1))
print("dP across streak (psi):", round(float(dp_streak_psi), 1))
print("dP across equal sand span (psi):", round(float(dp_sand_psi), 1))

lockCopying code is a Full Access feature.

arrow_back

11.1 Grid-Refinement Study - When Does the Answer Stop Moving?

11.3 Injection-Production Pair - Pressure Support & Mass Balance

arrow_forward