Decline Curve Analysis in Python: Fit the Arps Equation and Forecast EUR

By Saad Iqbal

It’s Friday afternoon and you’ve got eighteen wells to report on before the reserves committee meets Monday. You open the spreadsheet, filter to the first API number, paste the monthly rates into a scatter chart, and start dragging a trendline by eye. Does this one look more exponential or hyperbolic? You nudge the b-factor until the curve looks right, write down an EUR that feels plausible, and move to well number two. By well twelve your eyes have stopped trusting themselves, and you already know next month you’ll do this exact same thing again, by hand, from scratch.

Decline curve analysis python workflows exist precisely to kill that Friday. The Arps decline equation hasn’t changed since 1945, but the way engineers fit it has: instead of dragging a trendline, you fit it numerically with SciPy, get a repeatable EUR with a confidence interval attached, and re-run the whole portfolio in the time it takes to get coffee. This tutorial builds that workflow from a raw production export to an automated weekly refresh, using real equations, a real worked fit, and numbers you can check yourself.

What the Arps Equation Actually Says

J.J. Arps formalized three decline shapes in 1945, and upstream engineers still default to the hyperbolic case because it’s the one that fits shale and tight-rock wells best:

q(t) = q_i / (1 + b·D_i·t)^1/b

• q_i — initial production rate
• D_i — initial nominal decline rate (per year)
• b — the hyperbolic exponent; b=0 collapses to exponential decline, b=1 gives harmonic decline, and most unconventional wells land somewhere between 0.6 and 1.3 in their early flush-production months

Integrating that rate equation gives cumulative production, which is the number the reserves committee actually wants:

N_p(t) = [q_i^b / ((1−b)·D_i)] · [q_i^1−b − q(t)^1−b]

You almost never need to evaluate that by hand — that’s what Step 3 below is for — but every automated forecast is doing exactly this arithmetic underneath, so it’s worth being able to read it before you trust a tool’s output.

Prerequisites

• Python 3.9+ installed from python.org, with pandas, scipy, and matplotlib (pip install pandas scipy matplotlib)
• A Jupyter notebook or any Python editor to run the code blocks interactively
• A production history export with at least 12–24 months of monthly or daily rate data (date and rate columns, at minimum)
• Basic comfort with logarithms and the idea of a best-fit curve — no calculus required to use the code, though the equation above uses it

Step 1: Import Production History into Python

Start by getting the production export into a clean, time-indexed pandas DataFrame. Most operator reporting systems export a CSV with a date column and a rate column; the only real job here is making sure the date parses correctly and the series is sorted.

import pandas as pd
import numpy as np

df = pd.read_csv("well_production.csv", parse_dates=["date"])
df = df.sort_values("date").reset_index(drop=True)

# Convert calendar dates into decimal years since first production,
# which is the time basis the Arps equation expects
t0 = df["date"].iloc[0]
df["t_years"] = (df["date"] - t0).dt.days / 365.25

print(df.head())
print(f"History length: {df['t_years'].iloc[-1]:.2f} years, {len(df)} data points")

Two things to check before moving on: make sure shut-in days or workover gaps haven’t been silently zero-filled (that will bias the fit low), and make sure the rate column is a true daily or monthly average, not a cumulative total disguised as a rate.

Step 2: Fit the Arps Decline Curve with SciPy

This is the step that replaces the eyeballed trendline. scipy.optimize.curve_fit finds the q_i, D_i, and b that minimize squared error against your observed rates, using nonlinear least squares. Here’s a worked example built on 24 months of synthetic-but-realistic monthly rates (true parameters: q_i=850 bbl/d, D_i=0.70/yr, b=0.85, with 2.5% random noise added to mimic normal field reporting scatter):

from scipy.optimize import curve_fit

def arps_hyperbolic(t, qi, Di, b):
    return qi / (1 + b * Di * t) ** (1.0 / b)

t = df["t_years"].values
q = df["rate_bopd"].values

# Reasonable starting guesses and physical bounds keep the solver honest
p0 = [q[0], 0.6, 0.8]
bounds = ([1, 0.05, 0.05], [5000, 3.0, 1.95])

popt, pcov = curve_fit(arps_hyperbolic, t, q, p0=p0, bounds=bounds, maxfev=20000)
qi_fit, Di_fit, b_fit = popt
perr = np.sqrt(np.diag(pcov))

print(f"qi = {qi_fit:.1f} bbl/d (+/- {perr[0]:.1f})")
print(f"Di = {Di_fit:.3f} /yr (+/- {perr[1]:.3f})")
print(f"b  = {b_fit:.3f} (+/- {perr[2]:.3f})")

Run against that 24-month synthetic dataset, this actually produces q_i = 882.8 bbl/d, D_i = 0.878/yr, and b = 1.203 — close to, but not identical to, the true generating parameters. That gap is the lesson, not a bug: D_i and b trade off against each other in the fit, so a noisy or short history can land on a different (D_i, b) pair that still reproduces the observed rates almost as well. The residual root-mean-square error here is 2.3% of the mean rate, which is a reasonable fit — but a good residual doesn’t guarantee the “true” parameters, only a curve that tracks history well.

Step 3: Generate the Forecast and Cross-Check It

With fitted parameters in hand, forecasting forward and integrating to EUR is direct arithmetic:

def cumulative_production(t, qi, Di, b):
    """Np(t) in the same volume units as qi, t in years."""
    qt = arps_hyperbolic(t, qi, Di, b)
    np_bbl = (qi**b / ((1 - b) * Di)) * (qi**(1 - b) - qt**(1 - b))
    return np_bbl * 365.25  # convert daily-rate basis to annualized volume

from scipy.optimize import brentq

econ_limit_bopd = 10.0  # stripper economic limit for this example
t_limit = brentq(lambda tt: arps_hyperbolic(tt, qi_fit, Di_fit, b_fit) - econ_limit_bopd,
                  0.1, 200)
eur_bbl = cumulative_production(t_limit, qi_fit, Di_fit, b_fit)

print(f"Time to {econ_limit_bopd} bbl/d economic limit: {t_limit:.1f} years")
print(f"EUR at economic limit: {eur_bbl:,.0f} bbl")

On the original (true-parameter) model this well reaches a 10 bbl/d economic limit at 71.7 years and books an EUR of roughly 1.44 million barrels — the kind of long tail that makes b-factor sensitivity worth stress-testing before it goes into a reserves report. Don’t treat this forecast as gospel on its own: cross-check the implied flowing bottomhole pressure against a current nodal analysis operating point. If the decline curve says the well should still be flowing 80 bbl/d at a tubing pressure the nodal model says can’t lift that rate, the curve fit is extrapolating past where the well’s actual deliverability lives, and you need artificial lift or workover assumptions in the forecast, not just more curve.

Step 4: Automate the Weekly Refresh

A single well’s fit is a nice exercise; a portfolio of two hundred wells, refit automatically every week against the newest production data, is what actually gets this off your Friday afternoon. Wrap steps 1–3 in a function keyed by well ID, loop it over every API number in your production database, and write the results — qi, Di, b, EUR, and a fit-quality flag — to a summary table:

def fit_well(well_df, well_id):
    try:
        t = well_df["t_years"].values
        q = well_df["rate_bopd"].values
        p0 = [q[0], 0.6, 0.8]
        bounds = ([1, 0.05, 0.05], [q[0]*3, 3.0, 1.95])
        popt, pcov = curve_fit(arps_hyperbolic, t, q, p0=p0, bounds=bounds, maxfev=20000)
        resid = q - arps_hyperbolic(t, *popt)
        rmse_pct = np.sqrt(np.mean(resid**2)) / np.mean(q) * 100
        return {"well_id": well_id, "qi": popt[0], "Di": popt[1], "b": popt[2],
                "rmse_pct": rmse_pct, "status": "ok"}
    except Exception as e:
        return {"well_id": well_id, "status": f"failed: {e}"}

results = [fit_well(group, wid) for wid, group in df.groupby("well_id")]
summary = pd.DataFrame(results)
summary.to_csv("weekly_decline_summary.csv", index=False)

Save that script, then point a scheduler at it. If your production database already lives in Power BI or Looker Studio, a scheduled Python step in Power Automate or an open-source flow in n8n can trigger the refresh every Monday morning and drop the output straight into your dashboard. If you’d rather not maintain the correlation math yourself, the petropt package ships a ready-made arps_decline function built on the same exponential/hyperbolic/harmonic forms covered here, which is worth a look once you’ve built the fit by hand and trust what it’s doing.

Expected Result and Common Pitfalls

Run this end to end on a well with a clean 18–24 month history and you should see curve_fit converge in well under a second, a residual RMSE in the low single-digit percent of mean rate, and an EUR that lands within a reasonable range of whatever your reservoir engineers booked by hand last quarter. If it doesn’t, check these three things first:

• Short histories inflate b. With fewer than about 12 months of data, b and D_i become statistically hard to tell apart, and the solver will happily push b toward its upper bound to minimize error on too little information. Treat any fit from under a year of data as provisional.
• Un-flagged shut-ins corrupt the curve. A week of zero production from a compressor outage looks, to a least-squares solver, like a sudden catastrophic decline. Filter or flag non-producing days before fitting, don’t just feed in the raw zeros.
• An economic-limit forecast that never lands sanely is a red flag, not a feature. If brentq returns a time-to-limit of hundreds of years, your fitted b is likely pinned near its upper bound and the “long tail” is a fitting artifact, not a reservoir behavior — cap b at a value your basin’s offset wells actually support (commonly 1.0–1.3 for unconventional reservoirs) and refit.

Once this is running on a schedule, the Friday-afternoon trendline-dragging exercise becomes a five-minute sanity check on numbers the computer already produced — which is exactly the point. If you want the equivalent treatment for the deliverability side of the well, the nodal analysis walkthrough builds the IPR/VLP intersection in Python the same way, and if you’re weighing whether to buy rather than build this capability, the AI reservoir simulation and forecasting tools roundup covers where platforms like ComboCurve fit relative to a DIY script like this one.