Archie Parameter Sensitivity

Level 2

Chapter 7: Well Log Analysis

descriptionProblem

Water saturation from the Archie equation rides on four parameters: a, m, n, and Rw. None of them are measured on the logging truck; they come from core, from analogues, or from a regional default. So before you book a barrel you should know how much your answer moves when those numbers move.

This is the OD-001 log over OML 58. The OIL PAY sand sits at 7400–7600 ft. You'll build effective porosity once, then run Archie three ways and compare the average Sw in the pay zone.

Build the petrophysics chain on the full log:

Vshale (linear): IGR = (GR - gr_clean)/(gr_shale - gr_clean), clipped

0..1, with gr_clean = GR.quantile(0.05) and gr_shale = GR.quantile(0.95).

Density porosity: PHID = (2.65 - RHOB)/(2.65 - 1.0), clipped 0..0.45.
Average / effective: PHIA = (PHID + NPHI)/2,

PHIE = PHIA*(1 - Vsh), clipped 0..0.40.

Then write two functions:

archie_sw(phie, rt, a, m, n, Rw), vectorized,

Sw = ((a*Rw)/(PHIEm * RT))(1/n), clipped 0..1. Only trust it where PHIE > 0.01; elsewhere return Sw = 1.0 (no porosity, no hydrocarbon).

pay_zone_mean_sw(logs, a, m, n, Rw): build PHIE, take the rows with

7400 <= DEPTH < 7600, and return the mean Sw over the pay zone.

Run three cases on the pay zone and store each mean Sw:

Base: a=0.81, m=2.0, n=2.0, Rw=0.04
Higher cementation: m=2.3 (others base)
Higher water resistivity: Rw=0.08 (others base)

Put them in a dict sens = {"base": ..., "high_m": ..., "high_rw": ...}, and set moves_sw_most to the string "high_m" or "high_rw", whichever case shifts Sw furthest from the base case. Both raising m and raising Rw make the zone look less favourable (higher Sw); the question is which one your reserve number is more exposed to.

lightbulbHints (0/3)

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

codeYour solution

main.py

import numpy as np
import pandas as pd

np.random.seed(7)
depth = np.arange(7000.0, 7850.0, 0.5)
_cond = [depth < 7200, depth < 7400, depth < 7600, depth < 7700, depth >= 7700]
_gr   = np.select(_cond, [105.0, 40.0, 32.0, 110.0, 45.0])
_rt   = np.select(_cond, [4.0,   2.2,  35.0, 4.0,   2.5])
_rhob = np.select(_cond, [2.54,  2.31, 2.26, 2.55,  2.33])
_nphi = np.select(_cond, [0.33,  0.21, 0.17, 0.34,  0.22])
logs = pd.DataFrame({
    "DEPTH": depth,
    "GR":   _gr   + np.random.normal(0, 4.0,  len(depth)),
    "RT":   _rt   * np.exp(np.random.normal(0, 0.10, len(depth))),
    "RHOB": _rhob + np.random.normal(0, 0.02, len(depth)),
    "NPHI": np.clip(_nphi + np.random.normal(0, 0.015, len(depth)), 0.0, 1.0),
})
STEP_FT = 0.5

PAY_TOP_FT = 7400.0
PAY_BASE_FT = 7600.0

def build_phie(logs):
    """Effective porosity over the whole log (linear Vsh -> PHID -> PHIE)."""
    # TODO: gr_clean / gr_shale from the 5th / 95th GR percentiles
    # TODO: IGR -> Vsh (clip 0..1)
    # TODO: PHID = (2.65 - RHOB)/(2.65 - 1.0), clip 0..0.45
    # TODO: PHIA = (PHID + NPHI)/2 ; PHIE = PHIA*(1 - Vsh), clip 0..0.40
    pass

def archie_sw(phie, rt, a, m, n, Rw):
    """Vectorized Archie Sw, clipped 0..1; Sw = 1.0 where PHIE <= 0.01."""
    # TODO: Sw = ((a*Rw)/(phie**m * rt))**(1/n), guard PHIE <= 0.01 -> 1.0
    pass

def pay_zone_mean_sw(logs, a, m, n, Rw):
    """Mean Archie Sw over the 7400-7600 ft pay zone."""
    # TODO: build PHIE, mask 7400 <= DEPTH < 7600, return mean Sw
    pass

sens = {
    "base":    None,  # TODO a=.81, m=2.0, n=2.0, Rw=.04
    "high_m":  None,  # TODO m=2.3
    "high_rw": None,  # TODO Rw=.08
}
moves_sw_most = None  # TODO "high_m" or "high_rw"

print("pay-zone mean Sw:", sens)
print("parameter that moves Sw most:", moves_sw_most)

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
import pandas as pd

np.random.seed(7)
depth = np.arange(7000.0, 7850.0, 0.5)
_cond = [depth < 7200, depth < 7400, depth < 7600, depth < 7700, depth >= 7700]
_gr   = np.select(_cond, [105.0, 40.0, 32.0, 110.0, 45.0])
_rt   = np.select(_cond, [4.0,   2.2,  35.0, 4.0,   2.5])
_rhob = np.select(_cond, [2.54,  2.31, 2.26, 2.55,  2.33])
_nphi = np.select(_cond, [0.33,  0.21, 0.17, 0.34,  0.22])
logs = pd.DataFrame({
    "DEPTH": depth,
    "GR":   _gr   + np.random.normal(0, 4.0,  len(depth)),
    "RT":   _rt   * np.exp(np.random.normal(0, 0.10, len(depth))),
    "RHOB": _rhob + np.random.normal(0, 0.02, len(depth)),
    "NPHI": np.clip(_nphi + np.random.normal(0, 0.015, len(depth)), 0.0, 1.0),
})
STEP_FT = 0.5

PAY_TOP_FT = 7400.0
PAY_BASE_FT = 7600.0


def build_phie(logs):
    """Effective porosity over the whole log (linear Vsh -> PHID -> PHIE)."""
    gr_clean = logs["GR"].quantile(0.05)
    gr_shale = logs["GR"].quantile(0.95)
    igr = (logs["GR"] - gr_clean) / (gr_shale - gr_clean)
    vsh = igr.clip(0.0, 1.0)

    phid = ((2.65 - logs["RHOB"]) / (2.65 - 1.0)).clip(0.0, 0.45)
    phia = (phid + logs["NPHI"]) / 2.0
    phie = (phia * (1.0 - vsh)).clip(0.0, 0.40)
    return phie


def archie_sw(phie, rt, a, m, n, Rw):
    """Vectorized Archie Sw, clipped 0..1; Sw = 1.0 where PHIE <= 0.01."""
    phie = np.asarray(phie, dtype=float)
    rt = np.asarray(rt, dtype=float)
    with np.errstate(divide="ignore", invalid="ignore"):
        sw = ((a * Rw) / (phie ** m * rt)) ** (1.0 / n)
    sw = np.where(phie > 0.01, sw, 1.0)
    return np.clip(sw, 0.0, 1.0)


def pay_zone_mean_sw(logs, a, m, n, Rw):
    """Mean Archie Sw over the 7400-7600 ft pay zone."""
    phie = build_phie(logs)
    mask = (logs["DEPTH"] >= PAY_TOP_FT) & (logs["DEPTH"] < PAY_BASE_FT)
    sw = archie_sw(phie[mask].to_numpy(), logs["RT"][mask].to_numpy(),
                   a, m, n, Rw)
    return float(np.mean(sw))


sens = {
    "base":    pay_zone_mean_sw(logs, a=0.81, m=2.0, n=2.0, Rw=0.04),
    "high_m":  pay_zone_mean_sw(logs, a=0.81, m=2.3, n=2.0, Rw=0.04),
    "high_rw": pay_zone_mean_sw(logs, a=0.81, m=2.0, n=2.0, Rw=0.08),
}

_d_m = abs(sens["high_m"] - sens["base"])
_d_rw = abs(sens["high_rw"] - sens["base"])
moves_sw_most = "high_rw" if _d_rw > _d_m else "high_m"

print("pay-zone mean Sw:", sens)
print("parameter that moves Sw most:", moves_sw_most)

lockCopying code is a Full Access feature.

arrow_back

7.3 Density-Neutron Crossplot

7.5 Net Pay Cutoff Sensitivity

arrow_forward