Exerciseschevron_rightChapter 9chevron_right9.10
fitness_center

Exercise 9.10

Probabilistic P10/P50/P90 EUR

Level 3
Chapter 9: Decline Curve Analysis
descriptionProblem

Download production data for a real well from a public database (the North Dakota Industrial Commission publishes monthly production for every well in the state at dmr.nd.gov). Clean the data, fit Arps models, generate a probabilistic forecast, and report P10/P50/P90 EUR. Document your assumptions and present your results as a one-page technical summary suitable for a reserves review meeting.

---

A single decline fit gives you one EUR number. A reserves review wants the uncertainty: how good, how likely, how bad. You handle that by sweeping the uncertain Arps inputs (qi_bopd, di_per_month, b_factor) across plausible ranges, computing an EUR for each combination, and reading percentiles off the resulting distribution of EURs.

You are booking reserves for well OD-014 on OML 58. An ensemble of 200 deterministic EUR realizations (in barrels) has already been computed for you on a parameter grid and stored in EUR_SAMPLES; treat it as your raw sample of equally-likely outcomes.

## The exceedance convention (read this carefully)

Reserves people quote EUR with the exceedance convention, where the number is the probability the true EUR is at least that big:

LabelMeaningHow to read it off the sample
P9090% chance the EUR exceeds this: the conservative / proved estimatethe 10th percentile of the samples
P5050% chance: the medianthe 50th percentile
P10only 10% chance: the optimistic estimatethe 90th percentile

Because a higher exceedance probability means a lower EUR threshold, the ordering is always P90 < P50 < P10. The percentile index is 100 - P: P90 maps to np.percentile(samples, 10), P10 maps to np.percentile(samples, 90).

## What to build

Write p_values(eur_samples) → dict that returns the three reserves numbers using the exceedance convention above:

keyvalue
p90np.percentile(eur_samples, 10)
p50np.percentile(eur_samples, 50)
p10np.percentile(eur_samples, 90)

Each value should be a plain float. Then call it on EUR_SAMPLES and store the result in result.

Sanity check your output: result["p90"] < result["p50"] < result["p10"], and result["p50"] should sit right at the median of EUR_SAMPLES.

lightbulbHints (0/3)

Stuck? Reveal hints one at a time — they progress from nudge to near-solution.

codeYour solution
main.py
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


# ── Embedded Arps EUR function (do not edit) ─────────────────────────────
def eur_hyperbolic(qi, Di, b, q_abandon):
    if qi <= q_abandon or b >= 1:
        return float("nan")
    return (qi ** b / (Di * (1.0 - b))) * (qi ** (1.0 - b) - q_abandon ** (1.0 - b))


# ── OD-014 (OML 58): 200 deterministic EUR realizations ──────────────────
# We sweep the three uncertain Arps inputs across plausible ranges on a grid
# (5 qi x 5 Di x 8 b = 200 combos) and compute a hyperbolic EUR for each.
# Each realization is one equally-likely outcome for the well's reserves.
Q_ABANDON_BOPD = 20.0


def _build_eur_samples():
    qi_grid = np.linspace(1500.0, 2100.0, 5)      # initial rate (bopd)
    di_grid = np.linspace(0.09, 0.15, 5)          # initial decline (per month)
    b_grid = np.linspace(0.30, 0.90, 8)           # Arps b-factor
    samples = []
    for qi_bopd in qi_grid:
        for di_per_month in di_grid:
            for b_factor in b_grid:
                eur_bbl = eur_hyperbolic(qi_bopd, di_per_month, b_factor,
                                         Q_ABANDON_BOPD)
                samples.append(float(eur_bbl))
    return np.array(samples, dtype=float)


EUR_SAMPLES = _build_eur_samples()


def p_values(eur_samples):
    """Return {p90, p50, p10} EUR using the exceedance convention.

    p90 = 10th percentile (conservative), p50 = median,
    p10 = 90th percentile (optimistic). Each a plain float.
    """
    samples = np.asarray(eur_samples, dtype=float)
    p90 = float(np.percentile(samples, 10))
    p50 = float(np.percentile(samples, 50))
    p10 = float(np.percentile(samples, 90))
    return {"p90": p90, "p50": p50, "p10": p10}


result = p_values(EUR_SAMPLES)

print("n samples:", len(EUR_SAMPLES))
print("result:", result)

lockCopying code is a Full Access feature.