OOIP Sensitivity Tornado

Level 3

Chapter 10: Material Balance

descriptionProblem

The OOIP estimate depends on several input parameters. Using the undersaturated reservoir from this chapter, vary each of the following independently by $\pm$ 20% and record the effect on estimated OOIP:

(a) $c_f$ , (b) $c_w$ , (c) $S_{wi}$ , (d) $B_{oi}$

Create a tornado chart: a horizontal bar chart showing how each parameter's uncertainty range affects the result, with the widest-impact parameter at the top. What does this tell you about where to focus your data collection efforts?

---

A single OOIP number from a Havlena-Odeh straight line looks reassuringly precise, but it is only as good as the inputs you fed it. Before you carry ooip_stb = 50 MMSTB into the field-development plan for the undersaturated block on OML 58, you need to know which inputs the answer actually leans on. That is what a tornado chart is for: vary each uncertain parameter over its plausible band, re-solve for OOIP, and rank the parameters by how far they move the result. The widest bar is where your money should go on data acquisition.

You already have the chapter's self-consistent undersaturated dataset wired up (it solves to ooip_stb ~= 50e6 STB through the origin with R^2 ~= 1). The embedded helpers calculate_mbe_terms(...) and ooip_through_origin(x, F) do the material-balance arithmetic. Do not reimplement them. Above the bubble point the gas term cancels (Rp - Rs = 0), so the underground withdrawal is simply F = Np * Bo, and the expansion driving the depletion is Eo + Efw.

The four parameters in play are the rock compressibility cf, the connate-water compressibility cw, the initial water saturation swi, and the initial oil FVF boi. Write:

ooip_for_params(cf, cw, swi, boi) -> N (a float). Recompute the

Havlena-Odeh OOIP for the canonical pressure/production history, but with these four inputs substituted in. Build the MBE terms, drop the t = 0 row (where Eo + Efw = 0), and fit ooip_through_origin(Eo + Efw, F). At the base inputs this must return ~= 50e6 STB.

sensitivity(base_params) -> dict mapping each of "cf", "cw",

"swi", "boi" to a tuple (N_low, N_high, swing), where N_low is OOIP with that one parameter set to 0.8 base (all others at base), N_high is OOIP at 1.2 base, and swing = abs(N_high - N_low). base_params is a dict with keys "cf", "cw", "swi", "boi".

Then build the tornado chart: a horizontal bar (ax.barh) per parameter, sorted so the largest swing sits at the top. Plot each bar spanning N_low -> N_high (a left=/width= bar, or center it on the base OOIP) so the reader sees both the size and the direction of each parameter's pull.

Store the sensitivity result in tornado and set base_params = {"cf": CF, "cw": CW, "swi": SWI, "boi": BOI}.

On this dataset the rock compressibility cf carries the widest swing, edging out boi, with swi and cw well behind. The reading is blunt: a +/-20% error in cf moves the booked OOIP by several MMSTB, so a special core-analysis program to pin down rock compressibility buys you more certainty than refining the water properties. Note too that swinging boi by a full +/-20% drives the FVF outside its physical band and the straight line stops making sense: a reminder that a tornado chart shows mathematical leverage, not a license to perturb every input over the same naive range.

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

# ── Material-balance helpers (DO NOT EDIT - copy verbatim) ───────────────
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

# ── Canonical undersaturated dataset (OML 58, above the bubble point) ────
# Self-consistent: F = N * (Eo + Efw) with N = 50,000,000 STB, R^2 ~= 1.0.
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

# Base inputs the tornado fans out from (+/-20% each).
CF, CW, SWI, BOI = cf, cw, Swi, Boi
PCT = 0.20

def ooip_for_params(cf, cw, swi, boi):
    """Havlena-Odeh OOIP (STB) for the canonical history with these inputs.

Build the MBE terms, drop the t=0 row (Eo + Efw = 0 there), then fit
    ooip_through_origin(Eo + Efw, F). At the base inputs this is ~= 50e6.
    """
    # TODO: F, Eo, _, Efw, _, _ = calculate_mbe_terms(
    # TODO:     P_psi, Np_stb, Gp_scf, Wp_stb, Bo, Bg, Rs, Bw,
    # TODO:     boi, Bgi, Rsi, swi, cw, cf, Pi)
    # TODO: x = (Eo + Efw)[1:]; Fv = F[1:]
    # TODO: N, _r2 = ooip_through_origin(x, Fv)
    # TODO: return float(N)
    pass

def sensitivity(base_params):
    """Map each param -> (N_low, N_high, swing) for a +/-20% one-at-a-time sweep.

base_params has keys "cf", "cw", "swi", "boi". For each, hold the other
    three at base, set this one to 0.8*base then 1.2*base, recompute OOIP,
    and record swing = abs(N_high - N_low).
    """
    # TODO: out = {}
    # TODO: for name in ("cf", "cw", "swi", "boi"):
    # TODO:     lo = dict(base_params); hi = dict(base_params)
    # TODO:     lo[name] = base_params[name] * (1.0 - PCT)
    # TODO:     hi[name] = base_params[name] * (1.0 + PCT)
    # TODO:     n_low  = ooip_for_params(lo["cf"], lo["cw"], lo["swi"], lo["boi"])
    # TODO:     n_high = ooip_for_params(hi["cf"], hi["cw"], hi["swi"], hi["boi"])
    # TODO:     out[name] = (n_low, n_high, abs(n_high - n_low))
    # TODO: return out
    pass

base_params = {"cf": CF, "cw": CW, "swi": SWI, "boi": BOI}
tornado = None  # TODO: sensitivity(base_params)

if tornado is None:
    print("tornado: None")
else:
    base_ooip = ooip_for_params(CF, CW, SWI, BOI)
    print(f"base OOIP = {base_ooip:,.0f} STB")
    order = sorted(tornado, key=lambda k: tornado[k][2], reverse=True)
    for k in order:
        lo, hi, sw = tornado[k]
        print(f"  {k:4s}: low={lo:,.0f}  high={hi:,.0f}  swing={sw:,.0f}")

# TODO: build the tornado chart - barh per param, widest swing at top.
    # fig, ax = plt.subplots(figsize=(9, 5))
    # for i, k in enumerate(reversed(order)):   # widest last -> top after barh
    #     lo, hi, sw = tornado[k]
    #     ax.barh(i, width=hi - lo, left=lo, color="#4682B4")
    # ax.set_yticks(range(len(order)))
    # ax.set_yticklabels(list(reversed(order)))
    # ax.axvline(base_ooip, color="k", lw=1, ls="--")
    # ax.set_xlabel("OOIP (STB)")
    # ax.set_title("OOIP sensitivity tornado - OML 58")
    # plt.tight_layout(); plt.show()

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


# ── Material-balance helpers (DO NOT EDIT - copy verbatim) ───────────────
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


# ── Canonical undersaturated dataset (OML 58, above the bubble point) ────
# Self-consistent: F = N * (Eo + Efw) with N = 50,000,000 STB, R^2 ~= 1.0.
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

# Base inputs the tornado fans out from (+/-20% each).
CF, CW, SWI, BOI = cf, cw, Swi, Boi
PCT = 0.20


def ooip_for_params(cf, cw, swi, boi):
    """Havlena-Odeh OOIP (STB) for the canonical history with these inputs."""
    F, Eo, _, Efw, _, _ = 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)[1:]          # drop t=0 (Eo + Efw = 0 there)
    Fv = F[1:]
    N, _r2 = ooip_through_origin(x, Fv)
    return float(N)


def sensitivity(base_params):
    """Map each param -> (N_low, N_high, swing) for a +/-20% one-at-a-time sweep."""
    out = {}
    for name in ("cf", "cw", "swi", "boi"):
        lo = dict(base_params)
        hi = dict(base_params)
        lo[name] = base_params[name] * (1.0 - PCT)
        hi[name] = base_params[name] * (1.0 + PCT)
        n_low = ooip_for_params(lo["cf"], lo["cw"], lo["swi"], lo["boi"])
        n_high = ooip_for_params(hi["cf"], hi["cw"], hi["swi"], hi["boi"])
        out[name] = (n_low, n_high, abs(n_high - n_low))
    return out


base_params = {"cf": CF, "cw": CW, "swi": SWI, "boi": BOI}
tornado = sensitivity(base_params)

base_ooip = ooip_for_params(CF, CW, SWI, BOI)
print(f"base OOIP = {base_ooip:,.0f} STB")
order = sorted(tornado, key=lambda k: tornado[k][2], reverse=True)
for k in order:
    lo, hi, sw = tornado[k]
    print(f"  {k:4s}: low={lo:,.0f}  high={hi:,.0f}  swing={sw:,.0f}")

# ── Tornado chart: barh per param, widest swing at the top ───────────────
fig, ax = plt.subplots(figsize=(9, 5))
for i, k in enumerate(reversed(order)):          # widest plotted last -> top
    lo, hi, sw = tornado[k]
    left = min(lo, hi)
    width = abs(hi - lo)
    ax.barh(i, width=width, left=left, color="#4682B4", edgecolor="k")
ax.set_yticks(range(len(order)))
ax.set_yticklabels(list(reversed(order)))
ax.axvline(base_ooip, color="k", lw=1.0, ls="--", label="base OOIP")
ax.set_xlabel("OOIP (STB)")
ax.set_title("OOIP sensitivity tornado - OML 58")
ax.legend(loc="lower right")
plt.tight_layout()
plt.show()

lockCopying code is a Full Access feature.

arrow_back

10.6 Schilthuis Water Influx

10.8 Campbell Plot - Flat Means Depletion

arrow_forward