pycsamt.interp.uncertainty#

Monte Carlo uncertainty propagation for EM hydro-geophysical models.

Petrophysical parameters (ρ_w, Archie m, n, porosity prior) are never known exactly from the EM data alone. This module propagates that parametric uncertainty through the full chain

(ρ_w, m, n, φ) → EMHydroModel → K, Sw, φ, water-table depth, T, S

by drawing N random realisations of the free parameters from user-specified prior distributions, running EMHydroModel for each, and accumulating ensemble statistics (mean, std, P10, P50, P90).

EM method context#

  • TDEM — ρ_w uncertainty dominates (water quality unknown in advance). Typical range: 5–100 Ω·m for freshwater targets.

  • AMT — m uncertainty is significant in fractured basement (cementation exponent varies 1.3–2.5 in fractured vs porous media).

  • EMAP / MT — mixed uncertainty; wide ρ_w range for basin brines.

Typical use#

>>> from pycsamt.interp.uncertainty import (
...     UncertaintyBounds, MonteCarloHydro
... )
>>> bounds = UncertaintyBounds(
...     rho_w_range=(5.0, 80.0),      # fresh water range
...     m_range=(1.5, 2.2),
...     phi_prior_range=(0.15, 0.40),
... )
>>> mc = MonteCarloHydro(resistivity_model, config, bounds, n_samples=300)
>>> unc = mc.run()
>>> print(unc.p90_wt - unc.p10_wt)   # WT depth uncertainty (m)
>>> print(unc.cv_K)                    # coefficient of variation of K

References

Classes

MonteCarloHydro(resistivity_model, config, ...)

Monte Carlo uncertainty propagation for a quantitative EM hydro model.

UncertaintyBounds([rho_w_range, m_range, ...])

Prior distribution specification for Monte Carlo sampling.

UncertaintyResult(resistivity_model, config, ...)

Ensemble statistics from a MonteCarloHydro run.
class pycsamt.interp.uncertainty.UncertaintyBounds(rho_w_range=None, m_range=None, n_range=None, phi_prior_range=None, dist='uniform')[source]#

Bases: PyCSAMTObject

Prior distribution specification for Monte Carlo sampling.

Each *_range parameter activates that parameter as a free variable. Un-specified parameters are held fixed at the central PetrophysicalConfig values.

Parameters:

  • rho_w_range ((float, float), optional) – Pore-water resistivity range (Ω·m). uniform: (low, high); normal: (mean, std). Typical fresh-water range: (5, 100); brackish: (0.5, 20).

  • m_range ((float, float), optional) – Archie cementation exponent. Typical: (1.3, 2.5) for mixed lithology.

  • n_range ((float, float), optional) – Archie saturation exponent. Typical: (1.8, 2.5) for clean sand.

  • phi_prior_range ((float, float), optional) – Porosity prior. Typical: (0.10, 0.45).

  • dist (str) – Sampling distribution: 'uniform' (default, non-informative) or 'normal' (Gaussian — range is interpreted as (mean, std)).

Notes

At least one *_range must be set. rho_w_range alone covers the dominant source of uncertainty in most shallow EM surveys.

rho_w_range: tuple[float, float] | None = None#
m_range: tuple[float, float] | None = None#
n_range: tuple[float, float] | None = None#
phi_prior_range: tuple[float, float] | None = None#
dist: str = 'uniform'#
property n_free: int[source]#

Number of free (uncertain) parameters.

property free_names: list[str][source]#

Names of the free parameters.

sample(cfg, n, rng)[source]#

Draw n parameter samples and return a list of configs.

Parameters:

  • cfg (PetrophysicalConfig) – Central (best-estimate) configuration — used as fallback for fixed parameters.

  • n (int) – Number of Monte Carlo samples.

  • rng (np.random.Generator)

Return type:

list of PetrophysicalConfig, length n

class pycsamt.interp.uncertainty.UncertaintyResult(resistivity_model, config, bounds, n_samples, method_tag='', sampled_params=<factory>, mean_K=<factory>, std_K=<factory>, p10_K=<factory>, p50_K=<factory>, p90_K=<factory>, cv_K=<factory>, mean_Sw=<factory>, std_Sw=<factory>, p10_Sw=<factory>, p90_Sw=<factory>, mean_phi=<factory>, std_phi=<factory>, mean_wt=<factory>, std_wt=<factory>, p10_wt=<factory>, p50_wt=<factory>, p90_wt=<factory>, wt_detection_rate=<factory>, mean_T=<factory>, std_T=<factory>, p10_T=<factory>, p50_T=<factory>, p90_T=<factory>)[source]#

Bases: PyCSAMTObject

Ensemble statistics from a MonteCarloHydro run.

All 2-D arrays have shape (n_z, n_x); 1-D station arrays have shape (n_x,). K values are stored on the linear scale (m/s).

Variables:

  • resistivity_model (ResistivityModel)

  • config (PetrophysicalConfig) – Central configuration used to seed the Monte Carlo.

  • bounds (UncertaintyBounds)

  • n_samples (int)

  • method_tag (str)

  • sampled_params (ndarray (n_samples, n_free)) – Matrix of sampled parameter values, columns in bounds.free_names order.

  • (m/s) (Hydraulic conductivity K)

  • ------------------------------

  • p90_K (mean_K, std_K, p10_K, p50_K,)

  • cv_K (ndarray (n_z, n_x)) – Coefficient of variation std_K / mean_K. Dimensionless uncertainty map.

  • Sw (Water saturation)

  • -------------------

  • p90_Sw (mean_Sw, std_Sw, p10_Sw,)

  • φ (Porosity)

  • ----------

  • std_phi (mean_phi,)

  • (m) (Water-table depth)

  • ---------------------

  • p90_wt (mean_wt, std_wt, p10_wt, p50_wt,)

  • wt_detection_rate (ndarray (n_x,)) – Fraction of samples in which the water table was detected (0–1).

  • (m²/s) (Transmissivity T)

  • -----------------------

  • p90_T (mean_T, std_T, p10_T, p50_T,)

Parameters:
resistivity_model: ResistivityModel#
config: PetrophysicalConfig#
bounds: UncertaintyBounds#
n_samples: int#
method_tag: str = ''#
sampled_params: ndarray#
mean_K: ndarray#
std_K: ndarray#
p10_K: ndarray#
p50_K: ndarray#
p90_K: ndarray#
cv_K: ndarray#
mean_Sw: ndarray#
std_Sw: ndarray#
p10_Sw: ndarray#
p90_Sw: ndarray#
mean_phi: ndarray#
std_phi: ndarray#
mean_wt: ndarray#
std_wt: ndarray#
p10_wt: ndarray#
p50_wt: ndarray#
p90_wt: ndarray#
wt_detection_rate: ndarray#
mean_T: ndarray#
std_T: ndarray#
p10_T: ndarray#
p50_T: ndarray#
p90_T: ndarray#
prob_wt_shallower_than(depth_m, *, wt_ensemble=None)[source]#

P(water-table depth < depth_m) per station column.

Parameters:

  • depth_m (float) – Reference depth (m, positive downward).

  • wt_ensemble (ndarray (n_samples, n_x), optional) – Full water-table ensemble from MonteCarloHydro.run_ensemble(). If None, approximated from the stored P10/P50/P90 assuming a Gaussian distribution with mean=mean_wt, std=std_wt.

Return type:

ndarray (n_x,) — probability in [0, 1]

station_report()[source]#

Per-station uncertainty summary.

Return type:

list[dict]

to_csv(path)[source]#

Write per-station uncertainty summary to CSV.

Parameters:

path (str | Path)

Return type:

Path

class pycsamt.interp.uncertainty.MonteCarloHydro(resistivity_model, config, bounds, *, n_samples=200, seed=42, method_tag='', verbose=False)[source]#

Bases: PyCSAMTObject

Monte Carlo uncertainty propagation for a quantitative EM hydro model.

Parameters:

  • resistivity_model (ResistivityModel) – Source inversion model.

  • config (PetrophysicalConfig) – Central (best-estimate) petrophysical configuration.

  • bounds (UncertaintyBounds) – Prior distributions for the free parameters.

  • n_samples (int) – Number of Monte Carlo realisations (default 200). Typical: 100 for quick runs, 500–1000 for publication-quality.

  • seed (int) – Random seed for reproducibility (default 42).

  • method_tag (str, optional) – EM method label inherited by the output.

  • verbose (bool) – Print progress every 50 iterations (default False).

Examples

>>> bounds = UncertaintyBounds(rho_w_range=(5.0, 100.0), m_range=(1.4, 2.4))
>>> mc  = MonteCarloHydro(rm, cfg, bounds, n_samples=300)
>>> unc = mc.run()
>>> unc.cv_K.max()          # worst-case K uncertainty (fraction)
>>> unc.p90_wt - unc.p10_wt  # WT depth spread per station (m)
run()[source]#

Execute the Monte Carlo ensemble and return UncertaintyResult.

Return type:

UncertaintyResult

run_ensemble()[source]#

Like run() but also returns the raw WT and T ensembles.

Returns:

  • unc (UncertaintyResult)

  • wt_ensemble (ndarray (n_samples, n_x))

  • T_ensemble (ndarray (n_samples, n_x))

Return type:

tuple[UncertaintyResult, ndarray, ndarray]