PyHydroGeophysX package#

Subpackages#

Module contents#

WatershedGeo - A comprehensive package for geophysical modeling and inversion in watershed monitoring.

This package integrates MODFLOW hydrological model outputs with geophysical forward modeling and inversion, specializing in electrical resistivity tomography (ERT) and seismic velocity modeling.

class PyHydroGeophysX.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.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.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.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) → Tuple[ndarray, ndarray][source]#

Calculate P-wave velocity for porous rocks.

Parameters:

  • porosity – Porosity values (array)

  • saturation – Saturation values (array)

  • bulk_modulus – Bulk modulus of the solid matrix (GPa)

  • shear_modulus – Shear modulus of the solid matrix (GPa)

  • mineral_density – Density of the solid matrix (kg/m³)

  • depth – Depth for pressure estimation (m, default: 1.0)

Returns:

Tuple of high-bound P-wave velocity (m/s) and low-bound P-wave velocity (m/s)

class PyHydroGeophysX.MODFLOWWaterContent(model_directory: str, idomain: ndarray)[source]#

Bases: HydroModelOutput

Class for processing water content data from MODFLOW simulations.

get_timestep_info() → List[Tuple[int, int, float, float]][source]#

Get information about each timestep in the WaterContent file.

Returns:

List of tuples (kstp, kper, pertim, totim) for each timestep

load_time_range(start_idx: int = 0, end_idx: int | None = None, nlay: int = 3) → ndarray[source]#

Load water content for a range of timesteps.

Parameters:

  • start_idx – Starting timestep index (default: 0)

  • end_idx – Ending timestep index (exclusive, default: None loads all)

  • nlay – Number of layers in the model (default: 3)

Returns:

Water content array with shape (timesteps, nlay, nrows, ncols)

load_timestep(timestep_idx: int, nlay: int = 3) → ndarray[source]#

Load water content for a specific timestep.

Parameters:

  • timestep_idx – Index of the timestep to load

  • nlay – Number of layers in the model

Returns:

Water content array with shape (nlay, nrows, ncols)

class PyHydroGeophysX.ProfileInterpolator(point1: List[int], point2: List[int], surface_data: ndarray, origin_x: float = 0.0, origin_y: float = 0.0, pixel_width: float = 1.0, pixel_height: float = -1.0, num_points: int = 200)[source]#

Bases: object

Class for handling interpolation of data to/from profiles.

interpolate_3d_data(data: ndarray) → ndarray[source]#

Interpolate 3D data (n_layers, ny, nx) to profile.

Parameters:

data – 3D array of values

Returns:

Array of interpolated values (n_layers, n_profile_points)

interpolate_layer_data(layer_data: List[ndarray]) → ndarray[source]#

Interpolate multiple layer data to profile.

Parameters:

layer_data – List of 2D arrays for each layer

Returns:

Array of interpolated values (n_layers, n_profile_points)

interpolate_to_mesh(property_values: ndarray, depth_values: ndarray, mesh_x: ndarray, mesh_y: ndarray, mesh_markers: ndarray, ID: ndarray, layer_markers: list = [3, 0, 2]) → ndarray[source]#

Interpolate property values from profile to mesh with layer-specific handling.

Parameters:

  • property_values – Property values array (n_points or n_layers, n_points)

  • depth_values – Depth values array (n_layers, n_points)

  • mesh_x – Coordinates of mesh cells

  • mesh_y – Coordinates of mesh cells

  • mesh_markers – Markers indicating different layers in mesh

  • layer_markers – List of marker values for each layer

Returns:

Interpolated values for mesh cells

class PyHydroGeophysX.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.VRH_model(f=[0.35, 0.25, 0.2, 0.125, 0.075], K=[55.4, 36.6, 75.6, 46.7, 50.4], G=[28.1, 45, 25.6, 23.65, 27.4], rho=[2560, 2650, 2630, 2540, 3050])[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.create_difference_gif(mesh, models: Sequence[ndarray], reference: ndarray, filename: str, *, mode: str = 'difference', titles: Sequence[str] | None = None, cmap: str = 'RdBu_r', symmetric: bool = True, label: str = '', figsize: Tuple[float, float] = (8, 2.5), dpi: int = 150, duration: int = 100, coverage: ndarray | Sequence[ndarray] | None = None) → str[source]#

Create a GIF animation showing model changes relative to a reference.

Parameters:

  • mesh (pygimli.Mesh) –

  • models (sequence of array-like) – Model snapshots for each frame.

  • reference (array-like) – Baseline model to subtract / divide.

  • mode ('difference' | 'ratio' | 'percent_change') –

  • symmetric (bool) – Center colorbar at zero / one.

:param (other parameters same as create_timelapse_gif()):

Returns:

Path to the saved GIF file.

Return type:

str

PyHydroGeophysX.create_timelapse_gif(mesh, models: Sequence[ndarray], filename: str, *, titles: Sequence[str] | None = None, cmap=None, cmin: float | None = None, cmax: float | None = None, log_scale: bool = False, label: str = '', xlabel: str = 'Distance (m)', ylabel: str = 'Elevation (m)', coverage: ndarray | Sequence[ndarray] | None = None, figsize: Tuple[float, float] = (8, 2.5), dpi: int = 150, duration: int = 100, first_frame_duration: int = 500, loop: int = 0) → str[source]#

Create a GIF animation of time-lapse model snapshots.

Each frame renders the model on the given mesh using pg.show and captures it as a PIL image. Requires Pillow.

Parameters:

  • mesh (pygimli.Mesh) – Mesh shared by all models.

  • models (sequence of array-like) – Model values for each frame.

  • filename (str) – Output GIF path.

  • titles (sequence of str, optional) – Title per frame.

  • cmap (str or Colormap, optional) – Colormap. Defaults to BlueDarkRed18_18_r if available.

  • cmin (float, optional) – Fixed color limits.

  • cmax (float, optional) – Fixed color limits.

  • log_scale (bool) – Logarithmic color scale.

  • label (str) – Colorbar label.

  • coverage (array or list of arrays, optional) – Coverage mask(s).

  • figsize (tuple) – Figure size per frame.

  • dpi (int) – Resolution.

  • duration (int) – Milliseconds per frame.

  • first_frame_duration (int) – Duration of the first frame in ms (longer for visual pause).

  • loop (int) – Number of loops (0 = infinite).

Returns:

Path to the saved GIF file.

Return type:

str

PyHydroGeophysX.create_timelapse_mp4(mesh, models: Sequence[ndarray], filename: str, *, titles: Sequence[str] | None = None, cmap=None, cmin: float | None = None, cmax: float | None = None, log_scale: bool = False, label: str = '', xlabel: str = 'Distance (m)', ylabel: str = 'Elevation (m)', coverage: ndarray | Sequence[ndarray] | None = None, figsize: Tuple[float, float] = (8, 2.5), dpi: int = 150, fps: int = 10) → str[source]#

Create an MP4 video of time-lapse model snapshots using matplotlib.

Requires ffmpeg to be available on the system path.

Parameters:

  • mesh (pygimli.Mesh) –

  • models (sequence of array-like) –

  • filename (str) – Output .mp4 path.

  • fps (int) – Frames per second.

:param (other parameters same as create_timelapse_gif()):

Returns:

Path to the saved MP4 file.

Return type:

str

PyHydroGeophysX.export_mesh_to_vtk(mesh, filename: str, cell_data: Dict[str, ndarray] | None = None) → str[source]#

Export a PyGIMLi mesh to VTK unstructured grid format.

This uses PyGIMLi’s built-in mesh.exportVTK when available and then optionally injects extra cell data fields.

Parameters:

  • mesh (pygimli.Mesh) – The mesh to export.

  • filename (str) – Output .vtk path.

  • cell_data (dict of str -> array, optional) – Additional named cell-data arrays to write.

Returns:

Path to the written file.

Return type:

str

PyHydroGeophysX.export_points_to_vtk(points: ndarray, scalars: Dict[str, ndarray], filename: str) → str[source]#

Export point-cloud data (e.g. cell centers) to VTK PolyData.

Parameters:

  • points (array of shape (N, 3)) – XYZ coordinates.

  • scalars (dict of str -> array) – Named scalar values at each point.

  • filename (str) – Output .vtk path.

Returns:

Path to the written file.

Return type:

str

PyHydroGeophysX.export_structured_vtk(values: ndarray, filename: str, *, dx: float = 1.0, dy: float = 1.0, dz: float = 1.0, origin: tuple = (0.0, 0.0, 0.0), scalar_name: str = 'model') → str[source]#

Export a 3-D numpy array to VTK structured-points format.

Useful for MODFLOW/ParFlow grids that map directly to a regular grid.

Parameters:

  • values (3-D array) – Shape (nz, ny, nx) — layer, row, column ordering.

  • filename (str) – Output .vtk path.

  • dx (float) – Cell spacing in each direction.

  • dy (float) – Cell spacing in each direction.

  • dz (float) – Cell spacing in each direction.

  • origin (tuple of float) – (x0, y0, z0) origin of the grid.

  • scalar_name (str) – Name for the scalar dataset.

Returns:

Path to the written file.

Return type:

str

PyHydroGeophysX.export_structured_vtk_multi(scalars: Dict[str, ndarray], filename: str, *, dx: float = 1.0, dy: float = 1.0, dz: float = 1.0, origin: tuple = (0.0, 0.0, 0.0)) → str[source]#

Export multiple 3-D arrays as named scalars in a single VTK file.

All arrays must have the same shape (nz, ny, nx).

Parameters:

  • scalars (dict of str -> 3-D array) – Named scalar datasets (e.g. {"resistivity": res, "water_content": wc}).

  • filename (str) – Output .vtk path.

Returns:

Path to the written file.

Return type:

str

PyHydroGeophysX.export_timelapse_structured_vtk(models: Sequence[ndarray], output_dir: str, prefix: str = 'timelapse', *, dx: float = 1.0, dy: float = 1.0, dz: float = 1.0, origin: tuple = (0.0, 0.0, 0.0), scalar_name: str = 'model') → List[str][source]#

Export time-lapse 3-D arrays as a numbered VTK series.

Parameters:

  • models (sequence of 3-D arrays) – Each entry has shape (nz, ny, nx).

  • output_dir (str) – Output directory.

  • prefix (str) – File prefix.

  • dx (float) – Grid spacing.

  • dy (float) – Grid spacing.

  • dz (float) – Grid spacing.

  • origin (tuple) – Grid origin.

  • scalar_name (str) – Scalar field name.

Returns:

Paths to generated VTK files.

Return type:

list of str

PyHydroGeophysX.export_timelapse_vtk(mesh, models: Sequence[ndarray], output_dir: str, prefix: str = 'timelapse', *, scalar_name: str = 'model', extra_data: Dict[str, Sequence[ndarray]] | None = None) → List[str][source]#

Export time-lapse models as a numbered VTK series for ParaView.

Creates files prefix_0000.vtk, prefix_0001.vtk, … and a .pvd collection file that ParaView can load as a time series.

Parameters:

  • mesh (pygimli.Mesh) – Mesh shared by all timesteps.

  • models (sequence of array-like) – Model values per timestep.

  • output_dir (str) – Directory for output files.

  • prefix (str) – Filename prefix.

  • scalar_name (str) – Name for the primary scalar field.

  • extra_data (dict of str -> sequence of arrays, optional) – Additional scalar fields per timestep.

Returns:

Paths to all generated VTK files.

Return type:

list of str

PyHydroGeophysX.interpolate_structure_to_profile(structure_data: List[ndarray], X_grid: ndarray, Y_grid: ndarray, X_pro: ndarray, Y_pro: ndarray) → ndarray[source]#

Interpolate multiple structure layers onto profile

Parameters:

  • structure_data – List of 2D arrays for each layer

  • X_grid – X coordinates of original grid

  • Y_grid – Y coordinates of original grid

  • X_pro – X coordinates of profile points

  • Y_pro – Y coordinates of profile points

Returns:

Array of interpolated values with shape (n_layers, n_points)

PyHydroGeophysX.interpolate_to_mesh(property_values: ndarray, profile_distance: ndarray, depth_values: ndarray, mesh_x: ndarray, mesh_y: ndarray, mesh_markers: ndarray, ID, layer_markers: list = [3, 0, 2]) → ndarray[source]#

Interpolate property values from profile to mesh with layer-specific handling.

Parameters:

  • property_values – Property values array (n_points)

  • profile_distance – Distance along profile (n_points)

  • depth_values – Depth values array (n_layers, n_points)

  • mesh_x – X coordinates of mesh cells

  • mesh_y – Y coordinates of mesh cells

  • mesh_markers – Markers indicating different layers in mesh

  • layer_markers – List of marker values for each layer

Returns:

Interpolated values for mesh cells

PyHydroGeophysX.interpolate_to_profile(data: ndarray, X_grid: ndarray, Y_grid: ndarray, X_pro: ndarray, Y_pro: ndarray, method: str = 'linear') → ndarray[source]#

Interpolate 2D data onto a profile line

Parameters:

  • data – 2D array of values to interpolate

  • X_grid – X coordinates of original grid (meshgrid)

  • Y_grid – Y coordinates of original grid (meshgrid)

  • X_pro – X coordinates of profile points

  • Y_pro – Y coordinates of profile points

  • method – Interpolation method (‘linear’ or ‘nearest’)

Returns:

Interpolated values along profile

PyHydroGeophysX.plot_convergence(chi2_history: Sequence[float], *, ax=None, target_chi2: float = 1.0, ylabel: str = '$\\chi^2$', title: str = 'Inversion Convergence') → Tuple[source]#

Plot chi-squared convergence curve.

Parameters:

  • chi2_history (sequence of float) – Chi-squared value per iteration.

  • target_chi2 (float) – Target misfit (plotted as a dashed line).

Return type:

fig, ax

PyHydroGeophysX.plot_coverage(mesh, coverage: ndarray, *, ax=None, cmap: str = 'YlGn', threshold: float | None = None, title: str = 'Data Coverage') → Tuple[source]#

Plot a coverage / sensitivity map.

Parameters:

  • mesh (pygimli.Mesh) –

  • coverage (array-like) – Coverage values per cell.

  • threshold (float, optional) – If given, overlay a contour at this level.

Return type:

fig, ax, cbar

PyHydroGeophysX.plot_difference_map(mesh, model_a: ndarray, model_b: ndarray, *, mode: str = 'difference', ax=None, cmap: str = 'RdBu_r', symmetric: bool = True, label: str = '', title: str = '', coverage: ndarray | None = None) → Tuple[source]#

Plot the difference or ratio between two models.

Parameters:

  • mesh (pygimli.Mesh) –

  • model_a (array-like) – Two model arrays. result = model_b - model_a (difference) or model_b / model_a (ratio).

  • model_b (array-like) – Two model arrays. result = model_b - model_a (difference) or model_b / model_a (ratio).

  • mode ('difference' | 'ratio' | 'percent_change') –

  • symmetric (bool) – If True, center the colorbar at zero (difference) or one (ratio).

  • coverage (array-like, optional) – Coverage mask.

Return type:

fig, ax, cbar

PyHydroGeophysX.plot_electrode_layout(positions: Dict[str, ndarray], *, ax=None, color_by: str = 'z', cmap: str = 'terrain', title: str = 'Electrode Layout') → Tuple[source]#

Scatter-plot electrode positions colored by elevation.

Parameters:

positions (dict) – Must contain 'x' and 'y' keys; optionally 'z'.

Return type:

fig, ax

PyHydroGeophysX.plot_model_section(mesh, values: ndarray, *, ax=None, cmap=None, cmin: float | None = None, cmax: float | None = None, log_scale: bool = False, label: str = '', xlabel: str = 'Distance (m)', ylabel: str = 'Elevation (m)', title: str = '', coverage: ndarray | None = None, orientation: str = 'vertical') → Tuple[source]#

Plot a 2D model cross-section on a PyGIMLi mesh.

Parameters:

  • mesh (pygimli.Mesh) – The mesh to plot on.

  • values (array-like) – Cell values (resistivity, velocity, water content, etc.).

  • ax (matplotlib.axes.Axes, optional) – Axes to draw on. Created if None.

  • cmap (str or Colormap, optional) – Colormap. Defaults to BlueDarkRed18_18_r if available.

  • cmin (float, optional) – Color limits.

  • cmax (float, optional) – Color limits.

  • log_scale (bool) – Use logarithmic color scaling.

  • label (str) – Colorbar label.

  • coverage (array-like, optional) – Coverage array for masking low-sensitivity cells.

  • orientation (str) – Colorbar orientation ('vertical' or 'horizontal').

Return type:

fig, ax, cbar

PyHydroGeophysX.plot_monitoring_timeseries(times: ndarray, series: Dict[str, ndarray], *, true_series: Dict[str, ndarray] | None = None, uncertainties: Dict[str, Tuple[ndarray, ndarray]] | None = None, ax=None, ylabel: str = 'Value', title: str = 'Monitoring Point Time Series') → Tuple[source]#

Plot estimated (and optionally true) time-series at monitoring points.

Parameters:

  • times (array-like) – Time axis.

  • series (dict of str -> array) – Estimated values keyed by point name.

  • true_series (dict of str -> array, optional) – True / reference values for comparison (dashed lines).

  • uncertainties (dict of str -> (lower, upper), optional) – Uncertainty bounds per point for shading.

Return type:

fig, ax

PyHydroGeophysX.plot_pseudosection_matrix(data_matrix: ndarray, *, ax=None, cmap=None, vmin: float | None = None, vmax: float | None = None, xlabel: str = 'Time', ylabel: str = 'Measurement #', label: str = 'Apparent resistivity ($\\Omega\\cdot$m)', title: str = '') → Tuple[source]#

Plot a time-lapse apparent resistivity matrix as a heatmap.

Parameters:

data_matrix (2-D array) – Shape (n_times, n_measurements) or similar.

Return type:

fig, ax, im

PyHydroGeophysX.plot_timelapse_snapshots(mesh, models: Sequence[ndarray], *, titles: Sequence[str] | None = None, ncols: int = 4, cmap=None, cmin: float | None = None, cmax: float | None = None, log_scale: bool = False, label: str = '', coverage: ndarray | Sequence[ndarray] | None = None, figsize_per_panel: Tuple[float, float] = (4.0, 2.5)) → Tuple[source]#

Plot a grid of time-lapse model snapshots.

Parameters:

  • mesh (pygimli.Mesh) – Mesh shared by all snapshots.

  • models (sequence of array-like) – Model arrays for each timestep.

  • titles (sequence of str, optional) – Panel titles. Defaults to 'Timestep 1', 'Timestep 2', …

  • ncols (int) – Number of columns.

  • cmap – Passed to pg.show.

  • cmin – Passed to pg.show.

  • cmax – Passed to pg.show.

  • log_scale – Passed to pg.show.

  • label – Passed to pg.show.

  • coverage (array or sequence of arrays, optional) – Coverage mask(s). If a single 1-D array it is reused for all panels. If 2-D, coverage[i] is used for panel i.

  • figsize_per_panel (tuple) – (width, height) per subplot panel.

Return type:

fig, axes

PyHydroGeophysX.plot_topography(topo_grid: ndarray, *, profile_endpoints: List[Tuple[float, float]] | None = None, ax=None, cmap: str = 'terrain', title: str = 'Surface Topography') → Tuple[source]#

Plot a 2-D topography grid with optional profile line overlay.

Parameters:

  • topo_grid (2-D array) – Elevation raster.

  • profile_endpoints (list of (row, col) tuples, optional) – If two points are given, draw the profile line.

Return type:

fig, ax

PyHydroGeophysX.prepare_2D_profile_data(data: ndarray, XX: ndarray, YY: ndarray, X_pro: ndarray, Y_pro: ndarray) → ndarray[source]#

Interpolate multiple 2D gridded data layers onto a profile line.

Parameters:

  • data – 3D array of gridded data (n_layers, ny, nx)

  • XX – Coordinate grids from meshgrid

  • YY – Coordinate grids from meshgrid

  • X_pro – Profile line coordinates

  • Y_pro – Profile line coordinates

Returns:

Interpolated values along profile (n_layers, n_profile_points)

PyHydroGeophysX.resistivity_to_saturation(resistivity, porosity, m, rho_fluid, n, sigma_sur=0, a=1.0)[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)

Return type:

array

PyHydroGeophysX.resistivity_to_water_content(resistivity, rhos, n, porosity, sigma_sur=0)[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.satK(Keff, Km, phi, Sat)[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.setup_profile_coordinates(point1: List[int], point2: List[int], surface_data: ndarray, origin_x: float = 0.0, origin_y: float = 0.0, pixel_width: float = 1.0, pixel_height: float = -1.0, num_points: int = 200) → Tuple[ndarray, ndarray, ndarray, ndarray, ndarray][source]#

Set up profile coordinates based on surface elevation data between two points

Parameters:

  • point1 – Starting point indices [col, row]

  • point2 – Ending point indices [col, row]

  • surface_data – 2D array of surface elevation data

  • origin_x – X coordinate of origin

  • origin_y – Y coordinate of origin

  • pixel_width – Width of each pixel

  • pixel_height – Height of each pixel (negative for top-down)

  • num_points – Number of points along profile

Returns:

X coordinates along profile Y_pro: Y coordinates along profile L_profile: Distances along profile XX: X coordinate grid YY: Y coordinate grid

Return type:

X_pro

PyHydroGeophysX.velDEM(phi, Km, Gm, rho_b, Sat, alpha)[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 (float): Saturation level of the fluid in the cracks (0 to 1, where 1 is fully saturated). 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.vel_porous(phi, Km, Gm, rho_b, Sat, depth=1)[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).

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 (float): Saturation level of the fluid in the pores (0 to 1, where 1 is fully saturated). 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.water_content_to_resistivity(water_content, rhos, n, porosity, sigma_sur=0)[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