petrophysics package#

The relationships that connect a hydrological state to a geophysical property. Each one is a stated model with its own assumptions, so the reference for the model you use is given alongside it.

Submodules#

PyHydroGeophysX.petrophysics.resistivity_models module#

Water content and saturation to electrical resistivity, and back.

References

Archie, G. E. (1942). The electrical resistivity log as an aid in determining some reservoir characteristics. Transactions of the AIME, 146(1), 54-62. The saturation and resistivity relation used by water_content_to_resistivity and its inverses.
Waxman, M. H. and Smits, L. J. M. (1968). Electrical conductivities in oil-bearing shaly sands. Society of Petroleum Engineers Journal, 8(2), 107-122. The surface-conduction term in WS_Model, which matters where clay contributes to the measured conductivity.

Simplified Waxman-Smits model for converting between water content and resistivity.

This implementation follows the Waxman-Smits model that expresses conductivity as:

σ = σsat * S^n + σs * S^(n-1)

where: - σ is the electrical conductivity of the formation - σsat is the conductivity at full saturation without surface effects (1/rhos) - σs is the surface conductivity - S is the water saturation (S = θ/φ where θ is water content and φ is porosity) - n is the saturation exponent

The resistivity is the reciprocal of conductivity: ρ = 1/σ

PyHydroGeophysX.petrophysics.resistivity_models.WS_Model(saturation: Any, porosity: Any, sigma_w: Any, m: Any, n: Any, sigma_s: Any = 0) → Any[source]#

Convert saturation to resistivity using the implemented Waxman-Smits form.

Uses σ = (S_w^n/F) * (σ_w + σ_s/S_w), where F = φ^(-m). Here σ_s is inside the formation-factor scaling; it is not interchangeable with the unscaled sigma_sur parameter used by the other converters.

Parameters:

saturation (array) – Saturation fraction, clipped to [0.001, 1].
porosity (array) – Porosity fraction (φ), positive and broadcast-compatible.
sigma_w (float) – Pore water conductivity (S/m).
m (float) – Cementation exponent
n (float) – Saturation exponent
sigma_s (float) – Surface-conductivity coefficient (S/m). Default is 0.

Returns:

Resistivity in ohm-m, clipped to [0.1, 1e6].

Return type:

array

Warns:

RuntimeWarning – Once per call if saturation, the conductivity floor (1e-6 S/m), or resistivity is clipped. Bounds and numerical outputs are unchanged. NaN values propagate; they are not repaired.

PyHydroGeophysX.petrophysics.resistivity_models.resistivity_to_porosity(resistivity: Any, saturation: Any, m: Any, rho_fluid: Any, n: Any, sigma_sur: Any = 0, a: Any = 1.0) → Any[source]#

Convert resistivity to porosity using Waxman-Smits model, given known saturation.

The function solves the Waxman-Smits equation for porosity: 1/rho = sigma_sat * S^n + sigma_sur * S^(n-1) where sigma_sat = 1/rhos and rhos = a * rho_fluid * porosity^(-m)

Rearranging: porosity = [(1/rho - sigma_sur * S^(n-1)) * (a * rho_fluid) / S^n]^(1/m)

Porosity only adds conduction, so where 1/rho <= sigma_sur * S**(n-1) the surface term alone conducts at least as well as the rock is measured to and no porosity satisfies the equation. Those values are NaN, with a RuntimeWarning naming the parameters and the resistivity range a porosity can explain.

Parameters:

resistivity (array) – Resistivity values (ohm-m)
saturation (array) – Saturation values (fraction, 0-1)
m (float) – Cementation exponent (typically 1.3-2.5)
rho_fluid (float) – Fluid resistivity (ohm-m)
n (float) – Saturation exponent (typically 1.8-2.2)
sigma_sur (float) – Surface conductivity (S/m). Default is 0 (no surface effects)
a (float) – Tortuosity factor. Default is 1.0

Returns:

Porosity values (fraction), clipped to [0.001, 0.99], shaped like the inputs broadcast together; NaN where no porosity exists. A float when resistivity and saturation are scalars.

Return type:

array

PyHydroGeophysX.petrophysics.resistivity_models.resistivity_to_saturation(resistivity: Any, porosity: Any, m: Any, rho_fluid: Any, n: Any, sigma_sur: Any = 0, a: Any = 1.0) → Any[source]#

Convert resistivity to saturation using Waxman-Smits model.

The function calculates saturated resistivity using Archie’s law: rhos = a * rho_fluid * porosity^(-m)

Then solves the Waxman-Smits equation: 1/rho = sigma_sat * S^n + sigma_sur * S^(n-1) where sigma_sat = 1/rhos

Parameters:

resistivity (array) – Resistivity values (ohm-m)
porosity (array) – Porosity values (fraction, 0-1)
m (float) – Cementation exponent (typically 1.3-2.5)
rho_fluid (float) – Fluid resistivity (ohm-m)
n (float) – Saturation exponent (typically 1.8-2.2)
sigma_sur (float) – Surface conductivity (S/m). Default is 0 (no surface effects)
a (float) – Tortuosity factor. Default is 1.0

Returns:

Saturation values (fraction, 0-1), shaped like the inputs broadcast together (NumPy rules, so a (cells, times) resistivity takes a scalar, (cells, 1) or (cells, times) porosity); a float for one value.

Return type:

array

Raises:

ValueError – sigma_sur is negative.

PyHydroGeophysX.petrophysics.resistivity_models.resistivity_to_saturation2(resistivity: Any, rhos: Any, n: Any, sigma_sur: Any = 0) → Any[source]#

Convert resistivity to saturation using Waxman-Smits model.

Surface conduction is solved with a safeguarded root finder. The zero-surface case uses Archie’s analytic solution without a dry-end floor. Saturations outside the attainable range for n >= 1 are clipped to [0, 1]. For 0 < n < 1, the larger physical root is preferred; if no physical root exists, the result is NaN with a RuntimeWarning.

Parameters:

resistivity (array) – Resistivity values
rhos (float) – Saturated resistivity without surface effects
n (float) – Saturation exponent
sigma_sur (float) – Surface conductivity. Default is 0 (no surface effects).

Returns:

Saturation values

Return type:

array

PyHydroGeophysX.petrophysics.resistivity_models.resistivity_to_water_content(resistivity: Any, rhos: Any, n: Any, porosity: Any, sigma_sur: Any = 0) → Any[source]#

Convert resistivity to water content using Waxman-Smits model.

Parameters:

resistivity (array) – Resistivity values
rhos (float) – Saturated resistivity without surface effects
n (float) – Saturation exponent
porosity (array) – Porosity values
sigma_sur (float) – Surface conductivity. Default is 0 (no surface effects).

Returns:

Volumetric water content values

Return type:

array

PyHydroGeophysX.petrophysics.resistivity_models.water_content_to_resistivity(water_content: Any, rhos: Any, n: Any, porosity: Any, sigma_sur: Any = 0) → Any[source]#

Convert water content to resistivity using Waxman-Smits model.

Parameters:

water_content (array) – Volumetric water content (θ)
rhos (float) – Saturated resistivity without surface effects
n (float) – Saturation exponent
porosity (array) – Porosity values (φ)
sigma_sur (float) – Surface conductivity. Default is 0 (no surface effects).

Returns:

Resistivity values

Return type:

array

PyHydroGeophysX.petrophysics.velocity_models module#

Porosity, saturation and mineralogy to seismic velocity.

References

Mindlin, R. D. (1949). Compliance of elastic bodies in contact. Journal of Applied Mechanics, 16, 259-268. The grain-contact theory behind HertzMindlinModel.
Gassmann, F. (1951). Über die Elastizität poröser Medien. Vierteljahrsschrift der Naturforschenden Gesellschaft in Zürich, 96, 1-23. Fluid substitution, applied in satK.
Hashin, Z. and Shtrikman, S. (1963). A variational approach to the theory of the elastic behaviour of multiphase materials. Journal of the Mechanics and Physics of Solids, 11(2), 127-140. The bounds used by the effective-medium models, alongside the Voigt-Reuss-Hill average in VRHModel.
Brie, A., Pampuri, F., Marsala, A. F., and Meazza, O. (1995). Shear sonic interpretation in gas-bearing sands. SPE Annual Technical Conference, SPE 30595. The fluid-mixing law in BrieModel.

PetrophysicalCoupling.compare_inversions_to_hydro uses the same velocity_model as its forward helper. Select "linear", "wyllie" or "raymer" directly, or use velocity_model="empirical" with empirical_model. DEM and Hertz-Mindlin, including the forward helper’s default, require a calibrated velocity_inverse(velocity, porosity) callable in petro_params. They no longer silently use a linear inverse. Empirical velocities outside the physical saturation range or at the forward clipping limits raise ValueError instead of producing misleading clipped estimates. Without a usable velocity inverse, or when no cell converts, the SRT entry is {"error": reason} and a warning is issued; the ERT and EM comparisons are still returned. Cells whose velocity lies outside the model’s range are NaN and counted in stats['n_masked'].

Seismic velocity models for relating rock properties to elastic wave velocities.

class PyHydroGeophysX.petrophysics.velocity_models.BaseVelocityModel[source]#

Bases: object

Base class for seismic velocity models.

calculate_velocity(**kwargs) → ndarray[source]#

Calculate seismic velocity from rock properties.

Parameters:: **kwargs – Rock properties specific to each model
Returns:: Seismic velocity values (Vp, Vs, or both)

class PyHydroGeophysX.petrophysics.velocity_models.BrieModel(exponent: float = 3.0)[source]#

Bases: object

Brie’s model for calculating the effective bulk modulus of a partially saturated medium.

calculate_fluid_modulus(saturation: float, water_modulus: float = 2.0, gas_modulus: float = 0.01) → float[source]#

Calculate effective fluid bulk modulus using Brie’s equation.

Parameters:

saturation – Water saturation (0 to 1)
water_modulus – Bulk modulus of water (GPa, default: 2.0)
gas_modulus – Bulk modulus of gas (GPa, default: 0.01)

Returns:

Effective fluid bulk modulus (GPa)

calculate_saturated_modulus(dry_modulus: float, mineral_modulus: float, porosity: float, saturation: float, water_modulus: float = 2.0, gas_modulus: float = 0.01) → float[source]#

Calculate the saturated bulk modulus based on Brie’s equation.

Parameters:

dry_modulus – Bulk modulus of the dry rock (GPa)
mineral_modulus – Bulk modulus of the mineral matrix (GPa)
porosity – Porosity of the rock
saturation – Water saturation (0 to 1)
water_modulus – Bulk modulus of water (GPa, default: 2.0)
gas_modulus – Bulk modulus of gas (GPa, default: 0.01)

Returns:

Saturated bulk modulus (GPa)

class PyHydroGeophysX.petrophysics.velocity_models.DEMModel[source]#

Bases: BaseVelocityModel

Differential Effective Medium (DEM) model for calculating elastic properties and seismic velocities of porous rocks.

calculate_velocity(porosity: ndarray, saturation: ndarray, bulk_modulus: float, shear_modulus: float, mineral_density: float, aspect_ratio: float = 0.1) → Tuple[ndarray, ndarray, ndarray][source]#

Calculate P-wave velocity using the DEM model.

Parameters:

porosity – Porosity values (array)
saturation – Saturation values (array)
bulk_modulus – Initial bulk modulus of the solid matrix (GPa)
shear_modulus – Initial shear modulus of the solid matrix (GPa)
mineral_density – Density of the solid matrix (kg/m³)
aspect_ratio – Aspect ratio of pores (default: 0.1)

Returns:

Effective bulk modulus (GPa), effective shear modulus (GPa), and P-wave velocity (m/s)

class PyHydroGeophysX.petrophysics.velocity_models.HertzMindlinModel(critical_porosity: float = 0.4, coordination_number: float = 4.0)[source]#

Bases: BaseVelocityModel

Hertz-Mindlin model and Hashin-Shtrikman bounds for seismic velocity in porous rocks.

calculate_velocity(porosity: ndarray, saturation: ndarray, bulk_modulus: float, shear_modulus: float, mineral_density: float, depth: float = 1.0, *, effective_pressure=None, fluid_modulus=None, fluid_density=None, return_moduli: bool = False)[source]#

Calculate high/low P-wave bounds, optionally with elastic moduli.

Parameters:

porosity (array-like) – Fractions in [0, 1]. All numerical inputs broadcast together.
saturation (array-like) – Fractions in [0, 1]. All numerical inputs broadcast together.
bulk_modulus (array-like) – Positive mineral moduli in GPa.
shear_modulus (array-like) – Positive mineral moduli in GPa.
mineral_density (array-like) – Positive mineral density in kg/m3.
depth (array-like) – Depth in metres, default 1. Ignored when effective_pressure is set.
effective_pressure (array-like, optional) – Nonnegative grain-contact pressure in Pa, replacing the default (mineral_density - 1000) * 9.8 * depth estimate. This is distinct from the pore-fluid pressure passed to batzle_wang_brine.
fluid_modulus (array-like, optional) – Liquid bulk modulus in GPa, replacing water’s default 2.0 in Brie’s liquid/gas mixture (gas=0.01 GPa, exponent=3). An explicit value <= 1e-12 GPa selects the dry-frame limit at any saturation, bypassing fluid mixing and substitution.
fluid_density (array-like, optional) – Liquid density in kg/m3, default 1000; used in bulk density only. The default pressure estimate does not depend on this parameter.
return_moduli (bool) – False preserves (Vp_high, Vp_low). True adds a dictionary with Vs_low, Vs_high, K_dry_low, K_dry_high, G_dry_low, G_dry_high, K_sat_low, K_sat_high and density. Moduli are GPa, velocities m/s, density kg/m3. Output arrays have the broadcast input shape.

Notes

Coordination number remains 4. A value around 9 is commonly used for idealized random packings, but should be calibrated to the rock (Mavko, Mukerji and Dvorkin, 2009, The Rock Physics Handbook, ch. 5). Zero effective pressure is handled by the zero-contact-stiffness limit. No change is made to the caller’s saturation or porosity.

class PyHydroGeophysX.petrophysics.velocity_models.VRHModel[source]#

Bases: BaseVelocityModel

Voigt-Reuss-Hill (VRH) mixing model for effective elastic properties of composites.

calculate_properties(fractions: List[float], bulk_moduli: List[float], shear_moduli: List[float], densities: List[float]) → Tuple[float, float, float][source]#

Calculate effective elastic properties using the VRH model.

Parameters:

fractions – Volume fractions of each mineral (must sum to 1)
bulk_moduli – Bulk moduli of each mineral (GPa)
shear_moduli – Shear moduli of each mineral (GPa)
densities – Densities of each mineral (kg/m³)

Returns:

Effective bulk modulus (GPa), effective shear modulus (GPa), and effective density (kg/m³)

calculate_velocity(fractions: List[float], bulk_moduli: List[float], shear_moduli: List[float], densities: List[float]) → Tuple[float, float][source]#

Calculate P-wave and S-wave velocities using the VRH model.

Parameters:

fractions – Volume fractions of each mineral (must sum to 1)
bulk_moduli – Bulk moduli of each mineral (GPa)
shear_moduli – Shear moduli of each mineral (GPa)
densities – Densities of each mineral (kg/m³)

Returns:

P-wave velocity (m/s) and S-wave velocity (m/s)

PyHydroGeophysX.petrophysics.velocity_models.VRH_model(f: Any = [0.35, 0.25, 0.2, 0.125, 0.075], K: Any = [55.4, 36.6, 75.6, 46.7, 50.4], G: Any = [28.1, 45, 25.6, 23.65, 27.4], rho: Any = [2560, 2650, 2630, 2540, 3050]) → Any[source]#

Implements the Voigt-Reuss-Hill (VRH) mixing model to estimate the effective bulk modulus (Km), shear modulus (Gm), and density (rho_b) of a composite material made from various minerals.

Parameters: f (list): Fraction of each mineral in the composite (must sum to 1). K (list): Bulk modulus of each mineral (GPa). G (list): Shear modulus of each mineral (GPa). rho (list): Density of each mineral (kg/m^3).

Returns: Km (float): Effective bulk modulus of the composite material (GPa). Gm (float): Effective shear modulus of the composite material (GPa). rho_b (float): Effective density of the composite material (kg/m^3).

PyHydroGeophysX.petrophysics.velocity_models.satK(Keff: Any, Km: Any, phi: Any, Sat: Any) → Any[source]#

Calculate the saturated bulk modulus (K_sat) based on Brie’s equation.

Parameters: Keff (float): Effective bulk modulus of the dry rock (GPa). Km (float): Bulk modulus of the matrix (GPa). phi (float): Porosity of the rock. Sat (float): Saturation level of the fluid in the pores.

Returns: float: Saturated bulk modulus (GPa).

PyHydroGeophysX.petrophysics.velocity_models.velDEM(phi: Any, Km: Any, Gm: Any, rho_b: Any, Sat: Any, alpha: Any) → Any[source]#

Calculate effective bulk modulus (Keff), shear modulus (Geff), and P-wave velocity (Vp) for a rock with varying porosity (phi) based on the DEM model, taking into account the saturation (Sat) and the crack aspect ratio (alpha).

Parameters: phi (np.array): Array of porosities. Km (float): Initial bulk modulus of the material (GPa). Gm (float): Initial shear modulus of the material (GPa). rho_b (float): Density of the solid phase (kg/m^3). Sat (array): Saturation fractions with one entry per phi value; scalars are not broadcast here. alpha (float): Crack aspect ratio.

Returns: Keff1 (np.array): Effective bulk modulus for each porosity value (GPa). Geff1 (np.array): Effective shear modulus for each porosity value (GPa). Vp (np.array): P-wave velocity for each porosity value (m/s).

PyHydroGeophysX.petrophysics.velocity_models.vel_porous(phi: Any, Km: Any, Gm: Any, rho_b: Any, Sat: Any, depth: Any = 1) → Any[source]#

Calculate P-wave velocity (Vp) for a rock with varying porosity (phi) based on the Hertz-Mindlin model and Hashin-Shtrikman bounds, taking into account the saturation (Sat).

Compatibility wrapper over HertzMindlinModel.calculate_velocity with critical porosity 0.4 and coordination number 4. Both return modes of that method and this function share the contact, HS and Gassmann implementation.

Parameters: phi (np.array): Array of porosities. Km (float): Bulk modulus of the solid phase (GPa). Gm (float): Shear modulus of the solid phase (GPa). rho_b (float): Density of the solid phase (kg/m^3). Sat (array): Saturation fractions; scalars and arrays broadcast with phi. depth (float): depth for pressure estimation (m)

Returns: Vp_h (np.array): P-wave velocity for each porosity value (upper Hashin-Shtrikman bound) (m/s). Vp_l (np.array): P-wave velocity for each porosity value (lower Hashin-Shtrikman bound) (m/s).

PyHydroGeophysX.petrophysics.velocity_models.velocity_to_water_content(velocity: ndarray, v_dry: float = 3500.0, v_sat: float = 4500.0, porosity: float | ndarray = 0.3, model: str = 'linear') → ndarray[source]#

Invert the empirical relationship used by water_content_to_velocity.

Supports linear, Wyllie and Raymer models with the same fixed constituent velocities as the forward function. Porosity must be known. Observations outside the model’s saturation range, or at the forward function’s clipping limits (200/8000 m/s), are rejected rather than silently assigned 0/1 saturation. Hertz-Mindlin and DEM require a separately calibrated inverse.

PyHydroGeophysX.petrophysics.velocity_models.water_content_to_velocity(water_content: ndarray, v_dry: float = 3500.0, v_sat: float = 4500.0, porosity: float | ndarray = 0.3, model: str = 'linear') → ndarray[source]#

Convert water content to seismic P-wave velocity using empirical relationships.

These are simplified mixing rules, not calibrated site predictions. Saturation is clipped to [0, 1] and output velocity to [200, 8000] m/s.

Parameters:

water_content – Volumetric water content (θ), typically 0-porosity
v_dry – Dry velocity for the linear model only (m/s), default 3500.
v_sat – Saturated velocity for the linear model only (m/s), default 4500.
porosity – Porosity values (φ), can be scalar or array
model – Velocity model type: - ‘linear’: Linear interpolation between dry and saturated velocities - ‘wyllie’: Wyllie time-average equation - ‘raymer’: Raymer-Hunt-Gardner equation

Returns:

P-wave velocity array (m/s)

Pore-fluid properties#

For units, worked examples and calibration limits, see Petrophysics and pore-fluid properties.

Pore-fluid properties with explicit SI inputs and empirical-model units.

class PyHydroGeophysX.petrophysics.fluid_properties.BrineProperties(bulk_modulus: ndarray, density: ndarray)[source]#

Bases: NamedTuple

Bulk modulus in GPa and density in kg/m3, with broadcast input shape.

bulk_modulus: ndarray#: Alias for field number 0

density: ndarray#: Alias for field number 1

PyHydroGeophysX.petrophysics.fluid_properties.batzle_wang_brine(p, T, salinity)[source]#

Return NaCl-brine bulk modulus (GPa) and density (kg/m3).

p is pore-fluid pressure in Pa (not effective grain-contact pressure); T is temperature in deg C; salinity is NaCl mass fraction in kg/kg of solution. Scalars and arrays broadcast together. The empirical liquid model is restricted here to 0–100 MPa, 0–100 deg C and 0–0.35 kg/kg; the caller must ensure the fluid is liquid and salt remains dissolved.

Implements the density and acoustic-velocity polynomials of Batzle and Wang (1992), Geophysics 57, 1396–1408, doi:10.1190/1.1443207. The published polynomial uses MPa and g/cm3 internally.

PyHydroGeophysX.petrophysics.fluid_properties.bulk_cec_to_qv(cec, porosity)[source]#

Convert CEC in eq/m3 bulk to Qv in meq/cm3 pore.

Qv = CEC / porosity * 1e-3. This accepts bulk-volume CEC, not mass-specific laboratory CEC (e.g. meq/100 g). Inputs broadcast. Porosity must be in (0, 1]; CEC must be finite and nonnegative.

PyHydroGeophysX.petrophysics.fluid_properties.ion_charge(name)[source]#

Return signed charge from repeated trailing signs (Ca++=2, CO3–=-2).

Only the repeated-sign convention is accepted. Formula digits belong to the species (NH4+ has charge +1); numeric charge notation is not parsed. Missing or mixed trailing signs raise ValueError.

PyHydroGeophysX.petrophysics.fluid_properties.pore_fluid_conductivity(conc, equivalent_conductivities=None, return_contributions=True, *, ionic_strength_correction=False, kohlrausch_coefficient=None)[source]#

Calculate pore-water conductivity at 25 deg C in S/m.

Parameters:

conc (mapping) – Ion names to dissolved-species concentrations in mol/L. Values may be scalars or broadcast-compatible arrays; supply speciated ions, not total analytical concentrations or neutral salts.
equivalent_conductivities (mapping, optional) – Overrides/additions to the limiting table, in S m2/eq. Do not supply molar conductivities: sigma_i = abs(z_i) * lambda_i * c_i * 1000. The default table is Kanciruk et al. (1986), ORNL/TM-10153, Table 8, CONCAL (p. 30). Unknown ions require an explicit coefficient.
return_contributions (bool) – If True return (sigma_w, {ion: sigma_i}), otherwise only sigma_w. An empty composition returns zero (no implicit water autoionization).
ionic_strength_correction (bool) – Off by default. Apply a calibrated Kohlrausch-type factor 1 - a * sqrt(I), where I = 0.5 * sum(c_i * z_i**2) in mol/L, to every contribution. Intended only for dilute strong electrolytes of the composition used to calibrate a; no speciation is performed.
kohlrausch_coefficient (float, optional) – Required when correction is enabled, in (mol/L)**(-1/2). There is no universal default. For a fixed electrolyte, Kohlrausch’s Lambda = Lambda0 - A * sqrt(c_salt) gives a = A / (Lambda0 * sqrt(I/c_salt)). See Metrohm, Conductometry - Conductivity Measurement, monograph 8.109.5021, p. 14. A negative resulting factor raises ValueError rather than clipping.

Notes

Correction is an empirical fixed-composition approximation, not a general mixed-electrolyte transport model. Contributions under this approximation retain their limiting relative proportions.

PyHydroGeophysX.petrophysics.fluid_properties.waxman_smits_surface_conductivity(sigma_w, qv)[source]#

Return the Waxman-Thomas surface coefficient in S/m at 25 deg C.

sigma_s = 4.6 * (1 - 0.6 * exp(-sigma_w / 1.3)) * qv, with sigma_w in S/m and qv in meq/cm3 pore. Thus 1.3 has units S/m and 4.6 has units (S/m)/(meq/cm3). Inputs broadcast.

Pass the result as sigma_s to WS_Model, which applies the formation factor. It is NOT the unscaled sigma_sur in the other converters. Reference: Waxman and Thomas (1974), doi:10.2118/4094-PA.

Module contents#

Petrophysical models for resistivity and seismic velocity relationships.

class PyHydroGeophysX.petrophysics.BaseVelocityModel[source]#

Bases: object

Base class for seismic velocity models.

calculate_velocity(**kwargs) → ndarray[source]#

Calculate seismic velocity from rock properties.

Parameters:: **kwargs – Rock properties specific to each model
Returns:: Seismic velocity values (Vp, Vs, or both)

class PyHydroGeophysX.petrophysics.BrieModel(exponent: float = 3.0)[source]#

Bases: object

Brie’s model for calculating the effective bulk modulus of a partially saturated medium.

calculate_fluid_modulus(saturation: float, water_modulus: float = 2.0, gas_modulus: float = 0.01) → float[source]#

Calculate effective fluid bulk modulus using Brie’s equation.

Parameters:

saturation – Water saturation (0 to 1)
water_modulus – Bulk modulus of water (GPa, default: 2.0)
gas_modulus – Bulk modulus of gas (GPa, default: 0.01)

Returns:

Effective fluid bulk modulus (GPa)

calculate_saturated_modulus(dry_modulus: float, mineral_modulus: float, porosity: float, saturation: float, water_modulus: float = 2.0, gas_modulus: float = 0.01) → float[source]#

Calculate the saturated bulk modulus based on Brie’s equation.

Parameters:

dry_modulus – Bulk modulus of the dry rock (GPa)
mineral_modulus – Bulk modulus of the mineral matrix (GPa)
porosity – Porosity of the rock
saturation – Water saturation (0 to 1)
water_modulus – Bulk modulus of water (GPa, default: 2.0)
gas_modulus – Bulk modulus of gas (GPa, default: 0.01)

Returns:

Saturated bulk modulus (GPa)

class PyHydroGeophysX.petrophysics.DEMModel[source]#

Bases: BaseVelocityModel

Differential Effective Medium (DEM) model for calculating elastic properties and seismic velocities of porous rocks.

calculate_velocity(porosity: ndarray, saturation: ndarray, bulk_modulus: float, shear_modulus: float, mineral_density: float, aspect_ratio: float = 0.1) → Tuple[ndarray, ndarray, ndarray][source]#

Calculate P-wave velocity using the DEM model.

Parameters:

porosity – Porosity values (array)
saturation – Saturation values (array)
bulk_modulus – Initial bulk modulus of the solid matrix (GPa)
shear_modulus – Initial shear modulus of the solid matrix (GPa)
mineral_density – Density of the solid matrix (kg/m³)
aspect_ratio – Aspect ratio of pores (default: 0.1)

Returns:

Effective bulk modulus (GPa), effective shear modulus (GPa), and P-wave velocity (m/s)

class PyHydroGeophysX.petrophysics.HertzMindlinModel(critical_porosity: float = 0.4, coordination_number: float = 4.0)[source]#

Bases: BaseVelocityModel

Hertz-Mindlin model and Hashin-Shtrikman bounds for seismic velocity in porous rocks.

calculate_velocity(porosity: ndarray, saturation: ndarray, bulk_modulus: float, shear_modulus: float, mineral_density: float, depth: float = 1.0, *, effective_pressure=None, fluid_modulus=None, fluid_density=None, return_moduli: bool = False)[source]#

Calculate high/low P-wave bounds, optionally with elastic moduli.

Parameters:

porosity (array-like) – Fractions in [0, 1]. All numerical inputs broadcast together.
saturation (array-like) – Fractions in [0, 1]. All numerical inputs broadcast together.
bulk_modulus (array-like) – Positive mineral moduli in GPa.
shear_modulus (array-like) – Positive mineral moduli in GPa.
mineral_density (array-like) – Positive mineral density in kg/m3.
depth (array-like) – Depth in metres, default 1. Ignored when effective_pressure is set.
effective_pressure (array-like, optional) – Nonnegative grain-contact pressure in Pa, replacing the default (mineral_density - 1000) * 9.8 * depth estimate. This is distinct from the pore-fluid pressure passed to batzle_wang_brine.
fluid_modulus (array-like, optional) – Liquid bulk modulus in GPa, replacing water’s default 2.0 in Brie’s liquid/gas mixture (gas=0.01 GPa, exponent=3). An explicit value <= 1e-12 GPa selects the dry-frame limit at any saturation, bypassing fluid mixing and substitution.
fluid_density (array-like, optional) – Liquid density in kg/m3, default 1000; used in bulk density only. The default pressure estimate does not depend on this parameter.
return_moduli (bool) – False preserves (Vp_high, Vp_low). True adds a dictionary with Vs_low, Vs_high, K_dry_low, K_dry_high, G_dry_low, G_dry_high, K_sat_low, K_sat_high and density. Moduli are GPa, velocities m/s, density kg/m3. Output arrays have the broadcast input shape.

Notes

Coordination number remains 4. A value around 9 is commonly used for idealized random packings, but should be calibrated to the rock (Mavko, Mukerji and Dvorkin, 2009, The Rock Physics Handbook, ch. 5). Zero effective pressure is handled by the zero-contact-stiffness limit. No change is made to the caller’s saturation or porosity.

class PyHydroGeophysX.petrophysics.VRHModel[source]#

Bases: BaseVelocityModel

Voigt-Reuss-Hill (VRH) mixing model for effective elastic properties of composites.

calculate_properties(fractions: List[float], bulk_moduli: List[float], shear_moduli: List[float], densities: List[float]) → Tuple[float, float, float][source]#

Calculate effective elastic properties using the VRH model.

Parameters:

fractions – Volume fractions of each mineral (must sum to 1)
bulk_moduli – Bulk moduli of each mineral (GPa)
shear_moduli – Shear moduli of each mineral (GPa)
densities – Densities of each mineral (kg/m³)

Returns:

Effective bulk modulus (GPa), effective shear modulus (GPa), and effective density (kg/m³)

calculate_velocity(fractions: List[float], bulk_moduli: List[float], shear_moduli: List[float], densities: List[float]) → Tuple[float, float][source]#

Calculate P-wave and S-wave velocities using the VRH model.

Parameters:

fractions – Volume fractions of each mineral (must sum to 1)
bulk_moduli – Bulk moduli of each mineral (GPa)
shear_moduli – Shear moduli of each mineral (GPa)
densities – Densities of each mineral (kg/m³)

Returns:

P-wave velocity (m/s) and S-wave velocity (m/s)

PyHydroGeophysX.petrophysics.VRH_model(f: Any = [0.35, 0.25, 0.2, 0.125, 0.075], K: Any = [55.4, 36.6, 75.6, 46.7, 50.4], G: Any = [28.1, 45, 25.6, 23.65, 27.4], rho: Any = [2560, 2650, 2630, 2540, 3050]) → Any[source]#

Implements the Voigt-Reuss-Hill (VRH) mixing model to estimate the effective bulk modulus (Km), shear modulus (Gm), and density (rho_b) of a composite material made from various minerals.

Parameters: f (list): Fraction of each mineral in the composite (must sum to 1). K (list): Bulk modulus of each mineral (GPa). G (list): Shear modulus of each mineral (GPa). rho (list): Density of each mineral (kg/m^3).

Returns: Km (float): Effective bulk modulus of the composite material (GPa). Gm (float): Effective shear modulus of the composite material (GPa). rho_b (float): Effective density of the composite material (kg/m^3).

PyHydroGeophysX.petrophysics.WS_Model(saturation: Any, porosity: Any, sigma_w: Any, m: Any, n: Any, sigma_s: Any = 0) → Any[source]#

Convert saturation to resistivity using the implemented Waxman-Smits form.

Uses σ = (S_w^n/F) * (σ_w + σ_s/S_w), where F = φ^(-m). Here σ_s is inside the formation-factor scaling; it is not interchangeable with the unscaled sigma_sur parameter used by the other converters.

Parameters:

saturation (array) – Saturation fraction, clipped to [0.001, 1].
porosity (array) – Porosity fraction (φ), positive and broadcast-compatible.
sigma_w (float) – Pore water conductivity (S/m).
m (float) – Cementation exponent
n (float) – Saturation exponent
sigma_s (float) – Surface-conductivity coefficient (S/m). Default is 0.

Returns:

Resistivity in ohm-m, clipped to [0.1, 1e6].

Return type:

array

Warns:

RuntimeWarning – Once per call if saturation, the conductivity floor (1e-6 S/m), or resistivity is clipped. Bounds and numerical outputs are unchanged. NaN values propagate; they are not repaired.

PyHydroGeophysX.petrophysics.batzle_wang_brine(p, T, salinity)[source]#

Return NaCl-brine bulk modulus (GPa) and density (kg/m3).

p is pore-fluid pressure in Pa (not effective grain-contact pressure); T is temperature in deg C; salinity is NaCl mass fraction in kg/kg of solution. Scalars and arrays broadcast together. The empirical liquid model is restricted here to 0–100 MPa, 0–100 deg C and 0–0.35 kg/kg; the caller must ensure the fluid is liquid and salt remains dissolved.

Implements the density and acoustic-velocity polynomials of Batzle and Wang (1992), Geophysics 57, 1396–1408, doi:10.1190/1.1443207. The published polynomial uses MPa and g/cm3 internally.

PyHydroGeophysX.petrophysics.bulk_cec_to_qv(cec, porosity)[source]#

Convert CEC in eq/m3 bulk to Qv in meq/cm3 pore.

Qv = CEC / porosity * 1e-3. This accepts bulk-volume CEC, not mass-specific laboratory CEC (e.g. meq/100 g). Inputs broadcast. Porosity must be in (0, 1]; CEC must be finite and nonnegative.

PyHydroGeophysX.petrophysics.ion_charge(name)[source]#

Return signed charge from repeated trailing signs (Ca++=2, CO3–=-2).

Only the repeated-sign convention is accepted. Formula digits belong to the species (NH4+ has charge +1); numeric charge notation is not parsed. Missing or mixed trailing signs raise ValueError.

PyHydroGeophysX.petrophysics.pore_fluid_conductivity(conc, equivalent_conductivities=None, return_contributions=True, *, ionic_strength_correction=False, kohlrausch_coefficient=None)[source]#

Calculate pore-water conductivity at 25 deg C in S/m.

Parameters:

conc (mapping) – Ion names to dissolved-species concentrations in mol/L. Values may be scalars or broadcast-compatible arrays; supply speciated ions, not total analytical concentrations or neutral salts.
equivalent_conductivities (mapping, optional) – Overrides/additions to the limiting table, in S m2/eq. Do not supply molar conductivities: sigma_i = abs(z_i) * lambda_i * c_i * 1000. The default table is Kanciruk et al. (1986), ORNL/TM-10153, Table 8, CONCAL (p. 30). Unknown ions require an explicit coefficient.
return_contributions (bool) – If True return (sigma_w, {ion: sigma_i}), otherwise only sigma_w. An empty composition returns zero (no implicit water autoionization).
ionic_strength_correction (bool) – Off by default. Apply a calibrated Kohlrausch-type factor 1 - a * sqrt(I), where I = 0.5 * sum(c_i * z_i**2) in mol/L, to every contribution. Intended only for dilute strong electrolytes of the composition used to calibrate a; no speciation is performed.
kohlrausch_coefficient (float, optional) – Required when correction is enabled, in (mol/L)**(-1/2). There is no universal default. For a fixed electrolyte, Kohlrausch’s Lambda = Lambda0 - A * sqrt(c_salt) gives a = A / (Lambda0 * sqrt(I/c_salt)). See Metrohm, Conductometry - Conductivity Measurement, monograph 8.109.5021, p. 14. A negative resulting factor raises ValueError rather than clipping.

Notes

Correction is an empirical fixed-composition approximation, not a general mixed-electrolyte transport model. Contributions under this approximation retain their limiting relative proportions.

PyHydroGeophysX.petrophysics.resistivity_to_saturation(resistivity: Any, porosity: Any, m: Any, rho_fluid: Any, n: Any, sigma_sur: Any = 0, a: Any = 1.0) → Any[source]#

Convert resistivity to saturation using Waxman-Smits model.

The function calculates saturated resistivity using Archie’s law: rhos = a * rho_fluid * porosity^(-m)

Then solves the Waxman-Smits equation: 1/rho = sigma_sat * S^n + sigma_sur * S^(n-1) where sigma_sat = 1/rhos

Parameters:

resistivity (array) – Resistivity values (ohm-m)
porosity (array) – Porosity values (fraction, 0-1)
m (float) – Cementation exponent (typically 1.3-2.5)
rho_fluid (float) – Fluid resistivity (ohm-m)
n (float) – Saturation exponent (typically 1.8-2.2)
sigma_sur (float) – Surface conductivity (S/m). Default is 0 (no surface effects)
a (float) – Tortuosity factor. Default is 1.0

Returns:

Saturation values (fraction, 0-1), shaped like the inputs broadcast together (NumPy rules, so a (cells, times) resistivity takes a scalar, (cells, 1) or (cells, times) porosity); a float for one value.

Return type:

array

Raises:

ValueError – sigma_sur is negative.

PyHydroGeophysX.petrophysics.resistivity_to_water_content(resistivity: Any, rhos: Any, n: Any, porosity: Any, sigma_sur: Any = 0) → Any[source]#

Convert resistivity to water content using Waxman-Smits model.

Parameters:

resistivity (array) – Resistivity values
rhos (float) – Saturated resistivity without surface effects
n (float) – Saturation exponent
porosity (array) – Porosity values
sigma_sur (float) – Surface conductivity. Default is 0 (no surface effects).

Returns:

Volumetric water content values

Return type:

array

PyHydroGeophysX.petrophysics.run_petrophysics_monte_carlo(resistivity: ~numpy.ndarray, markers: ~numpy.ndarray, layers: ~typing.Sequence[~typing.Mapping[str, ~typing.Any]], *, products: ~typing.Sequence[str] = ('water_content',), n_realizations: int = 100, seed: int = 7, saturation_value: float = 1.0, tortuosity_a: float = 1.0, timestep_indices: ~typing.Sequence[int] | None = None, progress: ~typing.Callable[[str], None] = <function noop>, return_realizations: bool = False, cell_chunk_size: int = 1024, log_resistivity_std: ~numpy.ndarray | None = None) → Dict[str, Any][source]#

Run reproducible layer-wise petrophysical uncertainty propagation.

Parameters:

resistivity (array_like) – Positive, finite resistivity in ohm-m, shaped (cells, times) or (cells,). A (times, cells) input is transposed only when its first axis cannot match markers. Square arrays are always interpreted as (cells, times).
markers (array_like) – Integer cell labels. Every label must have exactly one layer definition.
layers (sequence of mappings) – Each mapping requires marker. Parameters accept a scalar or a normal distribution specified by mean and nonnegative std. Resistivities use ohm-m, sigma_sur uses S/m, porosity is a fraction. Positive parameters are floored at 1e-6, sigma_sur at zero, and sampled porosity is clipped to [0.01, 0.9]. One draw per layer and realization is shared by all that layer’s cells and timesteps.
products (sequence of str) – Nonempty subset of water_content and porosity (volume fractions).
n_realizations (int) – Positive number of independent parameter draws.
seed (int) – Local random seed (default 7); NumPy’s global RNG is never modified.
saturation_value (float) – Assumed saturation for porosity estimation, in (0, 1].
tortuosity_a (float) – Positive, dimensionless Archie tortuosity factor.
timestep_indices (sequence of int, optional) – Selected zero-based columns, in output order; defaults to all columns.
progress (callable) – Receives progress messages as strings.
return_realizations (bool) – Also return samples shaped (realizations, cells, selected times).
cell_chunk_size (int) – Maximum cells processed at once. Exact percentiles and the random draws are independent of this setting. Full sample storage is only allocated when return_realizations is True.
log_resistivity_std (array_like, optional) – The resistivity’s own uncertainty: a standard deviation of ln(resistivity) per cell, shaped (cells,) or (cells, times), such as an inversion’s linearised posterior. Each realization then multiplies the resistivity by exp(std * e), e standard normal, drawn independently per cell and time from a stream of its own (so the draws do not depend on cell_chunk_size), and the statistics carry both sources. Without it the resistivity is taken as exact.

Returns:

statistics maps products to mean, population std, p10, p50, p90 arrays shaped (cells, selected times). Porosity draws with no Waxman-Smits root (surface conduction alone exceeds the measured conductivity) are NaN, left out of the porosity statistics, and counted in statistics["porosity"]["no_root_fraction"]; a cell with no solved draw is NaN. warnings lists one message per layer that had such draws (each also issued as a RuntimeWarning). Also includes params_used, seed and timestep_indices.

Return type:

dict

Raises:

ValueError – Invalid shapes, non-finite/nonpositive resistivity, missing or duplicate layer labels, invalid distributions, or unsupported products.

PyHydroGeophysX.petrophysics.satK(Keff: Any, Km: Any, phi: Any, Sat: Any) → Any[source]#

Calculate the saturated bulk modulus (K_sat) based on Brie’s equation.

Parameters: Keff (float): Effective bulk modulus of the dry rock (GPa). Km (float): Bulk modulus of the matrix (GPa). phi (float): Porosity of the rock. Sat (float): Saturation level of the fluid in the pores.

Returns: float: Saturated bulk modulus (GPa).

PyHydroGeophysX.petrophysics.velDEM(phi: Any, Km: Any, Gm: Any, rho_b: Any, Sat: Any, alpha: Any) → Any[source]#

Calculate effective bulk modulus (Keff), shear modulus (Geff), and P-wave velocity (Vp) for a rock with varying porosity (phi) based on the DEM model, taking into account the saturation (Sat) and the crack aspect ratio (alpha).

Parameters: phi (np.array): Array of porosities. Km (float): Initial bulk modulus of the material (GPa). Gm (float): Initial shear modulus of the material (GPa). rho_b (float): Density of the solid phase (kg/m^3). Sat (array): Saturation fractions with one entry per phi value; scalars are not broadcast here. alpha (float): Crack aspect ratio.

Returns: Keff1 (np.array): Effective bulk modulus for each porosity value (GPa). Geff1 (np.array): Effective shear modulus for each porosity value (GPa). Vp (np.array): P-wave velocity for each porosity value (m/s).

PyHydroGeophysX.petrophysics.vel_porous(phi: Any, Km: Any, Gm: Any, rho_b: Any, Sat: Any, depth: Any = 1) → Any[source]#

Calculate P-wave velocity (Vp) for a rock with varying porosity (phi) based on the Hertz-Mindlin model and Hashin-Shtrikman bounds, taking into account the saturation (Sat).

Compatibility wrapper over HertzMindlinModel.calculate_velocity with critical porosity 0.4 and coordination number 4. Both return modes of that method and this function share the contact, HS and Gassmann implementation.

Parameters: phi (np.array): Array of porosities. Km (float): Bulk modulus of the solid phase (GPa). Gm (float): Shear modulus of the solid phase (GPa). rho_b (float): Density of the solid phase (kg/m^3). Sat (array): Saturation fractions; scalars and arrays broadcast with phi. depth (float): depth for pressure estimation (m)

Returns: Vp_h (np.array): P-wave velocity for each porosity value (upper Hashin-Shtrikman bound) (m/s). Vp_l (np.array): P-wave velocity for each porosity value (lower Hashin-Shtrikman bound) (m/s).

PyHydroGeophysX.petrophysics.velocity_to_water_content(velocity: ndarray, v_dry: float = 3500.0, v_sat: float = 4500.0, porosity: float | ndarray = 0.3, model: str = 'linear') → ndarray[source]#

Invert the empirical relationship used by water_content_to_velocity.

Supports linear, Wyllie and Raymer models with the same fixed constituent velocities as the forward function. Porosity must be known. Observations outside the model’s saturation range, or at the forward function’s clipping limits (200/8000 m/s), are rejected rather than silently assigned 0/1 saturation. Hertz-Mindlin and DEM require a separately calibrated inverse.

PyHydroGeophysX.petrophysics.water_content_to_resistivity(water_content: Any, rhos: Any, n: Any, porosity: Any, sigma_sur: Any = 0) → Any[source]#

Convert water content to resistivity using Waxman-Smits model.

Parameters:

water_content (array) – Volumetric water content (θ)
rhos (float) – Saturated resistivity without surface effects
n (float) – Saturation exponent
porosity (array) – Porosity values (φ)
sigma_sur (float) – Surface conductivity. Default is 0 (no surface effects).

Returns:

Resistivity values

Return type:

array

PyHydroGeophysX.petrophysics.water_content_to_velocity(water_content: ndarray, v_dry: float = 3500.0, v_sat: float = 4500.0, porosity: float | ndarray = 0.3, model: str = 'linear') → ndarray[source]#

Convert water content to seismic P-wave velocity using empirical relationships.

These are simplified mixing rules, not calibrated site predictions. Saturation is clipped to [0, 1] and output velocity to [200, 8000] m/s.

Parameters:

water_content – Volumetric water content (θ), typically 0-porosity
v_dry – Dry velocity for the linear model only (m/s), default 3500.
v_sat – Saturated velocity for the linear model only (m/s), default 4500.
porosity – Porosity values (φ), can be scalar or array
model – Velocity model type: - ‘linear’: Linear interpolation between dry and saturated velocities - ‘wyllie’: Wyllie time-average equation - ‘raymer’: Raymer-Hunt-Gardner equation

Returns:

P-wave velocity array (m/s)

PyHydroGeophysX.petrophysics.waxman_smits_surface_conductivity(sigma_w, qv)[source]#

Return the Waxman-Thomas surface coefficient in S/m at 25 deg C.

sigma_s = 4.6 * (1 - 0.6 * exp(-sigma_w / 1.3)) * qv, with sigma_w in S/m and qv in meq/cm3 pore. Thus 1.3 has units S/m and 4.6 has units (S/m)/(meq/cm3). Inputs broadcast.

Pass the result as sigma_s to WS_Model, which applies the formation factor. It is NOT the unscaled sigma_sur in the other converters. Reference: Waxman and Thomas (1974), doi:10.2118/4094-PA.