Wellbore Stability and Geomechanics Calculator

In field units (psi), with depth z in feet and density RHOB in g/cc:

S_v(psi) = integral[RHOB(z) * 0.4335 * dz]

Where 0.4335 converts g/cc * ft to psi (= 0.4335 psi per ft per g/cc).

Simplified (average density):

S_v = rho_avg * 0.4335 * z    (psi, z in ft, rho in g/cc)

Typical Overburden Gradients by Depth

Log-Based Overburden Calculation

import math

def overburden_from_log(depths_ft, rhob_gcc):
    """
    Calculate overburden stress from density log.
    depths_ft: list of depth values (measured from surface, ft)
    rhob_gcc: list of bulk density values at each depth (g/cc)
    Returns: list of S_v values in psi at each depth
    """
    conv = 0.4335  # psi/(ft * g/cc)
    Sv = [0.0]
    for i in range(1, len(depths_ft)):
        dz = depths_ft[i] - depths_ft[i-1]
        rho_avg = (rhob_gcc[i] + rhob_gcc[i-1]) / 2
        Sv.append(Sv[-1] + rho_avg * conv * dz)
    return Sv

# Example: 3-layer overburden for 8,000 ft well
depths = [0, 2000, 5000, 8000]
rhob   = [2.10, 2.45, 2.60, 2.65]  # g/cc (increases with compaction)

Sv = overburden_from_log(depths, rhob)
print("Depth (ft) | Sv (psi) | Sv gradient (psi/ft)")
print("-" * 50)
for i, (d, sv) in enumerate(zip(depths, Sv)):
    grad = sv / d if d > 0 else 0
    print(f"{d:10,} | {sv:8,.0f} | {grad:.4f}")

Effective Overburden Stress

The effective overburden (Terzaghi effective stress principle):

sigma_v = S_v - alpha * P_p

Where alpha = Biot coefficient (0.6-1.0 for porous rock; 1.0 is conservative).

For Marcellus shale: alpha = 0.7-0.9 (lower than clean sandstone due to low porosity).

Module 2 — Pore Pressure

Capabilities

• Hydrostatic pore pressure for normal-pressured formations
• Eaton method for pore pressure from sonic or resistivity logs
• Equivalent mud weight (EMW) conversion
• PNGE applications: effective stress calculation, mud weight floor, cap rock integrity

Normal Pore Pressure Gradients

Marcellus shale context (WV): The Marcellus at 6,000-8,500 ft TVD in West Virginia is generally normally pressured or slightly overpressured in the core of the basin. Use 0.46-0.47 psi/ft as a reasonable default.

Eaton Method (Resistivity-Based)

Overpressure from compaction disequilibrium can be detected from resistivity log trends:

P_p = S_v - (S_v - P_p_normal) * (Ro/Rn)^1.2

Where:

• Ro = observed resistivity at depth of interest (ohm-m)
• Rn = normal (expected) resistivity at same depth from trend line (ohm-m)
• Exponent 1.2 is Eaton's empirical value (commonly used; range 0.6-2.0)

Sonic (interval transit time) version:

P_p = S_v - (S_v - P_p_normal) * (DTCn/DTC)^3.0

Where DTC = observed transit time, DTCn = normal trend transit time.

def eaton_resistivity(Sv_psi, Pp_normal_psi, Ro, Rn, n=1.2):
    """
    Eaton pore pressure from resistivity log.
    Sv_psi: overburden stress (psi)
    Pp_normal_psi: normal pore pressure at this depth (psi)
    Ro: observed resistivity (ohm-m)
    Rn: normal trend resistivity (ohm-m)
    n: Eaton exponent (default 1.2)
    """
    Pp = Sv_psi - (Sv_psi - Pp_normal_psi) * (Ro / Rn)**n
    return Pp

# Example: 7,500 ft depth
Sv = 7875
Pp_n = 7500 * 0.46
Ro = 8.0   # observed resistivity
Rn = 12.0  # normal trend (higher = more compacted = normal)

Pp_calc = eaton_resistivity(Sv, Pp_n, Ro, Rn)
EMW = Pp_calc / (7500 * 0.052)
print(f"Eaton pore pressure: {Pp_calc:.0f} psi")
print(f"Equivalent mud weight: {EMW:.2f} ppg")

Equivalent Mud Weight (EMW)

Convert any pressure to equivalent mud weight at a given depth:

EMW (ppg) = P (psi) / (0.052 * depth_ft)

Module 3 — Horizontal Stress Estimation

Capabilities

• Minimum horizontal stress (S_hmin) from poroelastic model
• Maximum horizontal stress (S_Hmax) from stress polygon (Zoback frictional limit)
• Stress regime classification
• Appalachian reference stress gradients and stress orientation
• PNGE applications: hydraulic fracture orientation, wellbore stability design

Poroelastic Model for S_hmin

The minimum horizontal stress from poroelastic theory (for a laterally constrained basin with no tectonic strain):

S_hmin = [nu/(1-nu)] * (S_v - alpha*P_p) + alpha*P_p + S_tect

Where:

• nu = Poisson's ratio of the rock
• alpha = Biot coefficient
• P_p = pore pressure
• S_tect = tectonic stress component (can be positive or negative)

For a basin with no active tectonic loading (S_tect = 0):

S_hmin = [nu/(1-nu)] * (S_v - alpha*P_p) + alpha*P_p

Frictional Equilibrium (Zoback Stress Polygon)

The maximum horizontal stress is bounded by the frictional strength of pre-existing faults. For coefficient of friction mu_f:

(S1/S3)_max = [(mu_f^2 + 1)^0.5 + mu_f]^2

Typical rock friction coefficient: mu_f = 0.6-0.8 (Byerlee's law).

Appalachian/Marcellus Reference Stress State

The Appalachian basin is in a strike-slip faulting regime in the deep Devonian/Mississippian section where the Marcellus resides:

S_Hmax > S_v > S_hmin

Appalachian reference stress gradients (TVD depth range 6,000-9,000 ft):

Source: Multiple published studies on Appalachian/Marcellus geomechanics (Lash and Engelder, 2011; Gale et al., 2014; Zoback et al., various).

def horizontal_stress_poroelastic(Sv, Pp, nu, alpha, S_tect=0):
    """
    Estimate minimum horizontal stress using poroelastic model.
    All pressures in psi, same depth.
    S_tect: tectonic stress additive term (psi). Positive = compression.
    """
    sigma_v_eff = Sv - alpha * Pp
    Shmin = (nu / (1 - nu)) * sigma_v_eff + alpha * Pp + S_tect
    return Shmin

# Marcellus at 7,500 ft TVD
depth = 7500
Sv = depth * 1.02
Pp = depth * 0.465
nu = 0.25
alpha = 0.8

Shmin_pe = horizontal_stress_poroelastic(Sv, Pp, nu, alpha)
Shmin_ref_lo = depth * 0.65
Shmin_ref_hi = depth * 0.75
Shmax_est = depth * 0.87

print(f"Depth: {depth} ft")
print(f"S_v:   {Sv:,.0f} psi  ({Sv/depth:.3f} psi/ft)")
print(f"P_p:   {Pp:,.0f} psi  ({Pp/depth:.3f} psi/ft)")
print(f"\nPoroelastic S_hmin (no tectonic): {Shmin_pe:,.0f} psi ({Shmin_pe/depth:.3f} psi/ft)")
print(f"Appalachian reference range:      {Shmin_ref_lo:,.0f}-{Shmin_ref_hi:,.0f} psi")
print(f"Estimated S_Hmax:                 {Shmax_est:,.0f} psi ({Shmax_est/depth:.3f} psi/ft)")

if Sv > Shmax_est > Shmin_pe:
    regime = "NORMAL FAULTING (S_v > S_Hmax > S_hmin)"
elif Shmax_est > Sv > Shmin_pe:
    regime = "STRIKE-SLIP (S_Hmax > S_v > S_hmin)"
elif Shmax_est > Shmin_pe > Sv:
    regime = "REVERSE FAULTING (S_Hmax > S_hmin > S_v)"