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.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.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