Exercise 10.4
Solution Gas Drive Below the Bubble Point
# Solution Gas Drive Below the Bubble Point
Above the bubble point the material-balance arithmetic is forgiving: the gas stays dissolved, the produced gas-oil ratio Rp equals the solution GOR Rs, and the (Rp - Rs) Bg term in the underground withdrawal cancels. You can get away with F = Np Bo and still close the balance.
Cross below the bubble point and that shortcut quietly poisons your reserves estimate. On OML 58 the NEMBE-12 block has fallen below its bubble point of Pb = 2,800 psi. Free gas now evolves out of solution: Rs drops, the gas formation-volume factor Bg grows, and Rp climbs above Rs as that liberated gas is produced at surface. Ignore the gas terms and you understate the underground withdrawal F, then back out an OOIP that is too small.
## The below-bubble-point PVT and production history
| P (psi) | Np (STB) | Bo (RB/STB) | Rs (SCF/STB) | Bg (RB/SCF) | Rp (SCF/STB) |
|---|---|---|---|---|---|
| 2800 | 1,100,000 | 1.330 | 600 | 0.00100 | 600 |
| 2600 | 1,300,000 | 1.310 | 540 | 0.00120 | 650 |
| 2400 | 1,500,000 | 1.290 | 480 | 0.00150 | 720 |
Bo now decreases as pressure drops (gas leaving the oil shrinks it), Rs falls, Bg rises, and Rp rises above Rs. Cumulative produced gas is Gp = Rp * Np, water production Wp = 0. The initial properties are the same as the canonical undersaturated study: Boi = 1.31000, Bgi = 0.00087, Rsi = 600, Swi = 0.22, cw = 3.2e-6, cf = 5.0e-6, Pi = 4000.
## Your task
The helper calculate_mbe_terms(...) (already provided, do not rewrite it) is the same MBE engine from the canonical study. The key terms it returns:
F = Np * (Bo + (Rp - Rs) * Bg) + Wp * Bw # underground withdrawal (RB)
Eo = (Bo - Boi) + (Rsi - Rs) * Bg # oil + dissolved-gas expansionWrite two functions:
below_pb_terms(P, Np, Gp, Wp, Bo, Bg, Rs, Bw): callcalculate_mbe_terms
with the module-level initial properties and return the tuple (F, Eo) (just the first two outputs). Because Rp = Gp / Np > Rs below Pb, the (Rp - Rs) Bg term is now non-zero; F includes the produced free gas, and Eo includes the liberated-gas term (Rsi - Rs) Bg.
naive_F(Np, Bo): the wrong shortcut,Np * Bo. This drops the gas
term entirely.
Then, on the dataset above:
- Call
below_pb_terms(...)and storeF_fullandEo(NumPy arrays). - Compute
F_naive = naive_F(Np_stb, Bo). - Store the gap
F_gap = F_full - F_naive: the reserves the naive balance
throws away.
You should find that at the bubble point (2,800 psi) Rp = Rs = 600, so the gas term cancels and F_gap = 0. As gas evolves, the gap opens up: 171,600 RB at 2,600 psi and 540,000 RB at 2,400 psi. Meanwhile Eo climbs from 0.020 to 0.160 RB/STB as the liberated-gas term grows; every term stays positive. Ignoring the liberated and produced gas below Pb corrupts the balance: the missing 540,000 RB at 2,400 psi alone is more than a third of a million reservoir barrels of withdrawal the naive method never sees.
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
# ── MBE engine (provided - do NOT modify) ────────────────────────────────
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
# ── NEMBE-12 below-bubble-point history (OML 58), Pb = 2800 psi ───────────
P_psi = np.array([2800.0, 2600.0, 2400.0])
Np_stb = np.array([1_100_000.0, 1_300_000.0, 1_500_000.0])
Bo = np.array([1.330, 1.310, 1.290]) # SHRINKS as gas leaves the oil
Rs = np.array([600.0, 540.0, 480.0]) # solution GOR drops below Pb
Bg = np.array([0.00100, 0.00120, 0.00150]) # free-gas FVF grows
Rp_gor = np.array([600.0, 650.0, 720.0]) # produced GOR climbs above Rs
Gp_scf = Rp_gor * Np_stb # Gp = Rp * Np
Bw = 1.02
Wp_stb = np.zeros(3)
# Initial properties - same as the canonical undersaturated study.
Boi, Bgi, Rsi, Swi, cw, cf, Pi = 1.31000, 0.00087, 600.0, 0.22, 3.2e-6, 5.0e-6, 4000.0
def below_pb_terms(P, Np, Gp, Wp, Bo, Bg, Rs, Bw):
"""Return (F, Eo) from calculate_mbe_terms using the module initial props."""
F, Eo, Eg, Efw, Rp, dP = calculate_mbe_terms(
P, Np, Gp, Wp, Bo, Bg, Rs, Bw, Boi, Bgi, Rsi, Swi, cw, cf, Pi)
return F, Eo
def naive_F(Np, Bo):
"""The WRONG shortcut: drop the gas term and use F = Np * Bo."""
return np.asarray(Np, float) * np.asarray(Bo, float)
F_full, Eo = below_pb_terms(P_psi, Np_stb, Gp_scf, Wp_stb, Bo, Bg, Rs, Bw)
F_naive = naive_F(Np_stb, Bo)
F_gap = F_full - F_naive
for i, p in enumerate(P_psi):
print(f"P={p:.0f} F_full={F_full[i]:,.0f} F_naive={F_naive[i]:,.0f} "
f"gap={F_gap[i]:,.0f} Eo={Eo[i]:.4f}")
lockCopying code is a Full Access feature.