Campbell Plot - Flat Means Depletion

Level 2

Chapter 10: Material Balance

descriptionProblem

A Campbell plot graphs $F/E_t$ vs. cumulative production $N_p$ . If OOIP is correct, the plot should be flat (horizontal). A rising trend suggests water influx or gas cap expansion not accounted for. Generate the Campbell plot for the data in this chapter. Is the interpretation consistent with the Havlena-Odeh plot?

---

The Campbell plot is the material-balance engineer's lie detector. It graphs the underground withdrawal divided by the total expansion, F / Et, against cumulative oil produced Np. The whole trick rests on the general MBE for a reservoir with no gas cap and no aquifer:

F = N * Et        where   Et = Eo + Efw

If your OOIP N is right and the only energy in the tank is fluid and rock expansion (pure depletion drive), then F / Et must equal N at every pressure step: a dead-flat horizontal line. The plot is reading N off the y-axis for you.

A rising F / Et trend is the tell-tale of energy you have not modelled: water encroaching from an aquifer (We) or a gas cap expanding both add to F without showing up in Et, so the ratio climbs as Np grows. The Campbell plot turns "is my OOIP and drive assumption consistent?" into "is this line flat?"

Reservoir Otumara-4 on OML 58 is undersaturated (pressure above the bubble point). Six pressure surveys are tabulated in starter.py. The calculate_mbe_terms(...) helper returns F, Eo, Efw (and more); copy it verbatim.

Write campbell(P, Np, Gp, Wp, Bo, Bg, Rs, Bw, Boi, Bgi, Rsi, Swi, cw, cf, Pi) that:

calls calculate_mbe_terms(...) to get F, Eo, Efw,
forms Et = Eo + Efw,
drops the first survey (t = 0, where F = 0 and Et = 0, so F/Et is

undefined),

returns a tuple (np_array, f_over_et_array): the producing Np values and

the matching F / Et ratios.

Then build a matplotlib figure that plots f_over_et vs np_array (a scatter or line) plus a horizontal reference line drawn at the mean of f_over_et (use ax.axhline(...)), so a flat profile sits right on the reference.

Call campbell(...) on the Otumara-4 data and store the results in np_array and f_over_et. You should find F / Et essentially constant at about 50,000,000 STB; its spread (max minus min) is well under 1% of the mean. Flat line, one OOIP, no surprises: Otumara-4 is on depletion drive and the volumetric OOIP is consistent. If this line had ramped upward across Np, you would go hunting for the aquifer you forgot to put in the model.

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
import matplotlib.pyplot as plt

# ── Otumara-4 (OML 58) - undersaturated, six pressure surveys ────────────
# Pressure stays ABOVE the bubble point, so the oil keeps expanding as P
# drops: Bo RISES toward Pb, and the dissolved gas stays in solution
# (Rs constant). The gas term in F cancels because Rp - Rs = 0 here.
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
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 campbell(P, Np, Gp, Wp, Bo, Bg, Rs, Bw, Boi, Bgi, Rsi, Swi, cw, cf, Pi):
    """Campbell diagnostic: return (np_array, f_over_et) skipping t=0.

F / Et should be ~flat (≈ OOIP) for a pure-depletion reservoir.
    """
    # TODO: F, Eo, Eg, Efw, Rp, dP = calculate_mbe_terms(P, Np, Gp, Wp, Bo, Bg,
    #                                                     Rs, Bw, Boi, Bgi, Rsi,
    #                                                     Swi, cw, cf, Pi)
    # TODO: Et = Eo + Efw
    # TODO: drop the first survey (t = 0, where F = Et = 0)
    #       np_prod  = np.asarray(Np, float)[1:]
    #       f_over_et = (F / Et)[1:]
    # TODO: return np_prod, f_over_et
    pass

np_array = None       # TODO: campbell(...)[0]
f_over_et = None      # TODO: campbell(...)[1]

# TODO: build a figure - scatter/line of f_over_et vs np_array, plus a
#       horizontal reference at the mean via ax.axhline(...).
# fig, ax = plt.subplots(figsize=(9, 5))
# ax.plot(np_array, f_over_et, "o-", label="F / Et")
# ax.axhline(np.mean(f_over_et), color="firebrick", linestyle="--", label="mean (≈ OOIP)")
# ax.set_xlabel("Cumulative oil produced  Np (STB)")
# ax.set_ylabel("F / Et  (STB)")
# ax.set_title("Otumara-4 Campbell plot - flat = depletion drive")
# ax.legend(); ax.grid(True, alpha=0.3)
# plt.show()

if np_array is None or f_over_et is None:
    print("campbell() not implemented yet")
else:
    print("Np (STB)        F/Et (STB)")
    for n, r in zip(np.asarray(np_array), np.asarray(f_over_et)):
        print(f"{n:>12,.0f}  {r:>14,.0f}")
    mean_ooip = float(np.mean(f_over_et))
    spread = float(np.max(f_over_et) - np.min(f_over_et))
    print(f"\nmean F/Et   = {mean_ooip:,.0f} STB  (≈ OOIP)")
    print(f"spread      = {spread:,.0f} STB  ({100 * spread / mean_ooip:.3f}% of mean)")

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
import matplotlib.pyplot as plt


# ── Otumara-4 (OML 58) - undersaturated, six pressure surveys ────────────
# Pressure stays ABOVE the bubble point, so the oil keeps expanding as P
# drops: Bo RISES toward Pb, and the dissolved gas stays in solution
# (Rs constant). The gas term in F cancels because Rp - Rs = 0 here.
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
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 campbell(P, Np, Gp, Wp, Bo, Bg, Rs, Bw, Boi, Bgi, Rsi, Swi, cw, cf, Pi):
    """Campbell diagnostic: return (np_array, f_over_et) skipping t=0.

    F / Et should be ~flat (≈ OOIP) for a pure-depletion reservoir.
    """
    F, Eo, Eg, Efw, Rp, dP = calculate_mbe_terms(
        P, Np, Gp, Wp, Bo, Bg, Rs, Bw, Boi, Bgi, Rsi, Swi, cw, cf, Pi
    )
    Et = Eo + Efw
    np_prod = np.asarray(Np, float)[1:]          # drop t = 0 (F = Et = 0)
    f_over_et = (F / Et)[1:]
    return np_prod, f_over_et


np_array, f_over_et = campbell(
    P_psi, Np_stb, Gp_scf, Wp_stb, Bo, Bg, Rs, Bw,
    Boi, Bgi, Rsi, Swi, cw, cf, Pi,
)

fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(np_array, f_over_et, "o-", color="steelblue", label="F / Et")
ax.axhline(float(np.mean(f_over_et)), color="firebrick", linestyle="--",
           label="mean (≈ OOIP)")
ax.set_xlabel("Cumulative oil produced  Np (STB)")
ax.set_ylabel("F / Et  (STB)")
ax.set_title("Otumara-4 Campbell plot - flat = depletion drive")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

print("Np (STB)        F/Et (STB)")
for n, r in zip(np.asarray(np_array), np.asarray(f_over_et)):
    print(f"{n:>12,.0f}  {r:>14,.0f}")
mean_ooip = float(np.mean(f_over_et))
spread = float(np.max(f_over_et) - np.min(f_over_et))
print(f"\nmean F/Et   = {mean_ooip:,.0f} STB  (≈ OOIP)")
print(f"spread      = {spread:,.0f} STB  ({100 * spread / mean_ooip:.3f}% of mean)")

lockCopying code is a Full Access feature.

arrow_back

10.7 OOIP Sensitivity Tornado

10.9 Recovery Factor So Far

arrow_forward