gwtransport

Raises:

ValueError – If the cumulative flow over the entire flow_tedges window does not reach retardation_factor * aquifer_pore_volume, indicating the flow timeseries is too short to fully characterise the aquifer.

diffusion_fast#

Fast Diffusive Transport Corrections via Gaussian Smoothing.

This module provides a computationally efficient approximation of diffusion/dispersion using Gaussian smoothing. It is much faster than gwtransport.diffusion but less physically accurate, especially under variable flow conditions.

Both diffusion_fast and gwtransport.diffusion add microdispersion and molecular diffusion on top of macrodispersion captured by the APVD.

When to use diffusion_fast vs diffusion:

Use diffusion_fast when: Speed is critical, flow and time steps are relatively constant, or you need real-time processing
Use diffusion when: Physical accuracy is critical, flow varies significantly, or you’re analyzing periods with changing conditions

See Macrodispersion and Microdispersion for background on macrodispersion and microdispersion.

This module implements diffusion/dispersion processes that modify advective transport in aquifer systems. Diffusion causes spreading and smoothing of concentration or temperature fronts as they travel through the aquifer. While advection moves compounds with water flow, diffusion causes spreading due to molecular diffusion, mechanical dispersion, and thermal diffusion (for temperature).

Limitation: When flow_out is not provided, this fast approximation works best when flow and tedges are relatively constant. The underlying assumption is that dx (spatial step between cells) remains approximately constant, which holds for steady flow but breaks down under highly variable conditions. When input time bins are non-uniform (e.g., from signal compression), provide flow_out to apply Gaussian smoothing on the (typically uniform) output grid instead. For scenarios with significant flow variability, consider using gwtransport.diffusion instead.

Available functions:

infiltration_to_extraction() - Apply diffusion during infiltration to extraction transport. Combines advection (via residence time) with diffusion (via Gaussian smoothing). Computes position-dependent diffusion based on local residence time and returns concentration or temperature in extracted water on the extraction time grid.
extraction_to_infiltration() - Inverse diffusion via Tikhonov regularization. Builds the combined forward matrix (advection + Gaussian diffusion) and solves the inverse problem to reconstruct infiltration concentrations from extraction data.
gamma_infiltration_to_extraction() - Gamma-distributed pore volumes with fast Gaussian diffusion. Convenience wrapper that parameterizes the pore volume distribution as a gamma distribution.
gamma_extraction_to_infiltration() - Gamma-distributed pore volumes, inverse fast Gaussian diffusion. Symmetric inverse of gamma_infiltration_to_extraction.
convolve_diffusion() - Apply variable-sigma Gaussian filtering. Extends scipy.ndimage.gaussian_filter1d to position-dependent sigma using sparse matrix representation for efficiency. Handles boundary conditions via nearest-neighbor extrapolation.
create_example_data() - Generate test data for demonstrating diffusion effects with signals having varying time steps and corresponding sigma arrays. Useful for testing and validation.

This file is part of gwtransport which is released under AGPL-3.0 license. See the ./LICENSE file or go to gwtransport/gwtransport for full license details.

gwtransport.diffusion_fast.infiltration_to_extraction(*, cin, flow, tedges, cout_tedges, aquifer_pore_volumes, mean_streamline_length, mean_molecular_diffusivity, mean_longitudinal_dispersivity, retardation_factor=1.0, suppress_dispersion_warning=False, asymptotic_cutoff_sigma=3.0, flow_out=None, suppress_uniform_tedges_check=False)[source]#

Compute diffusion during 1D transport using fast Gaussian smoothing.

Combines advection (via pore volume distribution) with diffusion (via Gaussian smoothing) to produce extraction concentrations on the cout_tedges time grid. This is the fast approximate counterpart of gwtransport.diffusion.infiltration_to_extraction().

When flow_out is provided, diffusion is applied on the output grid (advect-then-smooth), which correctly handles non-uniform tedges (e.g., from signal compression). Without flow_out, diffusion is applied on the input grid and tedges must have constant time steps.

Parameters:

cin (TypeAliasType) – Concentration or temperature of the compound in the infiltrating water.
flow (TypeAliasType) – Flow rate of water in the aquifer [m3/day].
tedges (DatetimeIndex) – Time edges for cin and flow data. Has length len(cin) + 1.
cout_tedges (DatetimeIndex) – Time edges for output data bins. Has length of desired output + 1.
aquifer_pore_volumes (TypeAliasType) – Array of aquifer pore volumes [m3] representing the distribution of flow paths.
mean_streamline_length (float) – Mean travel distance [m]. Must be positive.
mean_molecular_diffusivity (float) – Mean effective molecular diffusivity [m2/day]. Must be non-negative.
mean_longitudinal_dispersivity (float) – Mean longitudinal dispersivity [m]. Must be non-negative.
retardation_factor (float, default: 1.0) – Retardation factor (default 1.0). Values > 1.0 indicate slower transport.
suppress_dispersion_warning (bool, default: False) – Suppress warning about combining multiple pore volumes with dispersivity. Default is False.
asymptotic_cutoff_sigma (float | None, default: 3.0) – Truncate the Gaussian kernel at this many standard deviations. Default is 3.0.
flow_out (TypeAliasType | None, default: None) – Flow rate [m3/day] on the output grid (aligned to cout_tedges). Required when tedges has non-uniform time steps. Length must equal len(cout_tedges) - 1. Default is None.
suppress_uniform_tedges_check (bool, default: False) – Skip the constant-dt check on tedges when flow_out is None. Default is False.

Returns:

Bin-averaged concentration in extracted water. Length equals len(cout_tedges) - 1. NaN values indicate time periods with no valid contributions from the infiltration data.

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

See also

gwtransport.diffusion.infiltration_to_extraction: Physically rigorous analytical solution that supports per-pore-volume arrays for streamline_length, molecular_diffusivity, and longitudinal_dispersivity.
extraction_to_infiltration: Inverse operation
Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity: Macrodispersion vs microdispersion

Notes

A single mean_streamline_length is shared across all pore volumes. This assumes streamline-length heterogeneity is captured solely through the pore-volume distribution: the effective cross-sectional area A_eff = V_p / L_mean varies with V_p while the path length L is held fixed. This is appropriate for many systems but breaks down for partially-penetrating wells or wedge-shaped capture zones, where path length itself varies between streamtubes. In those cases, use gwtransport.diffusion.infiltration_to_extraction() with a per-streamtube streamline_length array instead.

gwtransport.diffusion_fast.extraction_to_infiltration(*, cout, flow, tedges, cout_tedges, aquifer_pore_volumes, mean_streamline_length, mean_molecular_diffusivity, mean_longitudinal_dispersivity, retardation_factor=1.0, suppress_dispersion_warning=False, asymptotic_cutoff_sigma=3.0, regularization_strength=1e-10, flow_out=None, suppress_uniform_tedges_check=False)[source]#

Reconstruct infiltration concentration from extracted water via Tikhonov inversion.

Inverts the forward transport model (Gaussian diffusion + advection) by building the combined forward matrix and solving via Tikhonov regularization. This is the fast approximate counterpart of gwtransport.diffusion.extraction_to_infiltration().

When flow_out is provided, diffusion is applied on the output grid, which correctly handles non-uniform tedges. Without flow_out, diffusion is applied on the input grid and tedges must have constant time steps.

Parameters:

cout (TypeAliasType) – Concentration of the compound in extracted water.
flow (TypeAliasType) – Flow rate of water in the aquifer [m3/day].
tedges (DatetimeIndex) – Time edges for cin (output) and flow data. Has length of len(flow) + 1.
cout_tedges (DatetimeIndex) – Time edges for cout data bins. Has length of len(cout) + 1.
aquifer_pore_volumes (TypeAliasType) – Array of aquifer pore volumes [m3].
mean_streamline_length (float) – Mean travel distance [m]. Must be positive.
mean_molecular_diffusivity (float) – Mean effective molecular diffusivity [m2/day]. Must be non-negative.
mean_longitudinal_dispersivity (float) – Mean longitudinal dispersivity [m]. Must be non-negative.
retardation_factor (float, default: 1.0) – Retardation factor (default 1.0).
suppress_dispersion_warning (bool, default: False) – Suppress warning about combining multiple pore volumes with dispersivity. Default is False.
asymptotic_cutoff_sigma (float | None, default: 3.0) – Truncate the Gaussian kernel at this many standard deviations. Default is 3.0.
regularization_strength (float, default: 1e-10) – Tikhonov regularization parameter. Default is 1e-10.
flow_out (TypeAliasType | None, default: None) – Flow rate [m3/day] on the output grid (aligned to cout_tedges). Required when tedges has non-uniform time steps. Length must equal len(cout_tedges) - 1. Default is None.
suppress_uniform_tedges_check (bool, default: False) – Skip the constant-dt check on tedges when flow_out is None. Default is False.

Returns:

Bin-averaged concentration in infiltrating water. Length equals len(tedges) - 1. NaN values indicate time periods with no valid contributions from the extraction data.

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

Warns:

UserWarning – When the forward matrix is rank-deficient. This occurs with constant flow when the residence time is an integer multiple of the time step width. To fix, adjust aquifer_pore_volumes slightly (e.g., multiply by 1.001).

See also

infiltration_to_extraction: Forward operation
gwtransport.diffusion.extraction_to_infiltration: Analytically correct deconvolution that supports per-pore-volume arrays for streamline_length, molecular_diffusivity, and longitudinal_dispersivity.
Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity: Macrodispersion vs microdispersion

Notes

A single mean_streamline_length is shared across all pore volumes. This assumes streamline-length heterogeneity is captured solely through the pore-volume distribution: the effective cross-sectional area A_eff = V_p / L_mean varies with V_p while the path length L is held fixed. This is appropriate for many systems but breaks down for partially-penetrating wells or wedge-shaped capture zones, where path length itself varies between streamtubes. In those cases, use gwtransport.diffusion.extraction_to_infiltration() with a per-streamtube streamline_length array instead.

gwtransport.diffusion_fast.gamma_infiltration_to_extraction(*, cin, flow, tedges, cout_tedges, mean=None, std=None, loc=0.0, alpha=None, beta=None, n_bins=100, mean_streamline_length, mean_molecular_diffusivity, mean_longitudinal_dispersivity, retardation_factor=1.0, suppress_dispersion_warning=False, asymptotic_cutoff_sigma=3.0, flow_out=None, suppress_uniform_tedges_check=False)[source]#

Compute extracted concentration with fast Gaussian diffusion for gamma-distributed pore volumes.

Combines advective transport (based on gamma-distributed pore volumes) with fast Gaussian diffusion. This is a convenience wrapper around infiltration_to_extraction() that parameterizes the aquifer pore volume distribution as a (shifted) gamma distribution.

Provide either (mean, std) or (alpha, beta); loc is optional and defaults to 0.

Parameters:

cin (TypeAliasType) – Concentration of the compound in infiltrating water.
flow (TypeAliasType) – Flow rate of water in the aquifer [m3/day].
tedges (DatetimeIndex) – Time edges for cin and flow data. Has length len(cin) + 1.
cout_tedges (DatetimeIndex) – Time edges for output data bins. Has length of desired output + 1.
mean (float | None, default: None) – Mean of the gamma distribution of the aquifer pore volume. Must be strictly greater than loc.
std (float | None, default: None) – Standard deviation of the gamma distribution of the aquifer pore volume (invariant under the loc shift).
loc (float, default: 0.0) – Location (minimum pore volume) of the gamma distribution. Must satisfy 0 <= loc < mean. Default is 0.0.
alpha (float | None, default: None) – Shape parameter of gamma distribution of the aquifer pore volume (must be > 0).
beta (float | None, default: None) – Scale parameter of gamma distribution of the aquifer pore volume (must be > 0).
n_bins (int, default: 100) – Number of bins to discretize the gamma distribution. Default is 100.
mean_streamline_length (float) – Mean travel distance through the aquifer [m] averaged over all aquifer pore volumes. Must be positive. For per-pore-volume arrays, use gwtransport.diffusion.gamma_infiltration_to_extraction().
mean_molecular_diffusivity (float) – Mean effective molecular diffusivity [m2/day] averaged over all aquifer pore volumes. Must be non-negative. For per-pore-volume arrays, use gwtransport.diffusion.gamma_infiltration_to_extraction().
mean_longitudinal_dispersivity (float) – Mean longitudinal dispersivity [m] averaged over all aquifer pore volumes. Must be non-negative. For per-pore-volume arrays, use gwtransport.diffusion.gamma_infiltration_to_extraction().
retardation_factor (float, default: 1.0) – Retardation factor (default 1.0). Values > 1.0 indicate slower transport.
suppress_dispersion_warning (bool, default: False) – Suppress warning about combining multiple pore volumes with dispersivity. Default is False.
asymptotic_cutoff_sigma (float | None, default: 3.0) – Truncate the Gaussian kernel at this many standard deviations. Default is 3.0.
flow_out (TypeAliasType | None, default: None) – Flow rate [m3/day] on the output grid (aligned to cout_tedges). Required when tedges has non-uniform time steps. Length must equal len(cout_tedges) - 1. Default is None.
suppress_uniform_tedges_check (bool, default: False) – Skip the constant-dt check on tedges when flow_out is None. Default is False.

Returns:

Bin-averaged concentration in extracted water. Length equals len(cout_tedges) - 1. NaN values indicate time periods with no valid contributions from the infiltration data.

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

See also

infiltration_to_extraction: Transport with explicit pore volume distribution
gamma_extraction_to_infiltration: Reverse operation (deconvolution)
gwtransport.gamma.bins: Create gamma distribution bins
gwtransport.diffusion.gamma_infiltration_to_extraction: Physically rigorous analytical solution that supports per-pore-volume arrays for streamline_length, molecular_diffusivity, and longitudinal_dispersivity.
gwtransport.advection.gamma_infiltration_to_extraction: Pure advection (no dispersion)
Gamma Distribution Model: Two-parameter pore volume model
Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity: Macrodispersion vs microdispersion

Notes

The APVD is only time-invariant under the steady-streamlines assumption (see 2. Steady Streamlines).

The spreading from the gamma-distributed pore volumes represents macrodispersion (aquifer-scale heterogeneity). When std comes from calibration on measurements, it absorbs all mixing: macrodispersion, microdispersion, and an average molecular diffusion contribution. When std comes from streamline analysis, it represents macrodispersion only; microdispersion and molecular diffusion can be added via the dispersion parameters. See Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity for guidance on when to add microdispersion.

Examples

>>> import pandas as pd
>>> import numpy as np
>>> from gwtransport.diffusion_fast import gamma_infiltration_to_extraction
>>>
>>> tedges = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
>>> cout_tedges = pd.date_range(start="2020-01-05", end="2020-01-25", freq="D")
>>> cin = np.zeros(len(tedges) - 1)
>>> cin[5:10] = 1.0
>>> flow = np.ones(len(tedges) - 1) * 100.0
>>>
>>> cout = gamma_infiltration_to_extraction(
...     cin=cin,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     mean=500.0,
...     std=100.0,
...     n_bins=5,
...     mean_streamline_length=100.0,
...     mean_molecular_diffusivity=1e-4,
...     mean_longitudinal_dispersivity=1.0,
... )

gwtransport.diffusion_fast.gamma_extraction_to_infiltration(*, cout, flow, tedges, cout_tedges, mean=None, std=None, loc=0.0, alpha=None, beta=None, n_bins=100, mean_streamline_length, mean_molecular_diffusivity, mean_longitudinal_dispersivity, retardation_factor=1.0, suppress_dispersion_warning=False, asymptotic_cutoff_sigma=3.0, regularization_strength=1e-10, flow_out=None, suppress_uniform_tedges_check=False)[source]#

Reconstruct infiltration concentration from extracted water for gamma-distributed pore volumes.

Inverts the forward transport model (fast Gaussian diffusion + advection with gamma-distributed pore volumes) via Tikhonov regularization. This is a convenience wrapper around extraction_to_infiltration() that parameterizes the aquifer pore volume distribution as a (shifted) gamma distribution.

Provide either (mean, std) or (alpha, beta); loc is optional and defaults to 0.

Parameters:

cout (TypeAliasType) – Concentration of the compound in extracted water.
flow (TypeAliasType) – Flow rate of water in the aquifer [m3/day].
tedges (DatetimeIndex) – Time edges for cin (output) and flow data. Has length of len(flow) + 1.
cout_tedges (DatetimeIndex) – Time edges for cout data bins. Has length of len(cout) + 1.
mean (float | None, default: None) – Mean of the gamma distribution of the aquifer pore volume. Must be strictly greater than loc.
std (float | None, default: None) – Standard deviation of the gamma distribution of the aquifer pore volume (invariant under the loc shift).
loc (float, default: 0.0) – Location (minimum pore volume) of the gamma distribution. Must satisfy 0 <= loc < mean. Default is 0.0.
alpha (float | None, default: None) – Shape parameter of gamma distribution of the aquifer pore volume (must be > 0).
beta (float | None, default: None) – Scale parameter of gamma distribution of the aquifer pore volume (must be > 0).
n_bins (int, default: 100) – Number of bins to discretize the gamma distribution. Default is 100.
mean_streamline_length (float) – Mean travel distance through the aquifer [m] averaged over all aquifer pore volumes. Must be positive. For per-pore-volume arrays, use gwtransport.diffusion.gamma_extraction_to_infiltration().
mean_molecular_diffusivity (float) – Mean effective molecular diffusivity [m2/day] averaged over all aquifer pore volumes. Must be non-negative. For per-pore-volume arrays, use gwtransport.diffusion.gamma_extraction_to_infiltration().
mean_longitudinal_dispersivity (float) – Mean longitudinal dispersivity [m] averaged over all aquifer pore volumes. Must be non-negative. For per-pore-volume arrays, use gwtransport.diffusion.gamma_extraction_to_infiltration().
retardation_factor (float, default: 1.0) – Retardation factor (default 1.0). Values > 1.0 indicate slower transport.
suppress_dispersion_warning (bool, default: False) – Suppress warning about combining multiple pore volumes with dispersivity. Default is False.
asymptotic_cutoff_sigma (float | None, default: 3.0) – Truncate the Gaussian kernel at this many standard deviations. Default is 3.0.
regularization_strength (float, default: 1e-10) – Tikhonov regularization parameter. Default is 1e-10.
flow_out (TypeAliasType | None, default: None) – Flow rate [m3/day] aligned to cout_tedges. When provided, the Gaussian matrix operates on the output grid. Length must equal len(cout_tedges) - 1. Default is None.
suppress_uniform_tedges_check (bool, default: False) – When True, skip the check that tedges has constant time steps (only relevant when flow_out is None). Default is False.

Returns:

Bin-averaged concentration in infiltrating water. Length equals len(tedges) - 1. NaN values indicate time periods with no valid contributions from the extraction data.

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

See also

extraction_to_infiltration: Deconvolution with explicit pore volume distribution
gamma_infiltration_to_extraction: Forward operation (convolution)
gwtransport.gamma.bins: Create gamma distribution bins
gwtransport.diffusion.gamma_extraction_to_infiltration: Physically rigorous analytical solution that supports per-pore-volume arrays for streamline_length, molecular_diffusivity, and longitudinal_dispersivity.
gwtransport.advection.gamma_extraction_to_infiltration: Pure advection (no dispersion)
Gamma Distribution Model: Two-parameter pore volume model
Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity: Macrodispersion vs microdispersion

Notes

The APVD is only time-invariant under the steady-streamlines assumption (see 2. Steady Streamlines).

The spreading from the gamma-distributed pore volumes represents macrodispersion (aquifer-scale heterogeneity). When std comes from calibration on measurements, it absorbs all mixing: macrodispersion, microdispersion, and an average molecular diffusion contribution. When std comes from streamline analysis, it represents macrodispersion only; microdispersion and molecular diffusion can be added via the dispersion parameters. See Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity for guidance on when to add microdispersion.

Examples

>>> import pandas as pd
>>> import numpy as np
>>> from gwtransport.diffusion_fast import gamma_extraction_to_infiltration
>>>
>>> tedges = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
>>> cout_tedges = pd.date_range(start="2020-01-05", end="2020-01-25", freq="D")
>>> cout = np.zeros(len(cout_tedges) - 1)
>>> cout[5:10] = 1.0
>>> flow = np.ones(len(tedges) - 1) * 100.0
>>>
>>> cin = gamma_extraction_to_infiltration(
...     cout=cout,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     mean=500.0,
...     std=100.0,
...     n_bins=5,
...     mean_streamline_length=100.0,
...     mean_molecular_diffusivity=1e-4,
...     mean_longitudinal_dispersivity=1.0,
... )

gwtransport.diffusion_fast.convolve_diffusion(*, input_signal, sigma_array, asymptotic_cutoff_sigma=3.0)[source]#

Apply Gaussian filter with position-dependent sigma values.

This function extends scipy.ndimage.gaussian_filter1d by allowing the standard deviation (sigma) of the Gaussian kernel to vary at each point in the signal. It implements the filter using a sparse convolution matrix where each row represents a Gaussian kernel with a locally-appropriate standard deviation.

Kernel weights are computed by integrating the Gaussian CDF over each bin, which is more accurate than point-sampling the PDF for small sigma values.

Parameters:

input_signal (TypeAliasType) – One-dimensional input array to be filtered.
sigma_array (TypeAliasType) – One-dimensional array of standard deviation values, must have same length as input_signal.
asymptotic_cutoff_sigma (float | None, default: 3.0) – Truncate the filter at this many standard deviations. Set to None to disable truncation. Default is 3.0.

Returns:

The filtered input signal. Has the same shape as input_signal.

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

Raises:

ValueError – If input_signal and sigma_array do not have the same length.

See also

scipy.ndimage.gaussian_filter1d: Fixed-sigma Gaussian filtering

Notes

At the boundaries, the outer values are repeated to avoid edge effects (equivalent to mode=’nearest’ in scipy.ndimage.gaussian_filter1d).

Examples

>>> import numpy as np
>>> from gwtransport.diffusion_fast import convolve_diffusion
>>> # Create a sample signal
>>> x = np.linspace(0, 10, 1000)
>>> signal = np.exp(-((x - 3) ** 2)) + 0.5 * np.exp(-((x - 7) ** 2) / 0.5)

>>> # Create position-dependent sigma values
>>> diffusivity = 0.1
>>> dt = 0.001 * (1 + np.sin(2 * np.pi * x / 10))
>>> dx = x[1] - x[0]
>>> sigma_array = np.sqrt(2 * diffusivity * dt) / dx

>>> # Apply the filter
>>> filtered = convolve_diffusion(input_signal=signal, sigma_array=sigma_array)

gwtransport.diffusion_fast.create_example_data(*, nx=1000, domain_length=10.0, diffusivity=0.1, seed=None)[source]#

Create example data for demonstrating variable-sigma diffusion.

Parameters:

nx (int, default: 1000) – Number of spatial points. Default is 1000.
domain_length (float, default: 10.0) – Domain length. Default is 10.0.
diffusivity (float, default: 0.1) – Diffusivity. Default is 0.1.
seed (int | None, default: None) – Random seed for reproducibility. Default is None.

Return type:

tuple[ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]], ndarray[tuple[Any, ...], dtype[floating]]]

Returns:

x (numpy.ndarray) – Spatial coordinates.
signal (numpy.ndarray) – Initial signal (sum of two Gaussians with noise).
sigma_array (numpy.ndarray) – Array of sigma values varying in space.
dt (numpy.ndarray) – Array of time steps varying in space.

diffusion#

Analytical solutions for 1D advection-dispersion transport.

This module implements analytical solutions for solute transport in 1D aquifer systems, combining advection with longitudinal dispersion. The solutions are based on the error function (erf) and its integrals.

Key functions:

infiltration_to_extraction() - Main transport function combining advection and dispersion with explicit pore volume distribution and streamline lengths.
extraction_to_infiltration() - Inverse operation (deconvolution with dispersion).
gamma_infiltration_to_extraction() - Gamma-distributed pore volumes with dispersion. Models aquifer heterogeneity with 2-parameter gamma distribution. Parameterizable via (alpha, beta) or (mean, std). Discretizes gamma distribution into equal-probability bins.
gamma_extraction_to_infiltration() - Gamma-distributed pore volumes, deconvolution with dispersion. Symmetric inverse of gamma_infiltration_to_extraction.

The dispersion is characterized by the longitudinal dispersion coefficient D_L, which is computed internally from:

D_L = D_m + alpha_L * v_s

where: - D_m is the molecular diffusion coefficient [m^2/day] - alpha_L is the longitudinal dispersivity [m] - v_s is the retarded velocity [m/day], computed as v_s = L / (R * tau_mean)

The velocity v_s is computed per (pore_volume, output_bin) using the mean residence time (which includes retardation), correctly accounting for time-varying flow. This formulation ensures that: - Molecular diffusion spreading scales with residence time: sqrt(D_m * tau) - Mechanical dispersion spreading scales with travel distance: sqrt(alpha_L * L)

This module adds microdispersion (alpha_L) and molecular diffusion (D_m) on top of macrodispersion captured by the pore volume distribution (APVD). Both represent velocity heterogeneity at different scales. Microdispersion is an aquifer property; macrodispersion depends additionally on hydrological boundary conditions. See Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity for guidance on when to use each approach and how to avoid double-counting spreading effects.

gwtransport.diffusion.infiltration_to_extraction(*, cin, flow, tedges, cout_tedges, aquifer_pore_volumes, streamline_length, molecular_diffusivity, longitudinal_dispersivity, retardation_factor=1.0, suppress_dispersion_warning=False, asymptotic_cutoff_sigma=3.0)[source]#

Compute extracted concentration with advection and longitudinal dispersion.

This function models 1D solute transport through an aquifer system, combining advective transport (based on residence times) with longitudinal dispersion (diffusive spreading during transport).

The physical model assumes: 1. Water infiltrates with concentration cin at time t_in 2. Water travels distance L through aquifer with residence time tau = V_pore / Q 3. During transport, longitudinal dispersion causes spreading 4. At extraction, the concentration is a diffused breakthrough curve

The longitudinal dispersion coefficient D_L is computed internally as:

D_L = D_m + alpha_L * v_s

where v_s = L / (R * tau_mean) is the retarded velocity computed from the mean residence time (which includes retardation) for each (pore_volume, output_bin) combination. For the dispersivity term, the R factors cancel (alpha_L * v_s * R * tau = alpha_L * L), giving correct distance-dependent mechanical dispersion. See the molecular_diffusivity parameter for how the molecular term interacts with retardation.

Parameters:

cin (TypeAliasType) – Concentration of the compound in infiltrating water [concentration units]. Length must match the number of time bins defined by tedges. The model assumes this value is constant over each interval [tedges[i], tedges[i+1]).
flow (TypeAliasType) – Flow rate of water in the aquifer [m3/day]. Length must match cin and the number of time bins defined by tedges. The model assumes this value is constant over each interval [tedges[i], tedges[i+1]).
tedges (DatetimeIndex) – Time edges defining bins for both cin and flow data. Has length of len(cin) + 1.
cout_tedges (DatetimeIndex) – Time edges for output data bins. Has length of desired output + 1. The output concentration is averaged over each bin.
aquifer_pore_volumes (TypeAliasType) – Array of aquifer pore volumes [m3] representing the distribution of flow paths. Each pore volume determines the residence time for that flow path: tau = V_pore / Q.
streamline_length (TypeAliasType) – Array of travel distances [m] corresponding to each pore volume. Must have the same length as aquifer_pore_volumes.
molecular_diffusivity (TypeAliasType) –
Effective molecular diffusivity [m2/day]. Can be a scalar (same for all pore volumes) or an array with the same length as aquifer_pore_volumes. Must be non-negative. For solute transport, this is the molecular diffusion coefficient D_m [m2/day] — typically ~1e-5 m2/day, negligible compared to mechanical dispersion. For heat transport, pass the thermal diffusivity D_th = lambda / (rho*c)_eff [m2/day], typically 0.01-0.1 m2/day.

Internally, this value enters the dispersion formula as D_L = molecular_diffusivity + alpha_L * v_s, where v_s = L / (R * tau) is the retarded velocity. For the dispersivity term, the R factors cancel (alpha_L * L / (R * tau) * R * tau = alpha_L * L), giving correct distance-dependent spreading. For the molecular term, the spreading scales as molecular_diffusivity * R * tau. This is exact when molecular_diffusivity equals D_m/R — which is negligible for solutes (D_m ~ 1e-5) and correct for heat (where thermal diffusivity already represents D_eff, not D_m * R).
longitudinal_dispersivity (TypeAliasType) – Longitudinal dispersivity [m]. Can be a scalar (same for all pore volumes) or an array with the same length as aquifer_pore_volumes. Must be non-negative. Represents microdispersion (mechanical dispersion from pore-scale velocity variations). Set to 0 for pure molecular diffusion.
retardation_factor (float, default: 1.0) – Retardation factor of the compound in the aquifer (default 1.0). Values > 1.0 indicate slower transport due to sorption.
suppress_dispersion_warning (bool, default: False) – If True, suppress the warning when using multiple pore volumes with non-zero longitudinal_dispersivity. Default is False.
asymptotic_cutoff_sigma (float | None, default: 3.0) – Performance optimization. Cells where the erf argument magnitude exceeds this threshold are assigned asymptotic values (±1) instead of computing the expensive integral. Since erf(3) ≈ 0.99998, the default of 3.0 provides excellent accuracy with significant speedup. Set to None to disable the optimization. Default is 3.0.

Returns:

Bin-averaged concentration in the extracted water. Same units as cin. Length equals len(cout_tedges) - 1. NaN values indicate time periods with no valid contributions from the infiltration data.

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

Raises:

ValueError – If input dimensions are inconsistent, if diffusivity is negative, or if aquifer_pore_volumes and streamline_length have different lengths.

Warns:

UserWarning – If multiple pore volumes are used with non-zero longitudinal_dispersivity. This may lead to double-counting of spreading effects. Suppress with suppress_dispersion_warning=True if this is intentional.

See also

extraction_to_infiltration: Inverse operation (deconvolution)
gwtransport.advection.infiltration_to_extraction: Pure advection (no dispersion)
gwtransport.diffusion_fast.infiltration_to_extraction: Fast Gaussian approximation
Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity: Macrodispersion vs microdispersion

Notes

The algorithm constructs a coefficient matrix W where cout = W @ cin:

For each pore volume, compute the breakthrough curve contribution: - delta_volume: volume between infiltration event and extraction point - step_widths: convert volume to spatial distance x (erf coordinate) - time_active: diffusion time, limited by residence time
For each infiltration time edge, compute the erf response at all extraction time edges using analytical space-time averaging.
Convert erf response to breakthrough fraction: frac = 0.5 * (1 + erf)
Coefficient for bin: coeff[i,j] = frac_start - frac_end This represents the fraction of cin[j] that arrives in cout[i].
Average coefficients across all pore volumes.

The error function solution assumes an initial step function that diffuses over time. The position coordinate x represents the distance from the concentration front to the observation point.

Examples

Basic usage with constant flow:

>>> import pandas as pd
>>> import numpy as np
>>> from gwtransport.diffusion import infiltration_to_extraction
>>>
>>> # Create time edges
>>> tedges = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
>>> cout_tedges = pd.date_range(start="2020-01-05", end="2020-01-25", freq="D")
>>>
>>> # Input concentration (step function) and constant flow
>>> cin = np.zeros(len(tedges) - 1)
>>> cin[5:10] = 1.0  # Pulse of concentration
>>> flow = np.ones(len(tedges) - 1) * 100.0  # 100 m3/day
>>>
>>> # Single pore volume of 500 m3, travel distance 100 m
>>> aquifer_pore_volumes = np.array([500.0])
>>> streamline_length = np.array([100.0])
>>>
>>> # Compute with dispersion (molecular diffusion + dispersivity)
>>> # Scalar values broadcast to all pore volumes
>>> cout = infiltration_to_extraction(
...     cin=cin,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     aquifer_pore_volumes=aquifer_pore_volumes,
...     streamline_length=streamline_length,
...     molecular_diffusivity=1e-4,  # m2/day, same for all pore volumes
...     longitudinal_dispersivity=1.0,  # m, same for all pore volumes
... )

With multiple pore volumes (heterogeneous aquifer):

>>> # Distribution of pore volumes and corresponding travel distances
>>> aquifer_pore_volumes = np.array([400.0, 500.0, 600.0])
>>> streamline_length = np.array([80.0, 100.0, 120.0])
>>>
>>> # Scalar diffusion parameters broadcast to all pore volumes
>>> cout = infiltration_to_extraction(
...     cin=cin,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     aquifer_pore_volumes=aquifer_pore_volumes,
...     streamline_length=streamline_length,
...     molecular_diffusivity=1e-4,  # m2/day
...     longitudinal_dispersivity=1.0,  # m
... )

gwtransport.diffusion.extraction_to_infiltration(*, cout, flow, tedges, cout_tedges, aquifer_pore_volumes, streamline_length, molecular_diffusivity, longitudinal_dispersivity, retardation_factor=1.0, suppress_dispersion_warning=False, asymptotic_cutoff_sigma=3.0, regularization_strength=1e-10)[source]#

Compute infiltration concentration from extracted water (deconvolution with dispersion).

Inverts the forward transport model by building the forward coefficient matrix W_forward from infiltration_to_extraction() and solving W_forward @ cin = cout via Tikhonov regularization. Well-determined modes are dominated by the data; poorly-determined modes are pulled toward the physically motivated target (transpose-and-normalize of the forward matrix).

Parameters:

cout (TypeAliasType) – Concentration of the compound in extracted water [concentration units]. Length must match the number of time bins defined by cout_tedges.
flow (TypeAliasType) – Flow rate of water in the aquifer [m3/day]. Length must match the number of time bins defined by tedges.
tedges (DatetimeIndex) – Time edges defining bins for cin (output) and flow data. Has length of len(flow) + 1. Output cin has length len(tedges) - 1.
cout_tedges (DatetimeIndex) – Time edges for cout data bins. Has length of len(cout) + 1. Can have different time alignment and resolution than tedges.
aquifer_pore_volumes (TypeAliasType) – Array of aquifer pore volumes [m3] representing the distribution of flow paths. Each pore volume determines the residence time for that flow path: tau = V_pore / Q.
streamline_length (TypeAliasType) – Array of travel distances [m] corresponding to each pore volume. Must have the same length as aquifer_pore_volumes.
molecular_diffusivity (TypeAliasType) – Effective molecular diffusivity [m2/day]. Can be a scalar (same for all pore volumes) or an array with the same length as aquifer_pore_volumes. Must be non-negative. See infiltration_to_extraction() for details on the physical interpretation and the interaction with retardation_factor.
longitudinal_dispersivity (TypeAliasType) – Longitudinal dispersivity [m]. Can be a scalar (same for all pore volumes) or an array with the same length as aquifer_pore_volumes. Must be non-negative.
retardation_factor (float, default: 1.0) – Retardation factor of the compound in the aquifer (default 1.0). Values > 1.0 indicate slower transport due to sorption.
suppress_dispersion_warning (bool, default: False) – If True, suppress the warning when using multiple pore volumes with non-zero longitudinal_dispersivity. Default is False.
asymptotic_cutoff_sigma (float | None, default: 3.0) – Performance optimization for the forward matrix construction. Default is 3.0.
regularization_strength (float, default: 1e-10) – Tikhonov regularization parameter λ. See gwtransport.advection.extraction_to_infiltration() for details. Default is 1e-10.

Returns:

Bin-averaged concentration in the infiltrating water. Same units as cout. Length equals len(tedges) - 1. NaN values indicate time periods with no valid contributions from the extraction data.

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

Raises:

ValueError – If input dimensions are inconsistent, if diffusivity is negative, or if aquifer_pore_volumes and streamline_length have different lengths.

Warns:

UserWarning – If multiple pore volumes are used with non-zero longitudinal_dispersivity.

See also

infiltration_to_extraction: Forward operation (convolution)
gwtransport.advection.extraction_to_infiltration: Pure advection (no dispersion)
Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity: Macrodispersion vs microdispersion

Notes

The algorithm builds the forward coefficient matrix W_forward (same as used by infiltration_to_extraction()) and solves W_forward @ cin = cout using gwtransport.utils.solve_tikhonov(). This ensures mathematical consistency between forward and inverse operations.

Examples

Basic usage with constant flow:

>>> import pandas as pd
>>> import numpy as np
>>> from gwtransport.diffusion import extraction_to_infiltration
>>>
>>> # Create time edges: tedges for cin/flow, cout_tedges for cout
>>> tedges = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
>>> cout_tedges = pd.date_range(start="2020-01-05", end="2020-01-25", freq="D")
>>>
>>> # Extracted concentration and constant flow
>>> cout = np.zeros(len(cout_tedges) - 1)
>>> cout[5:10] = 1.0  # Observed pulse at extraction
>>> flow = np.ones(len(tedges) - 1) * 100.0  # 100 m3/day
>>>
>>> # Single pore volume of 500 m3, travel distance 100 m
>>> aquifer_pore_volumes = np.array([500.0])
>>> streamline_length = np.array([100.0])
>>>
>>> # Reconstruct infiltration concentration
>>> cin = extraction_to_infiltration(
...     cout=cout,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     aquifer_pore_volumes=aquifer_pore_volumes,
...     streamline_length=streamline_length,
...     molecular_diffusivity=1e-4,
...     longitudinal_dispersivity=1.0,
... )

gwtransport.diffusion.gamma_infiltration_to_extraction(*, cin, flow, tedges, cout_tedges, mean=None, std=None, loc=0.0, alpha=None, beta=None, n_bins=100, streamline_length, molecular_diffusivity, longitudinal_dispersivity, retardation_factor=1.0, suppress_dispersion_warning=False, asymptotic_cutoff_sigma=3.0)[source]#

Compute extracted concentration with advection-dispersion for gamma-distributed pore volumes.

Combines advective transport (based on gamma-distributed pore volumes) with longitudinal dispersion (diffusive spreading during transport). This is a convenience wrapper around infiltration_to_extraction() that parameterizes the aquifer pore volume distribution as a (shifted) gamma distribution.

Provide either (mean, std) or (alpha, beta); loc is optional and defaults to 0.

Parameters:

cin (TypeAliasType) – Concentration of the compound in infiltrating water.
flow (TypeAliasType) – Flow rate of water in the aquifer [m3/day].
tedges (DatetimeIndex) – Time edges for cin and flow data. Has length len(cin) + 1.
cout_tedges (DatetimeIndex) – Time edges for output data bins. Has length of desired output + 1.
mean (float | None, default: None) – Mean of the gamma distribution of the aquifer pore volume. Must be strictly greater than loc.
std (float | None, default: None) – Standard deviation of the gamma distribution of the aquifer pore volume (invariant under the loc shift).
loc (float, default: 0.0) – Location (minimum pore volume) of the gamma distribution. Must satisfy 0 <= loc < mean. Default is 0.0.
alpha (float | None, default: None) – Shape parameter of gamma distribution of the aquifer pore volume (must be > 0).
beta (float | None, default: None) – Scale parameter of gamma distribution of the aquifer pore volume (must be > 0).
n_bins (int, default: 100) – Number of bins to discretize the gamma distribution. Default is 100.
streamline_length (float) – Travel distance through the aquifer [m]. Applied uniformly to all gamma-discretized pore volumes.
molecular_diffusivity (float) – Effective molecular diffusivity [m2/day]. Must be non-negative. See infiltration_to_extraction() for details on the interaction with retardation_factor.
longitudinal_dispersivity (float) – Longitudinal dispersivity [m]. Must be non-negative.
retardation_factor (float, default: 1.0) – Retardation factor (default 1.0). Values > 1.0 indicate slower transport.
suppress_dispersion_warning (bool, default: False) – Suppress warning about combining multiple pore volumes with dispersivity. Default is False.
asymptotic_cutoff_sigma (float | None, default: 3.0) – Performance optimization. Cells where the erf argument magnitude exceeds this threshold are assigned asymptotic values. Default is 3.0.

Returns:

Bin-averaged concentration in the extracted water. Length equals len(cout_tedges) - 1. NaN values indicate time periods with no valid contributions from the infiltration data.

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

See also

infiltration_to_extraction: Transport with explicit pore volume distribution
gamma_extraction_to_infiltration: Reverse operation (deconvolution)
gwtransport.gamma.bins: Create gamma distribution bins
gwtransport.advection.gamma_infiltration_to_extraction: Pure advection (no dispersion)
Gamma Distribution Model: Two-parameter pore volume model
Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity: Macrodispersion vs microdispersion

Notes

The APVD is only time-invariant under the steady-streamlines assumption (see 2. Steady Streamlines).

The spreading from the gamma-distributed pore volumes represents macrodispersion (aquifer-scale heterogeneity). When std comes from calibration on measurements, it absorbs all mixing: macrodispersion, microdispersion, and an average molecular diffusion contribution. When std comes from streamline analysis, it represents macrodispersion only; microdispersion and molecular diffusion can be added via the dispersion parameters. See Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity for guidance on when to add microdispersion.

Examples

>>> import pandas as pd
>>> import numpy as np
>>> from gwtransport.diffusion import gamma_infiltration_to_extraction
>>>
>>> tedges = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
>>> cout_tedges = pd.date_range(start="2020-01-05", end="2020-01-25", freq="D")
>>> cin = np.zeros(len(tedges) - 1)
>>> cin[5:10] = 1.0
>>> flow = np.ones(len(tedges) - 1) * 100.0
>>>
>>> cout = gamma_infiltration_to_extraction(
...     cin=cin,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     mean=500.0,
...     std=100.0,
...     n_bins=5,
...     streamline_length=100.0,
...     molecular_diffusivity=1e-4,
...     longitudinal_dispersivity=1.0,
... )

gwtransport.diffusion.gamma_extraction_to_infiltration(*, cout, flow, tedges, cout_tedges, mean=None, std=None, loc=0.0, alpha=None, beta=None, n_bins=100, streamline_length, molecular_diffusivity, longitudinal_dispersivity, retardation_factor=1.0, suppress_dispersion_warning=False, asymptotic_cutoff_sigma=3.0, regularization_strength=1e-10)[source]#

Compute infiltration concentration from extracted water for gamma-distributed pore volumes.

Inverts the forward transport model (advection + dispersion with gamma-distributed pore volumes) via Tikhonov regularization. This is a convenience wrapper around extraction_to_infiltration() that parameterizes the aquifer pore volume distribution as a (shifted) gamma distribution.

Provide either (mean, std) or (alpha, beta); loc is optional and defaults to 0.

Parameters:

cout (TypeAliasType) – Concentration of the compound in extracted water.
flow (TypeAliasType) – Flow rate of water in the aquifer [m3/day].
tedges (DatetimeIndex) – Time edges for cin (output) and flow data. Has length of len(flow) + 1.
cout_tedges (DatetimeIndex) – Time edges for cout data bins. Has length of len(cout) + 1.
mean (float | None, default: None) – Mean of the gamma distribution of the aquifer pore volume. Must be strictly greater than loc.
std (float | None, default: None) – Standard deviation of the gamma distribution of the aquifer pore volume (invariant under the loc shift).
loc (float, default: 0.0) – Location (minimum pore volume) of the gamma distribution. Must satisfy 0 <= loc < mean. Default is 0.0.
alpha (float | None, default: None) – Shape parameter of gamma distribution of the aquifer pore volume (must be > 0).
beta (float | None, default: None) – Scale parameter of gamma distribution of the aquifer pore volume (must be > 0).
n_bins (int, default: 100) – Number of bins to discretize the gamma distribution. Default is 100.
streamline_length (float) – Travel distance through the aquifer [m]. Applied uniformly to all gamma-discretized pore volumes.
molecular_diffusivity (float) – Effective molecular diffusivity [m2/day]. Must be non-negative. See infiltration_to_extraction() for details on the interaction with retardation_factor.
longitudinal_dispersivity (float) – Longitudinal dispersivity [m]. Must be non-negative.
retardation_factor (float, default: 1.0) – Retardation factor (default 1.0). Values > 1.0 indicate slower transport.
suppress_dispersion_warning (bool, default: False) – Suppress warning about combining multiple pore volumes with dispersivity. Default is False.
asymptotic_cutoff_sigma (float | None, default: 3.0) – Performance optimization for the forward matrix construction. Default is 3.0.
regularization_strength (float, default: 1e-10) – Tikhonov regularization parameter. Default is 1e-10.

Returns:

Bin-averaged concentration in the infiltrating water. Length equals len(tedges) - 1. NaN values indicate time periods with no valid contributions from the extraction data.

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

See also

extraction_to_infiltration: Deconvolution with explicit pore volume distribution
gamma_infiltration_to_extraction: Forward operation (convolution)
gwtransport.gamma.bins: Create gamma distribution bins
gwtransport.advection.gamma_extraction_to_infiltration: Pure advection (no dispersion)
Gamma Distribution Model: Two-parameter pore volume model
Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity: Macrodispersion vs microdispersion

Notes

The APVD is only time-invariant under the steady-streamlines assumption (see 2. Steady Streamlines).

The spreading from the gamma-distributed pore volumes represents macrodispersion (aquifer-scale heterogeneity). When std comes from calibration on measurements, it absorbs all mixing: macrodispersion, microdispersion, and an average molecular diffusion contribution. When std comes from streamline analysis, it represents macrodispersion only; microdispersion and molecular diffusion can be added via the dispersion parameters. See Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity for guidance on when to add microdispersion.

Examples

>>> import pandas as pd
>>> import numpy as np
>>> from gwtransport.diffusion import gamma_extraction_to_infiltration
>>>
>>> tedges = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
>>> cout_tedges = pd.date_range(start="2020-01-05", end="2020-01-25", freq="D")
>>> cout = np.zeros(len(cout_tedges) - 1)
>>> cout[5:10] = 1.0
>>> flow = np.ones(len(tedges) - 1) * 100.0
>>>
>>> cin = gamma_extraction_to_infiltration(
...     cout=cout,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     mean=500.0,
...     std=100.0,
...     n_bins=5,
...     streamline_length=100.0,
...     molecular_diffusivity=1e-4,
...     longitudinal_dispersivity=1.0,
... )

fronttracking#

Front tracking module for exact nonlinear transport modeling.

fronttracking.events#

Event Detection for Front Tracking.

This module provides exact analytical event detection for the front tracking algorithm. All intersections are computed using closed-form formulas with no numerical iteration or tolerance-based checks.

Events include: - Characteristic-characteristic collisions - Shock-shock collisions - Shock-characteristic collisions - Rarefaction boundary interactions - Outlet crossings

All calculations return exact floating-point results with machine precision.

class gwtransport.fronttracking.events.EventType(*values)[source]#

Bases: Enum

All possible event types in front tracking simulation.

CHAR_CHAR_COLLISION = 'characteristic_collision'#: Two characteristics intersect (will form shock).

SHOCK_SHOCK_COLLISION = 'shock_collision'#: Two shocks collide (will merge).

SHOCK_CHAR_COLLISION = 'shock_characteristic_collision'#: Shock catches or is caught by characteristic.

RAREF_CHAR_COLLISION = 'rarefaction_characteristic_collision'#: Rarefaction boundary intersects with characteristic.

SHOCK_RAREF_COLLISION = 'shock_rarefaction_collision'#: Shock intersects with rarefaction boundary.

RAREF_RAREF_COLLISION = 'rarefaction_rarefaction_collision'#: Rarefaction boundary intersects with another rarefaction boundary.

OUTLET_CROSSING = 'outlet_crossing'#: Wave crosses outlet boundary.

INLET_CHANGE = 'inlet_concentration_change'#: Inlet concentration changes (creates new wave).

FLOW_CHANGE = 'flow_change'#: Flow rate changes (all waves get new velocities).

class gwtransport.fronttracking.events.Event(time, event_type, waves_involved, location, flow_new=None, boundary_type=None)[source]#

Bases: object

Represents a single event in the simulation.

Events are ordered by time for processing in chronological order.

Parameters:

time (float) – Time when event occurs [days]
event_type (EventType) – Type of event
waves_involved (list) – List of wave objects involved in this event
location (float) – Volumetric position where event occurs [m³]
flow_new (float | None, default: None) – New flow rate for FLOW_CHANGE events [m³/day]
boundary_type (str | None, default: None) – Which rarefaction boundary collided: 'head' or 'tail'. Set for rarefaction collision events.

Examples

event = Event(
    time=15.5,
    event_type=EventType.SHOCK_CHAR_COLLISION,
    waves_involved=[shock1, char1],
    location=250.0,
)
print(f"Event at t={event.time}: {event.event_type.value}")

time: float#

event_type: EventType#

waves_involved: list#

location: float#

__init__(time, event_type, waves_involved, location, flow_new=None, boundary_type=None)#

flow_new: float | None = None#

boundary_type: str | None = None#

__lt__(other)[source]#

Events ordered by time (for priority queue).

Returns:: True if this event occurs before other.
Return type:: bool

__repr__()[source]#

Return string representation of Event.

Returns:: Human-readable summary of the event.
Return type:: str

gwtransport.fronttracking.events.find_characteristic_intersection(char1, char2, t_current)[source]#

Find exact analytical intersection of two characteristics.

Solves the linear system:: V1 = v1_start + vel1*(t - t1_start) V2 = v2_start + vel2*(t - t2_start) V1 = V2
This reduces to:: t = (v2_start - v1_start - vel2*t2_start + vel1*t1_start) / (vel1 - vel2)

Parameters:

char1 (CharacteristicWave) – First characteristic
char2 (CharacteristicWave) – Second characteristic
t_current (float) – Current simulation time [days]

Returns:

(t_intersect, v_intersect) if intersection exists in future, None otherwise

Return type:

tuple[float, float] | None

Notes

Returns None if: - Characteristics are parallel (velocities equal within machine precision) - Intersection would occur in the past (t <= t_current) - Either characteristic is not yet active at intersection time

The algorithm uses exact floating-point arithmetic with no tolerance checks except for detecting parallel lines (abs(vel1 - vel2) < 1e-15).

Examples

result = find_characteristic_intersection(char1, char2, t_current=10.0)
if result:
    t_int, v_int = result
    print(f"Intersection at t={t_int:.6f}, V={v_int:.6f}")

gwtransport.fronttracking.events.find_shock_shock_intersection(shock1, shock2, t_current)[source]#

Find exact analytical intersection of two shocks.

Similar to characteristic intersection but uses shock velocities from Rankine-Hugoniot condition.

Parameters:

shock1 (ShockWave) – First shock
shock2 (ShockWave) – Second shock
t_current (float) – Current simulation time [days]

Returns:

(t_intersect, v_intersect) if intersection exists in future, None otherwise

Return type:

tuple[float, float] | None

Notes

Shock velocities are constant (already computed from Rankine-Hugoniot), so this is a simple linear intersection problem.

Examples

result = find_shock_shock_intersection(shock1, shock2, t_current=10.0)
if result:
    t_int, v_int = result
    print(f"Shocks collide at t={t_int:.6f}, V={v_int:.6f}")

gwtransport.fronttracking.events.find_shock_characteristic_intersection(shock, char, t_current)[source]#

Find exact analytical intersection of shock and characteristic.

Parameters:

shock (ShockWave) – Shock wave
char (CharacteristicWave) – Characteristic wave
t_current (float) – Current simulation time [days]

Returns:

(t_intersect, v_intersect) if intersection exists in future, None otherwise

Return type:

tuple[float, float] | None

Examples

result = find_shock_characteristic_intersection(shock, char, t_current=10.0)
if result:
    t_int, v_int = result
    print(f"Shock catches characteristic at t={t_int:.6f}, V={v_int:.6f}")

gwtransport.fronttracking.events.find_rarefaction_boundary_intersections(raref, other_wave, t_current)[source]#

Find intersections of rarefaction head/tail with another wave.

A rarefaction has two boundaries (head and tail), each traveling at characteristic velocities. This function finds intersections of both boundaries with the given wave.

Parameters:

raref (RarefactionWave) – Rarefaction wave
other_wave (Wave) – Any other wave (Characteristic, Shock, or Rarefaction)
t_current (float) – Current simulation time [days]

Returns:

List of (t_intersect, v_intersect, boundary_type) where boundary_type is either ‘head’ or ‘tail’

Return type:

list[tuple[float, float, str]]

Notes

The head travels at velocity corresponding to c_head, and the tail at velocity corresponding to c_tail. Both are treated as characteristics for intersection calculation.

Examples

intersections = find_rarefaction_boundary_intersections(
    raref, char, t_current=10.0
)
for t, v, boundary in intersections:
    print(f"{boundary} intersects at t={t:.3f}, V={v:.3f}")

gwtransport.fronttracking.events.find_outlet_crossing(wave, v_outlet, t_current)[source]#

Find exact analytical time when wave crosses outlet.

For characteristics and shocks, solves:: v_start + velocity*(t - t_start) = v_outlet

For rarefactions, finds when head (leading edge) crosses.

Parameters:

wave (Wave) – Any wave type (Characteristic, Shock, or Rarefaction)
v_outlet (float) – Outlet position [m³]
t_current (float) – Current simulation time [days]

Returns:

Time when wave crosses outlet, or None if: - Wave already past outlet - Wave moving away from outlet - Wave not yet active

Return type:

Notes

This function assumes waves always move in positive V direction (toward outlet). Negative velocities would indicate unphysical backward flow.

Examples

t_cross = find_outlet_crossing(shock, v_outlet=500.0, t_current=10.0)
if t_cross:
    print(f"Shock exits at t={t_cross:.3f} days")

fronttracking.handlers#

Event Handlers for Front Tracking.

This module provides handlers for all wave interaction events in the front tracking algorithm. Each handler receives waves involved in an event and returns new waves created by the interaction.

All handlers enforce physical correctness: - Mass conservation (Rankine-Hugoniot condition) - Entropy conditions (Lax condition for shocks) - Causality (no backward-traveling information)

Handlers modify wave states in-place by deactivating parent waves and creating new child waves.

gwtransport.fronttracking.handlers.handle_characteristic_collision(char1, char2, t_event, v_event)[source]#

Handle collision of two characteristics → create shock or rarefaction.

When two characteristics with different concentrations intersect, they form a shock discontinuity. The faster characteristic (lower concentration for n>1) catches the slower one from behind.

Parameters:

char1 (CharacteristicWave) – First characteristic
char2 (CharacteristicWave) – Second characteristic
t_event (float) – Time of collision [days]
v_event (float) – Position of collision [m³]

Returns:

Single shock or rarefaction wave created at collision point

Return type:

list[ShockWave | RarefactionWave]

Raises:

RuntimeError – If the characteristic collision creates a non-entropic shock, indicating a bug in the collision detection logic.

Notes

The shock has: - c_left: concentration from faster (upstream) characteristic - c_right: concentration from slower (downstream) characteristic - velocity: computed from Rankine-Hugoniot condition

The parent characteristics are deactivated.

Examples

shock = handle_characteristic_collision(char1, char2, t=15.0, v=100.0)
assert shock.satisfies_entropy()
assert not char1.is_active  # Parent deactivated

gwtransport.fronttracking.handlers.handle_shock_collision(shock1, shock2, t_event, v_event)[source]#

Handle collision of two shocks → merge into single shock.

When two shocks collide, they merge into a single shock that connects the left state of the upstream shock to the right state of the downstream shock.

Parameters:

shock1 (ShockWave) – First shock
shock2 (ShockWave) – Second shock
t_event (float) – Time of collision [days]
v_event (float) – Position of collision [m³]

Returns:

Single merged shock wave

Return type:

list[ShockWave]

Raises:

RuntimeError – If shock velocities are not set, or if the merged shock violates the entropy condition.

Notes

The merged shock has: - c_left: from the faster (upstream) shock - c_right: from the slower (downstream) shock - velocity: recomputed from Rankine-Hugoniot

The parent shocks are deactivated.

Examples

merged = handle_shock_collision(shock1, shock2, t=20.0, v=200.0)
assert merged.satisfies_entropy()
assert not shock1.is_active  # Parents deactivated

gwtransport.fronttracking.handlers.handle_shock_characteristic_collision(shock, char, t_event, v_event)[source]#

Handle shock catching or being caught by characteristic.

When the attempted shock would violate entropy (indicating expansion rather than compression), a rarefaction wave is created instead to preserve mass balance. This addresses High Priority #1 from FRONT_TRACKING_REBUILD_PLAN.md.

The outcome depends on which wave is faster: - If shock is faster: shock catches characteristic, absorbs it - If characteristic is faster: characteristic catches shock, modifies it

Parameters:

shock (ShockWave) – Shock wave
char (CharacteristicWave) – Characteristic wave
t_event (float) – Time of collision [days]
v_event (float) – Position of collision [m³]

Returns:

List containing new wave(s): ShockWave if compression, RarefactionWave if expansion, or empty list in edge cases

Return type:

Raises:

RuntimeError – If the shock velocity is not set (i.e., shock.velocity is None).

Notes

The characteristic concentration modifies one side of the shock: - If shock catches char: modifies c_right - If char catches shock: modifies c_left

If the new shock satisfies entropy → return shock (compression) If entropy violated → create rarefaction instead (expansion)

Examples

new_shock = handle_shock_characteristic_collision(shock, char, t=25.0, v=300.0)
if new_shock:
    assert new_shock[0].satisfies_entropy()

gwtransport.fronttracking.handlers.handle_shock_rarefaction_collision(shock, raref, t_event, v_event, boundary_type)[source]#

Handle shock interacting with rarefaction fan with wave splitting.

Implements proper wave splitting for shock-rarefaction interactions, addressing High Priority #2 from FRONT_TRACKING_REBUILD_PLAN.md.

This is the most complex interaction. A shock can:

Catch the rarefaction tail: shock penetrates into rarefaction fan, creating both a modified rarefaction and a continuing shock
Be caught by rarefaction head: creates compression wave

Parameters:

shock (ShockWave) – Shock wave
raref (RarefactionWave) – Rarefaction wave
t_event (float) – Time of collision [days]
v_event (float) – Position of collision [m³]
boundary_type (str | None) – Which boundary collided: ‘head’ or ‘tail’

Returns:

List of new waves created: may include shock and modified rarefaction for tail collision, or compression shock for head collision

Return type:

Raises:

RuntimeError – If the shock velocity is not set (i.e., shock.velocity is None) when processing a head collision.

Notes

Tail collision: Shock penetrates rarefaction, creating: - New shock continuing through rarefaction - Modified rarefaction with compressed tail (if rarefaction not fully overtaken)

Head collision: Rarefaction head catches shock, may create compression shock

Examples

waves = handle_shock_rarefaction_collision(
    shock, raref, t=30.0, v=400.0, boundary_type="tail"
)

gwtransport.fronttracking.handlers.handle_rarefaction_characteristic_collision(raref, char, t_event, v_event, boundary_type)[source]#

Handle rarefaction boundary intersecting with characteristic.

Implements the safe option (b) of the front tracking rebuild plan: when a characteristic intersects either boundary of a rarefaction, the characteristic is absorbed into the rarefaction provided that the characteristic’s concentration matches the boundary concentration to within a tight tolerance. If the concentrations differ by more than the tolerance, deactivating the characteristic would silently destroy mass, so an informative RuntimeError is raised instead.

The matching tolerance is based on the rarefaction’s own concentration range; characteristics that are tangent to (and therefore physically indistinguishable from) the rarefaction boundary pass the check, while truly distinct characteristics are flagged.

Parameters:

raref (RarefactionWave) – Rarefaction wave whose boundary the characteristic crosses.
char (CharacteristicWave) – Characteristic wave that intersects the rarefaction boundary.
t_event (float) – Time of collision [days].
v_event (float) – Position of collision [m^3].
boundary_type (str | None) – Which boundary collided: 'head' or 'tail'.

Returns:

Empty list – the characteristic is absorbed; no new waves are created.

Return type:

Raises:

RuntimeError – If the characteristic’s concentration does not match the colliding rarefaction boundary concentration within tolerance, indicating that absorption would silently violate mass balance.

Notes

Future enhancement: properly split the rarefaction at the collision point and create new wave(s) representing the post-interaction state.

gwtransport.fronttracking.handlers.handle_rarefaction_rarefaction_collision(raref1, raref2, t_event, v_event, boundary_type)[source]#

Handle collision between two rarefaction boundaries.

This handler is intentionally conservative: it records the fact that two rarefaction fans have intersected but does not yet modify the wave topology. Full entropic treatment of rarefaction-rarefaction interactions (potentially involving wave splitting) is reserved for a dedicated future enhancement.

Parameters:

raref1 (RarefactionWave) – First rarefaction wave in the collision.
raref2 (RarefactionWave) – Second rarefaction wave in the collision.
t_event (float) – Time of the boundary intersection [days].
v_event (float) – Position of the intersection [m³].
boundary_type (str | None) – Boundary of the first rarefaction that intersected: ‘head’ or ‘tail’.

Returns:

Empty list; no new waves are created at this stage.

Return type:

Notes

Waves remain active so that concentration queries remain valid.
The FrontTracker records the event in its diagnostics history.
This is consistent with the design goal of exact analytical computation while deferring complex topology changes.

gwtransport.fronttracking.handlers.handle_outlet_crossing(wave, t_event, v_outlet)[source]#

Handle wave crossing outlet boundary.

The wave exits the domain. It remains in the wave list for querying concentration at earlier times but is marked for different handling.

Parameters:

wave (Wave) – Any wave type (Characteristic, Shock, or Rarefaction)
t_event (float) – Time when wave exits [days]
v_outlet (float) – Outlet position [m³]

Returns:

Event record with details about the crossing

Return type:

dict

Notes

Waves are NOT deactivated when they cross the outlet. They remain active for concentration queries at points between their origin and outlet.

The event record includes: - time: crossing time - type: ‘outlet_crossing’ - wave: reference to wave object - concentration_left: upstream concentration - concentration_right: downstream concentration

Examples

event = handle_outlet_crossing(shock, t=50.0, v_outlet=500.0)
print(f"Wave exited at t={event['time']:.2f}")

gwtransport.fronttracking.handlers.recreate_characteristic_with_new_flow(char, t_change, flow_new)[source]#

Create new characteristic at current position with new flow.

When flow changes, existing characteristics must be recreated with updated velocities. The concentration remains constant, but velocity becomes flow_new / R(concentration).

Parameters:

char (CharacteristicWave) – Existing characteristic to recreate
t_change (float) – Time of flow change [days]
flow_new (float) – New flow rate [m³/day]

Returns:

New characteristic at current position with new flow

Return type:

CharacteristicWave

Raises:

ValueError – If the characteristic is not yet active at t_change.

Notes

The parent characteristic should be deactivated by the caller.

Examples

char_old = CharacteristicWave(
    t_start=0.0, v_start=0.0, flow=100.0, concentration=5.0, sorption=sorption
)
char_new = recreate_characteristic_with_new_flow(
    char_old, t_change=10.0, flow_new=200.0
)
assert char_new.flow == 200.0
assert char_new.concentration == 5.0  # Concentration unchanged

gwtransport.fronttracking.handlers.recreate_shock_with_new_flow(shock, t_change, flow_new)[source]#

Create new shock at current position with new flow.

When flow changes, shock velocity must be recomputed using Rankine-Hugoniot condition with the new flow: s = flow*(c_R - c_L) / (C_total(c_R) - C_total(c_L)).

Parameters:

shock (ShockWave) – Existing shock to recreate
t_change (float) – Time of flow change [days]
flow_new (float) – New flow rate [m³/day]

Returns:

New shock at current position with updated velocity

Return type:

ShockWave

Raises:

ValueError – If the shock is not yet active at t_change.

Notes

The parent shock should be deactivated by the caller. Shock velocity is automatically recomputed in ShockWave.__post_init__.

Examples

shock_old = ShockWave(
    t_start=0.0,
    v_start=0.0,
    flow=100.0,
    c_left=10.0,
    c_right=2.0,
    sorption=sorption,
)
shock_new = recreate_shock_with_new_flow(
    shock_old, t_change=10.0, flow_new=200.0
)
assert shock_new.flow == 200.0
assert (
    shock_new.velocity == 2 * shock_old.velocity
)  # Velocity scales linearly with flow

gwtransport.fronttracking.handlers.recreate_rarefaction_with_new_flow(raref, t_change, flow_new)[source]#

Create new rarefaction at current position with new flow.

When flow changes, rarefaction head and tail velocities are updated. The fan structure (c_head, c_tail) is preserved, but the self-similar solution pivots at the flow change point.

Parameters:

raref (RarefactionWave) – Existing rarefaction to recreate
t_change (float) – Time of flow change [days]
flow_new (float) – New flow rate [m³/day]

Returns:

New rarefaction at current position with updated velocities

Return type:

RarefactionWave

Raises:

ValueError – If the rarefaction is not yet active at t_change.

Notes

The parent rarefaction should be deactivated by the caller. The rarefaction fan “pivots” at (v_at_change, t_change).

Before: R(C) = flow_old * (t - t_start_old) / (v - v_start_old) After: R(C) = flow_new * (t - t_change) / (v - v_at_change)

Examples

raref_old = RarefactionWave(
    t_start=0.0,
    v_start=0.0,
    flow=100.0,
    c_head=10.0,
    c_tail=2.0,
    sorption=sorption,
)
raref_new = recreate_rarefaction_with_new_flow(
    raref_old, t_change=10.0, flow_new=200.0
)
assert raref_new.flow == 200.0
assert raref_new.c_head == 10.0  # Concentrations unchanged

gwtransport.fronttracking.handlers.handle_flow_change(t_change, flow_new, active_waves)[source]#

Handle flow change event by recreating all active waves with new flow.

When flow changes, all existing waves must be recreated at their current positions with updated velocities. This maintains exact analytical computation while correctly handling time-varying flow.

Parameters:

t_change (float) – Time of flow change [days]
flow_new (float) – New flow rate [m³/day]
active_waves (list) – All currently active waves

Returns:

New waves created at current positions with new flow

Return type:

Raises:

TypeError – If an unrecognized wave type is encountered in active_waves.

Notes

Parent waves are deactivated by this handler.

Physical interpretation: - Characteristics: velocity changes from flow_old/R(c) to flow_new/R(c) - Shocks: Rankine-Hugoniot velocity recomputed with new flow - Rarefactions: fan pivots at (v_change, t_change)

Examples

new_waves = handle_flow_change(
    t_change=10.0, flow_new=200.0, active_waves=[char1, shock1, raref1]
)
assert len(new_waves) == 3
assert all(w.flow == 200.0 for w in new_waves)

gwtransport.fronttracking.handlers.create_inlet_waves_at_time(c_prev, c_new, t, flow, sorption, v_inlet=0.0)[source]#

Create appropriate waves when inlet concentration changes.

Analyzes the concentration change and creates the physically correct wave type based on characteristic velocities.

Parameters:

c_prev (float) – Previous concentration [mass/volume]
c_new (float) – New concentration [mass/volume]
t (float) – Time of concentration change [days]
flow (float) – Flow rate [m³/day]
sorption (NonlinearSorption | ConstantRetardation) – Sorption parameters
v_inlet (float, default: 0.0) – Inlet position [m³], default 0.0

Returns:

List of newly created waves (typically 1 wave per concentration change)

Return type:

Notes

Wave type logic: - vel_new > vel_prev: Compression → create ShockWave - vel_new < vel_prev: Expansion → create RarefactionWave - vel_new == vel_prev: Contact discontinuity → create CharacteristicWave

For shocks, verifies entropy condition before creation.

Examples

# Step increase from zero creates characteristic
waves = create_inlet_waves_at_time(
    c_prev=0.0, c_new=10.0, t=10.0, flow=100.0, sorption=sorption
)
assert isinstance(waves[0], CharacteristicWave)
# Step between nonzero values creates shock for n>1 (compression)
waves = create_inlet_waves_at_time(
    c_prev=2.0, c_new=10.0, t=10.0, flow=100.0, sorption=sorption
)
assert isinstance(waves[0], ShockWave)

fronttracking.math#

Mathematical Foundation for Front Tracking with Nonlinear Sorption.

This module provides exact analytical computations for: - Freundlich, Langmuir, and constant retardation models - Shock velocities via Rankine-Hugoniot condition - Characteristic velocities and positions - First arrival time calculations - Entropy condition verification

All computations are exact analytical formulas with no numerical tolerances.

This file is part of gwtransport which is released under AGPL-3.0 license. See the ./LICENSE file or go to gwtransport/gwtransport for full license details.

class gwtransport.fronttracking.math.NonlinearSorption[source]#

Bases: ABC

Abstract base for concentration-dependent sorption models.

Subclasses must implement retardation, total_concentration, and concentration_from_retardation. Shock velocity and entropy checking are provided generically via the Rankine-Hugoniot and Lax conditions.

See also

FreundlichSorption: Freundlich isotherm implementation.
LangmuirSorption: Langmuir isotherm implementation.
ConstantRetardation: Linear (constant R) retardation model.

abstractmethod retardation(c)[source]#

Compute retardation factor R(C).

Return type:: float | ndarray[tuple[Any, ...], dtype[double]]

abstractmethod total_concentration(c)[source]#

Compute total concentration (dissolved + sorbed per unit pore volume).

Return type:: float | ndarray[tuple[Any, ...], dtype[double]]

abstractmethod concentration_from_retardation(r)[source]#

Invert retardation factor to obtain concentration.

Return type:: float | ndarray[tuple[Any, ...], dtype[double]]

shock_velocity(c_left, c_right, flow)[source]#

Compute shock velocity via Rankine-Hugoniot condition.

The Rankine-Hugoniot condition ensures mass conservation across the shock:: s_shock = [flux(C_R) - flux(C_L)] / [C_total(C_R) - C_total(C_L)]

where flux(C) = flow * C (only dissolved species are transported).

Parameters:

c_left (float) – Concentration upstream (behind) shock [mass/volume].
c_right (float) – Concentration downstream (ahead of) shock [mass/volume].
flow (float) – Flow rate [volume/time]. Must be positive.

Returns:

s_shock – Shock velocity [volume/time].

Return type:

Notes

The Rankine-Hugoniot condition is derived from integrating the conservation law across the shock discontinuity. It ensures that the total mass flux (advective transport) is conserved.

check_entropy_condition(c_left, c_right, shock_vel, flow)[source]#

Verify Lax entropy condition for physical admissibility of shock.

The Lax entropy condition ensures that characteristics flow INTO the shock from both sides, which is required for physical shocks:

λ(C_L) > s_shock > λ(C_R)

where λ(C) = flow / R(C) is the characteristic velocity.

Parameters:

c_left (float) – Concentration upstream of shock [mass/volume].
c_right (float) – Concentration downstream of shock [mass/volume].
shock_vel (float) – Shock velocity [volume/time].
flow (float) – Flow rate [volume/time].

Returns:

satisfies – True if shock satisfies entropy condition (is physical).

Return type:

bool

class gwtransport.fronttracking.math.FreundlichSorption(k_f, n, bulk_density, porosity, c_min=1e-12)[source]#

Bases: NonlinearSorption

Freundlich sorption isotherm with exact analytical methods.

The Freundlich isotherm is: s(C) = k_f * C^(1/n)

where: - s is sorbed concentration [mass/mass of solid] - C is dissolved concentration [mass/volume of water] - k_f is Freundlich coefficient [(volume/mass)^(1/n)] - n is Freundlich exponent (dimensionless)

For n > 1: Higher C travels faster For n < 1: Higher C travels slower For n = 1: linear (not supported, use ConstantRetardation instead)

Parameters:

k_f (float) – Freundlich coefficient [(m³/kg)^(1/n)]. Must be positive.
n (float) – Freundlich exponent [-]. Must be positive and != 1.
bulk_density (float) – Bulk density of porous medium [kg/m³]. Must be positive.
porosity (float) – Porosity [-]. Must be in (0, 1).
c_min (float, default: 1e-12) – Minimum concentration threshold. For n>1, prevents infinite retardation as C→0. Default: 0.1 for n>1, 0.0 for n<1 (set automatically if not provided).

Notes

The retardation factor is defined as:

R(C) = 1 + (rho_b/n_por) * ds/dC: = 1 + (rho_b*k_f)/(n_por*n) * C^((1/n)-1)

For Freundlich sorption, R depends on C, which creates nonlinear wave behavior.

For n>1 (higher C travels faster), R(C)→∞ as C→0, which can cause extremely slow wave propagation. The c_min parameter prevents this by enforcing a minimum concentration, making R(C) finite for all C≥0.

Examples

>>> sorption = FreundlichSorption(
...     k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
... )
>>> r = sorption.retardation(5.0)
>>> c_back = sorption.concentration_from_retardation(r)
>>> bool(np.isclose(c_back, 5.0))
True

k_f: float#: Freundlich coefficient [(m³/kg)^(1/n)].

n: float#: Freundlich exponent [-].

bulk_density: float#: Bulk density of porous medium [kg/m³].

porosity: float#: Porosity [-].

c_min: float = 1e-12#: Minimum concentration threshold to prevent infinite retardation.

__post_init__()[source]#

Validate parameters after initialization.

Raises:: ValueError – If any parameter is outside its valid range: k_f <= 0, n <= 0, n == 1, bulk_density <= 0, porosity outside (0, 1), or c_min < 0.

retardation(c)[source]#

Compute retardation factor R(C).

The retardation factor relates concentration velocity to pore water velocity:: v_C = flow / R(C)
For Freundlich sorption:: R(C) = 1 + (rho_b*k_f)/(n_por*n) * C^((1/n)-1)

Parameters:: c (float | ndarray[tuple[Any, ...], dtype[double]]) – Dissolved concentration [mass/volume]. Non-negative.
Returns:: r – Retardation factor [-]. Always >= 1.0.
Return type:: float | ndarray[tuple[Any, ...], dtype[double]]

Notes

For n > 1: R decreases with increasing C (higher C travels faster)
For n < 1: R increases with increasing C (higher C travels slower)
Concentrations at or below c_min return R=1 if c_min=0, else are clipped to c_min

total_concentration(c)[source]#

Compute total concentration (dissolved + sorbed per unit pore volume).

Total concentration includes both dissolved and sorbed mass:

C_total = C + (rho_b/n_por) * s(C): = C + (rho_b/n_por) * k_f * C^(1/n)

Parameters:: c (float | ndarray[tuple[Any, ...], dtype[double]]) – Dissolved concentration [mass/volume]. Non-negative.
Returns:: c_total – Total concentration [mass/volume]. Always >= c.
Return type:: float | ndarray[tuple[Any, ...], dtype[double]]

Notes

This is the conserved quantity in the transport equation:: ∂C_total/∂t + ∂(flow*C)/∂v = 0

The flux term only includes dissolved concentration because sorbed mass is immobile. Concentrations at or below c_min return C if c_min=0, else use c_min.

concentration_from_retardation(r)[source]#

Invert retardation factor to obtain concentration analytically.

Given R, solves R = retardation(C) for C. This is used in rarefaction waves where the self-similar solution gives R as a function of position and time.

Parameters:: r (float | ndarray[tuple[Any, ...], dtype[double]]) – Retardation factor [-]. Must be >= 1.0.
Returns:: c – Dissolved concentration [mass/volume]. Non-negative.
Return type:: float | ndarray[tuple[Any, ...], dtype[double]]
Raises:: ValueError – If the retardation exponent is effectively zero (i.e., n ≈ 1), making inversion undefined.

Notes

This inverts the relation:: R = 1 + (rho_b*k_f)/(n_por*n) * C^((1/n)-1)
The analytical solution is:: C = [(R-1) * n_por*n / (rho_b*k_f)]^(n/(1-n))

For n = 1 (linear sorption), the exponent n/(1-n) is undefined, which is why linear sorption must use ConstantRetardation class instead.

Examples

>>> sorption = FreundlichSorption(
...     k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
... )
>>> r = sorption.retardation(5.0)
>>> c = sorption.concentration_from_retardation(r)
>>> bool(np.isclose(c, 5.0, rtol=1e-14))
True

__init__(k_f, n, bulk_density, porosity, c_min=1e-12)#

class gwtransport.fronttracking.math.ConstantRetardation(retardation_factor)[source]#

Bases: object

Constant (linear) retardation model.

For linear sorption: s(C) = K_d * C This gives constant retardation: R(C) = 1 + (rho_b/n_por) * K_d = constant

This is a special case where concentration-dependent behavior disappears. Used for conservative tracers or as approximation for weak sorption.

Parameters:: retardation_factor (float) – Constant retardation factor [-]. Must be >= 1.0. R = 1.0 means no retardation (conservative tracer).

Notes

With constant retardation: - All concentrations travel at same velocity: flow / R - No rarefaction waves form (all concentrations travel together) - Shocks occur only at concentration discontinuities at inlet - Solution reduces to simple time-shifting

This is equivalent to using infiltration_to_extraction_series in the gwtransport package.

Examples

>>> sorption = ConstantRetardation(retardation_factor=2.0)
>>> sorption.retardation(5.0)
2.0
>>> sorption.retardation(10.0)
2.0

retardation_factor: float#: Constant retardation factor [-]. Must be >= 1.0.

__post_init__()[source]#

Validate parameters after initialization.

Raises:: ValueError – If retardation_factor is less than 1.0.

retardation(c)[source]#

Return constant retardation factor (independent of concentration).

Parameters:: c (float) – Dissolved concentration (not used for constant retardation).
Returns:: r – Constant retardation factor.
Return type:: float

total_concentration(c)[source]#

Compute total concentration for linear sorption.

For constant retardation:: C_total = C * R

Parameters:: c (float) – Dissolved concentration [mass/volume].
Returns:: c_total – Total concentration [mass/volume].
Return type:: float

concentration_from_retardation(r)[source]#

Not applicable for constant retardation.

With constant R, all concentrations have the same retardation, so inversion is not meaningful. This method raises an error.

Raises:: NotImplementedError – Always raised for constant retardation.
Return type:: float

shock_velocity(c_left, c_right, flow)[source]#

Compute shock velocity for constant retardation.

With constant R, the shock velocity simplifies to:: s_shock = flow / R

This is the same as all characteristic velocities.

Parameters:

c_left (float) – Concentration upstream of shock (not used for constant R).
c_right (float) – Concentration downstream of shock (not used for constant R).
flow (float) – Flow rate [volume/time].

Returns:

s_shock – Shock velocity [volume/time].

Return type:

check_entropy_condition(c_left, c_right, shock_vel, flow)[source]#

Check entropy condition for constant retardation.

With constant R, all characteristic velocities are equal, so the entropy condition is trivially satisfied for any shock (or rather, shocks don’t really exist - they’re just contact discontinuities).

Parameters:

c_left (float) – Concentration upstream.
c_right (float) – Concentration downstream.
shock_vel (float) – Shock velocity.
flow (float) – Flow rate.

Returns:

satisfies – Always True for constant retardation.

Return type:

bool

__init__(retardation_factor)#

class gwtransport.fronttracking.math.LangmuirSorption(s_max, k_l, bulk_density, porosity)[source]#

Bases: NonlinearSorption

Langmuir sorption isotherm with exact analytical methods.

The Langmuir isotherm is: s(C) = s_max * C / (K_L + C)

where: - s is sorbed concentration [mass/mass of solid] - C is dissolved concentration [mass/volume of water] - s_max is maximum sorption capacity [mass/mass of solid] - K_L is half-saturation constant [mass/volume]

Retardation always decreases with C (favorable isotherm), and R(0) is finite — unlike Freundlich with n > 1, no minimum concentration threshold is needed.

Parameters:

s_max (float) – Maximum sorption capacity [mass/mass of solid]. Must be positive.
k_l (float) – Half-saturation constant [mass/volume]. Concentration at which s = s_max / 2. Must be positive.
bulk_density (float) – Bulk density of porous medium [kg/m³]. Must be positive.
porosity (float) – Porosity [-]. Must be in (0, 1).

See also

FreundlichSorption: Freundlich isotherm (unbounded sorption).
ConstantRetardation: Linear (constant R) retardation model.
Non-Linear Sorption: Exact Solutions: Background on nonlinear sorption.

Notes

The retardation factor is defined as:: R(C) = 1 + (rho_b * s_max * K_L) / (n_por * (K_L + C)^2)

Key properties:

R(0) = 1 + rho_b * s_max / (n_por * K_L) – finite for all parameters
R -> 1 as C -> infinity (all sorption sites saturated)
R always decreases with increasing C (higher C travels faster)
Shocks form on concentration increases, rarefaction fans on decreases

Examples

>>> sorption = LangmuirSorption(
...     s_max=0.1, k_l=5.0, bulk_density=1500.0, porosity=0.3
... )
>>> r = sorption.retardation(5.0)
>>> c_back = sorption.concentration_from_retardation(r)
>>> bool(np.isclose(c_back, 5.0))
True

s_max: float#: Maximum sorption capacity [mass/mass of solid].

k_l: float#: Half-saturation constant [mass/volume].

bulk_density: float#: Bulk density of porous medium [kg/m³].

__init__(s_max, k_l, bulk_density, porosity)#

porosity: float#: Porosity [-].

__post_init__()[source]#

Validate parameters after initialization.

Raises:: ValueError – If any parameter is outside its valid range: s_max <= 0, k_l <= 0, bulk_density <= 0, or porosity outside (0, 1).

retardation(c)[source]#

Compute retardation factor R(C).

For Langmuir sorption:: R(C) = 1 + A / (K_L + C)²

where A = rho_b * s_max * K_L / n_por.

Parameters:: c (float | ndarray[tuple[Any, ...], dtype[double]]) – Dissolved concentration [mass/volume]. Non-negative.
Returns:: r – Retardation factor [-]. Always >= 1.0.
Return type:: float | ndarray[tuple[Any, ...], dtype[double]]

Notes

R(0) = 1 + rho_b * s_max / (n_por * K_L) — always finite
R decreases with increasing C (higher C travels faster)
R → 1 as C → ∞ (all sorption sites saturated)

total_concentration(c)[source]#

Compute total concentration (dissolved + sorbed per unit pore volume).

For Langmuir sorption:: C_total = C + (rho_b / n_por) * s_max * C / (K_L + C)

Parameters:: c (float | ndarray[tuple[Any, ...], dtype[double]]) – Dissolved concentration [mass/volume]. Non-negative.
Returns:: c_total – Total concentration [mass/volume]. Always >= c.
Return type:: float | ndarray[tuple[Any, ...], dtype[double]]

concentration_from_retardation(r)[source]#

Invert retardation factor to obtain concentration analytically.

Given R, solves R = 1 + A / (K_L + C)² for C:: C = sqrt(A / (R - 1)) - K_L

Parameters:: r (float | ndarray[tuple[Any, ...], dtype[double]]) – Retardation factor [-]. Must be >= 1.0.
Returns:: c – Dissolved concentration [mass/volume]. Non-negative.
Return type:: float | ndarray[tuple[Any, ...], dtype[double]]

Notes

For R <= 1, returns 0.0 (unphysical region). For R >= R(0) = 1 + A/K_L², returns 0.0 (at or below zero concentration).

Examples

>>> sorption = LangmuirSorption(
...     s_max=0.1, k_l=5.0, bulk_density=1500.0, porosity=0.3
... )
>>> r = sorption.retardation(5.0)
>>> c = sorption.concentration_from_retardation(r)
>>> bool(np.isclose(c, 5.0, rtol=1e-14))
True

gwtransport.fronttracking.math.SorptionModel = gwtransport.fronttracking.math.NonlinearSorption | gwtransport.fronttracking.math.ConstantRetardation#: Type alias for all sorption models accepted by the front-tracking solver.

gwtransport.fronttracking.math.characteristic_velocity(c, flow, sorption)[source]#

Compute characteristic velocity for given concentration.

In smooth regions of the solution, concentration travels at velocity:: v = flow / R(C)

Parameters:

c (float) – Dissolved concentration [mass/volume].
flow (float) – Flow rate [volume/time].
sorption (NonlinearSorption | ConstantRetardation) – Sorption model.

Returns:

velocity – Characteristic velocity [volume/time].

Return type:

Examples

>>> sorption = FreundlichSorption(
...     k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
... )
>>> v = characteristic_velocity(c=5.0, flow=100.0, sorption=sorption)
>>> v > 0
True

gwtransport.fronttracking.math.characteristic_position(c, flow, sorption, t_start, v_start, t)[source]#

Compute exact position of characteristic at time t.

Characteristics propagate linearly in time:: V(t) = v_start + velocity * (t - t_start)

where velocity = flow / R(C) is constant along the characteristic.

Parameters:

c (float) – Concentration carried by characteristic [mass/volume].
flow (float) – Flow rate [volume/time].
sorption (NonlinearSorption | ConstantRetardation) – Sorption model.
t_start (float) – Time when characteristic starts [days].
v_start (float) – Starting position [volume].
t (float) – Time at which to compute position [days].

Returns:

position – Position at time t [volume], or None if t < t_start.

Return type:

Examples

>>> sorption = ConstantRetardation(retardation_factor=2.0)
>>> v = characteristic_position(
...     c=5.0, flow=100.0, sorption=sorption, t_start=0.0, v_start=0.0, t=10.0
... )
>>> bool(np.isclose(v, 500.0))  # v = 100/2 * 10 = 500
True

gwtransport.fronttracking.math.compute_first_front_arrival_time(cin, flow, tedges, aquifer_pore_volume, sorption)[source]#

Compute exact time when first wave reaches outlet (v_max).

This function returns the precise moment when the first non-zero concentration wave from the inlet arrives at the outlet. This marks the end of the spin-up period.

For the typical case where the first inlet change creates a characteristic (e.g., 0→C transition), this is when that characteristic reaches v_max.

For cases with rarefaction waves:

n>1 (higher C travels faster): The head of a rarefaction (higher C) arrives first.
n<1 (lower C travels faster): The head of a rarefaction (lower C) arrives first.

Algorithm: 1. Find first index where cin > 0 2. Compute residence time for this concentration from inlet to outlet 3. Account for piecewise constant flow during transit 4. Return arrival time in days from tedges[0]

Parameters:

cin (ndarray[tuple[Any, ...], dtype[floating]]) – Inlet concentration [mass/volume]. Length = len(tedges) - 1.
flow (ndarray[tuple[Any, ...], dtype[floating]]) – Flow rate [volume/time]. Length = len(tedges) - 1.
tedges (DatetimeIndex) – Time bin edges. Length = len(cin) + 1. Expected to be DatetimeIndex.
aquifer_pore_volume (float) – Total pore volume [volume]. Must be positive.
sorption (NonlinearSorption | ConstantRetardation) – Sorption model.

Returns:

t_first_arrival – Time when first wave reaches outlet, measured in days from tedges[0]. Returns np.inf if no concentration ever arrives.

Return type:

See also

compute_first_fully_informed_bin_edge: Get first valid output bin edge

Notes

The residence time accounts for retardation:: residence_time = aquifer_pore_volume * R(C) / flow_avg
For piecewise constant flow, we integrate:: ∫₀^residence_time flow(t) dt = aquifer_pore_volume * R(C)

This function computes the EXACT crossing time in days, not a bin edge. Use compute_first_fully_informed_bin_edge() to get the corresponding output bin edge for masking purposes.

Examples

>>> import pandas as pd
>>> cin = np.array([0.0, 10.0] + [10.0] * 10)  # First bin zero, then nonzero
>>> flow = np.array([100.0] * 12)  # Constant flow
>>> tedges = pd.date_range("2020-01-01", periods=13, freq="D")
>>> sorption = ConstantRetardation(retardation_factor=2.0)
>>> t_first = compute_first_front_arrival_time(cin, flow, tedges, 500.0, sorption)
>>> # Result is in days from tedges[0]
>>> bool(np.isclose(t_first, 11.0))  # 1 day (offset) + 10 days (travel time)
True

gwtransport.fronttracking.math.compute_first_fully_informed_bin_edge(cin, flow, tedges, aquifer_pore_volume, sorption, output_tedges)[source]#

Compute left edge of first output bin that is fully informed.

A bin [t_i, t_{i+1}] is “fully informed” if all water exiting during that interval originated from known inlet conditions (not unknown initial conditions). This occurs when t_i >= first front arrival time.

This function is useful for: - Masking output bins during spin-up period - Determining which output times are valid for analysis - Plotting valid vs spin-up regions

Rarefaction handling:

For n>1: Rarefaction tail (lower C, slower) arrives after head. Once the first wave (head) arrives, subsequent bins are informed.
For n<1: Rarefaction tail (higher C, slower) arrives after head. Same principle applies.

In both cases, once the leading edge of the inlet-generated wave structure reaches the outlet, all subsequent output is determined by inlet history.

Parameters:

cin (ndarray[tuple[Any, ...], dtype[floating]]) – Inlet concentration [mass/volume]. Length = len(tedges) - 1.
flow (ndarray[tuple[Any, ...], dtype[floating]]) – Flow rate [volume/time]. Length = len(tedges) - 1.
tedges (DatetimeIndex) – Inlet time bin edges. Length = len(cin) + 1. Expected to be DatetimeIndex.
aquifer_pore_volume (float) – Total pore volume [volume]. Must be positive.
sorption (NonlinearSorption | ConstantRetardation) – Sorption model.
output_tedges (DatetimeIndex) – Output time bin edges. These are the bins for which we want to determine the first fully-informed bin. Expected to be DatetimeIndex.

Returns:

t_first_bin – Left edge of first output bin that is fully informed, measured in days from output_tedges[0]. Returns last edge in days if no bin is fully informed. Returns np.inf if output_tedges is empty.

Return type:

See also

compute_first_front_arrival_time: Get exact arrival time

Notes

This differs from compute_first_front_arrival_time in that it returns a BIN EDGE (from output_tedges), not the exact crossing time.

Both functions return time in days, but measured from different reference points: - compute_first_front_arrival_time: days from tedges[0] - compute_first_fully_informed_bin_edge: days from output_tedges[0]

For masking output arrays:

import pandas as pd

t_first_bin = compute_first_fully_informed_bin_edge(...)
# Convert output_tedges to days from output_tedges[0]
tedges_days = (output_tedges - output_tedges[0]) / pd.Timedelta(days=1)
mask = tedges_days[:-1] >= t_first_bin
cout_valid = cout[mask]

Examples

>>> import pandas as pd
>>> # Exact arrival at ~11 days from tedges[0]
>>> cin = np.array([0.0, 10.0, 10.0])
>>> flow = np.array([100.0, 100.0, 100.0])
>>> tedges = pd.date_range("2020-01-01", periods=4, freq="D")
>>> output_tedges = pd.date_range("2020-01-01", periods=20, freq="D")
>>> sorption = ConstantRetardation(retardation_factor=2.0)
>>> t_bin = compute_first_fully_informed_bin_edge(
...     cin, flow, tedges, 500.0, sorption, output_tedges
... )
>>> # Result is in days from output_tedges[0]
>>> t_bin >= 11.0  # First bin edge >= arrival time
True

fronttracking.output#

Concentration Extraction from Front Tracking Solutions.

This module computes outlet concentrations from front tracking wave solutions using exact analytical integration. All calculations maintain machine precision with no numerical dispersion.

Functions#

concentration_at_point(v, t, waves, sorption): Compute concentration at any point (v, t) in domain
compute_breakthrough_curve(t_array, v_outlet, waves, sorption): Compute concentration at outlet over time array
compute_bin_averaged_concentration_exact(t_edges, v_outlet, waves, sorption): Compute bin-averaged concentrations using exact analytical integration
compute_domain_mass(t, v_outlet, waves, sorption): Compute total mass in domain [0, v_outlet] at time t using exact analytical integration
compute_cumulative_inlet_mass(t, cin, flow, tedges_days): Compute cumulative mass entering domain from t=0 to t
compute_cumulative_outlet_mass(t, v_outlet, waves, sorption, flow, tedges_days): Compute cumulative mass exiting domain from t=0 to t
find_last_rarefaction_start_time(v_outlet, waves): Find time when last rarefaction head reaches outlet
integrate_rarefaction_total_mass(raref, v_outlet, t_start, sorption, flow): Compute total mass exiting through rarefaction to infinity
compute_total_outlet_mass(v_outlet, waves, sorption, flow, tedges_days): Compute total integrated outlet mass until all mass has exited

This file is part of gwtransport which is released under AGPL-3.0 license. See the ./LICENSE file or go to gwtransport/gwtransport for full license details.

gwtransport.fronttracking.output.concentration_at_point(v, t, waves, sorption)[source]#

Compute concentration at point (v, t) with exact analytical value.

Searches through all waves to find which wave controls the concentration at the given point in space-time. Uses exact analytical solutions for characteristics, shocks, and rarefaction fans.

Parameters:

v (float) – Position [m³].
t (float) – Time [days].
waves (Sequence[Wave]) – All waves in the simulation (active and inactive).
sorption (NonlinearSorption | ConstantRetardation) – Sorption model.

Returns:

concentration – Concentration at point (v, t) [mass/volume].

Return type:

Notes

Wave Priority: The algorithm checks waves in this order: 1. Rarefaction waves (if point is inside rarefaction fan) 2. Shocks (discontinuities) 3. Rarefaction tails (if rarefaction tail has passed point) 4. Characteristics (smooth regions)

If no active wave controls the point, returns 0.0 (initial condition).

Rarefaction Tail Behavior: After a rarefaction tail passes a query point, that point maintains the tail concentration as a plateau. This ensures proper concentration behavior after rarefaction waves pass through.

Machine Precision: All position and concentration calculations use exact analytical formulas. Numerical tolerance is only used for equality checks (v == v_shock).

Examples

sorption = FreundlichSorption(
    k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
)
# After running simulation with waves...
c = concentration_at_point(v=250.0, t=5.0, waves=all_waves, sorption=sorption)
c >= 0.0

gwtransport.fronttracking.output.compute_breakthrough_curve(t_array, v_outlet, waves, sorption)[source]#

Compute concentration at outlet over time array.

This is the breakthrough curve: C(v_outlet, t) for all t in t_array.

Parameters:

t_array (ndarray[tuple[Any, ...], dtype[floating]]) – Time points [days]. Must be sorted in ascending order.
v_outlet (float) – Outlet position [m³].
waves (Sequence[Wave]) – All waves in the simulation.
sorption (NonlinearSorption | ConstantRetardation) – Sorption model.

Returns:

c_out – Array of concentrations matching t_array [mass/volume].

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

See also

concentration_at_point: Point-wise concentration
compute_bin_averaged_concentration_exact: Bin-averaged concentrations

Examples

t_array = np.linspace(0, 100, 1000)
c_out = compute_breakthrough_curve(
    t_array, v_outlet=500.0, waves=all_waves, sorption=sorption
)
len(c_out) == len(t_array)

gwtransport.fronttracking.output.identify_outlet_segments(t_start, t_end, v_outlet, waves, sorption)[source]#

Identify which waves control outlet concentration in time interval [t_start, t_end].

Finds all wave crossing events at the outlet and constructs segments where concentration is constant or varying (rarefaction).

Parameters:

t_start (float) – Start of time interval [days].
t_end (float) – End of time interval [days].
v_outlet (float) – Outlet position [m³].
waves (Sequence[Wave]) – All waves in the simulation.
sorption (NonlinearSorption | ConstantRetardation) – Sorption model.

Returns:

segments – List of segment dictionaries, each containing:

’t_start’float
Segment start time
’t_end’float
Segment end time
’type’str
’constant’ or ‘rarefaction’
’concentration’float
For constant segments
’wave’Wave
For rarefaction segments
’c_start’float
Concentration at segment start
’c_end’float
Concentration at segment end

Return type:

list[dict]

Notes

Segments are constructed by: 1. Finding all wave crossing events at outlet in [t_start, t_end] 2. Sorting events chronologically 3. Creating constant-concentration segments between events 4. Handling rarefaction fans with time-varying concentration

The segments completely partition the time interval [t_start, t_end].

gwtransport.fronttracking.output.integrate_rarefaction_exact(raref, v_outlet, t_start, t_end, sorption)[source]#

Exact analytical integral of rarefaction concentration over time at fixed position.

Computes integral of C(t) dt from t_start to t_end where C(t) is the self-similar rarefaction solution at the outlet. Dispatches to isotherm-specific closed-form formulas.

Parameters:

raref (RarefactionWave) – Rarefaction wave controlling the outlet.
v_outlet (float) – Outlet position [m3].
t_start (float) – Integration start time [days]. Can be -np.inf.
t_end (float) – Integration end time [days]. Can be np.inf.
sorption (NonlinearSorption | ConstantRetardation) – Sorption model (FreundlichSorption or LangmuirSorption).

Returns:

integral – Exact integral value [concentration * time].

Return type:

Raises:

TypeError – If sorption model does not support exact rarefaction integration.

gwtransport.fronttracking.output.compute_bin_averaged_concentration_exact(t_edges, v_outlet, waves, sorption)[source]#

Compute bin-averaged concentration using EXACT analytical integration.

For each time bin [t_i, t_{i+1}], computes:: C_avg = (1/(t_{i+1} - t_i)) * ∫_{t_i}^{t_{i+1}} C(v_outlet, t) dt

This is the critical function for maintaining machine precision in output. All integrations use exact analytical formulas with no numerical quadrature.

Parameters:

t_edges (ndarray[tuple[Any, ...], dtype[floating]]) – Time bin edges [days]. Length N+1 for N bins.
v_outlet (float) – Outlet position [m³].
waves (Sequence[Wave]) – All waves from front tracking simulation.
sorption (NonlinearSorption | ConstantRetardation) – Sorption model.

Returns:

c_avg – Bin-averaged concentrations [mass/volume]. Length N.

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

Raises:

ValueError – If any time bin has non-positive width (t_edges[i] >= t_edges[i+1]), or if an unknown segment type is encountered during integration.

See also

concentration_at_point: Point-wise concentration
compute_breakthrough_curve: Breakthrough curve
integrate_rarefaction_exact: Exact rarefaction integration

Notes

Algorithm:

For each bin [t_i, t_{i+1}]:
1. Identify which wave segments control outlet during this period
2. For each segment, compute: Constant C gives integral = C * Δt, Rarefaction C(t) uses exact analytical integral formula
3. Sum segment integrals and divide by bin width

Machine Precision:

Constant segments: exact to floating-point precision
Rarefaction segments: uses closed-form antiderivative
No numerical quadrature or interpolation
Maintains mass balance to ~1e-14 relative error

Rarefaction Integration:

For Freundlich sorption, rarefaction concentration at outlet varies as:

C(t) = [(kappa*t + mu - 1)/alpha]^(1/beta)

The exact integral is:

∫ C dt = (1/(alpha^(1/beta)*kappa*exponent)) * (kappa*t + mu - 1)^exponent

where exponent = 1/beta + 1.

Examples

# After running front tracking simulation
t_edges = np.array([0.0, 10.0, 20.0, 30.0])
c_avg = compute_bin_averaged_concentration_exact(
    t_edges=t_edges,
    v_outlet=500.0,
    waves=tracker.state.waves,
    sorption=sorption,
)
len(c_avg) == len(t_edges) - 1
np.all(c_avg >= 0)

gwtransport.fronttracking.output.compute_domain_mass(t, v_outlet, waves, sorption)[source]#

Compute total mass in domain [0, v_outlet] at time t using exact analytical integration.

Implements runtime mass balance verification as described in High Priority #3 of FRONT_TRACKING_REBUILD_PLAN.md. Integrates concentration over space:

M(t) = ∫₀^v_outlet C(v, t) dv

using exact analytical formulas for each wave type.

Parameters:

t (float) – Time at which to compute domain mass [days].
v_outlet (float) – Outlet position (domain extent) [m³].
waves (Sequence[Wave]) – All waves in the simulation.
sorption (NonlinearSorption | ConstantRetardation) – Sorption model.

Returns:

mass – Total mass in domain [mass]. Computed to machine precision (~1e-14).

Return type:

See also

compute_cumulative_inlet_mass: Cumulative inlet mass
compute_cumulative_outlet_mass: Cumulative outlet mass
concentration_at_point: Point-wise concentration

Notes

Algorithm:

Partition spatial domain [0, v_outlet] into segments where concentration is controlled by a single wave or is constant.
For each segment, compute mass analytically: - Constant C: mass = C * Δv - Rarefaction C(v): use exact spatial integration formula
Sum all segment masses.

Wave Priority (same as concentration_at_point):

Rarefactions (if position is inside rarefaction fan)
Shocks and rarefaction tails (most recent to pass)
Characteristics (smooth regions)

Rarefaction Spatial Integration:

For a rarefaction fan at fixed time t, concentration varies with position v according to the self-similar solution:

R(C) = flow*(t - t_origin)/(v - v_origin)

The spatial integral ∫ C dv is computed analytically using the inverse retardation relation.

Integration Precision:

Constant concentration regions: Exact analytical (C_total * dv)
Rarefaction regions: High-precision trapezoidal quadrature (10+ points)
Overall accuracy: ~1e-10 to 1e-12 relative error
Sufficient for runtime verification; primary outputs remain exact

Examples

# After running simulation to time t=10.0
mass = compute_domain_mass(
    t=10.0, v_outlet=500.0, waves=tracker.state.waves, sorption=sorption
)
mass >= 0.0

gwtransport.fronttracking.output.compute_cumulative_inlet_mass(t, cin, flow, tedges_days)[source]#

Compute cumulative mass entering domain from t=0 to t.

Integrates inlet mass flux over time:: M_in(t) = ∫₀^t cin(τ) * flow(τ) dτ

using exact analytical integration of piecewise-constant functions.

Parameters:

t (float) – Time up to which to integrate [days].
cin (TypeAliasType) – Inlet concentration time series [mass/volume]. Piecewise constant within bins defined by tedges_days.
flow (TypeAliasType) – Flow rate time series [m³/day]. Piecewise constant within bins defined by tedges_days.
tedges_days (TypeAliasType) – Time bin edges [days]. Length len(cin) + 1.

Returns:

mass_in – Cumulative mass entered [mass].

Return type:

Notes

For piecewise-constant cin and flow:: M_in = Σ cin[i] * flow[i] * dt[i]

where the sum is over all bins from tedges_days[0] to t. Partial bins are handled exactly.

Examples

mass_in = compute_cumulative_inlet_mass(
    t=50.0, cin=cin, flow=flow, tedges_days=tedges_days
)
mass_in >= 0.0

gwtransport.fronttracking.output.find_last_rarefaction_start_time(v_outlet, waves)[source]#

Find the time when the last rarefaction head reaches the outlet.

For rarefactions, we integrate analytically so we only need to know when the rarefaction starts at the outlet (head arrival).

Parameters:

v_outlet (float) – Outlet position [m³].
waves (Sequence[Wave]) – All waves in the simulation.

Returns:

t_last – Time when last rarefaction head reaches outlet [days]. For non-rarefaction waves, uses their arrival time. Returns 0.0 if no waves reach the outlet.

Return type:

Notes

This function finds when the last wave structure starts at the outlet. For rarefactions, this is the head arrival time. The tail may arrive much later (or at infinite time for rarefactions to C=0), but the total mass in the rarefaction is computed analytically.

Examples

t_last = find_last_rarefaction_start_time(v_outlet=500.0, waves=all_waves)
t_last >= 0.0

gwtransport.fronttracking.output.compute_cumulative_outlet_mass(t, v_outlet, waves, sorption, flow, tedges_days)[source]#

Compute cumulative mass exiting domain from t=0 to t.

Integrates outlet mass flux over time:: M_out(t) = ∫₀^t cout(τ) * flow(τ) dτ

using exact analytical integration. Outlet concentration cout(τ) is obtained from the wave solution, and flow(τ) is piecewise constant.

Parameters:

t (float) – Time up to which to integrate [days].
v_outlet (float) – Outlet position [m³].
waves (Sequence[Wave]) – All waves in the simulation.
sorption (NonlinearSorption | ConstantRetardation) – Sorption model.
flow (TypeAliasType) – Flow rate time series [m³/day]. Piecewise constant within bins defined by tedges_days.
tedges_days (TypeAliasType) – Time bin edges [days]. Length len(flow) + 1.

Returns:

mass_out – Cumulative mass exited [mass].

Return type:

Notes

The outlet concentration is obtained from wave solution via concentration_at_point(v_outlet, τ, waves, sorption).

For each flow bin [t_i, t_{i+1}], the mass flux integral is computed exactly using identify_outlet_segments and analytical integration (constant segments and rarefaction segments).

Examples

mass_out = compute_cumulative_outlet_mass(
    t=50.0,
    v_outlet=500.0,
    waves=tracker.state.waves,
    sorption=sorption,
    flow=flow,
    tedges_days=tedges_days,
)
mass_out >= 0.0

gwtransport.fronttracking.output.integrate_rarefaction_total_mass(raref, v_outlet, t_start, sorption, flow)[source]#

Compute total mass exiting through a rarefaction.

For a rarefaction that starts at the outlet at time t_start, compute the total mass that will exit through the rarefaction. Integration endpoint depends on the rarefaction tail concentration:

If c_tail ≈ 0: Integrate to infinity (tail extends infinitely)
If c_tail > 0: Integrate to finite tail arrival time

Parameters:

raref (RarefactionWave) – Rarefaction wave.
v_outlet (float) – Outlet position [m³].
t_start (float) – Time when rarefaction head reaches outlet [days].
sorption (NonlinearSorption | ConstantRetardation) – Sorption model.
flow (float) – Flow rate [m³/day] (assumed constant).

Returns:

total_mass – Total mass that exits through rarefaction [mass].

Return type:

matplotlib.figure.Figure

Notes

Uses the exact analytical integral:: M_total = ∫_{t_start}^{t_end} Q * C(t) dt

where C(t) is the concentration at the outlet from the rarefaction wave.

For n > 1 (favorable sorption), rarefactions typically decrease to C=0, so t_end = ∞ and the integral converges.

For n < 1 (unfavorable sorption), rarefactions typically increase from low C to high C, so c_tail > 0 and the tail arrives at finite time t_end = t_start + (v_outlet - v_start) / tail_velocity.

Examples

mass = integrate_rarefaction_total_mass(
    raref=raref_wave,
    v_outlet=500.0,
    t_start=40.0,
    sorption=sorption,
    flow=100.0,
)
mass >= 0.0

gwtransport.fronttracking.output.compute_total_outlet_mass(v_outlet, waves, sorption, flow, tedges_days)[source]#

Compute total integrated outlet mass until all mass has exited.

Automatically determines when the last wave passes the outlet and integrates the outlet mass flux until that time, regardless of the provided tedges extent.

Parameters:

v_outlet (float) – Outlet position [m³].
waves (Sequence[Wave]) – All waves in the simulation.
sorption (NonlinearSorption | ConstantRetardation) – Sorption model.
flow (TypeAliasType) – Flow rate time series [m³/day]. Piecewise constant within bins defined by tedges_days.
tedges_days (TypeAliasType) – Time bin edges [days]. Length len(flow) + 1.

Return type:

tuple[float, float]

Returns:

total_mass_out (float) – Total mass that exits through outlet [mass].
t_integration_end (float) – Time until which integration was performed [days]. This is the time when the last wave passes the outlet.

See also

compute_cumulative_outlet_mass: Cumulative outlet mass up to time t
find_last_rarefaction_start_time: Find when last rarefaction starts
integrate_rarefaction_total_mass: Total mass in rarefaction to infinity

Notes

This function: 1. Finds when the last rarefaction starts at the outlet (head arrival) 2. Integrates outlet mass flux until that time 3. Adds analytical integral of rarefaction mass from start to infinity

For rarefactions to C=0, the tail has infinite arrival time but the total mass is finite and computed analytically.

Examples

total_mass, t_end = compute_total_outlet_mass(
    v_outlet=500.0,
    waves=tracker.state.waves,
    sorption=sorption,
    flow=flow,
    tedges_days=tedges_days,
)
total_mass >= 0.0
t_end >= tedges_days[0]

fronttracking.plot#

Visualization functions for front tracking.

This module provides plotting utilities for visualizing front-tracking simulations: - V-t diagrams showing wave propagation in space-time - Breakthrough curves showing concentration at outlet over time

This file is part of gwtransport which is released under AGPL-3.0 license. See the ./LICENSE file or go to gwtransport/gwtransport for full license details.

gwtransport.fronttracking.plot.plot_vt_diagram(state, ax=None, t_max=None, figsize=(14, 10), *, show_inactive=False, show_events=False)[source]#

Create V-t diagram showing all waves in space-time.

Plots characteristics (blue lines), shocks (red lines), and rarefactions (green fans) in the (time, position) plane. This visualization shows how waves propagate and interact throughout the simulation.

Parameters:

state (FrontTrackerState) – Complete simulation state containing all waves.
ax (matplotlib.axes.Axes, optional) – Existing axes to plot into. If None, a new figure and axes are created using figsize.
t_max (float | None, default: None) – Maximum time to plot [days]. If None, uses final simulation time * 1.2.
figsize (tuple[float, float], default: (14, 10)) – Figure size in inches (width, height). Default (14, 10).
show_inactive (bool, default: False) – Whether to show inactive waves (deactivated by interactions). Default False.
show_events (bool, default: False) – Whether to show wave interaction events as markers. Default False.

Returns:

fig – Figure object containing the V-t diagram.

Return type:

Notes

Characteristics appear as blue lines (constant velocity).
Shocks appear as thick red lines (jump discontinuities).
Rarefactions appear as green fans (smooth transition regions).
Outlet position is shown as a horizontal dashed line.
Only waves within domain [0, v_outlet] are plotted.

Examples

from gwtransport.fronttracking.solver import FrontTracker

tracker = FrontTracker(cin, flow, tedges, aquifer_pore_volume, sorption)
tracker.run()
fig = plot_vt_diagram(tracker.state)
fig.savefig("vt_diagram.png")

gwtransport.fronttracking.plot.plot_breakthrough_curve(state, ax=None, t_max=None, n_rarefaction_points=50, figsize=(12, 6), t_first_arrival=None)[source]#

Plot exact analytical concentration breakthrough curve at outlet.

Uses wave segment information to plot the exact analytical solution without discretization. Constant concentration regions are plotted as horizontal lines, and rarefaction regions are plotted using their exact self-similar solutions.

Parameters:

state (FrontTrackerState) – Complete simulation state containing all waves.
ax (matplotlib.axes.Axes, optional) – Existing axes to plot into. If None, a new figure and axes are created using figsize.
t_max (float | None, default: None) – Maximum time to plot [days]. If None, uses final simulation time * 1.1.
n_rarefaction_points (int, default: 50) – Number of points to use for plotting rarefaction segments (analytical curves). Default 50 per rarefaction segment.
figsize (tuple[float, float], default: (12, 6)) – Figure size in inches (width, height). Default (12, 6).
t_first_arrival (float | None, default: None) – First arrival time for marking spin-up period [days]. If None, spin-up period is not plotted.

Returns:

fig – Figure object containing the breakthrough curve.

Return type:

matplotlib.figure.Figure

Notes

Uses identify_outlet_segments to get exact analytical segment boundaries
Constant concentration segments plotted as horizontal lines (no discretization)
Rarefaction segments plotted using exact self-similar solution
Shocks appear as instantaneous jumps at exact crossing times
No bin averaging or discretization artifacts

Examples

from gwtransport.fronttracking.solver import FrontTracker

tracker = FrontTracker(cin, flow, tedges, aquifer_pore_volume, sorption)
tracker.run()
fig = plot_breakthrough_curve(tracker.state)
fig.savefig("exact_breakthrough.png")

gwtransport.fronttracking.plot.plot_wave_interactions(state, figsize=(14, 8), ax=None)[source]#

Plot event timeline showing wave interactions.

Creates a scatter plot showing when and where different types of wave interactions occur during the simulation.

Parameters:

state (FrontTrackerState) – Complete simulation state containing all events.
ax (matplotlib.axes.Axes, optional) – Existing axes to plot into. If None, a new figure and axes are created using figsize.
figsize (tuple[float, float], default: (14, 8)) – Figure size in inches (width, height). Default (14, 8).

Returns:

fig – Figure object containing the event timeline.

Return type:

matplotlib.figure.Figure

Notes

Each event type is shown with a different color and marker.
Outlet crossings are shown separately from internal collisions.
Event locations are plotted in the (time, position) plane.

Examples

from gwtransport.fronttracking.solver import FrontTracker

tracker = FrontTracker(cin, flow, tedges, aquifer_pore_volume, sorption)
tracker.run()
fig = plot_wave_interactions(tracker.state)
fig.savefig("wave_interactions.png")

gwtransport.fronttracking.plot.plot_inlet_concentration(tedges, cin, ax=None, *, t_first_arrival=None, event_markers=None, color='blue', t_max=None, xlabel='Time [days]', ylabel='Concentration', title='Inlet Concentration', figsize=(8, 5), **step_kwargs)[source]#

Plot inlet concentration as step function with optional markers.

Parameters:

tedges (pandas.DatetimeIndex) – Time bin edges for inlet concentration. Length = len(cin) + 1.
cin (ArrayLike) – Inlet concentration values. Length = len(tedges) - 1.
ax (matplotlib.axes.Axes, optional) – Existing axes to plot into. If None, creates new figure.
t_first_arrival (float, optional) – First arrival time to mark with vertical line [days].
event_markers (list of dict, optional) – Event markers to add. Each dict should have keys: ‘time’, ‘label’, ‘color’.
color (str, optional) – Color for inlet concentration line. Default ‘blue’.
t_max (float, optional) – Maximum time for x-axis [days]. If None, uses full range.
xlabel (str, optional) – Label for x-axis. Default ‘Time [days]’.
ylabel (str, optional) – Label for y-axis. Default ‘Concentration’.
title (str, optional) – Plot title. Default ‘Inlet Concentration’.
figsize (tuple, optional) – Figure size if creating new figure. Default (8, 5).
**step_kwargs – Additional arguments passed to ax.plot().

Returns:

fig (matplotlib.figure.Figure) – Figure object.
ax (matplotlib.axes.Axes) – Axes object.

gwtransport.fronttracking.plot.plot_front_tracking_summary(structure, tedges, cin, cout_tedges, cout, *, figsize=(16, 10), show_exact=True, show_bin_averaged=True, show_events=True, show_inactive=False, t_max=None, title=None, inlet_color='blue', outlet_exact_color='blue', outlet_binned_color='red', first_arrival_color='green')[source]#

Create comprehensive 3-panel summary figure for front tracking simulation.

Creates a multi-panel visualization with: - Top-left: V-t diagram showing wave propagation - Top-right: Inlet concentration time series - Bottom: Outlet concentration (exact and/or bin-averaged)

Parameters:

structure (dict) – Structure returned from infiltration_to_extraction_front_tracking_detailed. Must contain keys: ‘tracker_state’, ‘t_first_arrival’.
tedges (pandas.DatetimeIndex) – Time bin edges for inlet concentration. Length = len(cin) + 1.
cin (ArrayLike) – Inlet concentration values. Length = len(tedges) - 1.
cout_tedges (pandas.DatetimeIndex) – Output time bin edges for bin-averaged concentration. Length = len(cout) + 1.
cout (ArrayLike) – Bin-averaged output concentration values. Length = len(cout_tedges) - 1.
figsize (tuple, optional) – Figure size (width, height). Default (16, 10).
show_exact (bool, optional) – Whether to show exact analytical breakthrough curve. Default True.
show_bin_averaged (bool, optional) – Whether to show bin-averaged concentration. Default True.
show_events (bool, optional) – Whether to show wave interaction events on V-t diagram. Default True.
show_inactive (bool, optional) – Whether to show inactive waves on V-t diagram. Default False.
t_max (float, optional) – Maximum time for plots [days]. If None, uses input data range.
title (str, optional) – Overall figure title. If None, uses generic title.
inlet_color (str, optional) – Color for inlet concentration. Default ‘blue’.
outlet_exact_color (str, optional) – Color for exact outlet curve. Default ‘blue’.
outlet_binned_color (str, optional) – Color for bin-averaged outlet. Default ‘red’.
first_arrival_color (str, optional) – Color for first arrival marker. Default ‘green’.

Returns:

fig (matplotlib.figure.Figure) – Figure object.
axes (dict) – Dictionary with keys ‘vt’, ‘inlet’, ‘outlet’ containing axes objects.

gwtransport.fronttracking.plot.plot_sorption_comparison(pulse_favorable_structure, pulse_unfavorable_structure, pulse_tedges, pulse_cin, dip_favorable_structure, dip_unfavorable_structure, dip_tedges, dip_cin, *, figsize=(16, 12), t_max_pulse=None, t_max_dip=None)[source]#

Compare how each inlet produces different outputs with n>1 vs n<1 sorption.

Creates a 2x3 grid: - Row 1: Pulse inlet and its outputs with n>1 and n<1 sorption - Row 2: Dip inlet and its outputs with n>1 and n<1 sorption

This demonstrates how the SAME inlet timeseries produces DIFFERENT breakthrough curves depending on the sorption isotherm.

Parameters:

pulse_favorable_structure (dict) – Structure from pulse inlet with n>1 (higher C travels faster).
pulse_unfavorable_structure (dict) – Structure from pulse inlet with n<1 (lower C travels faster).
pulse_tedges (pandas.DatetimeIndex) – Time bin edges for pulse inlet. Length = len(pulse_cin) + 1.
pulse_cin (ArrayLike) – Pulse inlet concentration (e.g., 0→10→0). Length = len(pulse_tedges) - 1.
dip_favorable_structure (dict) – Structure from dip inlet with n>1 (higher C travels faster).
dip_unfavorable_structure (dict) – Structure from dip inlet with n<1 (lower C travels faster).
dip_tedges (pandas.DatetimeIndex) – Time bin edges for dip inlet. Length = len(dip_cin) + 1.
dip_cin (ArrayLike) – Dip inlet concentration (e.g., 10→2→10). Length = len(dip_tedges) - 1.
figsize (tuple, optional) – Figure size (width, height). Default (16, 12).
t_max_pulse (float, optional) – Max time for pulse plots [days]. If None, auto-computed.
t_max_dip (float, optional) – Max time for dip plots [days]. If None, auto-computed.

Returns:

fig (matplotlib.figure.Figure) – Figure object.
axes (numpy.ndarray) – 2x3 array of axes objects.

fronttracking.solver#

Front Tracking Solver - Event-Driven Simulation Engine.

This module implements the main event-driven front tracking solver for nonlinear sorption transport. The solver maintains a list of waves and processes collision events chronologically using exact analytical calculations.

The algorithm: 1. Initialize waves from inlet boundary conditions 2. Find next event (earliest collision or outlet crossing) 3. Advance time to event 4. Handle event (create new waves, deactivate old ones) 5. Repeat until no more events

All calculations are exact analytical with machine precision.

class gwtransport.fronttracking.solver.FrontTrackerState(waves, events, t_current, v_outlet, sorption, cin, flow, tedges)[source]#

Bases: object

Complete state of the front tracking simulation.

This dataclass holds all information about the current simulation state, including all waves (active and inactive), event history, and simulation parameters.

Parameters:

waves (list[Wave]) – All waves created during simulation (includes inactive waves)
events (list[dict]) – Event history with details about each event
t_current (float) – Current simulation time [days from tedges[0]]
v_outlet (float) – Outlet position [m³]
sorption (NonlinearSorption | ConstantRetardation) – Sorption parameters
cin (ndarray) – Inlet concentration time series [mass/volume]
flow (ndarray) – Flow rate time series [m³/day]
tedges (DatetimeIndex) – Time bin edges

Examples

state = FrontTrackerState(
    waves=[],
    events=[],
    t_current=0.0,
    v_outlet=500.0,
    sorption=sorption,
    cin=cin,
    flow=flow,
    tedges=tedges,
)

waves: list[Wave]#

events: list[dict]#

t_current: float#

v_outlet: float#

sorption: NonlinearSorption | ConstantRetardation#

cin: ndarray#

flow: ndarray#

tedges: DatetimeIndex#

__init__(waves, events, t_current, v_outlet, sorption, cin, flow, tedges)#

class gwtransport.fronttracking.solver.FrontTracker(cin, flow, tedges, aquifer_pore_volume, sorption)[source]#

Bases: object

Event-driven front tracking solver for nonlinear sorption transport.

This is the main simulation engine that orchestrates wave propagation, event detection, and event handling. The solver maintains a list of waves and processes collision events chronologically.

Parameters:

cin (ndarray) – Inlet concentration time series [mass/volume]
flow (ndarray) – Flow rate time series [m³/day]
tedges (DatetimeIndex) – Time bin edges [days]
aquifer_pore_volume (float) – Total pore volume [m³]
sorption (NonlinearSorption | ConstantRetardation) – Sorption parameters

state#

Complete simulation state

Type:: FrontTrackerState

t_first_arrival#

First arrival time (end of spin-up period) [days]

Type:: float

Notes

The solver uses exact analytical calculations throughout with no numerical tolerances or iterative methods. All wave interactions are detected and handled with machine precision.

The spin-up period (t < t_first_arrival) is affected by unknown initial conditions. Results are only valid for t >= t_first_arrival.

Examples

tracker = FrontTracker(
    cin=cin,
    flow=flow,
    tedges=tedges,
    aquifer_pore_volume=500.0,
    sorption=sorption,
)
tracker.run(max_iterations=1000)
# Access results
print(f"Total events: {len(tracker.state.events)}")
print(f"Active waves: {sum(1 for w in tracker.state.waves if w.is_active)}")

__init__(cin, flow, tedges, aquifer_pore_volume, sorption)[source]#

Initialize tracker with inlet conditions and physical parameters.

Parameters:

cin (ndarray) – Inlet concentration time series [mass/volume]
flow (ndarray) – Flow rate time series [m³/day]
tedges (DatetimeIndex) – Time bin edges [days]
aquifer_pore_volume (float) – Total pore volume [m³]
sorption (NonlinearSorption | ConstantRetardation) – Sorption parameters

Raises:

ValueError – If input arrays have incompatible lengths or invalid values

find_next_event()[source]#

Find the next event (earliest in time).

Searches all possible wave interactions and returns the earliest event. Uses a priority queue (min-heap) to efficiently find the minimum time.

Returns:: Next event to process, or None if no future events
Return type:: Event | None

Notes

Checks for: - Characteristic-characteristic collisions - Shock-shock collisions - Shock-characteristic collisions - Rarefaction-characteristic collisions (head/tail) - Shock-rarefaction collisions (head/tail) - Outlet crossings for all wave types

All collision times are computed using exact analytical formulas.

handle_event(event)[source]#

Handle an event by calling appropriate handler and updating state.

Dispatches to the correct event handler based on event type, then updates the simulation state with any new waves created.

Parameters:: event (Event) – Event to handle
Raises:: RuntimeError – If the FLOW_CHANGE event does not have flow_new set.

Notes

Event handlers may: - Deactivate parent waves - Create new child waves - Record event details in history - Verify physical correctness (entropy, mass balance)

run(max_iterations=10000, *, verbose=False)[source]#

Run simulation until no more events or max_iterations reached.

Processes events chronologically by repeatedly finding the next event, advancing time, and handling the event. Continues until no more events exist or the iteration limit is reached.

Parameters:

max_iterations (int, default: 10000) – Maximum number of events to process. Default 10000. Prevents infinite loops in case of bugs.
verbose (bool, default: False) – Print progress messages. Default False.

Notes

The simulation stops when: - No more events exist (all waves have exited or become inactive) - max_iterations is reached (safety limit)

After completion, results are available in: - self.state.waves: All waves (active and inactive) - self.state.events: Complete event history - self.t_first_arrival: End of spin-up period

verify_physics(*, check_mass_balance=False, mass_balance_rtol=1e-12)[source]#

Verify physical correctness of current state.

Implements High Priority #3 from FRONT_TRACKING_REBUILD_PLAN.md by adding runtime mass balance verification using exact analytical integration.

Checks: - All shocks satisfy Lax entropy condition - All rarefactions have proper head/tail ordering - Mass balance: mass_in_domain + mass_out = mass_in (to specified tolerance)

Parameters:

check_mass_balance (bool, default: False) – Enable mass balance verification. Default False (opt-in for now).
mass_balance_rtol (float, default: 1e-12) – Relative tolerance for mass balance check. Default 1e-6. This tolerance accounts for: - Midpoint approximation in spatial integration of rarefactions - Numerical precision in wave position calculations - Piecewise-constant approximations in domain partitioning

Raises:

RuntimeError – If physics violation is detected

Notes

Mass balance equation:: mass_in_domain(t) + mass_out_cumulative(t) = mass_in_cumulative(t)

All mass calculations use exact analytical integration where possible: - Inlet/outlet temporal integrals: exact for piecewise-constant functions - Domain spatial integrals: exact for constants, midpoint rule for rarefactions - Overall precision: ~1e-10 to 1e-12 relative error

fronttracking.validation#

Physics validation utilities for front tracking.

This module provides functions to verify physical correctness of front-tracking simulations, including entropy conditions, concentration bounds, and mass conservation.

This file is part of gwtransport which is released under AGPL-3.0 license. See the ./LICENSE file or go to gwtransport/gwtransport for full license details.

gwtransport.fronttracking.validation.verify_physics(structure, cout, cout_tedges, cin, *, verbose=True, rtol=1e-10)[source]#

Run comprehensive physics verification checks on front tracking results.

Performs the following checks: 1. Entropy condition for all shocks 2. No negative concentrations (within tolerance) 3. Output concentration ≤ input maximum 4. Finite first arrival time 5. No NaN values after spin-up period 6. Events chronologically ordered 7. Rarefaction concentration ordering (head vs tail) 8. Total integrated outlet mass (until all mass exits)

Parameters:

structure (dict) – Structure returned from infiltration_to_extraction_front_tracking_detailed. Must contain keys: ‘waves’, ‘t_first_arrival’, ‘events’, ‘n_shocks’, ‘n_rarefactions’, etc.
cout (ArrayLike) – Bin-averaged output concentrations.
cout_tedges (pandas.DatetimeIndex) – Output time edges for bins.
cin (ArrayLike) – Input concentrations.
verbose (bool, optional) – If True, print detailed results. If False, only return summary. Default True.
rtol (float, optional) – Relative tolerance for numerical checks. Default 1e-10.

Returns:

results – Dictionary containing: - ‘all_passed’: bool - True if all checks passed - ‘n_checks’: int - Total number of checks performed - ‘n_passed’: int - Number of checks that passed - ‘failures’: list of str - Description of failed checks (empty if all passed) - ‘summary’: str - One-line summary (e.g., “✓ All 8 checks passed”)

Return type:

dict

Examples

results = verify_physics(structure, cout, cout_tedges, cin, verbose=False)
print(results["summary"])
# ✓ All 8 checks passed
assert results["all_passed"]

fronttracking.waves#

Wave Representation for Front Tracking.

This module implements wave classes for representing characteristics, shocks, and rarefaction waves in the front tracking algorithm. Each wave type knows how to compute its position and concentration at any point in space-time.

This file is part of gwtransport which is released under AGPL-3.0 license. See the ./LICENSE file or go to gwtransport/gwtransport for full license details.

class gwtransport.fronttracking.waves.Wave(t_start, v_start, flow, *, is_active=True)[source]#

Bases: ABC

Abstract base class for all wave types in front tracking.

All waves share common attributes and must implement methods for computing position and concentration. Waves can be active or inactive (deactivated waves are preserved for history but don’t participate in future interactions).

Parameters:

t_start (float) – Time when wave forms [days].
v_start (float) – Position where wave forms [m³].
flow (float) – Flow rate at formation time [m³/day].
is_active (bool, default: True) – Whether wave is currently active. Default True.

t_start: float#: Time when wave forms [days].

v_start: float#: Position where wave forms [m³].

flow: float#: Flow rate at formation time [m³/day].

is_active: bool = True#: Whether wave is currently active.

abstractmethod position_at_time(t)[source]#

Compute wave position at time t.

Parameters:: t (float) – Time [days].
Returns:: position – Position [m³], or None if t < t_start or wave is inactive.
Return type:: float | None

abstractmethod concentration_left()[source]#

Get concentration on left (upstream) side of wave.

Returns:: c_left – Upstream concentration [mass/volume].
Return type:: float

abstractmethod concentration_right()[source]#

Get concentration on right (downstream) side of wave.

Returns:: c_right – Downstream concentration [mass/volume].
Return type:: float

abstractmethod concentration_at_point(v, t)[source]#

Compute concentration at point (v, t) if wave controls it.

Parameters:

v (float) – Position [m³].
t (float) – Time [days].

Returns:

concentration – Concentration [mass/volume] if wave controls this point, None otherwise.

Return type:

__init__(t_start, v_start, flow, *, is_active=True)#

class gwtransport.fronttracking.waves.CharacteristicWave(t_start, v_start, flow, concentration, sorption, *, is_active=True)[source]#

Bases: Wave

Characteristic line along which concentration is constant.

In smooth regions, concentration travels at velocity flow/R(C). Along each characteristic line, the concentration value is constant. This is the fundamental solution element for hyperbolic conservation laws.

Parameters:

t_start (float) – Formation time [days].
v_start (float) – Starting position [m³].
flow (float) – Flow rate [m³/day].
concentration (float) – Constant concentration carried [mass/volume].
sorption (NonlinearSorption | ConstantRetardation) – Sorption model determining velocity.
is_active (bool, default: True) – Activity status. Default True.

Examples

>>> sorption = ConstantRetardation(retardation_factor=2.0)
>>> char = CharacteristicWave(
...     t_start=0.0, v_start=0.0, flow=100.0, concentration=5.0, sorption=sorption
... )
>>> char.velocity()
50.0
>>> char.position_at_time(10.0)
500.0

concentration: float#: Constant concentration carried [mass/volume].

sorption: NonlinearSorption | ConstantRetardation#: Sorption model determining velocity.

velocity()[source]#

Compute characteristic velocity.

The velocity is v = flow / R(C), where R is the retardation factor.

Returns:: velocity – Characteristic velocity [m³/day].
Return type:: float

position_at_time(t)[source]#

Compute position at time t.

Characteristics propagate linearly: V(t) = v_start + velocity*(t - t_start).

Parameters:: t (float) – Time [days].
Returns:: position – Position [m³], or None if t < t_start or inactive.
Return type:: float | None

concentration_left()[source]#

Get upstream concentration (same as concentration for characteristics).

Returns:: c_left – Upstream concentration [mass/volume].
Return type:: float

concentration_right()[source]#

Get downstream concentration (same as concentration for characteristics).

Returns:: c_right – Downstream concentration [mass/volume].
Return type:: float

concentration_at_point(v, t)[source]#

Get concentration if point is on this characteristic.

For practical purposes, we check if the characteristic has reached position v by time t.

Parameters:

v (float) – Position [m³].
t (float) – Time [days].

Returns:

concentration – Concentration if point is on characteristic, None otherwise.

Return type:

Notes

In practice, this method is used by higher-level algorithms to determine which wave controls a given point. The exact point-on-line check is handled by the solver.

__init__(t_start, v_start, flow, concentration, sorption, *, is_active=True)#

class gwtransport.fronttracking.waves.ShockWave(t_start, v_start, flow, c_left, c_right, sorption, velocity=None, *, is_active=True)[source]#

Bases: Wave

Shock wave (discontinuity) with jump in concentration.

Shocks form when faster water overtakes slower water, creating a sharp front. The shock velocity is determined by the Rankine-Hugoniot condition to ensure mass conservation across the discontinuity.

Parameters:

t_start (float) – Formation time [days].
v_start (float) – Formation position [m³].
flow (float) – Flow rate [m³/day].
c_left (float) – Concentration upstream (behind) shock [mass/volume].
c_right (float) – Concentration downstream (ahead of) shock [mass/volume].
sorption (NonlinearSorption | ConstantRetardation) – Sorption model.
is_active (bool, default: True) – Activity status. Default True.
velocity (float | None, default: None) – Shock velocity computed from Rankine-Hugoniot. Computed automatically if not provided.

Examples

>>> sorption = FreundlichSorption(
...     k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
... )
>>> shock = ShockWave(
...     t_start=0.0,
...     v_start=0.0,
...     flow=100.0,
...     c_left=10.0,
...     c_right=2.0,
...     sorption=sorption,
... )
>>> shock.velocity > 0
True
>>> shock.satisfies_entropy()
True

c_left: float#: Concentration upstream (behind) shock [mass/volume].

c_right: float#: Concentration downstream (ahead of) shock [mass/volume].

sorption: NonlinearSorption | ConstantRetardation#: Sorption model.

velocity: float | None = None#: Shock velocity computed from Rankine-Hugoniot [m³/day].

__post_init__()[source]#

Compute shock velocity from Rankine-Hugoniot condition.

Return type:: None

position_at_time(t)[source]#

Compute shock position at time t.

Shock propagates at constant velocity: V(t) = v_start + velocity*(t - t_start).

Parameters:: t (float) – Time [days].
Returns:: position – Position [m³], or None if t < t_start or inactive.
Return type:: float | None
Raises:: RuntimeError – If shock velocity is None (should have been set in __post_init__).

concentration_left()[source]#

Get upstream concentration.

Returns:: c_left – Upstream concentration [mass/volume].
Return type:: float

concentration_right()[source]#

Get downstream concentration.

Returns:: c_right – Downstream concentration [mass/volume].
Return type:: float

concentration_at_point(v, t)[source]#

Get concentration at point based on which side of shock.

Parameters:

v (float) – Position [m³].
t (float) – Time [days].

Returns:

concentration – c_left if upstream of shock, c_right if downstream, None if shock hasn’t formed yet.

Return type:

Notes

At the exact shock position, returns the average of left and right values. This is a convention for the singular point; in practice, the shock is infinitesimally thin.

satisfies_entropy()[source]#

Check if shock satisfies Lax entropy condition.

The entropy condition ensures characteristics flow INTO the shock from both sides, which is required for physical admissibility.

Returns:: satisfies – True if shock satisfies entropy condition.
Return type:: bool
Raises:: RuntimeError – If shock velocity is None (should have been set in __post_init__).

Notes

The condition is: lambda(c_left) > shock_velocity > lambda(c_right), where lambda(C) = flow / R(C) is the characteristic velocity.

Shocks that violate this condition are unphysical and should be replaced by rarefaction waves.

__init__(t_start, v_start, flow, c_left, c_right, sorption, velocity=None, *, is_active=True)#

class gwtransport.fronttracking.waves.RarefactionWave(t_start, v_start, flow, c_head, c_tail, sorption, *, is_active=True)[source]#

Bases: Wave

Rarefaction (expansion fan) with smooth concentration gradient.

Rarefactions form when slower water follows faster water, creating an expanding region where concentration varies smoothly. The solution is self-similar: C = C(V/t).

Parameters:

t_start (float) – Formation time [days].
v_start (float) – Formation position [m³].
flow (float) – Flow rate [m³/day].
c_head (float) – Concentration at leading edge (faster) [mass/volume].
c_tail (float) – Concentration at trailing edge (slower) [mass/volume].
sorption (NonlinearSorption | ConstantRetardation) – Sorption model (must be concentration-dependent).
is_active (bool, default: True) – Activity status. Default True.

Raises:

ValueError – If head velocity <= tail velocity (would be compression, not rarefaction).

Examples

>>> sorption = FreundlichSorption(
...     k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
... )
>>> raref = RarefactionWave(
...     t_start=0.0,
...     v_start=0.0,
...     flow=100.0,
...     c_head=10.0,
...     c_tail=2.0,
...     sorption=sorption,
... )
>>> raref.head_velocity() > raref.tail_velocity()
True
>>> raref.contains_point(v=150.0, t=20.0)
True

c_head: float#: Concentration at leading edge (faster) [mass/volume].

__init__(t_start, v_start, flow, c_head, c_tail, sorption, *, is_active=True)#

c_tail: float#: Concentration at trailing edge (slower) [mass/volume].

sorption: NonlinearSorption | ConstantRetardation#: Sorption model (must be concentration-dependent).

__post_init__()[source]#

Verify this is actually a rarefaction (head faster than tail).

Raises:: ValueError – If head velocity is not greater than tail velocity, indicating a compression wave (shock) rather than a rarefaction.

head_velocity()[source]#

Compute velocity of rarefaction head (leading edge).

Returns:: velocity – Head velocity [m³/day].
Return type:: float

tail_velocity()[source]#

Compute velocity of rarefaction tail (trailing edge).

Returns:: velocity – Tail velocity [m³/day].
Return type:: float

head_position_at_time(t)[source]#

Compute position of rarefaction head at time t.

Parameters:: t (float) – Time [days].
Returns:: position – Head position [m³], or None if t < t_start or inactive.
Return type:: float | None

tail_position_at_time(t)[source]#

Compute position of rarefaction tail at time t.

Parameters:: t (float) – Time [days].
Returns:: position – Tail position [m³], or None if t < t_start or inactive.
Return type:: float | None

position_at_time(t)[source]#

Return head position (leading edge of rarefaction).

This implements the abstract Wave method.

Parameters:: t (float) – Time [days].
Returns:: position – Head position [m³].
Return type:: float | None

contains_point(v, t)[source]#

Check if point (v, t) is inside the rarefaction fan.

Parameters:

v (float) – Position [m³].
t (float) – Time [days].

Returns:

contains – True if point is between tail and head.

Return type:

bool

concentration_left()[source]#

Get upstream concentration (tail).

Returns:: c_left – Upstream concentration at the trailing edge [mass/volume].
Return type:: float

concentration_right()[source]#

Get downstream concentration (head).

Returns:: c_right – Downstream concentration at the leading edge [mass/volume].
Return type:: float

concentration_at_point(v, t)[source]#

Compute concentration at point (v, t) via self-similar solution.

Within the rarefaction fan, concentration varies smoothly according to:: R(C) = flow * (t - t_start) / (v - v_start)

This is inverted to find C using the sorption model.

Parameters:

v (float) – Position [m³].
t (float) – Time [days].

Returns:

concentration – Concentration if point is in rarefaction, None otherwise.

Return type:

Notes

The self-similar solution automatically maintains mass balance and provides the exact analytical form of the concentration profile.

For ConstantRetardation, rarefactions don’t form (all concentrations travel at same speed), so this returns None.

Examples

>>> sorption = FreundlichSorption(
...     k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
... )
>>> raref = RarefactionWave(0.0, 0.0, 100.0, 10.0, 2.0, sorption)
>>> c = raref.concentration_at_point(v=150.0, t=20.0)
>>> c is not None
True
>>> 2.0 <= c <= 10.0
True

gamma#

Gamma Distribution Utilities for Aquifer Pore Volume Heterogeneity.

This module provides utilities for working with gamma distributions to model heterogeneous aquifer pore volumes in groundwater transport analysis. The gamma distribution offers a flexible three-parameter model (shape, scale, location) for representing the natural variability in flow path lengths and residence times within aquifer systems. In heterogeneous aquifers, water travels through multiple flow paths with different pore volumes; the location parameter additionally represents a guaranteed minimum pore volume (for example, immobile porosity or a geometric minimum travel distance).

Parameterizations#

Two equivalent parameterizations are supported, each optionally with a location shift:

(mean, std, loc) — physically intuitive. mean is the total expected value, std is the spread (invariant under shift), and loc is the lower bound of support. Constraint: 0 <= loc < mean.
(alpha, beta, loc) — scipy-style. alpha is shape, beta is scale, and loc is the lower bound of support. Constraint: alpha > 0, beta > 0, loc >= 0.

Conversion formulas (with constraint mean > loc):

alpha = ((mean - loc) / std) ** 2 beta = std ** 2 / (mean - loc) mean = alpha * beta + loc std = sqrt(alpha) * beta

When loc == 0 the three-parameter model reduces to the standard two-parameter gamma distribution.

Available functions:

parse_parameters() - Parse and validate gamma distribution parameters from either (mean, std, loc) or (alpha, beta, loc). Ensures exactly one parameter pair is provided and validates positivity and ordering constraints.
mean_std_loc_to_alpha_beta() - Convert physically intuitive (mean, std, loc) parameters to gamma shape/scale parameters.
alpha_beta_loc_to_mean_std() - Convert gamma (alpha, beta, loc) parameters back to (mean, std) for physical interpretation.
bins() - Primary function for transport modeling. Creates discrete probability bins from the (optionally shifted) gamma distribution with equal-probability bins (default) or custom quantile edges. Returns bin edges, expected values (mean pore volume within each bin), and probability masses (weight in transport calculations).
bin_masses() - Calculate probability mass for custom bin edges using the incomplete gamma function. Lower-level function used internally by bins().

This file is part of gwtransport which is released under AGPL-3.0 license. See the ./LICENSE file or go to gwtransport/gwtransport for full license details.

gwtransport.gamma.parse_parameters(*, mean=None, std=None, loc=0.0, alpha=None, beta=None)[source]#

Parse parameters for gamma distribution.

Either (mean, std) or (alpha, beta) must be provided. loc is optional and defaults to 0, which recovers the standard two-parameter gamma distribution.

Parameters:

mean (float | None, default: None) – Mean of the gamma distribution. Must be strictly greater than loc.
std (float | None, default: None) – Standard deviation of the gamma distribution. Must be positive. See Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity for what std represents depending on APVD source. std is invariant under the loc shift.
loc (float, default: 0.0) – Location (horizontal shift) of the gamma distribution; the lower bound of support. Must satisfy loc >= 0 and, when mean is supplied, loc < mean. Default is 0.0.
alpha (float | None, default: None) – Shape parameter of gamma distribution (must be > 0).
beta (float | None, default: None) – Scale parameter of gamma distribution (must be > 0).

Return type:

tuple[float, float, float]

Returns:

alpha (float) – Shape parameter of gamma distribution.
beta (float) – Scale parameter of gamma distribution.
loc (float) – Location parameter of gamma distribution.

Raises:

ValueError – If neither (mean, std) nor (alpha, beta) is provided, if only one of a pair is provided, if alpha or beta are not positive, if loc is negative, or if mean <= loc.

gwtransport.gamma.mean_std_loc_to_alpha_beta(*, mean, std, loc=0.0)[source]#

Convert mean, standard deviation, and location of gamma distribution to shape/scale.

The two-parameter shape/scale representation (alpha, beta) is derived from the excess-over-loc moments: mean_excess = mean - loc, std_excess = std.

Parameters:

mean (float) – Mean of the gamma distribution. Must be strictly greater than loc.
std (float) – Standard deviation of the gamma distribution. Must be positive. See Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity for what std represents depending on APVD source. std is invariant under the loc shift.
loc (float, default: 0.0) – Location (horizontal shift) of the gamma distribution. Must satisfy 0 <= loc < mean. Default is 0.0.

Return type:

tuple[float, float]

Returns:

alpha (float) – Shape parameter of gamma distribution.
beta (float) – Scale parameter of gamma distribution.

Raises:

ValueError – If std is not positive, if loc is negative, or if mean <= loc.

See also

alpha_beta_loc_to_mean_std: Convert shape/scale/loc parameters to mean and std.
parse_parameters: Parse and validate gamma distribution parameters.

Examples

>>> from gwtransport.gamma import mean_std_loc_to_alpha_beta
>>> mean_pore_volume = 30000.0  # m³
>>> std_pore_volume = 8100.0  # m³
>>> alpha, beta = mean_std_loc_to_alpha_beta(
...     mean=mean_pore_volume, std=std_pore_volume
... )
>>> print(f"Shape parameter (alpha): {alpha:.2f}")
Shape parameter (alpha): 13.72
>>> print(f"Scale parameter (beta): {beta:.2f}")
Scale parameter (beta): 2187.00

With a 5000 m³ minimum pore volume:

>>> alpha, beta = mean_std_loc_to_alpha_beta(mean=30000.0, std=8100.0, loc=5000.0)
>>> print(f"Shape parameter (alpha): {alpha:.2f}")
Shape parameter (alpha): 9.45
>>> print(f"Scale parameter (beta): {beta:.2f}")
Scale parameter (beta): 2646.00

gwtransport.gamma.alpha_beta_loc_to_mean_std(*, alpha, beta, loc=0.0)[source]#

Convert shape, scale, and location of gamma distribution to mean and standard deviation.

Parameters:

alpha (float) – Shape parameter of the gamma distribution. Must be positive.
beta (float) – Scale parameter of the gamma distribution. Must be positive.
loc (float, default: 0.0) – Location (horizontal shift) of the gamma distribution. Must be non-negative. Default is 0.0.

Return type:

tuple[float, float]

Returns:

mean (float) – Mean of the gamma distribution, equal to alpha * beta + loc.
std (float) – Standard deviation of the gamma distribution, equal to sqrt(alpha) * beta. std is invariant under the loc shift.

Raises:

ValueError – If loc is negative.

See also

mean_std_loc_to_alpha_beta: Convert mean/std/loc to shape and scale parameters.
parse_parameters: Parse and validate gamma distribution parameters.

Examples

>>> from gwtransport.gamma import alpha_beta_loc_to_mean_std
>>> alpha = 13.72  # shape parameter
>>> beta = 2187.0  # scale parameter
>>> mean, std = alpha_beta_loc_to_mean_std(alpha=alpha, beta=beta)
>>> print(f"Mean pore volume: {mean:.0f} m³")
Mean pore volume: 3000... m³
>>> print(f"Std pore volume: {std:.0f} m³")
Std pore volume: 810... m³

gwtransport.gamma.bins(*, mean=None, std=None, loc=0.0, alpha=None, beta=None, n_bins=100, quantile_edges=None)[source]#

Divide a (shifted) gamma distribution into bins and compute bin properties.

If n_bins is provided, the gamma distribution is divided into n_bins equal-mass bins. If quantile_edges is provided, the distribution is divided into bins defined by those quantile edges. The quantile edges must lie in [0, 1] with size n_bins + 1; the first and last entries must be 0 and 1.

Parameters:

mean (float | None, default: None) – Mean of the gamma distribution. Must be strictly greater than loc.
std (float | None, default: None) – Standard deviation of the gamma distribution. Must be positive.
loc (float, default: 0.0) – Location (horizontal shift) of the gamma distribution; the lower bound of support. Must satisfy 0 <= loc < mean (or loc >= 0 when using alpha/beta). Default is 0.0.
alpha (float | None, default: None) – Shape parameter of gamma distribution (must be > 0).
beta (float | None, default: None) – Scale parameter of gamma distribution (must be > 0).
n_bins (int, default: 100) – Number of bins to divide the gamma distribution (must be > 1). Default is 100.
quantile_edges (ndarray | None, default: None) – Quantile edges for binning. Must be in the range [0, 1] and of size n_bins + 1. The first and last quantile edges must be 0 and 1, respectively. If provided, n_bins is ignored.

Returns:

Dictionary with keys of type str and values of type numpy.ndarray:

lower_bound: lower bounds of bins (first one equals loc)
upper_bound: upper bounds of bins (last one is inf)
edges: bin edges (lower_bound[0], upper_bound[0], …, upper_bound[-1])
expected_values: expected values in bins. Is what you would expect to observe if you repeatedly sampled from the probability distribution, but only considered samples that fall within that particular bin.
probability_mass: probability mass in bins (invariant under loc shift).

Return type:

dict[str, ndarray[tuple[Any, ...], dtype[floating]]]

Raises:

ValueError – If n_bins is not greater than 1, or if parameter validation in parse_parameters() fails.

See also

bin_masses: Calculate probability mass for bins.
mean_std_loc_to_alpha_beta: Convert mean/std/loc to alpha/beta parameters.
gwtransport.advection.gamma_infiltration_to_extraction: Use bins for transport modeling.
Gamma Distribution Model: Two-parameter pore volume model.
Shifted Gamma Distribution (Minimum Pore Volume): Shifted gamma with minimum pore volume.
Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity: What std represents (macrodispersion vs total spreading).
8. Gamma Distribution Adequacy: When gamma distribution is adequate.

Examples

Create equal-mass bins for a gamma distribution:

>>> from gwtransport.gamma import bins
>>> result = bins(mean=30000.0, std=8100.0, n_bins=5)

With a location parameter representing a minimum pore volume:

>>> result = bins(mean=30000.0, std=8100.0, loc=5000.0, n_bins=5)
>>> float(result["edges"][0])
5000.0

Create bins with custom quantile edges:

>>> import numpy as np
>>> quantiles = np.array([0.0, 0.25, 0.5, 0.75, 1.0])
>>> result = bins(mean=30000.0, std=8100.0, quantile_edges=quantiles)
>>> print(f"Number of bins: {len(result['probability_mass'])}")
Number of bins: 4

gwtransport.gamma.bin_masses(*, alpha, beta, bin_edges)[source]#

Calculate probability mass for each bin in a standard (unshifted) gamma distribution.

Is the area under the Gamma(alpha, beta) PDF between the bin edges. This lower-level function operates on the unshifted gamma distribution; if a location shift is needed, callers should subtract loc from their physical bin edges before passing them in. Because probability mass is invariant under a location shift, the result is identical to that of the shifted distribution.

Parameters:

alpha (float) – Shape parameter of gamma distribution (must be > 0).
beta (float) – Scale parameter of gamma distribution (must be > 0).
bin_edges (TypeAliasType) – Bin edges of the unshifted distribution. Array of increasing values of size len(bins) + 1. Must be non-negative.

Returns:

Probability mass for each bin.

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

Raises:

ValueError – If alpha or beta are not positive, if bin_edges contains fewer than two values, if bin_edges are not monotonically increasing, or if any bin_edges are negative.

logremoval#

Log Removal Calculations for First-Order Decay Processes.

This module provides utilities to calculate log removal values from first-order decay processes, including pathogen inactivation and radioactive decay. The module supports basic log removal calculations and parallel flow arrangements where multiple flow paths operate simultaneously.

First-Order Decay Model#

The log removal from any first-order decay process is:

Log Removal = log10_decay_rate * residence_time

where log10_decay_rate has units [log10/day] and residence_time has units [days]. This is equivalent to exponential decay C_out/C_in = 10^(-mu * t), where mu is the log10 decay rate and t is residence time. The natural-log decay rate constant lambda [1/day] is related to mu by lambda = mu * ln(10).

This model applies to any process that follows first-order kinetics:

Pathogen inactivation: viruses, bacteria, and protozoa lose infectivity over time
Radioactive decay: isotopes used for groundwater dating (tritium, CFC, SF6)
Chemical degradation: first-order breakdown of contaminants

Pathogen Removal in Bank Filtration#

For pathogen removal during soil passage, total removal consists of two distinct mechanisms (Schijven and Hassanizadeh, 2000):

Inactivation (time-dependent): Pathogens lose infectivity over time through biological decay. This follows first-order kinetics and is modeled by this module as LR_decay = log10_decay_rate * residence_time. The inactivation rate depends strongly on temperature and pathogen type.
Attachment (geometry-dependent): Pathogens are physically removed by adsorption to soil grains and straining. This depends on aquifer geometry, distance, soil properties, and pH, and is NOT modeled by this module. Users should add this component separately based on site-specific data.

Total log removal = LR_decay (this module) + LR_attachment (user-specified).

At the Castricum dune recharge site, Schijven et al. (1999) found that attachment contributed approximately 97% of total MS2 removal, with inactivation contributing only 3%. Inactivation rates for common model viruses at 10 degrees C are typically 0.02-0.11 log10/day (Schijven and Hassanizadeh, 2000, Table 7).

Available functions:

residence_time_to_log_removal() - Calculate log removal from residence times and decay rate coefficient. Uses formula: Log Removal = log10_decay_rate * residence_time. Handles single values, 1D arrays, or multi-dimensional arrays of residence times. Returns log removal values with same shape as input.
decay_rate_to_log10_decay_rate() - Convert a natural-log decay rate constant lambda [1/day] to a log10 decay rate mu [log10/day].
log10_decay_rate_to_decay_rate() - Convert a log10 decay rate mu [log10/day] to a natural-log decay rate constant lambda [1/day].
parallel_mean() - Calculate weighted average log removal for parallel flow systems. Computes overall efficiency when multiple treatment paths operate in parallel with different log removal values and flow fractions. Uses formula: Total Log Removal = -log10(sum(F_i * 10^(-LR_i))) where F_i is flow fraction and LR_i is log removal for path i. Supports multi-dimensional arrays via axis parameter for batch processing. Assumes equal flow distribution if flow_fractions not provided.
gamma_pdf() - Compute probability density function (PDF) of log removal given gamma-distributed residence time. Since R = mu*T and T ~ Gamma(alpha, beta), R follows a Gamma(alpha, mu*beta) distribution.
gamma_cdf() - Compute cumulative distribution function (CDF) of log removal given gamma-distributed residence time. Returns probability that log removal is less than or equal to specified values.
gamma_mean() - Compute effective (parallel) mean log removal for gamma-distributed residence time. Uses the moment generating function of the gamma distribution to compute the log-weighted average: LR_eff = alpha * log10(1 + beta * mu * ln(10)).
gamma_find_flow_for_target_mean() - Find flow rate that produces specified target effective mean log removal given gamma-distributed aquifer pore volume. Solves inverse problem: flow = apv_beta * mu * ln(10) / (10^(target_mean / apv_alpha) - 1).

This file is part of gwtransport which is released under AGPL-3.0 license. See the ./LICENSE file or go to gwtransport/gwtransport for full license details.

gwtransport.logremoval.residence_time_to_log_removal(*, residence_times, log10_decay_rate)[source]#

Compute log removal given residence times and a log10 decay rate.

This function calculates the log removal based on residence times and a log10 decay rate coefficient using first-order decay:

Log Removal = log10_decay_rate * residence_time

This corresponds to exponential decay of pathogen concentration: C_out/C_in = 10^(-log10_decay_rate * residence_time).

Parameters:

residence_times (TypeAliasType) – Array of residence times in days. Must be positive values.
log10_decay_rate (float) – Log10 decay rate coefficient (log10/day). Relates residence time to log removal efficiency via first-order decay.

Returns:

log_removals – Array of log removal values corresponding to the input residence times. Same shape as input residence_times.

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

See also

decay_rate_to_log10_decay_rate: Convert natural-log decay rate to log10 decay rate
log10_decay_rate_to_decay_rate: Convert log10 decay rate to natural-log decay rate
gamma_mean: Compute mean log removal for gamma-distributed residence times
gamma_find_flow_for_target_mean: Find flow rate to achieve target log removal
parallel_mean: Calculate weighted average for parallel flow systems
gwtransport.residence_time.residence_time: Compute residence times from flow and pore volume
Residence Time: Time in aquifer determines pathogen contact time

Notes

Log removal is a logarithmic measure of pathogen reduction: - Log 1 = 90% reduction - Log 2 = 99% reduction - Log 3 = 99.9% reduction

The first-order decay model is mathematically identical to radioactive decay used in tracer dating. To convert a published natural-log decay rate lambda [1/day] to log10_decay_rate mu [log10/day], use decay_rate_to_log10_decay_rate().

Examples

>>> import numpy as np
>>> from gwtransport.logremoval import residence_time_to_log_removal
>>> residence_times = np.array([10.0, 20.0, 50.0])
>>> log10_decay_rate = 0.2
>>> residence_time_to_log_removal(
...     residence_times=residence_times, log10_decay_rate=log10_decay_rate
... )
array([ 2.,  4., 10.])

>>> # Single residence time
>>> residence_time_to_log_removal(residence_times=5.0, log10_decay_rate=0.3)
np.float64(1.5)

>>> # 2D array of residence times
>>> residence_times_2d = np.array([[10.0, 20.0], [30.0, 40.0]])
>>> residence_time_to_log_removal(
...     residence_times=residence_times_2d, log10_decay_rate=0.1
... )
array([[1., 2.],
       [3., 4.]])

gwtransport.logremoval.decay_rate_to_log10_decay_rate(decay_rate)[source]#

Convert a natural-log decay rate constant to a log10 decay rate.

Converts lambda [1/day] to mu [log10/day] using the relationship mu = lambda / ln(10).

Parameters:: decay_rate (float) – Natural-log first-order decay rate constant lambda (1/day). For example, from tracer dating: lambda = ln(2) / half_life.
Returns:: log10_decay_rate – Log10 decay rate mu (log10/day).
Return type:: float

See also

log10_decay_rate_to_decay_rate: Inverse conversion
residence_time_to_log_removal: Apply the log10 decay rate

Examples

>>> from gwtransport.logremoval import decay_rate_to_log10_decay_rate
>>> import numpy as np
>>> # Convert a decay rate of ln(2)/30 (half-life of 30 days)
>>> decay_rate = np.log(2) / 30
>>> decay_rate_to_log10_decay_rate(decay_rate)
0.01003...

gwtransport.logremoval.log10_decay_rate_to_decay_rate(log10_decay_rate)[source]#

Convert a log10 decay rate to a natural-log decay rate constant.

Converts mu [log10/day] to lambda [1/day] using the relationship lambda = mu * ln(10).

Parameters:: log10_decay_rate (float) – Log10 decay rate mu (log10/day).
Returns:: decay_rate – Natural-log first-order decay rate constant lambda (1/day).
Return type:: float

See also

decay_rate_to_log10_decay_rate: Inverse conversion

Examples

>>> from gwtransport.logremoval import log10_decay_rate_to_decay_rate
>>> log10_decay_rate_to_decay_rate(0.2)
0.4605...

gwtransport.logremoval.parallel_mean(*, log_removals, flow_fractions=None, axis=None)[source]#

Calculate the weighted average log removal for a system with parallel flows.

This function computes the overall log removal efficiency of a parallel filtration system. If flow_fractions is not provided, it assumes equal distribution of flow across all paths.

The calculation uses the formula:

Total Log Removal = -log10(sum(F_i * 10^(-LR_i)))

Where: - F_i = fraction of flow through system i (decimal, sum to 1.0) - LR_i = log removal of system i

Parameters:

log_removals (TypeAliasType) – Array of log removal values for each parallel flow. Each value represents the log10 reduction of pathogens. For multi-dimensional arrays, the parallel mean is computed along the specified axis.
flow_fractions (TypeAliasType | None, default: None) – Array of flow fractions for each parallel flow. Must sum to 1.0 along the specified axis and have compatible shape with log_removals. If None, equal flow distribution is assumed (default is None).
axis (int | None, default: None) – Axis along which to compute the parallel mean for multi-dimensional arrays. If None, the array is treated as 1-dimensional (default is None).

Returns:

The combined log removal value for the parallel system. If log_removals is multi-dimensional and axis is specified, returns an array with the specified axis removed.

Return type:

floating | ndarray[tuple[Any, ...], dtype[floating]]

Raises:

ValueError – If flow_fractions does not sum to 1.0 along the specified axis.

See also

residence_time_to_log_removal: Compute log removal from residence times

Notes

Log removal is a logarithmic measure of pathogen reduction:

Log 1 = 90% reduction
Log 2 = 99% reduction
Log 3 = 99.9% reduction

For parallel flows, the combined removal is typically less effective than the best individual removal but better than the worst. For systems in series, log removals would be summed directly.

Examples

>>> import numpy as np
>>> from gwtransport.logremoval import parallel_mean
>>> # Three parallel streams with equal flow and log removals of 3, 4, and 5
>>> log_removals = np.array([3, 4, 5])
>>> parallel_mean(log_removals=log_removals)
np.float64(3.431798275933005)

>>> # Two parallel streams with weighted flow
>>> log_removals = np.array([3, 5])
>>> flow_fractions = np.array([0.7, 0.3])
>>> parallel_mean(log_removals=log_removals, flow_fractions=flow_fractions)
np.float64(3.153044674980176)

>>> # Multi-dimensional array: parallel mean along axis 1
>>> log_removals_2d = np.array([[3, 4, 5], [2, 3, 4]])
>>> parallel_mean(log_removals=log_removals_2d, axis=1)
array([3.43179828, 2.43179828])

gwtransport.logremoval.gamma_pdf(*, r, rt_alpha, rt_beta, rt_loc=0.0, log10_decay_rate)[source]#

Compute the PDF of log removal given (shifted) gamma-distributed residence time.

With residence time T = T0 + rt_loc where T0 ~ Gamma(rt_alpha, rt_beta), the log removal R = mu * T follows a shifted gamma distribution with shape rt_alpha, scale mu * rt_beta, and location mu * rt_loc.

Parameters:

r (TypeAliasType) – Log removal values at which to compute the PDF.
rt_alpha (float) – Shape parameter of the gamma distribution for residence time.
rt_beta (float) – Scale parameter of the gamma distribution for residence time (days).
rt_loc (float, default: 0.0) – Location (minimum residence time, days) of the residence time distribution. Must be non-negative. Default is 0.0.
log10_decay_rate (float) – Log10 decay rate mu (log10/day). Relates residence time to log removal via R = mu * T.

Returns:

pdf – PDF values corresponding to the input r values.

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

Raises:

ValueError – If rt_loc is negative.

See also

gamma_cdf: Cumulative distribution function of log removal
gamma_mean: Mean of the log removal distribution

gwtransport.logremoval.gamma_cdf(*, r, rt_alpha, rt_beta, rt_loc=0.0, log10_decay_rate)[source]#

Compute the CDF of log removal given (shifted) gamma-distributed residence time.

With residence time T = T0 + rt_loc where T0 ~ Gamma(rt_alpha, rt_beta), the CDF is P(R <= r) = P(mu*(T0 + rt_loc) <= r) = P(T0 <= (r - mu*rt_loc)/mu) which is the CDF of a shifted gamma distribution with location mu * rt_loc.

Parameters:

r (TypeAliasType) – Log removal values at which to compute the CDF.
rt_alpha (float) – Shape parameter of the gamma distribution for residence time.
rt_beta (float) – Scale parameter of the gamma distribution for residence time (days).
rt_loc (float, default: 0.0) – Location (minimum residence time, days) of the residence time distribution. Must be non-negative. Default is 0.0.
log10_decay_rate (float) – Log10 decay rate mu (log10/day). Relates residence time to log removal via R = mu * T.

Returns:

cdf – CDF values corresponding to the input r values.

Return type:

ndarray[tuple[Any, ...], dtype[floating]]

Raises:

ValueError – If rt_loc is negative.

See also

gamma_pdf: Probability density function of log removal
gamma_mean: Mean of the log removal distribution

gwtransport.logremoval.gamma_mean(*, rt_alpha, rt_beta, rt_loc=0.0, log10_decay_rate)[source]#

Compute the effective (parallel) mean log removal for (shifted) gamma-distributed residence time.

When water travels through multiple flow paths with gamma-distributed residence times, the effective log removal is determined by mixing the output concentrations (not by averaging individual log removals). For a shifted gamma distribution T = T0 + rt_loc with T0 ~ Gamma(alpha, beta), factoring the moment generating function gives:

LR_eff = -log10(E[10^(-mu*T)]): = -log10(10^(-mu*rt_loc) * E[10^(-mu*T0)]) = mu * rt_loc + alpha * log10(1 + beta * mu * ln(10))

The rt_loc term shifts the whole log-removal distribution by a constant mu * rt_loc; the alpha/beta term is unchanged. This is always less than the arithmetic mean mu * (alpha * beta + rt_loc) because short residence time paths contribute disproportionately to the output concentration.

Parameters:

rt_alpha (float) – Shape parameter of the gamma distribution for residence time.
rt_beta (float) – Scale parameter of the gamma distribution for residence time (days).
rt_loc (float, default: 0.0) – Location (minimum residence time, days) of the residence time distribution. Must be non-negative. Default is 0.0.
log10_decay_rate (float) – Log10 decay rate mu (log10/day).

Returns:

mean – Effective (parallel) mean log removal value.

Return type:

Raises:

ValueError – If rt_loc is negative.

See also

gamma_find_flow_for_target_mean: Find flow for target mean log removal
parallel_mean: Discrete version of this calculation
gamma_pdf: PDF of the log removal distribution
gamma_cdf: CDF of the log removal distribution
The Central Concept: Pore Volume Distribution: Why residence times are distributed

gwtransport.logremoval.gamma_find_flow_for_target_mean(*, target_mean, apv_alpha, apv_beta, apv_loc=0.0, log10_decay_rate)[source]#

Find the flow rate that produces a target effective mean log removal.

Given a (shifted) gamma-distributed aquifer pore volume with parameters (apv_alpha, apv_beta, apv_loc), the residence time distribution at flow Q is a shifted gamma with shape apv_alpha, scale apv_beta/Q, and location apv_loc/Q. From gamma_mean():

LR_eff = mu * apv_loc / Q + apv_alpha * log10(1 + (apv_beta/Q) * mu * ln(10))

For apv_loc == 0 this is closed-form:

Q = apv_beta * mu * ln(10) / (10^(target_mean / apv_alpha) - 1)

For apv_loc > 0 the equation is transcendental and solved numerically with scipy.optimize.brentq() by bracketing the root in 1/Q.

Parameters:

target_mean (float) – Target effective mean log removal value. Must be positive.
apv_alpha (float) – Shape parameter of the gamma distribution for aquifer pore volume.
apv_beta (float) – Scale parameter of the gamma distribution for aquifer pore volume.
apv_loc (float, default: 0.0) – Location (minimum aquifer pore volume) of the gamma distribution. Must be non-negative. Default is 0.0.
log10_decay_rate (float) – Log10 decay rate mu (log10/day).

Returns:

flow – Flow rate (same units as apv_beta per day) that produces the target mean log removal.

Return type: