Petrophysics and Formation Evaluation

def vshale(igr, method="larionov_older"):
    if method == "linear":
        return igr
    elif method == "larionov_older":
        return 0.33 * (2**(2.0 * igr) - 1)
    elif method == "larionov_clastic":
        return 0.083 * (2**(3.7 * igr) - 1)
    elif method == "steiber":
        return igr / (3.0 - 2.0 * igr)
    else:
        raise ValueError(f"Unknown method: {method}")

Appalachian guidance:

• Marcellus Shale: use larionov_older (Devonian). GR_clean ≈ 20 API (Onondaga limestone), GR_shale ≈ 200–280 API
• Utica/Point Pleasant: GR_clean ≈ 15–25, GR_shale ≈ 150–200 API

Module 2 — Porosity

Density Porosity

def phi_density(rhob, rho_matrix=2.65, rho_fluid=1.0):
    """
    phi_D = (rho_matrix - rhob) / (rho_matrix - rho_fluid)
    rhob: bulk density from RHOB log (g/cc)
    rho_matrix: 2.65 sandstone, 2.71 limestone, 2.87 dolomite
    rho_fluid: 1.0 freshwater mud, 1.1 saltwater mud
    """
    return (rho_matrix - rhob) / (rho_matrix - rho_fluid)

Sonic Porosity (Wyllie Time-Average)

def phi_sonic(dt, dt_matrix=55.5, dt_fluid=189.0, cp=1.0):
    """
    phi_S = cp * (dt - dt_matrix) / (dt_fluid - dt_matrix)
    dt: interval transit time from DT log (us/ft)
    dt_matrix: 55.5 sandstone, 47.5 limestone, 43.5 dolomite
    dt_fluid: 189 freshwater, 185 saltwater
    cp: compaction correction (1.0 normally consolidated)
    """
    return cp * (dt - dt_matrix) / (dt_fluid - dt_matrix)

def cp_correction(dt_shale, dt_shale_normal=100.0):
    """Hilchie compaction correction: cp = dt_shale / 100"""
    return dt_shale / dt_shale_normal

Neutron-Density Crossplot Porosity

import math

def phi_nd(phi_n, phi_d):
    """
    Gaymard-Poupon crossplot: phi_ND = sqrt((phi_N^2 + phi_D^2) / 2)
    phi_n: neutron porosity (decimal, NPHI log)
    phi_d: density porosity (decimal)
    In gas zones: phi_N decreases, phi_D increases — crossplot minimizes bias.
    """
    return math.sqrt((phi_n**2 + phi_d**2) / 2.0)

Matrix Parameters Reference

Module 3 — Water Saturation

Archie Equation (clean formations)

def sw_archie(rt, phi, rw, a=1.0, m=2.0, n=2.0):
    """
    Sw^n = a * Rw / (phi^m * Rt)
    rt:  true resistivity (ohm-m) from deep log
    phi: porosity (decimal)
    rw:  formation water resistivity (ohm-m)
    a=1.0, m=2.0, n=2.0: tortuosity, cementation, saturation exponents
    """
    return (a * rw / (phi**m * rt))**(1.0 / n)

def rw_from_salinity(salinity_ppm, temp_f):
    """Approximate Rw from NaCl-equivalent salinity and temperature (Arps)."""
    rw_75 = (400000.0 / salinity_ppm)**0.88
    return rw_75 * (75.0 + 6.77) / (temp_f + 6.77)

Simandoux (shaly sand)

def sw_simandoux(rt, phi, rw, vsh, rsh, m=2.0):
    """
    Solves: phi^m * Sw^2 / (0.4*Rw) + Vsh*Sw / Rsh = 1 / Rt
    rsh: shale resistivity (ohm-m)
    vsh: shale volume (decimal)
    Returns the lower root (Sw solution).
    """
    A = phi**m / (0.4 * rw)
    B = vsh / rsh
    C = 1.0 / rt
    discriminant = B**2 + 4 * A * C
    return (-B + discriminant**0.5) / (2.0 * A)

Appalachian Formation Water Resistivity Reference

Module 4 — Irreducible Sw and Moveable Hydrocarbons

Bulk Volume Water (Buckles Plot)

def bvw(sw, phi):
    """
    BVW = Sw * phi
    Constant BVW with depth -> formation at irreducible Sw (no water production).
    Varying BVW -> transition zone, expect water production.
    """
    return sw * phi

Timur Equation (permeability from porosity and Swirr)

def permeability_timur(phi, swirr):
    """
    Timur (1968): k = 0.136 * phi^4.4 / swirr^2  (mD)
    Rearranged for Swirr: swirr = sqrt(0.136 * phi^4.4 / k)
    """
    return 0.136 * phi**4.4 / swirr**2

def swirr_timur(phi, k):
    """Irreducible water saturation from porosity and permeability."""
    return (0.136 * phi**4.4 / k)**0.5

Moveable Hydrocarbon Index

def mhi(sw, sxo):
    """
    MHI = Sw / Sxo
    sxo: flushed zone saturation from shallow resistivity (Rxo)
    MHI < 0.7: good mobility   MHI > 0.9: tight, little movability
    """
    return sw / sxo

Module 5 — Brittleness Index (Completion Quality)

Rickman et al. (2008) Brittleness

def brittleness_rickman(ym_gpa, pr,
                         ym_brit=70.0, ym_duct=15.0,
                         pr_brit=0.15, pr_duct=0.40):
    """
    Brittleness (%) = 50 * (YM_norm + PR_norm) * 100
    ym_gpa: dynamic Young's modulus (GPa) from sonic
    pr: Poisson's ratio from sonic
    High brittleness (>40%): good frac candidate
    Low brittleness (<20%): ductile, poor frac response
    """
    ym_n = max(0.0, min(1.0, (ym_gpa - ym_duct) / (ym_brit - ym_duct)))
    pr_n = max(0.0, min(1.0, (pr - pr_duct) / (pr_brit - pr_duct)))
    return 0.5 * (ym_n + pr_n) * 100.0

Dynamic Elastic Moduli from Logs

def dynamic_moduli(dtc, dts, rhob):
    """
    dtc: compressional slowness (us/ft)
    dts: shear slowness (us/ft)
    rhob: bulk density (g/cc)
    Returns: (E_gpa, nu) dynamic Young's modulus and Poisson's ratio
    """
    vp = 304800.0 / dtc   # ft/s
    vs = 304800.0 / dts   # ft/s
    rho = rhob * 1000.0   # kg/m3
    vp_m = vp * 0.3048    # m/s
    vs_m = vs * 0.3048    # m/s
    mu = rho * vs_m**2    # shear modulus (Pa)
    r2 = (vp_m / vs_m)**2
    E = mu * (3*r2 - 4) / (r2 - 1) / 1e9   # GPa
    nu = (r2 - 2.0) / (2.0 * (r2 - 1.0))
    return E, nu

Note: Dynamic-to-static correction: E_static ≈ 0.4–0.7 × E_dynamic for shales. Use core measurements to calibrate before geomechanical modeling.

Marcellus brittleness reference:

Module 6 — Net Pay Cutoffs

def net_pay(depths, vsh_log, phi_log, sw_log,
            vsh_cut=0.35, phi_cut=0.05, sw_cut=0.60,
            dz=0.5):
    """
    Returns gross thickness, net pay thickness, and N/G ratio.
    dz: sample interval (ft) — default 0.5 ft for digital logs
    """
    gross = len(depths) * dz
    pay_intervals = [(d, vsh, phi, sw)
                     for d, vsh, phi, sw
                     in zip(depths, vsh_log, phi_log, sw_log)
                     if vsh <= vsh_cut and phi >= phi_cut and sw <= sw_cut]
    net = len(pay_intervals) * dz
    return {
        "gross_ft": gross,
        "net_ft": net,
        "net_gross": net / gross if gross > 0 else 0,
        "pay_intervals": pay_intervals
    }

Typical Appalachian cutoffs:

Output Format

## Formation Evaluation — [Formation], [Well Name/API]
**Interval:** X,XXX – X,XXX ft MD | **Logged:** [curves available]

### Log Quality Check
| Curve | Min | Max | Mean | Issues |
|-------|-----|-----|------|--------|
| GR (API) | | | | |
| RHOB (g/cc) | | | | |
| NPHI (frac) | | | | |
| RT (ohm-m) | | | | |

### Petrophysical Summary (Net Pay)
| Property | Value | Method |
|----------|-------|--------|
| Net pay (ft) | | Vsh/phi/Sw cutoffs |
| Average phi_ND | | N-D crossplot |
| Average Sw | | Archie |
| Swirr estimate | | Buckles / Timur |
| Brittleness avg (%) | | Rickman |

### Parameter Assumptions
- GR_clean / GR_shale: XX / XX API
- Matrix: [sandstone / limestone / mixed]
- Rw: X.XX ohm-m (source: [DST / database / salinity gradient])
- a / m / n: 1.0 / 2.0 / 2.0 (note if shaly-sand correction applied)

### Confidence: [HIGH / MEDIUM / LOW]
[Note any data quality issues, missing curves, or calibration gaps]

Error Handling

Caveats

• Archie's equation assumes shale-free formation. In organic shales (Marcellus, Utica), TOC contributes independent conductivity — Archie Sw is unreliable. Use as a relative indicator only; calibrate to core.
• Dynamic Young's modulus from sonic logs exceeds static lab values by 20–100%. Apply dynamic-to-static correction before using in geomechanical or frac models.
• Larionov corrections reduce IGR-to-Vsh conversion; linear method overestimates shale content and underestimates net pay.
• All cutoff values are formation-specific. Values here are starting points — always calibrate to production data and core measurements.
• Log quality in horizontal wells can be poor due to invasion, tool standoff, and borehole geometry; evaluate LQC before interpretation.