pycsamt.interp.petrophysics#

Petrophysical transforms for EM hydrogeophysics.

Converts EM-derived resistivity (from AMT, MT, EMAP, TDEM inversions) into quantitative hydrogeological parameters: saturation, porosity, hydraulic conductivity, transmissivity, and water chemistry indicators.

All transforms are bidirectional — forward (hydro → ρ) and inverse (ρ → hydro) — and accept numpy arrays so they can be applied cell-by-cell to a full 2-D resistivity section from ResistivityModel.

Background#

Two petrophysical frameworks are implemented:

Archie (1942) — clean sandstones and crystalline rocks:

\[\rho = a \cdot \rho_w \cdot \phi^{-m} \cdot S_w^{-n}\]

where a is the tortuosity factor (≈ 1), m the cementation exponent (1.3–2.5), and n the saturation exponent (1.8–2.5).
Waxman-Smits (1968) — clay-bearing sediments (EMAP, shallow AMT):

\[\sigma = \frac{S_w^n}{F}(\sigma_w + \sigma_s / S_w)\]

where σ_s is the surface (clay) conductivity.

Additional module-level functions cover:

Kozeny-Carman hydraulic conductivity from porosity
Transmissivity and storativity from layer geometry
Water chemistry: ρ_w ↔ TDS ↔ EC
Hashin-Shtrikman effective-medium bounds
EM-specific helpers: Bostick depth (TDEM/AMT), skin depth, fracture-zone K

References

Functions

`aquifer_top_from_profile`(rho_log10, z_centers, *)	Detect the top of an aquifer/conductor in a 1-D resistivity profile.
`bostick_depth`(rho_a, freq)	Bostick (1977) pseudo-depth for MT/AMT/TDEM data.
`ec_mscm_to_rho`(ec, *[, temp_c])	Electrical conductivity (mS/cm) → resistivity (Ω·m).
`fractured_zone_K`(rho, *[, rho_matrix, ...])	Hydraulic conductivity of a fractured zone from EM resistivity.
`kozeny_carman_K`(phi, *[, d50_m, C, T, ...])	Hydraulic conductivity from porosity via Kozeny-Carman equation.
`rho_to_ec_mscm`(rho, *[, temp_c])	Formation/water resistivity (Ω·m) → EC (mS/cm).
`rho_to_hydraulic_conductivity`(rho, archie, *)	Estimate hydraulic conductivity K from EM resistivity.
`rho_w_to_tds`(rho_w, *[, temp_c])	Pore-water resistivity → Total Dissolved Solids (mg/L).
`skin_depth`(rho, freq)	EM skin depth (penetration depth) in metres.
`storativity`(phi, thickness, *[, ...])	Aquifer storativity \(S = S_s \cdot b\) (confined) or φ (unconfined).
`tds_to_rho_w`(tds, *[, temp_c])	TDS (mg/L) → pore-water resistivity (Ω·m).
`transmissivity`(K, thickness)	Aquifer transmissivity \(T = K \cdot b\) (m²/s).
`water_table_from_profile`(rho_log10, ...[, ...])	Estimate water table depth from a 1-D EM resistivity column.

Classes

`ArchieModel`([m, n, a])	Archie (1942) petrophysical model.
`HashinShtrikmanBounds`(rho_matrix, rho_fluid)	Hashin-Shtrikman (1962) bounds on effective resistivity.
`WaxmanSmitsModel`([m, n, a, sigma_s])	Waxman-Smits (1968) model for clay-bearing formations.

class pycsamt.interp.petrophysics.ArchieModel(m=1.8, n=2.0, a=1.0)[source]#

Bases: PyCSAMTObject

Archie (1942) petrophysical model.

Relates formation resistivity to porosity, water saturation, and pore-water resistivity via:

\[\rho = a \cdot \rho_w \cdot \phi^{-m} \cdot S_w^{-n}\]

Typical parameter ranges#

param m:: Cementation exponent.
type m:: float
param n:: Saturation exponent.
type n:: float
param a:: Tortuosity factor (default 1.0).
type a:: float

Examples

>>> archie = ArchieModel(m=1.8, n=2.0, a=1.0)
>>> archie.forward(phi=0.30, Sw=1.0, rho_w=0.025)   # fully saturated
0.309...
>>> archie.saturation(rho=1.2, phi=0.30, rho_w=0.025)
0.508...

m: float = 1.8#

n: float = 2.0#

a: float = 1.0#

formation_factor(phi)[source]#

Formation factor \(F = a \phi^{-m}\).

Parameters:: phi (array-like) – Porosity (fraction, 0–1).
Returns:: F
Return type:: ndarray

forward(phi, Sw, rho_w)[source]#

Formation resistivity from Archie’s law.

Parameters:

phi (array-like) – Porosity (fraction, 0–1).
Sw (array-like) – Water saturation (fraction, 0–1).
rho_w (float or array-like) – Pore-water resistivity (Ω·m).

Returns:

rho

Return type:

ndarray — formation resistivity (Ω·m)

saturation(rho, phi, rho_w)[source]#

Water saturation from measured resistivity (Archie inverse).

\[S_w = \left(\frac{a \rho_w}{\rho \phi^m}\right)^{1/n}\]

Parameters:

rho (array-like) – Formation resistivity (Ω·m).
phi (array-like) – Porosity (fraction, 0–1).
rho_w (float or array-like) – Pore-water resistivity (Ω·m).

Returns:

Sw

Return type:

ndarray — water saturation (0–1)

porosity(rho, Sw, rho_w)[source]#

Porosity from resistivity and saturation (Archie inverse).

\[\phi = \left(\frac{a \rho_w S_w^{-n}}{\rho}\right)^{1/m}\]

Parameters:

rho (array-like) – Formation resistivity (Ω·m).
Sw (array-like) – Water saturation (fraction, 0–1).
rho_w (float or array-like) – Pore-water resistivity (Ω·m).

Returns:

phi

Return type:

ndarray — porosity (0–1)

fluid_resistivity(rho, phi, Sw)[source]#

Pore-water resistivity from formation resistivity (Archie inverse).

\[\rho_w = \frac{\rho \phi^m S_w^n}{a}\]

Parameters:

rho (array-like) – Formation resistivity (Ω·m).
phi (array-like) – Porosity (fraction, 0–1).
Sw (array-like) – Water saturation (fraction, 0–1).

Returns:

rho_w

Return type:

ndarray — pore-water resistivity (Ω·m)

water_content(phi, Sw)[source]#

Volumetric water content \(\theta = \phi \cdot S_w\).

Parameters:

phi (Any)
Sw (Any)

Return type:

Any

Parameters:

m (float)
n (float)
a (float)

class pycsamt.interp.petrophysics.WaxmanSmitsModel(m=1.8, n=2.0, a=1.0, sigma_s=0.0)[source]#

Bases: PyCSAMTObject

Waxman-Smits (1968) model for clay-bearing formations.

Extends Archie’s law with a surface conductivity term σ_s that represents the excess conductance associated with clay minerals (cation exchange capacity). Use this model when clay content > ~10 % or when Archie overestimates saturation in shaly formations.

\[\sigma = \frac{S_w^n}{F}\bigl(\sigma_w + \sigma_s / S_w\bigr)\]

where F = a φ^{-m} is the formation factor.

Parameters:

m (float) – Cementation exponent (default 1.8).
n (float) – Saturation exponent (default 2.0).
a (float) – Tortuosity factor (default 1.0).
sigma_s (float) – Surface (clay) conductivity S/m (default 0 → reduces to Archie).

Notes

σ_s can be estimated from CEC (cation exchange capacity) measurements or calibrated against borehole resistivity logs. Values for AMT/MT targets typically range 0.001–0.05 S/m for moderate clay content.

Examples

>>> ws = WaxmanSmitsModel(m=1.8, n=2.0, sigma_s=0.01)
>>> ws.forward(phi=0.30, Sw=0.70, sigma_w=40.0)   # sigma_w in mS/m
8.31...  # Ω·m

m: float = 1.8#

n: float = 2.0#

a: float = 1.0#

sigma_s: float = 0.0#

forward(phi, Sw, sigma_w)[source]#

Formation resistivity from Waxman-Smits equation.

Parameters:

phi (array-like) – Porosity (fraction, 0–1).
Sw (array-like) – Water saturation (fraction, 0–1).
sigma_w (float or array-like) – Pore-water conductivity in mS/m (≈ 40 mS/m for fresh water). Internally converted to S/m.

Returns:

rho

Return type:

ndarray — formation resistivity (Ω·m)

saturation(rho, phi, sigma_w, *, tol=1e-08, max_iter=50)[source]#

Water saturation by numerical inversion of Waxman-Smits equation.

Parameters:

rho (array-like) – Formation resistivity (Ω·m).
phi (array-like) – Porosity (fraction, 0–1).
sigma_w (float or array-like) – Pore-water conductivity in mS/m.
tol (float) – Convergence tolerance for brentq (default 1e-8).
max_iter (int) – Maximum root-finding iterations (default 50).

Returns:

Sw

Return type:

ndarray — water saturation (0–1)

porosity(rho, Sw, sigma_w)[source]#

Porosity by direct inversion of Waxman-Smits (known saturation).

Rearranging WS:

\[F = \frac{S_w^n \,(\sigma_w + \sigma_s / S_w)}{\sigma_{obs}} \Rightarrow \phi = \left(F / a\right)^{-1/m}\]

Parameters:

rho (array-like) – Formation resistivity (Ω·m).
Sw (array-like) – Water saturation (fraction, 0–1).
sigma_w (float or array-like) – Pore-water conductivity in mS/m.

Returns:

phi

Return type:

ndarray — porosity (0–1)

class pycsamt.interp.petrophysics.HashinShtrikmanBounds(rho_matrix, rho_fluid)[source]#

Bases: PyCSAMTObject

Hashin-Shtrikman (1962) bounds on effective resistivity.

Provides theoretical upper and lower bounds for a two-phase medium (rock matrix + fluid) without assumptions about microstructure. The bounds are tighter than Voigt-Reuss and are used to validate whether an Archie-derived porosity is physically plausible.

For conductivity σ (= 1/ρ):

\[\sigma^{HS-} \leq \sigma_{eff} \leq \sigma^{HS+}\]

Parameters:

rho_matrix (float) – Resistivity of the rock matrix (Ω·m). Typical: granite 10⁴–10⁶, sandstone 100–1000, limestone 500–5000.
rho_fluid (float) – Resistivity of the pore fluid (Ω·m). Typical: fresh water 20–100, saline water 0.1–5.

Examples

>>> hs = HashinShtrikmanBounds(rho_matrix=1000.0, rho_fluid=25.0)
>>> lower, upper = hs.bounds(phi=0.25)
>>> print(f"{lower:.1f}  {upper:.1f}")

bounds(phi)[source]#

Compute Hashin-Shtrikman lower and upper conductivity bounds.

Parameters:

phi (array-like) – Porosity (volume fraction of fluid, 0–1).

Returns:

rho_lower (ndarray — lower resistivity bound (HS+, fluid-connected))
rho_upper (ndarray — upper resistivity bound (HS-, isolated fluid))

Return type:

tuple[ndarray, ndarray]

in_bounds(rho, phi, *, margin=0.0)[source]#

Return boolean mask: True where rho lies within HS bounds ± margin (log₁₀).

Parameters:

rho (Any)
phi (Any)
margin (float)

Return type:

ndarray

pycsamt.interp.petrophysics.kozeny_carman_K(phi, *, d50_m=0.00025, C=180.0, T=0.5, gravity=9.81, kinematic_viscosity=1e-06)[source]#

Hydraulic conductivity from porosity via Kozeny-Carman equation.

\[K = \frac{g}{\nu} \cdot \frac{d_{50}^2}{C} \cdot \frac{\phi^3}{(1-\phi)^2} \cdot T\]

where \(g/\nu\) converts intrinsic permeability (m²) to hydraulic conductivity (m/s): \(g = 9.81\) m/s², \(\nu \approx 10^{-6}\) m²/s at 20 °C.

Parameters:

phi (array-like) – Porosity (fraction, 0–1).
d50_m (float) – Median grain size in metres (default 2.5×10⁻⁴ m = 0.25 mm, upper fine sand / lower medium sand).
C (float) – Kozeny constant × shape factor (default 180 for spherical grains).
T (float) – Tortuosity correction factor (default 0.5; range 0.3–0.7).
gravity (float) – Gravitational acceleration (m/s²; default 9.81).
kinematic_viscosity (float) – Kinematic viscosity of water (m²/s; default 10⁻⁶ at 20 °C). Use 1.31×10⁻⁶ at 10 °C or 0.80×10⁻⁶ at 30 °C if temperature matters.

Returns:

K

Return type:

ndarray — hydraulic conductivity (m/s)

Notes

Typical K values at φ ≈ 0.30 with matching d50:

Clean gravel (d50 = 5 mm) : ~10⁻² – 10⁻¹ m/s
Coarse sand (d50 = 1 mm) : ~10⁻³ – 10⁻² m/s
Medium sand (d50 = 0.5 mm) : ~10⁻⁴ – 10⁻³ m/s
Fine sand (d50 = 0.1 mm) : ~10⁻⁶ – 10⁻⁵ m/s
Silty sand (d50 = 0.05 mm): ~10⁻⁷ – 10⁻⁶ m/s
Clay / silt : < 10⁻⁹ m/s

Examples

>>> kozeny_carman_K(0.30, d50_m=1e-3)   # coarse sand, φ = 0.30
1.16e-03  # m/s  (within typical coarse-sand range)

pycsamt.interp.petrophysics.rho_to_hydraulic_conductivity(rho, archie, *, rho_w=20.0, phi_prior=0.25, Sw=1.0, d50_m=0.00025, C=180.0, T=0.5, gravity=9.81, kinematic_viscosity=1e-06)[source]#

Estimate hydraulic conductivity K from EM resistivity.

Chain: ρ → φ (Archie inverse) → K (Kozeny-Carman).

This is valid only in fully or near-fully saturated zones (S_w ≈ 1). For unsaturated zones use the forward Archie model to estimate S_w first, then re-run on the saturated resistivity.

Parameters:

rho (array-like) – Formation resistivity (Ω·m) from EM inversion.
archie (ArchieModel) – Calibrated Archie model (m, n, a).
rho_w (float) – Pore-water resistivity (Ω·m). Default 20 Ω·m (fresh water at 25 °C, EC ≈ 0.5 mS/cm). Use ~0.2 Ω·m for seawater, ~2–5 Ω·m for brackish.
phi_prior (float) – Prior (reference) porosity used to clip unrealistic values.
Sw (float) – Water saturation assumed for the inversion (default 1.0 = saturated).
d50_m (float) – Median grain size (m) for Kozeny-Carman.
C (float) – Kozeny constant and tortuosity (see kozeny_carman_K()).
T (float) – Kozeny constant and tortuosity (see kozeny_carman_K()).
gravity (float) – Gravitational acceleration (m/s²; default 9.81).
kinematic_viscosity (float) – Kinematic viscosity of water (m²/s; default 10⁻⁶ at 20 °C).

Returns:

K

Return type:

ndarray — hydraulic conductivity estimate (m/s)

pycsamt.interp.petrophysics.transmissivity(K, thickness)[source]#

Aquifer transmissivity \(T = K \cdot b\) (m²/s).

Parameters:

K (array-like) – Hydraulic conductivity (m/s).
thickness (array-like) – Saturated thickness of the aquifer layer (m).

Returns:

T

Return type:

ndarray — transmissivity (m²/s)

pycsamt.interp.petrophysics.storativity(phi, thickness, *, specific_storage=0.0001)[source]#

Aquifer storativity \(S = S_s \cdot b\) (confined) or φ (unconfined).

For a confined aquifer: S = Ss × b (dimensionless). For an unconfined aquifer: S ≈ φ (specific yield).

Parameters:

phi (array-like) – Porosity (fraction, 0–1). Used as specific yield for unconfined.
thickness (array-like) – Saturated thickness (m).
specific_storage (float) – Specific storage Ss (m⁻¹) for confined conditions (default 10⁻⁴).

Returns:

S_confined (ndarray — confined storativity (dimensionless))
S_unconfined (ndarray — unconfined storativity ≈ specific yield (fraction))

Return type:

Any

pycsamt.interp.petrophysics.rho_w_to_tds(rho_w, *, temp_c=25.0)[source]#

Pore-water resistivity → Total Dissolved Solids (mg/L).

Uses the empirical Keller (1966) relation:

\[\mathrm{TDS} \approx \frac{640}{\mathrm{EC}_{(dS/m)}} = \frac{6400}{\sigma_{w(mS/cm)}}\]

With a temperature correction from the reference 25 °C:

\[\rho_{w,25} = \rho_{w,T} / [1 + 0.02(T - 25)]\]

Parameters:

rho_w (array-like) – Pore-water resistivity (Ω·m) at temperature temp_c.
temp_c (float) – Measurement temperature (°C; default 25).

Returns:

tds

Return type:

ndarray — TDS in mg/L (drinking limit: 500 mg/L; saline: > 3 000)

pycsamt.interp.petrophysics.tds_to_rho_w(tds, *, temp_c=25.0)[source]#

TDS (mg/L) → pore-water resistivity (Ω·m).

Inverse of rho_w_to_tds().

Parameters:

tds (array-like) – Total dissolved solids (mg/L).
temp_c (float) – Temperature (°C; default 25).

Returns:

rho_w

Return type:

ndarray — pore-water resistivity (Ω·m)

pycsamt.interp.petrophysics.ec_mscm_to_rho(ec, *, temp_c=25.0)[source]#

Electrical conductivity (mS/cm) → resistivity (Ω·m).

Parameters:

ec (Any)
temp_c (float)

Return type:

Any

pycsamt.interp.petrophysics.rho_to_ec_mscm(rho, *, temp_c=25.0)[source]#

Formation/water resistivity (Ω·m) → EC (mS/cm).

Parameters:

rho (Any)
temp_c (float)

Return type:

Any

pycsamt.interp.petrophysics.skin_depth(rho, freq)[source]#

EM skin depth (penetration depth) in metres.

\[\delta = 503 \sqrt{\rho / f}\]

This gives the 1/e amplitude depth for a plane EM wave, useful for estimating the sensitivity depth of AMT/MT measurements.

Parameters:

rho (array-like) – Apparent (or formation) resistivity (Ω·m).
freq (array-like) – Frequency (Hz).

Returns:

delta

Return type:

ndarray — skin depth (m)

pycsamt.interp.petrophysics.bostick_depth(rho_a, freq)[source]#

Bostick (1977) pseudo-depth for MT/AMT/TDEM data.

\[d_B = 503 \sqrt{\rho_a / f}\]

Identical numerically to skin depth but interpreted as the depth of maximum sensitivity for a given apparent resistivity and frequency. Use the full sequence {(ρ_a, f)} to build a pseudo-1D depth profile.

Parameters:

rho_a (array-like) – Apparent resistivity (Ω·m).
freq (array-like) – Frequency (Hz). For TDEM use equivalent frequency 1/(2πt).

Returns:

d

Return type:

ndarray — Bostick pseudo-depth (m)

pycsamt.interp.petrophysics.aquifer_top_from_profile(rho_log10, z_centers, *, rho_threshold_ohm_m=300.0, direction='low', min_depth=0.0)[source]#

Detect the top of an aquifer/conductor in a 1-D resistivity profile.

Scans the profile from the surface down and returns the depth of the first cell that crosses the threshold.

Parameters:

rho_log10 (array-like (n_z,)) – Log₁₀ resistivity profile (one station column).
z_centers (array-like (n_z,)) – Depth of each cell centre (m, positive downward).
rho_threshold_ohm_m (float) – Resistivity threshold. For aquifer detection (low resistivity): crossing downward indicates aquifer top.
direction ({'low', 'high'}) – 'low' — detect transition to low ρ (aquifer). 'high' — detect transition to high ρ (basement / caprock).
min_depth (float) – Ignore cells shallower than min_depth metres (default 0).

Returns:

depth

Return type:

float or None — top depth (m), or None if threshold not crossed.

pycsamt.interp.petrophysics.water_table_from_profile(rho_log10, z_centers, archie, *, rho_w=0.025, Sw_threshold=0.85, min_depth=0.5)[source]#

Estimate water table depth from a 1-D EM resistivity column.

Converts each cell’s log₁₀ρ to water saturation S_w via Archie’s inverse, then finds the shallowest depth where S_w ≥ Sw_threshold.

Parameters:

rho_log10 (array-like (n_z,)) – Log₁₀ resistivity column from EM inversion.
z_centers (array-like (n_z,)) – Depth of each cell (m).
archie (ArchieModel) – Calibrated Archie model.
rho_w (float) – Pore-water resistivity (Ω·m).
Sw_threshold (float) – Saturation value defining the water table (default 0.85).
min_depth (float) – Minimum search depth (m) to skip near-surface noise.

Returns:

depth

Return type:

float or None — estimated water table depth (m).

pycsamt.interp.petrophysics.fractured_zone_K(rho, *, rho_matrix=5000.0, aperture_m=0.001, cubic_factor=0.08333333333333333, kinematic_viscosity=1e-06, gravity=9.81)[source]#

Hydraulic conductivity of a fractured zone from EM resistivity.

Uses the parallel-plate cubic law for fracture flow:

\[K = \frac{g \, b^3}{12 \nu} \cdot f_\mathrm{vol}\]

where the volumetric fracture fraction f_vol is estimated from the resistivity contrast against the matrix:

\[f_\mathrm{vol} \approx 1 - \frac{\rho}{\rho_\mathrm{matrix}}\]

This is a rough first-pass estimate for fractured basement in AMT surveys. For rigorous fracture characterisation use the Bahr or distortion framework in pycsamt.emtools.

Parameters:

rho (array-like) – Measured formation resistivity (Ω·m).
rho_matrix (float) – Background (intact rock) resistivity (Ω·m).
aperture_m (float) – Representative fracture aperture (m).
cubic_factor (float) – Cubic-law pre-factor (1/12 for parallel plates).
kinematic_viscosity (float) – Kinematic viscosity of water at ~20 °C (m²/s, default 10⁻⁶).
gravity (float) – Gravitational acceleration (m/s²).

Returns:

K

Return type:

ndarray — effective hydraulic conductivity (m/s)