Exercise 10.2
Undersaturated Reservoir OOIP (Havlena-Odeh)
# Undersaturated Reservoir OOIP: Havlena-Odeh straight line
Reservoir D-Block on OML 58 is a clean undersaturated oil reservoir: the bubble point sits at 2,800 psi and the whole pressure-production history runs above it (4,000 -> 3,000 psi). No free gas evolves, so the solution GOR holds at Rsi = 600 SCF/STB throughout and Gp = Np * Rsi. There is no water production and no aquifer support; this is a candidate depletion-drive reservoir, and management wants an independent OOIP before booking reserves.
The Havlena-Odeh method rearranges the general material-balance equation into a straight line. With no gas cap and no water influx, the only energy is fluid + rock expansion, so
F = N * (Eo + Efw)where F is the underground withdrawal (RB), Eo is the oil + dissolved-gas expansion, and Efw is the rock + connate-water expansion. Plot F (y) against Eo + Efw (x) and the points fall on a line through the origin whose slope is the OOIP N in STB. A tight fit (R^2 near 1) is your evidence that the model, pure depletion, is the right one.
> Watch the Bo trend. Above the bubble point the oil has no free gas to > lose, so as pressure drops the oil expands and Bo increases (it peaks > at the bubble point). That makes Eo = (Bo - Boi) + (Rsi - Rs) * Bg positive > and rising. If your Eo + Efw comes out negative, your Bo trend is upside > down.
## Initial conditions (D-Block)
| quantity | value |
|---|---|
Boi | 1.31000 RB/STB |
Bgi | 0.00087 RB/SCF |
Rsi | 600 SCF/STB |
Swi | 0.22 |
cw | 3.2e-6 psi⁻¹ |
cf | 5.0e-6 psi⁻¹ |
Bw | 1.02 RB/STB (constant) |
Pi | 4000 psi |
The pressure-production table (Wp = 0 everywhere, Bg = 0.00090 constant above Pb, Rs = 600 constant above Pb) is given in starter.py. Note that the t = 0 row has Np = 0, so its F and Eo + Efw are both zero; it sits exactly on the origin and carries no information for the fit. Skip it.
Two helpers are already provided in starter.py (copy them verbatim, do not rewrite them):
calculate_mbe_terms(...)->(F, Eo, Eg, Efw, Rp, dP)as arrays over all
pressure steps.
ooip_through_origin(x, F)->(N, r2), the least-squares slope through the
origin and its R^2.
Write:
estimate_ooip(P, Np, Gp, Wp, Bo, Bg, Rs, Bw, Boi, Bgi, Rsi, Swi, cw, cf, Pi)
-> (ooip_stb, r2)Inside it: call calculate_mbe_terms(...), build x = Eo + Efw, drop the first (t = 0) row from both x and F, then call ooip_through_origin(x_skip, F_skip) and return (float(N), float(r2)).
Call estimate_ooip(...) on the D-Block arrays and store the tuple in result.
You should find OOIP ≈ 50,000,000 STB with R² > 0.999, a near-perfect straight line through the origin. Because Rp - Rs = 0 above the bubble point, the gas term in F cancels and F collapses to Np * Bo; the slight upward expansion of the oil plus the rock/water squeeze is exactly enough to push out the produced barrels from a 50-MMSTB tank. The drive ratio F / (Eo + Efw) is essentially flat at N across every step: the textbook signature of a pure depletion-drive reservoir with no aquifer and no gas cap.
Stuck? Reveal hints one at a time — they progress from nudge to near-solution.
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 estimate_ooip(P, Np, Gp, Wp, Bo, Bg, Rs, Bw, Boi, Bgi, Rsi, Swi, cw, cf, Pi):
"""Havlena-Odeh OOIP for an undersaturated reservoir.
Returns (ooip_stb, r2): the slope of F vs (Eo + Efw) through the origin,
fitted on every pressure step EXCEPT the t = 0 row (which sits on the origin).
"""
F, Eo, Eg, Efw, Rp, dP = calculate_mbe_terms(
P, Np, Gp, Wp, Bo, Bg, Rs, Bw, Boi, Bgi, Rsi, Swi, cw, cf, Pi
)
x = Eo + Efw
N, r2 = ooip_through_origin(x[1:], F[1:]) # drop the t = 0 origin row
return float(N), float(r2)
result = estimate_ooip(
P_psi, Np_stb, Gp_scf, Wp_stb, Bo, Bg, Rs, Bw, Boi, Bgi, Rsi, Swi, cw, cf, Pi
)
ooip_stb, r2 = result
print(f"OOIP = {ooip_stb:,.0f} STB R^2 = {r2:.6f}")
lockCopying code is a Full Access feature.