Part II: Petroleum Data Engineering

Chapter 9

Production Decline Curve Analysis

schedule15 min readfitness_center10 exercises

Why This Chapter Exists

Every well ever drilled will eventually stop producing. The question is not if production will decline, but how fast, by how much, and how much oil or gas remains to be recovered. These questions drive some of the largest financial decisions in the industry: whether to drill additional wells, when to invest in artificial lift, how to value an asset for acquisition, and when to abandon a field.

Decline Curve Analysis (DCA) is the primary tool for answering these questions. It fits mathematical models to a well's production history and extrapolates into the future. The method is deceptively simple --- you are fitting a curve to data points --- but the engineering judgment required to do it well is substantial. A poorly chosen model or a careless fit can overestimate reserves by millions of barrels, leading to investment decisions that destroy value.

This chapter teaches you to do DCA properly: understanding the physics behind each model, fitting them rigorously with Python, quantifying uncertainty, and knowing when each model applies and when it breaks down.

infoWhat You'll Learn

The physical basis for production decline in oil and gas wells
Arps' three decline models: exponential, hyperbolic, and harmonic
Fitting decline parameters to production data using scipy.optimize
Calculating Estimated Ultimate Recovery (EUR) and remaining reserves
Modified Hyperbolic and Duong models for unconventional wells
Probabilistic DCA: generating P10, P50, and P90 forecasts
Automating DCA across a portfolio of wells

The Physics of Decline

Before touching any equations, you need to understand why production declines in the first place.

When a well first begins producing, the reservoir pressure near the wellbore is high. Hydrocarbons flow toward the well driven by the pressure difference between the reservoir and the wellbore. Over time, as fluid is withdrawn, the pressure around the well drops. The driving force weakens. Flow rates decrease.

How quickly this happens depends on several factors:

Reservoir permeability --- how easily fluid moves through the rock. High-permeability reservoirs deliver fluid quickly but can also deplete quickly.
Drive mechanism --- whether the reservoir has natural energy support (water drive, gas cap expansion) or relies solely on fluid expansion.
Well completion --- the number and placement of perforations, the length of the horizontal section, whether hydraulic fracturing was performed.
Fluid properties --- gas wells decline differently from oil wells because gas is compressible and its flow behavior changes with pressure.

These physical factors determine the shape of the decline curve. Arps' empirical models capture that shape mathematically.

Arps' Decline Models

In 1945, J.J. Arps published a framework that remains the foundation of DCA to this day. He defined decline in terms of the loss ratio, which describes how quickly the production rate changes relative to its current value. From this single concept, three models emerge depending on how the loss ratio behaves over time.

The Loss Ratio

The instantaneous decline rate $D$ is defined as:

D = -\frac{1}{q}\frac{dq}{dt}

This is the fractional rate of change of production. If a well produces 1,000 bbl/d today and loses 10 bbl/d per day, $D = 0.01$ per day, or about 3.65 per year.

What matters is how $D$ changes over time:

If $D$ is constant --- the decline rate never changes --- you get exponential decline.
If $D$ decreases over time at a rate proportional to $D$ itself --- you get hyperbolic decline.
If $D$ decreases such that $1/D$ increases linearly --- you get harmonic decline (a special case of hyperbolic where $b = 1$ ).

Exponential Decline ($b = 0$)

The simplest case. The production rate decays as:

q(t) = q_i \cdot e^{-D_i t}

where $q_i$ is the initial rate and $D_i$ is the constant decline rate.

This model applies when the reservoir is in boundary-dominated flow with a strong drive mechanism --- typically a well producing under constant bottomhole pressure from a closed reservoir with good permeability. In practice, pure exponential decline is uncommon in the early life of a well but is often a reasonable assumption for mature wells in late-stage depletion.

The EUR (Estimated Ultimate Recovery) for exponential decline has a clean closed-form solution:

EUR = \frac{q_i}{D_i}

This is one of the most useful results in petroleum engineering. If you know the initial rate and the decline rate, you know the total recoverable volume.

Hyperbolic Decline ($0 < b < 1$)

Most real wells follow hyperbolic decline, where the decline rate itself slows over time:

q(t) = \frac{q_i}{(1 + b D_i t)^{1/b}}

The parameter $b$ controls the curvature. When $b = 0$ , the equation reduces to exponential. When $b = 1$ , it becomes harmonic. Values between 0 and 1 produce curves that decline steeply at first, then flatten.

The physical interpretation: early in a well's life, the near-wellbore region depletes rapidly, causing steep decline. As the drainage area expands, the decline slows because the well is drawing from a larger volume of rock.

The EUR for hyperbolic decline is:

EUR = \frac{q_i}{D_i(1-b)} \left[ 1 - \left(\frac{q_{ab}}{q_i}\right)^{1-b} \right]

where $q_{ab}$ is the abandonment rate --- the rate below which the well is no longer economic to operate.

warningThe $b > 1$ Problem

When $b > 1$ , the EUR integral diverges --- it predicts infinite recovery. This is physically impossible. If your curve fit returns $b > 1$ , something is wrong: the well may still be in transient flow, the data may contain artifacts, or the hyperbolic model may not be appropriate. This is one of the most common mistakes in DCA.

Harmonic Decline ($b = 1$)

A special case where:

q(t) = \frac{q_i}{1 + D_i t}

Harmonic decline produces the slowest rate of decline among Arps' models. It is occasionally observed in heavy oil reservoirs with strong water drive, where the displacement process maintains relatively high reservoir pressure.

Fitting Decline Curves in Python

With the theory established, let's build the tools to apply it.

main.py

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

# --- Define the Arps decline models ---

def exponential_decline(t, qi, Di):
    """Exponential decline: constant decline rate."""
    return qi * np.exp(-Di * t)

def hyperbolic_decline(t, qi, Di, b):
    """Hyperbolic decline: decline rate decreases over time."""
    return qi / (1 + b * Di * t) ** (1 / b)

def harmonic_decline(t, qi, Di):
    """Harmonic decline: special case of hyperbolic with b=1."""
    return qi / (1 + Di * t)

These three functions are the core of every DCA analysis. Each takes time as input and returns the predicted production rate. The parameters ( $q_i$ , $D_i$ , and optionally $b$ ) are what we need to estimate from data.

Loading and Preparing Production Data

Production data in the real world arrives as monthly or daily records, often with gaps, unit inconsistencies, and noise. Before fitting any model, the data must be cleaned.

main.py

import pandas as pd

# Simulated monthly production data for a Niger Delta well
# In practice, this comes from a production database or CSV export
np.random.seed(42)

months = np.arange(1, 61)  # 5 years of monthly data
qi_true = 2400  # True initial rate, bbl/d
Di_true = 0.08  # True initial decline rate, 1/month
b_true = 0.65   # True Arps b-factor

# Generate synthetic data with realistic noise
q_true = hyperbolic_decline(months, qi_true, Di_true, b_true)
noise = np.random.normal(0, 50, len(months))  # ±50 bbl/d measurement noise
q_observed = q_true + noise

# Simulate two months of shut-in (zero production) --- common in real data
q_observed[18] = 0   # Month 19: planned shutdown
q_observed[19] = 0   # Month 20: planned shutdown

# Build a production DataFrame
production = pd.DataFrame({
    "Month": months,
    "Oil_Rate_bpd": q_observed.round(1)
})

# Data cleaning: remove shut-in periods (zero or negative rates)
# These are operational events, not decline behavior
production_clean = production[production["Oil_Rate_bpd"] > 0].copy()
production_clean = production_clean.reset_index(drop=True)

print(f"Raw records:     {len(production)}")
print(f"After cleaning:  {len(production_clean)}")
print(f"Removed:         {len(production) - len(production_clean)} months (shut-in periods)")
print()
print(production_clean.head(10).to_string(index=False))

The cleaning step is critical. Shut-in periods, equipment failures, and well interventions create zeros or anomalous values that distort curve fits. Removing them is not "cheating" --- it is separating operational events from reservoir behavior.

Fitting the Models

Now we fit each of the three Arps models to the cleaned data and compare their quality of fit.

main.py

t_data = production_clean["Month"].values
q_data = production_clean["Oil_Rate_bpd"].values

# --- Fit Exponential ---
try:
    popt_exp, pcov_exp = curve_fit(
        exponential_decline, t_data, q_data,
        p0=[2000, 0.05],           # Initial guesses
        bounds=([0, 0], [5000, 1]) # Physical constraints
    )
    qi_exp, Di_exp = popt_exp
    q_pred_exp = exponential_decline(t_data, *popt_exp)
    residuals_exp = q_data - q_pred_exp
    rmse_exp = np.sqrt(np.mean(residuals_exp**2))
    print(f"Exponential:  qi = {qi_exp:,.0f} bbl/d,  Di = {Di_exp:.4f} /month,  RMSE = {rmse_exp:.1f} bbl/d")
except RuntimeError:
    print("Exponential fit failed to converge.")

# --- Fit Hyperbolic ---
try:
    popt_hyp, pcov_hyp = curve_fit(
        hyperbolic_decline, t_data, q_data,
        p0=[2000, 0.05, 0.5],
        bounds=([0, 0, 0], [5000, 1, 1])  # b constrained to [0, 1]
    )
    qi_hyp, Di_hyp, b_hyp = popt_hyp
    q_pred_hyp = hyperbolic_decline(t_data, *popt_hyp)
    residuals_hyp = q_data - q_pred_hyp
    rmse_hyp = np.sqrt(np.mean(residuals_hyp**2))
    print(f"Hyperbolic:   qi = {qi_hyp:,.0f} bbl/d,  Di = {Di_hyp:.4f} /month,  b = {b_hyp:.3f},  RMSE = {rmse_hyp:.1f} bbl/d")
except RuntimeError:
    print("Hyperbolic fit failed to converge.")

# --- Fit Harmonic ---
try:
    popt_har, pcov_har = curve_fit(
        harmonic_decline, t_data, q_data,
        p0=[2000, 0.05],
        bounds=([0, 0], [5000, 1])
    )
    qi_har, Di_har = popt_har
    q_pred_har = harmonic_decline(t_data, *popt_har)
    residuals_har = q_data - q_pred_har
    rmse_har = np.sqrt(np.mean(residuals_har**2))
    print(f"Harmonic:     qi = {qi_har:,.0f} bbl/d,  Di = {Di_har:.4f} /month,  RMSE = {rmse_har:.1f} bbl/d")
except RuntimeError:
    print("Harmonic fit failed to converge.")

The RMSE (Root Mean Square Error) tells you how well each model reproduces the observed data. Lower is better. But RMSE alone is not sufficient for model selection --- a model can fit historical data well and still produce unreasonable forecasts. Engineering judgment matters.

Visualizing the Fits

main.py

t_forecast = np.arange(1, 121)  # Forecast 10 years

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# --- Cartesian plot ---
ax1 = axes[0]
ax1.scatter(t_data, q_data, s=15, c="black", alpha=0.6, label="Observed", zorder=5)
ax1.plot(t_forecast, exponential_decline(t_forecast, *popt_exp),
         "--", color="#2196F3", label=f"Exponential (RMSE={rmse_exp:.0f})")
ax1.plot(t_forecast, hyperbolic_decline(t_forecast, *popt_hyp),
         "-", color="#E53935", linewidth=2, label=f"Hyperbolic (RMSE={rmse_hyp:.0f})")
ax1.plot(t_forecast, harmonic_decline(t_forecast, *popt_har),
         "-.", color="#43A047", label=f"Harmonic (RMSE={rmse_har:.0f})")
ax1.axvline(x=t_data[-1], color="gray", linestyle=":", alpha=0.5, label="End of history")
ax1.set_xlabel("Time (months)")
ax1.set_ylabel("Oil Rate (bbl/d)")
ax1.set_title("Cartesian Scale")
ax1.legend(fontsize=8)
ax1.set_ylim(bottom=0)

# --- Semi-log plot ---
ax2 = axes[1]
ax2.scatter(t_data, q_data, s=15, c="black", alpha=0.6, zorder=5)
ax2.plot(t_forecast, exponential_decline(t_forecast, *popt_exp), "--", color="#2196F3")
ax2.plot(t_forecast, hyperbolic_decline(t_forecast, *popt_hyp), "-", color="#E53935", linewidth=2)
ax2.plot(t_forecast, harmonic_decline(t_forecast, *popt_har), "-.", color="#43A047")
ax2.axvline(x=t_data[-1], color="gray", linestyle=":", alpha=0.5)
ax2.set_xlabel("Time (months)")
ax2.set_ylabel("Oil Rate (bbl/d)")
ax2.set_title("Semi-Log Scale")
ax2.set_yscale("log")

plt.tight_layout()
plt.show()

The semi-log plot is standard practice in DCA. On a semi-log scale, exponential decline appears as a straight line. Deviation from a straight line indicates hyperbolic or harmonic behavior. This visual check is one of the fastest ways to identify the appropriate model.

Calculating EUR and Remaining Reserves

The primary deliverable of any DCA analysis is a reserves estimate. EUR (Estimated Ultimate Recovery) represents the total volume of oil or gas a well will produce over its entire economic life.

main.py

def eur_exponential(qi, Di, q_abandon):
    """EUR for exponential decline, in same units as qi × time."""
    if qi <= q_abandon:
        return 0
    return (qi - q_abandon) / Di

def eur_hyperbolic(qi, Di, b, q_abandon):
    """EUR for hyperbolic decline."""
    if qi <= q_abandon or b >= 1:
        return np.nan
    return (qi**b / (Di * (1 - b))) * (qi**(1-b) - q_abandon**(1-b))

def eur_harmonic(qi, Di, q_abandon):
    """EUR for harmonic decline (b=1)."""
    if qi <= q_abandon:
        return 0
    return (qi / Di) * np.log(qi / q_abandon)

# Economic limit: the rate below which operating costs exceed revenue
# This varies by well, but 20 bbl/d is a reasonable assumption for
# a conventional onshore well in Nigeria
q_abandon = 20  # bbl/d

# Calculate cumulative production already recovered
months_produced = len(t_data)
cum_prod_actual = np.trapezoid(q_data, t_data)  # bbl-months
cum_prod_bbl = cum_prod_actual * 30.44       # Convert to bbl (avg days/month)

# EUR for each model (in bbl-months, converted to barrels)
eur_exp = eur_exponential(qi_exp, Di_exp, q_abandon) * 30.44
eur_hyp = eur_hyperbolic(qi_hyp, Di_hyp, b_hyp, q_abandon) * 30.44
eur_har = eur_harmonic(qi_har, Di_har, q_abandon) * 30.44

# Remaining reserves = EUR - cumulative production
remaining_exp = eur_exp - cum_prod_bbl
remaining_hyp = eur_hyp - cum_prod_bbl
remaining_har = eur_har - cum_prod_bbl

print("=" * 65)
print("              RESERVES ESTIMATION SUMMARY")
print("=" * 65)
print(f"  Abandonment rate:      {q_abandon} bbl/d")
print(f"  Cumulative production: {cum_prod_bbl:,.0f} bbl")
print()
print(f"  {'Model':<14} {'EUR (bbl)':>14} {'Remaining (bbl)':>16}")
print(f"  {'-'*14} {'-'*14} {'-'*16}")
print(f"  {'Exponential':<14} {eur_exp:>14,.0f} {remaining_exp:>16,.0f}")
print(f"  {'Hyperbolic':<14} {eur_hyp:>14,.0f} {remaining_hyp:>16,.0f}")
print(f"  {'Harmonic':<14} {eur_har:>14,.0f} {remaining_har:>16,.0f}")
print("=" * 65)
print()
print(f"  Range of EUR estimates: {min(eur_exp, eur_hyp, eur_har):,.0f}"
      f" to {max(eur_exp, eur_hyp, eur_har):,.0f} bbl")
print(f"  Spread: {max(eur_exp, eur_hyp, eur_har) - min(eur_exp, eur_hyp, eur_har):,.0f} bbl")

The spread between model estimates is important. A narrow spread means all models agree, and you can report EUR with reasonable confidence. A wide spread means the data does not clearly distinguish between decline behaviors, and you need either more production history or additional engineering judgment to select the appropriate model.

warningEUR Is Not Reserves

EUR is a technical estimate of total recoverable volume under current operating conditions. Proved reserves (1P), probable reserves (2P), and possible reserves (3P) are formal classifications defined by the SPE Petroleum Resources Management System (PRMS) and involve additional considerations including economics, regulatory approval, and development plans. DCA informs reserves estimation but does not replace it.

Modified Hyperbolic Decline for Unconventional Wells

Unconventional wells --- horizontal wells with multi-stage hydraulic fractures in shale or tight rock --- present a challenge for standard Arps models. These wells often exhibit $b$ values greater than 1 during early transient flow, which causes the hyperbolic model to predict infinite EUR.

The Modified Hyperbolic (also called the Robertson method) addresses this by switching from hyperbolic decline to exponential decline at a specified minimum decline rate $D_{min}$ :

q(t) = \begin{cases} \frac{q_i}{(1 + b D_i t)^{1/b}} & \text{if } D(t) > D_{min} \\[6pt] q_{switch} \cdot e^{-D_{min}(t - t_{switch})} & \text{if } D(t) \leq D_{min} \end{cases}

The switch point occurs when the instantaneous decline rate of the hyperbolic model drops to $D_{min}$ . After that, the well declines exponentially at a constant rate.

main.py

def modified_hyperbolic(t, qi, Di, b, Dmin):
    """Modified Hyperbolic: switches to exponential at Dmin."""
    q = np.zeros_like(t, dtype=float)

    for i, ti in enumerate(t):
        # Instantaneous decline rate at time ti (hyperbolic)
        D_t = Di / (1 + b * Di * ti)

        if D_t > Dmin:
            # Still in hyperbolic regime
            q[i] = qi / (1 + b * Di * ti) ** (1 / b)
        else:
            # Switch to exponential
            # Find the switch time
            t_switch = (Di / Dmin - 1) / (b * Di)
            q_switch = qi / (1 + b * Di * t_switch) ** (1 / b)
            q[i] = q_switch * np.exp(-Dmin * (ti - t_switch))

    return q

# Example: fit to the same data but allow b > 1
# and impose a minimum terminal decline of 0.5%/month
popt_mod, _ = curve_fit(
    lambda t, qi, Di, b: modified_hyperbolic(t, qi, Di, b, Dmin=0.005),
    t_data, q_data,
    p0=[2000, 0.1, 0.8],
    bounds=([0, 0, 0], [5000, 2, 2])  # b can exceed 1 safely now
)

qi_mod, Di_mod, b_mod = popt_mod
print(f"Modified Hyperbolic: qi = {qi_mod:,.0f} bbl/d, Di = {Di_mod:.4f}, b = {b_mod:.3f}")
print(f"Terminal decline rate (Dmin): 0.50%/month")

Probabilistic DCA: P10, P50, P90

A single best-fit curve gives a single EUR estimate. But all curve fits carry uncertainty --- the data has noise, the model is an approximation, and the future is unknown. Reporting a single number without quantifying uncertainty is misleading.

Probabilistic DCA generates a distribution of possible outcomes. The industry standard is to report three percentiles:

P90 --- the conservative estimate. There is a 90% probability that actual recovery will exceed this value.
P50 --- the most likely estimate. Equal probability of exceeding or falling short.
P10 --- the optimistic estimate. Only a 10% probability of exceeding this value.

main.py

def monte_carlo_dca(t_data, q_data, n_iterations=1000):
    """
    Generate probabilistic decline forecasts using
    bootstrap resampling of residuals.
    """
    # Fit the base model
    popt, pcov = curve_fit(
        hyperbolic_decline, t_data, q_data,
        p0=[2000, 0.05, 0.5],
        bounds=([0, 0, 0.01], [5000, 1, 0.99])
    )

    # Calculate residuals
    residuals = q_data - hyperbolic_decline(t_data, *popt)

    # Generate forecast timeline
    t_forecast = np.arange(1, 241)  # 20 years
    forecasts = np.zeros((n_iterations, len(t_forecast)))

    for i in range(n_iterations):
        # Resample residuals with replacement
        resampled_residuals = np.random.choice(residuals, size=len(t_data), replace=True)
        q_resampled = hyperbolic_decline(t_data, *popt) + resampled_residuals
        q_resampled = np.maximum(q_resampled, 1)  # No negative rates

        # Refit to perturbed data
        try:
            popt_i, _ = curve_fit(
                hyperbolic_decline, t_data, q_resampled,
                p0=popt,
                bounds=([0, 0, 0.01], [5000, 1, 0.99]),
                maxfev=5000
            )
            forecasts[i] = hyperbolic_decline(t_forecast, *popt_i)
        except RuntimeError:
            forecasts[i] = hyperbolic_decline(t_forecast, *popt)

    # Calculate percentiles
    p10 = np.percentile(forecasts, 90, axis=0)  # P10 = 90th percentile of rate
    p50 = np.percentile(forecasts, 50, axis=0)
    p90 = np.percentile(forecasts, 10, axis=0)  # P90 = 10th percentile of rate

    return t_forecast, p10, p50, p90

t_fc, p10, p50, p90 = monte_carlo_dca(t_data, q_data, n_iterations=500)

# Calculate EUR for each case
q_ab = 20  # abandonment rate
eur_p10 = np.trapezoid(np.maximum(p10, q_ab), t_fc) * 30.44
eur_p50 = np.trapezoid(np.maximum(p50, q_ab), t_fc) * 30.44
eur_p90 = np.trapezoid(np.maximum(p90, q_ab), t_fc) * 30.44

print("PROBABILISTIC RESERVES ESTIMATE")
print(f"  P90 (conservative):  {eur_p90:>12,.0f} bbl")
print(f"  P50 (most likely):   {eur_p50:>12,.0f} bbl")
print(f"  P10 (optimistic):    {eur_p10:>12,.0f} bbl")

main.py

fig, ax = plt.subplots(figsize=(10, 5))

ax.scatter(t_data, q_data, s=15, c="black", alpha=0.6, label="Observed", zorder=5)
ax.plot(t_fc, p50, "-", color="#E53935", linewidth=2, label="P50 (most likely)")
ax.plot(t_fc, p10, "--", color="#43A047", linewidth=1, label="P10 (optimistic)")
ax.plot(t_fc, p90, "--", color="#2196F3", linewidth=1, label="P90 (conservative)")
ax.fill_between(t_fc, p90, p10, alpha=0.08, color="gray", label="P10–P90 range")
ax.axvline(x=t_data[-1], color="gray", linestyle=":", alpha=0.5)
ax.axhline(y=q_ab, color="red", linestyle=":", alpha=0.3, label=f"Abandonment ({q_ab} bbl/d)")
ax.set_xlabel("Time (months)")
ax.set_ylabel("Oil Rate (bbl/d)")
ax.set_title("Probabilistic Decline Forecast")
ax.legend(fontsize=8)
ax.set_ylim(bottom=0)
plt.tight_layout()
plt.show()

The shaded region between P10 and P90 represents the range of plausible outcomes. A narrow band means the forecast is confident. A wide band means more data or additional engineering analysis is needed before committing capital.

Multi-Well DCA Automation

In practice, you rarely analyze a single well. A field evaluation requires DCA for every producing well --- often dozens or hundreds. Automating this process is where Python's value becomes undeniable.

main.py

def auto_dca(well_name, t, q, q_abandon=20):
    """
    Perform automated DCA for a single well.
    Returns fitted parameters and EUR.
    """
    # Remove zero/negative rates
    mask = q > 0
    t_clean = t[mask]
    q_clean = q[mask]

    if len(t_clean) < 6:
        return {"well": well_name, "status": "Insufficient data"}

    try:
        popt, _ = curve_fit(
            hyperbolic_decline, t_clean, q_clean,
            p0=[q_clean[0], 0.05, 0.5],
            bounds=([0, 0, 0.01], [q_clean[0]*2, 1, 0.99]),
            maxfev=10000
        )
        qi, Di, b = popt
        q_pred = hyperbolic_decline(t_clean, *popt)
        rmse = np.sqrt(np.mean((q_clean - q_pred)**2))

        eur = eur_hyperbolic(qi, Di, b, q_abandon) * 30.44
        cum = np.trapezoid(q_clean, t_clean) * 30.44
        remaining = eur - cum

        return {
            "well": well_name,
            "qi_bpd": round(qi),
            "Di_per_month": round(Di, 4),
            "b_factor": round(b, 3),
            "rmse_bpd": round(rmse, 1),
            "eur_bbl": round(eur),
            "cum_bbl": round(cum),
            "remaining_bbl": round(remaining),
            "status": "OK"
        }
    except RuntimeError:
        return {"well": well_name, "status": "Fit failed"}

# Generate synthetic data for a portfolio of 8 wells
np.random.seed(123)
portfolio = {}
for i in range(8):
    qi = np.random.uniform(800, 3000)
    Di = np.random.uniform(0.03, 0.12)
    b = np.random.uniform(0.3, 0.9)
    n_months = np.random.randint(24, 72)
    t = np.arange(1, n_months + 1)
    q = hyperbolic_decline(t, qi, Di, b) + np.random.normal(0, 40, n_months)
    q = np.maximum(q, 0)
    portfolio[f"OD-{i+1:03d}"] = (t, q)

# Run DCA on all wells
results = []
for name, (t, q) in portfolio.items():
    results.append(auto_dca(name, t, q))

results_df = pd.DataFrame(results)
print(results_df.to_string(index=False))
print()
print(f"Total portfolio EUR: {results_df['eur_bbl'].sum():,.0f} bbl")
print(f"Total remaining:     {results_df['remaining_bbl'].sum():,.0f} bbl")

This is the kind of workflow that takes days in a spreadsheet and minutes in Python. Each well gets a rigorous fit, quality metrics, and reserves estimates. The results are in a single table ready for reporting or further analysis.

Summary

Production decline is caused by reservoir pressure depletion as fluids are withdrawn. The rate and shape of decline depend on reservoir properties, drive mechanism, and well completion.
Arps' framework defines three models based on how the decline rate changes over time: exponential ( $b=0$ ), hyperbolic ( $0 < b < 1$ ), and harmonic ( $b=1$ ).
Exponential decline assumes a constant decline rate. It is the simplest model and applies best to mature wells in boundary-dominated flow. EUR = $q_i / D_i$ .
Hyperbolic decline is the most commonly used model. The decline rate itself decreases over time, producing a curve that flattens. The $b$ -factor controls the curvature.
EUR is calculated by integrating the decline curve to an economic abandonment rate. Different models can produce significantly different EUR estimates from the same data.
Modified Hyperbolic addresses the problem of $b > 1$ in unconventional wells by switching to exponential decline at a minimum decline rate.
Probabilistic DCA quantifies uncertainty by generating P10, P50, and P90 forecasts using bootstrap resampling. Single-point estimates without uncertainty ranges are incomplete.
Multi-well automation is where Python transforms DCA from a well-by-well exercise into a portfolio-scale analysis.

Exercises

fitness_center

Exercise 9.1Practice

-- Exponential vs. Hyperbolic

A well has the following monthly production data (bbl/d): 1800, 1650, 1520, 1400, 1295, 1200, 1110, 1030, 958, 892, 832, 778. Fit both exponential and...

arrow_forward