Exercise 10.6
Schilthuis Water Influx
Implement the Schilthuis water influx model for an aquifer with influx constant = 150 RB/psi/year. Using the pressure history from Exercise 10.2 and assuming annual time steps, calculate the cumulative water influx at each pressure step. Then re-estimate OOIP accounting for this influx.
---
The OML 58 block from this chapter looked like a textbook depletion drive: plot the underground withdrawal F against the expansion term Eo + Efw and the points fall on a straight line through the origin whose slope is the OOIP; about 50 million STB, with an R² of essentially 1.0. Clean. Flat Campbell plot. No drama.
But suppose the block is not fully closed. A modest aquifer is leaning on the reservoir, pushing water in as pressure drops. That water influx We is energy the reservoir did not have to spend from its own expansion, and if you leave it out of the balance, the aquifer's support masquerades as extra oil. The straight line through the origin tilts, and the slope you read off as "OOIP" is too big. You book barrels that aren't there.
Schilthuis gave the simplest aquifer model: steady-state influx, where the water enters at a rate proportional to the pressure drawdown (Pi - P). Over annual time steps (dt = 1 yr) the cumulative influx is
We(t) = k * Σ (Pi - P_i) * dt # summed over all steps up to time twith influx constant k = 150 RB/psi/yr. Pressure has fallen 200 psi each year from the initial 4,000 psi, so the drawdown (and therefore the influx rate) grows year over year.
Write two functions.
1. cumulative_we(P_psi, Pi, k, dt=1.0) -> a NumPy array of cumulative water influx We (RB), one entry per pressure step. Each entry is k times the running cumulative sum of (Pi - P) * dt. At the initial pressure the drawdown is zero, so We[0] = 0, and because every later drawdown is positive, We is monotonically increasing.
2. ooip_with_influx(F, Eo, Efw, We) -> the OOIP (STB) when influx is accounted for. Subtract the influx from the withdrawal first, then regress the corrected withdrawal (F - We) against (Eo + Efw) through the origin using the provided ooip_through_origin helper. (Skip the t = 0 row, where every term is zero.) Return just the slope N.
The MBE terms (F, Eo, Efw) and the ooip_through_origin helper are given to you, already computed from the canonical undersaturated dataset. Above the bubble point the gas term cancels, so F reduces to Np * Bo, Eo is positive and rises as pressure drops, and the no-influx OOIP comes out near 50 MMSTB.
Compute both numbers and store them:
we_rb=cumulative_we(P_psi, Pi, K_RB_PSI_YR)ooip_no_influx_stb=ooip_through_origin((Eo + Efw)[1:], F[1:])[0]ooip_with_influx_stb=ooip_with_influx(F, Eo, Efw, we_rb)
You should find the no-influx OOIP near 50.0 MMSTB and the influx-corrected OOIP near 34.9 MMSTB, roughly 30% lower. Subtracting a positive We lowers the apparent withdrawal at every step, so the line through the origin is less steep and the honest OOIP is smaller. The lesson: a reservoir with hidden aquifer support will over-state its oil in place if you treat it as a closed tank. The 15 million "extra" barrels were really the aquifer doing the work.
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 machinery (given - do not edit) ──────────────────────────────────
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
# ── OML 58 undersaturated reservoir - annual pressure history ────────────
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
# Schilthuis steady-state aquifer: influx constant and annual time step.
K_RB_PSI_YR = 150.0
DT_YR = 1.0
# MBE terms, already computed for you.
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,
)
def cumulative_we(P_psi, Pi, k, dt=1.0):
"""Schilthuis cumulative water influx We (RB) at each pressure step.
We(t) = k * running cumulative sum of (Pi - P) * dt.
We[0] = 0 (zero drawdown at Pi); strictly increasing thereafter.
"""
drawdown = Pi - np.asarray(P_psi, float)
return k * np.cumsum(drawdown * dt)
def ooip_with_influx(F, Eo, Efw, We):
"""OOIP (STB) with aquifer support removed from the withdrawal.
Regress the influx-corrected withdrawal (F - We) against (Eo + Efw)
through the origin (skipping the t=0 row) and return the slope N.
"""
x = (np.asarray(Eo, float) + np.asarray(Efw, float))[1:]
f_corr = (np.asarray(F, float) - np.asarray(We, float))[1:]
N, _ = ooip_through_origin(x, f_corr)
return N
we_rb = cumulative_we(P_psi, Pi, K_RB_PSI_YR)
ooip_no_influx_stb = ooip_through_origin((Eo + Efw)[1:], F[1:])[0]
ooip_with_influx_stb = ooip_with_influx(F, Eo, Efw, we_rb)
print("We (RB) by year:", np.round(we_rb, 0))
print(f"OOIP, no influx = {ooip_no_influx_stb:,.0f} STB")
print(f"OOIP, with influx = {ooip_with_influx_stb:,.0f} STB")
drop = 100.0 * (1.0 - ooip_with_influx_stb / ooip_no_influx_stb)
print(f"ignoring We overstates OOIP by {drop:,.1f}%")
lockCopying code is a Full Access feature.