Complete Material Balance Study

Level 3

Chapter 10: Material Balance

descriptionProblem

Combine everything from this chapter into a complete material balance study:

Load the data and calculate all GMBE terms
Perform Havlena-Odeh analysis and estimate OOIP
Diagnose the drive mechanism
Quantify uncertainty by estimating OOIP with $\pm$ 10% variations in $B_o$ , $c_f$ , and $S_{wi}$
Generate a one-page summary report with plots

Present your results as you would in a reserves review meeting.

---

This is the capstone of the chapter: the one call a reservoir engineer makes when fresh pressure-production data lands on the desk. Everything you built term by term: underground withdrawal F, the expansion drives Eo/Efw/Eg, the Havlena-Odeh straight line, the drive diagnosis, collapses into a single study function that returns the numbers you would read out loud in a reserves review.

The D-Block on OML 58 has produced for six pressure steps, all above the bubble point (an undersaturated reservoir). The pressure-production history and PVT are wired into starter.py. Because the reservoir is undersaturated, the oil expands as pressure drops, so Bo rises from Boi toward the bubble point. Eo = (Bo - Boi) + (Rsi - Rs) Bg is positive and growing. With no gas cap and no water influx, the underground withdrawal is balanced purely by oil-plus-rock expansion: F ≈ N (Eo + Efw).

Write one function:

material_balance_study(P_psi, Np_stb, Bo, Rs, Bg, Bw, Wp_stb, Gp_scf,
                       Boi, Bgi, Rsi, Swi, cw, cf, Pi)

It should call calculate_mbe_terms(...) (already provided), fit the Havlena-Odeh line through the origin with ooip_through_origin(...) on Eo + Efw vs F (dropping the t = 0 row, which sits exactly on the origin) and return a dict:

key	value
`ooip_stb`	OOIP, the slope `N` of `F` vs `(Eo + Efw)` through the origin
`r2`	the fit `R²` from `ooip_through_origin`
`drive`	`'depletion'`, `'water drive'`, or `'gas cap'`
`rf`	recovery factor so far: `Np_stb[-1] / ooip_stb`
`withdrawal_rb`	total reservoir-barrel withdrawal to date: `F[-1]`

Diagnosing the drive. Form the drive plot y = F / (Eo + Efw) over the producing steps (skip t = 0). For a pure depletion (volumetric) drive this ratio is flat; it just equals N at every step, because there is no extra energy from a gas cap or aquifer. A rising trend means an external drive is feeding the reservoir. A simple, robust rule: compare the last value to the first.

If y[-1] / y[0] is within a few percent of 1 (flat) -> 'depletion'.
If it rises strongly -> 'water drive'; a milder rise -> 'gas cap'.

The threshold is yours to pick, but on the D-Block the ratio is essentially 1.0 (a flat line) so the study must report drive == 'depletion'.

For the D-Block you should find OOIP ≈ 50,000,000 STB with an essentially perfect fit (R² > 0.999), a recovery factor near 1.9 % after this much depletion, and a depletion-drive diagnosis. Call material_balance_study(...) on the wired-in data and store the result in study.

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
import pandas as pd

# ── D-Block pressure-production history (OML 58), all above the bubble point ──
P_psi  = np.array([4000.0, 3800.0, 3600.0, 3400.0, 3200.0, 3000.0])
Np_stb = np.array([0.0, 192666.0, 384411.0, 575243.0, 765167.0, 954191.0])
Bo     = np.array([1.31000, 1.31314, 1.31629, 1.31943, 1.32258, 1.32572])  # RISES as P drops
Rs     = np.full(6, 600.0)            # constant above the bubble point
Bg     = np.full(6, 0.00090)          # cancels above Pb (Rp - Rs = 0); present for the signature
Bw     = 1.02
Wp_stb = np.zeros(6)
Gp_scf = Np_stb * 600.0               # Gp = Np * Rsi above Pb

# Initial conditions.
Boi, Bgi, Rsi, Swi, cw, cf, Pi = 1.31000, 0.00087, 600.0, 0.22, 3.2e-6, 5.0e-6, 4000.0

def calculate_mbe_terms(P, Np, Gp, Wp, Bo, Bg, Rs, Bw, Boi, Bgi, Rsi, Swi, cw, cf, Pi):
    P = np.asarray(P, float); Np = np.asarray(Np, float); Gp = np.asarray(Gp, float)
    Wp = np.asarray(Wp, float); Bo = np.asarray(Bo, float); Bg = np.asarray(Bg, float); Rs = np.asarray(Rs, float)
    dP = Pi - P
    Rp = np.where(Np > 0, Gp / Np, 0.0)
    F   = Np * (Bo + (Rp - Rs) * Bg) + Wp * Bw           # underground withdrawal (RB)
    Eo  = (Bo - Boi) + (Rsi - Rs) * Bg                   # oil + dissolved-gas expansion
    Efw = Boi * (cw * Swi + cf) * dP / (1.0 - Swi)       # rock + connate-water expansion
    Eg  = Boi / Bgi * (Bg - Bgi)                          # gas-cap expansion (per unit m)
    return F, Eo, Eg, Efw, Rp, dP

def ooip_through_origin(x, F):
    x = np.asarray(x, float); F = np.asarray(F, float)
    N = np.sum(x * F) / np.sum(x * x)                     # least-squares slope through origin
    ss_res = np.sum((F - N * x) ** 2)
    ss_tot = np.sum((F - np.mean(F)) ** 2)
    r2 = 1.0 - ss_res / ss_tot
    return N, r2

def material_balance_study(P_psi, Np_stb, Bo, Rs, Bg, Bw, Wp_stb, Gp_scf,
                           Boi, Bgi, Rsi, Swi, cw, cf, Pi):
    """One-call reservoir study on fresh pressure-production data.

Returns a dict with ooip_stb, r2, drive, rf, withdrawal_rb.
    """
    # TODO: F, Eo, Eg, Efw, Rp, dP = calculate_mbe_terms(
    #           P_psi, Np_stb, Gp_scf, Wp_stb, Bo, Bg, Rs, Bw,
    #           Boi, Bgi, Rsi, Swi, cw, cf, Pi)
    # TODO: x = Eo + Efw
    # TODO: fit Havlena-Odeh through the origin on the producing steps (drop t = 0):
    #       ooip_stb, r2 = ooip_through_origin(x[1:], F[1:])
    # TODO: drive plot y = F[1:] / x[1:]; if y[-1] / y[0] is within ~3% of 1 -> 'depletion'
    # TODO: rf = Np_stb[-1] / ooip_stb
    # TODO: withdrawal_rb = F[-1]
    # TODO: return {"ooip_stb": float(ooip_stb), "r2": float(r2), "drive": drive,
    #               "rf": float(rf), "withdrawal_rb": float(withdrawal_rb)}
    pass

study = None  # TODO: material_balance_study(P_psi, Np_stb, Bo, Rs, Bg, Bw, Wp_stb, Gp_scf,
#                                            Boi, Bgi, Rsi, Swi, cw, cf, Pi)

if study is None:
    print("study: None")
else:
    print(f"OOIP        = {study['ooip_stb']:,.0f} STB")
    print(f"R^2         = {study['r2']:.6f}")
    print(f"drive       = {study['drive']}")
    print(f"RF          = {study['rf'] * 100:.2f} %")
    print(f"withdrawal  = {study['withdrawal_rb']:,.0f} RB")

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
import pandas as pd


# ── D-Block pressure-production history (OML 58), all above the bubble point ──
P_psi  = np.array([4000.0, 3800.0, 3600.0, 3400.0, 3200.0, 3000.0])
Np_stb = np.array([0.0, 192666.0, 384411.0, 575243.0, 765167.0, 954191.0])
Bo     = np.array([1.31000, 1.31314, 1.31629, 1.31943, 1.32258, 1.32572])  # RISES as P drops
Rs     = np.full(6, 600.0)            # constant above the bubble point
Bg     = np.full(6, 0.00090)          # cancels above Pb (Rp - Rs = 0); present for the signature
Bw     = 1.02
Wp_stb = np.zeros(6)
Gp_scf = Np_stb * 600.0               # Gp = Np * Rsi above Pb

# Initial conditions.
Boi, Bgi, Rsi, Swi, cw, cf, Pi = 1.31000, 0.00087, 600.0, 0.22, 3.2e-6, 5.0e-6, 4000.0


def calculate_mbe_terms(P, Np, Gp, Wp, Bo, Bg, Rs, Bw, Boi, Bgi, Rsi, Swi, cw, cf, Pi):
    P = np.asarray(P, float); Np = np.asarray(Np, float); Gp = np.asarray(Gp, float)
    Wp = np.asarray(Wp, float); Bo = np.asarray(Bo, float); Bg = np.asarray(Bg, float); Rs = np.asarray(Rs, float)
    dP = Pi - P
    Rp = np.where(Np > 0, Gp / Np, 0.0)
    F   = Np * (Bo + (Rp - Rs) * Bg) + Wp * Bw           # underground withdrawal (RB)
    Eo  = (Bo - Boi) + (Rsi - Rs) * Bg                   # oil + dissolved-gas expansion
    Efw = Boi * (cw * Swi + cf) * dP / (1.0 - Swi)       # rock + connate-water expansion
    Eg  = Boi / Bgi * (Bg - Bgi)                          # gas-cap expansion (per unit m)
    return F, Eo, Eg, Efw, Rp, dP


def ooip_through_origin(x, F):
    x = np.asarray(x, float); F = np.asarray(F, float)
    N = np.sum(x * F) / np.sum(x * x)                     # least-squares slope through origin
    ss_res = np.sum((F - N * x) ** 2)
    ss_tot = np.sum((F - np.mean(F)) ** 2)
    r2 = 1.0 - ss_res / ss_tot
    return N, r2


def material_balance_study(P_psi, Np_stb, Bo, Rs, Bg, Bw, Wp_stb, Gp_scf,
                           Boi, Bgi, Rsi, Swi, cw, cf, Pi):
    """One-call reservoir study on fresh pressure-production data.

    Returns a dict with ooip_stb, r2, drive, rf, withdrawal_rb.
    """
    F, Eo, Eg, Efw, Rp, dP = calculate_mbe_terms(
        P_psi, Np_stb, Gp_scf, Wp_stb, Bo, Bg, Rs, Bw,
        Boi, Bgi, Rsi, Swi, cw, cf, Pi
    )
    x = Eo + Efw

    # Havlena-Odeh straight line through the origin, dropping the t = 0 row.
    ooip_stb, r2 = ooip_through_origin(x[1:], F[1:])

    # Drive plot F/(Eo+Efw): flat -> depletion; rising -> external drive.
    y = F[1:] / x[1:]
    rise = y[-1] / y[0]
    if rise <= 1.03:
        drive = "depletion"
    elif rise > 1.15:
        drive = "water drive"
    else:
        drive = "gas cap"

    Np_last = np.asarray(Np_stb, float)[-1]
    rf = Np_last / ooip_stb
    withdrawal_rb = float(F[-1])

    return {
        "ooip_stb": float(ooip_stb),
        "r2": float(r2),
        "drive": drive,
        "rf": float(rf),
        "withdrawal_rb": float(withdrawal_rb),
    }


study = material_balance_study(
    P_psi, Np_stb, Bo, Rs, Bg, Bw, Wp_stb, Gp_scf,
    Boi, Bgi, Rsi, Swi, cw, cf, Pi
)

print(f"OOIP        = {study['ooip_stb']:,.0f} STB")
print(f"R^2         = {study['r2']:.6f}")
print(f"drive       = {study['drive']}")
print(f"RF          = {study['rf'] * 100:.2f} %")
print(f"withdrawal  = {study['withdrawal_rb']:,.0f} RB")

lockCopying code is a Full Access feature.

arrow_back

9.9 DCA vs Volumetric Reserves

10.1 MBE Term Identification

arrow_forward