Exercise 7.9
Shaly-Sand Sw: Simandoux vs Archie
Archie's equation is the field standard for water saturation, but it assumes clean sand: it treats shale-bound water as if it were nothing. In a shaly sand the clay carries its own conductivity, which Archie blames on extra formation water, so Archie reads optimistically high Sw and can bury a pay zone under a fake "wet" reading.
The Simandoux equation adds the shale conductivity term:
For n = 2 this is quadratic in Sw. Write A = phie**m / (a*Rw) and B = vsh / Rsh; then A·Sw² + B·Sw − 1/Rt = 0, whose physical root is
You are evaluating the OIL PAY sand in OD-007 (OML 58). This particular sand is fairly clean, so Archie and Simandoux nearly agree on it; but the moment a sand carries real dispersed clay the gap opens fast, and that is exactly where Simandoux earns its keep. Implement both and measure the gap.
Implement three functions:
archie_sw(phie, rt, a, m, Rw): standard Archie withn = 2:
Sw = sqrt((a*Rw) / (phie*m rt)), clipped to 0..1. Vectorized over numpy arrays.
simandoux_sw(phie, rt, vsh, a, m, Rw, Rsh): the quadratic root above,
clipped to 0..1. Vectorized.
pay_zone(logs, gr_clean, gr_shale, a, m, Rw, Rsh): over the OIL PAY
sand (DEPTH in [7400, 7600)) compute Vsh (linear), PHIE, then both saturations, and return a dict with sw_archie_mean, sw_simandoux_mean, vsh_mean, and sw_uplift (= sw_archie_mean − sw_simandoux_mean).
The endpoints, porosity recipe, and parameters are defined for you in the starter. Simandoux always lands at or below Archie: the gap is the shale conductivity Archie ignores. In a clean sand it's a fraction of a saturation unit; in a shaly sand (Vsh ≈ 0.3) it can be 10+ saturation units. The difference between booking the sand and walking away from it.
Stuck? Reveal hints one at a time — they progress from nudge to near-solution.
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
# --- Petrophysical parameters (OD-007, OML 58) -------------------------------
GR_CLEAN = float(logs["GR"].quantile(0.05)) # clean-sand GR endpoint
GR_SHALE = float(logs["GR"].quantile(0.95)) # shale GR endpoint
RHO_MATRIX = 2.65
RHO_FLUID = 1.0
A = 0.81 # Archie tortuosity
M = 2.0 # cementation exponent
RW = 0.04 # formation-water resistivity (ohm-m)
RSH = 2.0 # shale resistivity (ohm-m)
def vshale_linear(gr, gr_clean, gr_shale):
"""Linear Vshale from the gamma-ray index, clipped 0..1."""
igr = (gr - gr_clean) / (gr_shale - gr_clean)
return np.clip(igr, 0.0, 1.0)
def porosity_effective(rhob, nphi, vsh):
"""PHIE = mean(density, neutron) porosity de-rated by shale, clipped 0..0.40."""
phid = (RHO_MATRIX - rhob) / (RHO_MATRIX - RHO_FLUID)
phid = np.clip(phid, 0.0, 0.45)
phia = (phid + nphi) / 2.0
phie = phia * (1.0 - vsh)
return np.clip(phie, 0.0, 0.40)
def archie_sw(phie, rt, a, m, Rw):
"""Archie Sw with n=2: sqrt((a*Rw)/(phie**m * rt)), clipped 0..1."""
phie = np.asarray(phie, dtype=float)
rt = np.asarray(rt, dtype=float)
sw = np.sqrt((a * Rw) / (phie ** m * rt))
return np.clip(sw, 0.0, 1.0)
def simandoux_sw(phie, rt, vsh, a, m, Rw, Rsh):
"""Simandoux Sw (n=2), the physical root of A*Sw**2 + B*Sw - 1/rt = 0,
with A = phie**m/(a*Rw) and B = vsh/Rsh. Clip 0..1."""
phie = np.asarray(phie, dtype=float)
rt = np.asarray(rt, dtype=float)
vsh = np.asarray(vsh, dtype=float)
A_coef = phie ** m / (a * Rw)
B_coef = vsh / Rsh
sw = (-B_coef + np.sqrt(B_coef ** 2 + 4.0 * A_coef / rt)) / (2.0 * A_coef)
return np.clip(sw, 0.0, 1.0)
def pay_zone(logs, gr_clean, gr_shale, a, m, Rw, Rsh):
"""Over the OIL PAY sand (7400 <= DEPTH < 7600) return a summary dict with
sw_archie_mean, sw_simandoux_mean, vsh_mean, sw_uplift."""
sand = logs[(logs["DEPTH"] >= 7400.0) & (logs["DEPTH"] < 7600.0)]
vsh = vshale_linear(sand["GR"].to_numpy(), gr_clean, gr_shale)
phie = porosity_effective(sand["RHOB"].to_numpy(), sand["NPHI"].to_numpy(), vsh)
rt = sand["RT"].to_numpy()
sw_a = archie_sw(phie, rt, a, m, Rw)
sw_s = simandoux_sw(phie, rt, vsh, a, m, Rw, Rsh)
sw_archie_mean = float(np.mean(sw_a))
sw_simandoux_mean = float(np.mean(sw_s))
return {
"sw_archie_mean": sw_archie_mean,
"sw_simandoux_mean": sw_simandoux_mean,
"vsh_mean": float(np.mean(vsh)),
"sw_uplift": sw_archie_mean - sw_simandoux_mean,
}
result = pay_zone(logs, GR_CLEAN, GR_SHALE, A, M, RW, RSH)
print(f"Vsh mean : {result['vsh_mean']:.3f}")
print(f"Sw (Archie) : {result['sw_archie_mean']:.3f}")
print(f"Sw (Simandoux) : {result['sw_simandoux_mean']:.3f}")
print(f"Archie uplift : {result['sw_uplift']:.3f}")
lockCopying code is a Full Access feature.