Decline Curve Analysis → EUR

Fit Arps decline models to a production history and forecast the well's estimated ultimate recovery.

What it does

Fits exponential, hyperbolic and harmonic Arps models to a rate–time history with SciPy, picks the best fit, forecasts production to an economic limit, and integrates the curve to report EUR and remaining reserves.

How it works

1Fit all three Arps models with SciPy's curve_fit.
2Score each by RMSE and keep the best fit.
3Forecast forward to your economic-limit rate.
4Integrate the rate–time curve for EUR and remaining reserves.
5Plot history vs. forecast — change the economic limit and re-run.

play_circle

Run it live

Runs entirely in your browser — no install. Edit any value and press Run to see the result change. This is the whole project; read it top to bottom.

decline-curve-eur.pystarting Python…

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# ── A monthly oil-rate history, standing in for a production CSV ─────────
# In Ch 9 you load real field data; here we make a noisy declining well.
t = np.arange(0, 36)                                   # months on production
q_true = 900 / (1 + 0.5 * 0.08 * t) ** (1 / 0.5)       # a hyperbolic truth
rng = np.random.default_rng(3)
q = q_true * (1 + rng.normal(0, 0.06, len(t)))         # noisy observations

# ── Arps decline models ─────────────────────────────────────────────────
def exponential(t, qi, Di):   return qi * np.exp(-Di * t)
def hyperbolic(t, qi, Di, b): return qi / (1 + b * Di * t) ** (1 / b)
def harmonic(t, qi, Di):      return qi / (1 + Di * t)

# ── Fit all three, keep the one with the lowest error ───────────────────
fits = {}
fits["Exponential"] = (exponential, curve_fit(exponential, t, q, p0=[900, .08], maxfev=10000)[0])
fits["Hyperbolic"]  = (hyperbolic,  curve_fit(hyperbolic,  t, q, p0=[900, .08, .5],
                       bounds=([0, 0, .01], [5000, 1, 1]), maxfev=10000)[0])
fits["Harmonic"]    = (harmonic,    curve_fit(harmonic,    t, q, p0=[900, .08], maxfev=10000)[0])

rmse = lambda f, p: np.sqrt(np.mean((f(t, *p) - q) ** 2))
best = min(fits, key=lambda k: rmse(*fits[k]))
f_best, p_best = fits[best]

# ── Forecast to an economic limit, integrate for EUR ────────────────────
Q_ECON = 30.0                                          # bopd economic limit (tweak!)
DAYS   = 30.4
t_fore = np.arange(0, 240)
q_fore = f_best(t_fore, *p_best)
t_end  = t_fore[q_fore >= Q_ECON][-1]
eur    = np.trapz(f_best(np.arange(0, t_end + 1), *p_best), dx=DAYS)  # bbl
cum    = np.trapz(q, dx=DAYS)

print(f"Best-fit model : {best}")
extra = f",  b = {p_best[2]:.2f}" if len(p_best) > 2 else ""
print(f"  qi = {p_best[0]:5.0f} bopd,  Di = {p_best[1]:.3f}/mo{extra}")
print(f"Economic limit : {Q_ECON:.0f} bopd reached at month {t_end}")
print(f"Cum to date    : {cum/1e3:6.1f} Mbbl")
print(f"EUR            : {eur/1e3:6.1f} Mbbl   (remaining {max(eur-cum,0)/1e3:.1f} Mbbl)")

# ── Plot history + fitted forecast ──────────────────────────────────────
plt.figure(figsize=(8, 5))
plt.scatter(t, q, s=18, c="0.35", label="History")
plt.plot(t_fore, q_fore, "orange", lw=2, label=f"{best} forecast")
plt.axhline(Q_ECON, ls="--", c="red", label=f"Economic limit ({Q_ECON:.0f} bopd)")
plt.axvline(t_end, ls=":", c="0.6")
plt.xlabel("Time (months)"); plt.ylabel("Oil rate (bopd)")
plt.title(f"Decline Curve Analysis  -  EUR {eur/1e3:.1f} Mbbl ({best})", fontweight="bold")
plt.legend(); plt.grid(alpha=0.3); plt.ylim(0, None); plt.tight_layout(); plt.show()

school

Want the full theory and the real-data version? Read the chapter this project is built from.

Read the chapter Companion repo