Exercise 15.3
Gas Material Balance - P/Z Plot & Original Gas In Place
For a volumetric gas reservoir (no water influx), the P/Z plot is linear: P/Z = (Pi/Zi)(1 - Gp/G). Given initial pressure 4,500 psia, temperature 220°F, and the following production data, estimate the original gas in place (G) and remaining reserves:
| Gp (Bscf) | P (psia) |
|---|---|
| 0 | 4,500 |
| 5 | 4,100 |
| 12 | 3,600 |
| 20 | 3,050 |
| 30 | 2,400 |
---
The classic depletion diagnostic for a volumetric gas reservoir (no water influx, no rock compaction worth modelling) is the P/Z plot. Plot P/Z against cumulative gas produced Gp and the points fall on a straight line:
The intercept is Pi/Zi and the line hits zero P/Z exactly when Gp = G, the original gas in place (OGIP). So OGIP is just -intercept / slope of the fitted line: extrapolate the depletion trend to abandonment and read off how much gas the reservoir ever held.
The verified Sutton pseudo-critical correlation pseudo_critical_properties and the Dranchuk-Abou-Kassem Z-factor solver z_factor_DAK are embedded for you; do not modify them or re-derive the physics.
Write gas_ogip(Gp_bscf, P_psia, T_F, gamma_g) returning a dict:
Ppc, Tpc = pseudo_critical_properties(gamma_g)and
Tpr = (T_F + 459.67) / Tpc (reservoir temperature is isothermal, so Tpr is the same for every point).
- For each pressure
PinP_psia, compute the reduced pressure
Ppr = P / Ppc, get Z = z_factor_DAK(Ppr, Tpr), and form P / Z.
- Fit a straight line of
P/ZversusGpwithnumpy.polyfit(Gp_bscf, PZ, 1),
which returns (slope, intercept). The slope is negative (depletion); the intercept is Pi/Zi.
- OGIP is
G = -intercept / slope(Bscf).
Return a dict with EXACTLY these keys:
'G_bscf': OGIPG(float, Bscf)'PiZi': the fitted interceptPi/Zi(float)'slope': the fitted slope (float, negative)'PZ': the list ofP/Zvalues, one per input pressure
Embedded constants and the call the tests read:
GP = [0.0, 5.0, 12.0, 20.0, 30.0] # cumulative produced, Bscf
P = [4500.0, 4100.0, 3600.0, 3050.0, 2400.0] # reservoir pressure, psia
T_F = 220.0 # reservoir temperature, deg F
GAMMA_G = 0.65 # gas specific gravity (air = 1)
mb = gas_ogip(GP, P, T_F, GAMMA_G)
ogip_bscf = mb['G_bscf']
remaining_bscf = ogip_bscf - GP[-1] # reserves left after the last Gp> Think about it: the line must slope downward; produce gas, pressure > drops, so P/Z falls. OGIP has to come out larger than anything you have > already produced (GP[-1] = 30 Bscf), and the remaining reserves are simply > what is left between today's Gp and the extrapolated G. Why does using > P/Z instead of raw P make the depletion trend a straight line you can > safely extrapolate?
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
from scipy.optimize import brentq
# ── Verified Sutton pseudo-critical correlation (do not edit) ────────────
def pseudo_critical_properties(gas_sg):
"""
Estimate pseudo-critical properties from gas specific gravity.
Uses the Sutton (1985) correlations.
Parameters
----------
gas_sg : float
Gas specific gravity (air = 1.0). Typical range: 0.55 to 0.90.
Returns
-------
tuple
(Ppc in psia, Tpc in deg R)
"""
Ppc = 756.8 - 131.0 * gas_sg - 3.6 * gas_sg**2
Tpc = 169.2 + 349.5 * gas_sg - 74.0 * gas_sg**2
return Ppc, Tpc
# ── Verified Dranchuk-Abou-Kassem Z-factor solver (do not edit) ──────────
def z_factor_DAK(Ppr, Tpr):
"""
Calculate gas Z-factor using the Dranchuk-Abou-Kassem (1975)
equation of state.
This correlation fits the Standing-Katz chart with 11 empirical
coefficients. Valid for 1.0 ≤ Tpr ≤ 3.0 and 0.2 ≤ Ppr ≤ 15.
Parameters
----------
Ppr : float
Pseudo-reduced pressure (dimensionless).
Tpr : float
Pseudo-reduced temperature (dimensionless).
Returns
-------
float
Compressibility factor Z.
"""
# DAK (1975) coefficients
A1, A2, A3 = 0.3265, -1.0700, -0.5339
A4, A5, A6 = 0.01569, -0.05165, 0.5475
A7, A8 = -0.7361, 0.1844
A9, A10, A11 = 0.1056, 0.6134, 0.7210
def residual(rho_r):
"""DAK residual function. Z = 0.27 * Ppr / (rho_r * Tpr)"""
Z = 0.27 * Ppr / (rho_r * Tpr)
F = (
A1 + A2/Tpr + A3/Tpr**3 + A4/Tpr**4 + A5/Tpr**5
) * rho_r + (
A6 + A7/Tpr + A8/Tpr**2
) * rho_r**2 - A9 * (
A7/Tpr + A8/Tpr**2
) * rho_r**5 + A10 * (
1 + A11*rho_r**2
) * (rho_r**2 / Tpr**3) * np.exp(-A11*rho_r**2) + 1 - Z
return F
# Solve for reduced density
rho_r = brentq(residual, 0.001, 5.0)
Z = 0.27 * Ppr / (rho_r * Tpr)
return Z
# ── Volumetric gas reservoir data (do not edit) ──────────────────────────
GP = [0.0, 5.0, 12.0, 20.0, 30.0] # cumulative produced, Bscf
P = [4500.0, 4100.0, 3600.0, 3050.0, 2400.0] # reservoir pressure, psia
T_F = 220.0 # reservoir temperature, deg F
GAMMA_G = 0.65 # gas specific gravity (air = 1)
def gas_ogip(Gp_bscf, P_psia, T_F, gamma_g):
"""Estimate OGIP from a P/Z material-balance plot (volumetric reservoir).
P/Z = (Pi/Zi)(1 - Gp/G) is a straight line; OGIP is where it hits zero.
Returns a dict with keys:
'G_bscf' = OGIP G (Bscf) = -intercept / slope
'PiZi' = fitted intercept Pi/Zi
'slope' = fitted slope (negative, depletion)
'PZ' = list of P/Z values, one per input pressure
"""
Ppc, Tpc = pseudo_critical_properties(gamma_g)
Tpr = (T_F + 459.67) / Tpc
PZ = []
for P in P_psia:
Ppr = P / Ppc
Z = z_factor_DAK(Ppr, Tpr)
PZ.append(P / Z)
slope, intercept = np.polyfit(Gp_bscf, PZ, 1)
G = -intercept / slope
return {
'G_bscf': float(G),
'PiZi': float(intercept),
'slope': float(slope),
'PZ': [float(x) for x in PZ],
}
mb = gas_ogip(GP, P, T_F, GAMMA_G)
ogip_bscf = mb['G_bscf']
remaining_bscf = ogip_bscf - GP[-1]
print("P/Z values:", [round(x, 1) for x in mb['PZ']])
print("Pi/Zi (intercept):", round(mb['PiZi'], 1))
print("slope (Bscf^-1):", round(mb['slope'], 3))
print("OGIP G (Bscf):", round(ogip_bscf, 2))
print("Remaining reserves (Bscf):", round(remaining_bscf, 2))
lockCopying code is a Full Access feature.