gwtransport

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.residence_time.residence_time: Compute residence times
gwtransport.diffusion.infiltration_to_extraction: Add microdispersion and molecular diffusion
Gamma Distribution Model: Two-parameter pore volume model
8. Gamma Distribution Adequacy: When gamma distribution is adequate

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 calibrating with the diffusion module, these three components are taken into account separately. When std comes from streamline analysis, it represents macrodispersion only; microdispersion and molecular diffusion can be added via the diffusion module or by combining variances (see Macrodispersion and Microdispersion in Solute Transport).

If calibrating with one compound (e.g., temperature) and predicting for another (e.g., a solute), the baked-in molecular diffusion contribution may need correction — see Macrodispersion and Microdispersion in Solute Transport. See Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity for guidance on when to add microdispersion using the diffusion module.

Examples

Basic usage with alpha and beta parameters:

>>> import pandas as pd
>>> import numpy as np
>>> from gwtransport.utils import compute_time_edges
>>> from gwtransport.advection import gamma_infiltration_to_extraction
>>>
>>> # Create input data with aligned time edges
>>> dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
>>> tedges = compute_time_edges(
...     tedges=None, tstart=None, tend=dates, number_of_bins=len(dates)
... )
>>>
>>> # Create output time edges (can be different alignment)
>>> cout_dates = pd.date_range(start="2020-01-05", end="2020-01-15", freq="D")
>>> cout_tedges = compute_time_edges(
...     tedges=None, tstart=None, tend=cout_dates, number_of_bins=len(cout_dates)
... )
>>>
>>> # Input concentration and flow (same length, aligned with tedges)
>>> cin = pd.Series(np.ones(len(dates)), index=dates)
>>> flow = pd.Series(np.ones(len(dates)) * 100, index=dates)  # 100 m3/day
>>>
>>> # Run gamma_infiltration_to_extraction with alpha/beta parameters
>>> cout = gamma_infiltration_to_extraction(
...     cin=cin,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     alpha=10.0,
...     beta=10.0,
...     n_bins=5,
... )
>>> cout.shape
(11,)

Using mean and std parameters instead:

>>> cout = gamma_infiltration_to_extraction(
...     cin=cin,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     mean=100.0,
...     std=20.0,
...     n_bins=5,
... )

With retardation factor:

>>> cout = gamma_infiltration_to_extraction(
...     cin=cin,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     alpha=10.0,
...     beta=10.0,
...     retardation_factor=2.0,  # Doubles residence time
... )

gwtransport.advection.gamma_extraction_to_infiltration(*, cout, flow, tedges, cin_tedges, alpha=None, beta=None, mean=None, std=None, n_bins=100, retardation_factor=1.0)[source]#

Compute the concentration of the infiltrating water from extracted water (deconvolution).

This function represents extraction to infiltration modeling (equivalent to deconvolution). It is symmetric to gamma_infiltration_to_extraction.

Provide either alpha and beta or mean and std.

Parameters:

cout (ArrayLike) – Concentration of the compound in extracted water or temperature of extracted water.
tedges (pandas.DatetimeIndex) – Time edges for the cout and flow data. Used to compute the cumulative concentration. Has a length of one more than cout and flow.
cin_tedges (pandas.DatetimeIndex) – Time edges for the output (infiltration) data. Used to compute the cumulative concentration. Has a length of one more than the desired output length.
flow (ArrayLike) – Flow rate of water in the aquifer [m3/day].
alpha (float, optional) – Shape parameter of gamma distribution of the aquifer pore volume (must be > 0)
beta (float, optional) – Scale parameter of gamma distribution of the aquifer pore volume (must be > 0)
mean (float, optional) – Mean of the gamma distribution.
std (float, optional) – Standard deviation of the gamma distribution.
n_bins (int) – Number of bins to discretize the gamma distribution.
retardation_factor (float, optional) – Retardation factor of the compound in the aquifer (default 1.0). Values > 1.0 indicate slower transport due to sorption/interaction.

Returns:

Concentration of the compound in the infiltrating water [ng/m3] or temperature.

Return type:

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.extraction_to_infiltration: Deconvolution with microdispersion and molecular diffusion
Gamma Distribution Model: Two-parameter pore volume model
8. Gamma Distribution Adequacy: When gamma distribution is adequate

Notes

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

Examples

Basic usage with alpha and beta parameters:

>>> import pandas as pd
>>> import numpy as np
>>> from gwtransport.utils import compute_time_edges
>>> from gwtransport.advection import gamma_extraction_to_infiltration
>>>
>>> # Create input data with aligned time edges
>>> dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
>>> tedges = compute_time_edges(
...     tedges=None, tstart=None, tend=dates, number_of_bins=len(dates)
... )
>>>
>>> # Create output time edges (can be different alignment)
>>> cin_dates = pd.date_range(start="2019-12-25", end="2020-01-15", freq="D")
>>> cin_tedges = compute_time_edges(
...     tedges=None, tstart=None, tend=cin_dates, number_of_bins=len(cin_dates)
... )
>>>
>>> # Input concentration and flow (same length, aligned with tedges)
>>> cout = pd.Series(np.ones(len(dates)), index=dates)
>>> flow = pd.Series(np.ones(len(dates)) * 100, index=dates)  # 100 m3/day
>>>
>>> # Run gamma_extraction_to_infiltration with alpha/beta parameters
>>> cin = gamma_extraction_to_infiltration(
...     cout=cout,
...     tedges=tedges,
...     cin_tedges=cin_tedges,
...     flow=flow,
...     alpha=10.0,
...     beta=10.0,
...     n_bins=5,
... )
>>> cin.shape
(22,)

Using mean and std parameters instead:

>>> cin = gamma_extraction_to_infiltration(
...     cout=cout,
...     tedges=tedges,
...     cin_tedges=cin_tedges,
...     flow=flow,
...     mean=100.0,
...     std=20.0,
...     n_bins=5,
... )

With retardation factor:

>>> cin = gamma_extraction_to_infiltration(
...     cout=cout,
...     tedges=tedges,
...     cin_tedges=cin_tedges,
...     flow=flow,
...     alpha=10.0,
...     beta=10.0,
...     retardation_factor=2.0,  # Doubles residence time
... )

gwtransport.advection.infiltration_to_extraction(*, cin, flow, tedges, cout_tedges, aquifer_pore_volumes, retardation_factor=1.0)[source]#

Compute the concentration of the extracted water using flow-weighted advection.

This function implements an infiltration to extraction advection model where cin and flow values correspond to the same aligned time bins defined by tedges.

The algorithm: 1. Computes residence times for each pore volume at cout time edges 2. Calculates infiltration time edges by subtracting residence times 3. Determines temporal overlaps between infiltration and cin time windows 4. Creates flow-weighted overlap matrices normalized by total weights 5. Computes weighted contributions and averages across pore volumes

Parameters:

cin (ArrayLike) – Concentration values of infiltrating water or temperature [concentration units]. Length must match the number of time bins defined by tedges.
flow (ArrayLike) – Flow rate values in the aquifer [m3/day]. Length must match cin and the number of time bins defined by tedges.
tedges (pandas.DatetimeIndex) – Time edges defining bins for both cin and flow data. Has length of len(cin) + 1 and len(flow) + 1.
cout_tedges (pandas.DatetimeIndex) – Time edges for output data bins. Has length of desired output + 1. Can have different time alignment and resolution than tedges.
aquifer_pore_volumes (ArrayLike) – Array of aquifer pore volumes [m3] representing the distribution of residence times in the aquifer system.
retardation_factor (float, optional) – Retardation factor of the compound in the aquifer (default 1.0). Values > 1.0 indicate slower transport due to sorption/interaction.

Returns:

Flow-weighted 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:

Raises:

ValueError – If tedges length doesn’t match cin/flow arrays plus one, or if infiltration time edges become non-monotonic (invalid input conditions).

See also

gamma_infiltration_to_extraction: Transport with gamma-distributed pore volumes
extraction_to_infiltration: Reverse operation (deconvolution)
infiltration_to_extraction_series: Simple time-shift for single pore volume
gwtransport.residence_time.residence_time: Compute residence times from flow and pore volume
gwtransport.residence_time.freundlich_retardation: Compute concentration-dependent retardation
The Central Concept: Pore Volume Distribution: Background on aquifer heterogeneity modeling
Core Transport Equation: Flow-weighted averaging approach

Examples

Basic usage with pandas Series:

>>> import pandas as pd
>>> import numpy as np
>>> from gwtransport.utils import compute_time_edges
>>> from gwtransport.advection import infiltration_to_extraction
>>>
>>> # Create input data
>>> dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
>>> tedges = compute_time_edges(
...     tedges=None, tstart=None, tend=dates, number_of_bins=len(dates)
... )
>>>
>>> # Create output time edges (different alignment)
>>> cout_dates = pd.date_range(start="2020-01-05", end="2020-01-15", freq="D")
>>> cout_tedges = compute_time_edges(
...     tedges=None, tstart=None, tend=cout_dates, number_of_bins=len(cout_dates)
... )
>>>
>>> # Input concentration and flow
>>> cin = pd.Series(np.ones(len(dates)), index=dates)
>>> flow = pd.Series(np.ones(len(dates)) * 100, index=dates)  # 100 m3/day
>>>
>>> # Define distribution of aquifer pore volumes
>>> aquifer_pore_volumes = np.array([50, 100, 200])  # m3
>>>
>>> # Run infiltration_to_extraction
>>> cout = infiltration_to_extraction(
...     cin=cin,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     aquifer_pore_volumes=aquifer_pore_volumes,
... )
>>> cout.shape
(11,)

Using array inputs instead of pandas Series:

>>> # Convert to arrays
>>> cin_values = cin.values
>>> flow_values = flow.values
>>>
>>> cout = infiltration_to_extraction(
...     cin=cin_values,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     aquifer_pore_volumes=aquifer_pore_volumes,
... )

With constant retardation factor (linear sorption):

>>> cout = infiltration_to_extraction(
...     cin=cin,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     aquifer_pore_volumes=aquifer_pore_volumes,
...     retardation_factor=2.0,  # Compound moves twice as slowly
... )

Note: For concentration-dependent retardation (nonlinear sorption), use infiltration_to_extraction_front_tracking_detailed instead, as this function only supports constant (float) retardation factors.

Using single pore volume:

>>> single_volume = np.array([100])  # Single 100 m3 pore volume
>>> cout = infiltration_to_extraction(
...     cin=cin,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     aquifer_pore_volumes=single_volume,
... )

gwtransport.advection.extraction_to_infiltration(*, cout, flow, tedges, cin_tedges, aquifer_pore_volumes, retardation_factor=1.0)[source]#

Compute the concentration of the infiltrating water from extracted water (deconvolution).

This function implements an extraction to infiltration advection model (inverse of infiltration_to_extraction) where cout and flow values correspond to the same aligned time bins defined by tedges.

SYMMETRIC RELATIONSHIP: - infiltration_to_extraction: cin + tedges → cout + cout_tedges - extraction_to_infiltration: cout + tedges → cin + cin_tedges

The algorithm (symmetric to infiltration_to_extraction): 1. Computes residence times for each pore volume at cint time edges 2. Calculates extraction time edges by adding residence times (reverse of infiltration_to_extraction) 3. Determines temporal overlaps between extraction and cout time windows 4. Creates flow-weighted overlap matrices normalized by total weights 5. Computes weighted contributions and averages across pore volumes

Parameters:

cout (ArrayLike) – Concentration values of extracted water [concentration units]. Length must match the number of time bins defined by tedges.
flow (ArrayLike) – Flow rate values in the aquifer [m3/day]. Length must match cout and the number of time bins defined by tedges.
tedges (pandas.DatetimeIndex) – Time edges defining bins for both cout and flow data. Has length of len(cout) + 1 and len(flow) + 1.
cin_tedges (pandas.DatetimeIndex) – Time edges for output (infiltration) data bins. Has length of desired output + 1. Can have different time alignment and resolution than tedges.
aquifer_pore_volumes (ArrayLike) – Array of aquifer pore volumes [m3] representing the distribution of residence times in the aquifer system.
retardation_factor (float, optional) – Retardation factor of the compound in the aquifer (default 1.0). Values > 1.0 indicate slower transport due to sorption/interaction.

Returns:

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

Return type:

Raises:

ValueError – If tedges length doesn’t match cout/flow arrays plus one, or if extraction time edges become non-monotonic (invalid input conditions).

See also

gamma_extraction_to_infiltration: Deconvolution with gamma-distributed pore volumes
infiltration_to_extraction: Forward operation (convolution)
extraction_to_infiltration_series: Simple time-shift for single pore volume
gwtransport.residence_time.residence_time: Compute residence times from flow and pore volume
gwtransport.residence_time.freundlich_retardation: Compute concentration-dependent retardation
The Central Concept: Pore Volume Distribution: Background on aquifer heterogeneity modeling
Core Transport Equation: Flow-weighted averaging approach

Examples

Basic usage with pandas Series:

>>> import pandas as pd
>>> import numpy as np
>>> from gwtransport.utils import compute_time_edges
>>> from gwtransport.advection import extraction_to_infiltration
>>>
>>> # Create input data
>>> dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
>>> tedges = compute_time_edges(
...     tedges=None, tstart=None, tend=dates, number_of_bins=len(dates)
... )
>>>
>>> # Create output time edges (different alignment)
>>> cint_dates = pd.date_range(start="2019-12-25", end="2020-01-15", freq="D")
>>> cin_tedges = compute_time_edges(
...     tedges=None, tstart=None, tend=cint_dates, number_of_bins=len(cint_dates)
... )
>>>
>>> # Input concentration and flow
>>> cout = pd.Series(np.ones(len(dates)), index=dates)
>>> flow = pd.Series(np.ones(len(dates)) * 100, index=dates)  # 100 m3/day
>>>
>>> # Define distribution of aquifer pore volumes
>>> aquifer_pore_volumes = np.array([50, 100, 200])  # m3
>>>
>>> # Run extraction_to_infiltration
>>> cin = extraction_to_infiltration(
...     cout=cout,
...     flow=flow,
...     tedges=tedges,
...     cin_tedges=cin_tedges,
...     aquifer_pore_volumes=aquifer_pore_volumes,
... )
>>> cin.shape
(22,)

Using array inputs instead of pandas Series:

>>> # Convert to arrays
>>> cout = cout.values
>>> flow = flow.values
>>>
>>> cin = extraction_to_infiltration(
...     cout=cout,
...     flow=flow,
...     tedges=tedges,
...     cin_tedges=cin_tedges,
...     aquifer_pore_volumes=aquifer_pore_volumes,
... )

With retardation factor:

>>> cin = extraction_to_infiltration(
...     cout=cout,
...     flow=flow,
...     tedges=tedges,
...     cin_tedges=cin_tedges,
...     aquifer_pore_volumes=aquifer_pore_volumes,
...     retardation_factor=2.0,  # Compound moves twice as slowly
... )

Using single pore volume:

>>> single_volume = np.array([100])  # Single 100 m3 pore volume
>>> cin = extraction_to_infiltration(
...     cout=cout,
...     flow=flow,
...     tedges=tedges,
...     cin_tedges=cin_tedges,
...     aquifer_pore_volumes=single_volume,
... )

gwtransport.advection.infiltration_to_extraction_front_tracking(*, cin, flow, tedges, cout_tedges, aquifer_pore_volumes, freundlich_k=None, freundlich_n=None, bulk_density=None, porosity=None, retardation_factor=None, max_iterations=10000)[source]#

Compute extracted concentration using exact front tracking with nonlinear sorption.

Uses event-driven analytical algorithm that tracks shock waves, rarefaction waves, and characteristics with machine precision. No numerical dispersion, exact mass balance to floating-point precision.

Parameters:

cin (ArrayLike) – Infiltration concentration [mg/L or any units]. Length = len(tedges) - 1.
flow (ArrayLike) – Flow rate [m³/day]. Must be positive. Length = len(tedges) - 1.
tedges (pandas.DatetimeIndex) – Time bin edges. Length = len(cin) + 1.
cout_tedges (pandas.DatetimeIndex) – Output time bin edges. Can be different from tedges. Length determines output array size.
aquifer_pore_volumes (ArrayLike) – Array of aquifer pore volumes [m³] representing the distribution of residence times in the aquifer system. Each pore volume must be positive.
freundlich_k (float, optional) – Freundlich coefficient [(m³/kg)^(1/n)]. Must be positive. Used if retardation_factor is None.
freundlich_n (float, optional) – Freundlich exponent [-]. Must be positive and != 1. Used if retardation_factor is None.
bulk_density (float, optional) – Bulk density [kg/m³]. Must be positive. Used if retardation_factor is None.
porosity (float, optional) – Porosity [-]. Must be in (0, 1). Used if retardation_factor is None.
retardation_factor (float, optional) – Constant retardation factor [-]. If provided, uses linear retardation instead of Freundlich sorption. Must be >= 1.0.
max_iterations (int, optional) – Maximum number of events. Default 10000.

Returns:

cout – Bin-averaged extraction concentration averaged across all pore volumes. Length = len(cout_tedges) - 1.

Return type:

Notes

Spin-up Period: The function computes the first arrival time t_first. Concentrations before t_first are affected by unknown initial conditions and should not be used for analysis. Use infiltration_to_extraction_front_tracking_detailed to access t_first.

Machine Precision: All calculations use exact analytical formulas. Mass balance is conserved to floating-point precision (~1e-14 relative error). No numerical tolerances are used for time/position calculations.

Physical Correctness: - All shocks satisfy Lax entropy condition - Rarefaction waves use self-similar solutions - Causality is strictly enforced

Examples

>>> import numpy as np
>>> import pandas as pd
>>>
>>> # Pulse injection with single pore volume
>>> tedges = pd.date_range("2020-01-01", periods=4, freq="10D")
>>> cin = np.array([0.0, 10.0, 0.0])
>>> flow = np.array([100.0, 100.0, 100.0])
>>> cout_tedges = pd.date_range("2020-01-01", periods=10, freq="5D")
>>>
>>> cout = infiltration_to_extraction_front_tracking(
...     cin=cin,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     aquifer_pore_volumes=np.array([500.0]),
...     freundlich_k=0.01,
...     freundlich_n=2.0,
...     bulk_density=1500.0,
...     porosity=0.3,
... )

With multiple pore volumes (distribution):

>>> aquifer_pore_volumes = np.array([400.0, 500.0, 600.0])
>>> cout = infiltration_to_extraction_front_tracking(
...     cin=cin,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     aquifer_pore_volumes=aquifer_pore_volumes,
...     freundlich_k=0.01,
...     freundlich_n=2.0,
...     bulk_density=1500.0,
...     porosity=0.3,
... )

See also

infiltration_to_extraction_front_tracking_detailed: Returns detailed structure
infiltration_to_extraction: Convolution-based approach for linear case
gamma_infiltration_to_extraction: For distributions of pore volumes
Non-Linear Sorption: Exact Solutions: Freundlich isotherm and front-tracking theory
1. Advection-Dominated Transport: When diffusion/dispersion is negligible

gwtransport.advection.infiltration_to_extraction_front_tracking_detailed(*, cin, flow, tedges, cout_tedges, aquifer_pore_volumes, freundlich_k=None, freundlich_n=None, bulk_density=None, porosity=None, retardation_factor=None, max_iterations=10000)[source]#

Compute extracted concentration with complete diagnostic information.

Returns both bin-averaged concentrations and detailed simulation structure for each pore volume.

Parameters:

cin (ArrayLike) – Infiltration concentration [mg/L or any units]. Length = len(tedges) - 1.
flow (ArrayLike) – Flow rate [m³/day]. Must be positive. Length = len(tedges) - 1.
tedges (pandas.DatetimeIndex) – Time bin edges. Length = len(cin) + 1.
cout_tedges (pandas.DatetimeIndex) – Output time bin edges. Can be different from tedges. Length determines output array size.
aquifer_pore_volumes (ArrayLike) – Array of aquifer pore volumes [m³] representing the distribution of residence times in the aquifer system. Each pore volume must be positive.
freundlich_k (float, optional) – Freundlich coefficient [(m³/kg)^(1/n)]. Must be positive. Used if retardation_factor is None.
freundlich_n (float, optional) – Freundlich exponent [-]. Must be positive and != 1. Used if retardation_factor is None.
bulk_density (float, optional) – Bulk density [kg/m³]. Must be positive. Used if retardation_factor is None.
porosity (float, optional) – Porosity [-]. Must be in (0, 1). Used if retardation_factor is None.
retardation_factor (float, optional) – Constant retardation factor [-]. If provided, uses linear retardation instead of Freundlich sorption. Must be >= 1.0.
max_iterations (int, optional) – Maximum number of events. Default 10000.

Return type:

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

Returns:

cout (numpy.ndarray) – Bin-averaged concentrations averaged across all pore volumes.
structures (list of dict) – List of detailed simulation structures, one for each pore volume, with keys:
- ’waves’: List[Wave] - All wave objects created during simulation
- ’events’: List[dict] - All events with times, types, and details
- ’t_first_arrival’: float - First arrival time (end of spin-up period)
- ’n_events’: int - Total number of events
- ’n_shocks’: int - Number of shocks created
- ’n_rarefactions’: int - Number of rarefactions created
- ’n_characteristics’: int - Number of characteristics created
- ’final_time’: float - Final simulation time
- ’sorption’: FreundlichSorption | ConstantRetardation - Sorption object
- ’tracker_state’: FrontTrackerState - Complete simulation state
- ’aquifer_pore_volume’: float - Pore volume for this simulation

Examples

cout, structures = infiltration_to_extraction_front_tracking_detailed(
    cin=cin,
    flow=flow,
    tedges=tedges,
    cout_tedges=cout_tedges,
    aquifer_pore_volumes=np.array([500.0]),
    freundlich_k=0.01,
    freundlich_n=2.0,
    bulk_density=1500.0,
    porosity=0.3,
)

# Access spin-up period for first pore volume
print(f"First arrival: {structures[0]['t_first_arrival']:.2f} days")

# Analyze events for first pore volume
for event in structures[0]["events"]:
    print(f"t={event['time']:.2f}: {event['type']}")

See also

infiltration_to_extraction_front_tracking: Returns concentrations only (simpler interface)
Non-Linear Sorption: Exact Solutions: Freundlich isotherm and front-tracking theory
1. Advection-Dominated Transport: When diffusion/dispersion is negligible

deposition#

Deposition Analysis for 1D Aquifer Systems.

This module analyzes compound transport by deposition in aquifer systems with tools for computing concentrations and deposition rates based on aquifer properties. The model assumes 1D groundwater flow where compound deposition occurs along the flow path, enriching the water. Deposition processes include pathogen attachment to aquifer matrix, particle filtration, or chemical precipitation. The module follows advection module patterns for consistency in forward (deposition to extraction) and reverse (extraction to deposition) calculations.

Available functions:

deposition_to_extraction() - Compute concentrations from deposition rates (convolution). Given deposition rate time series [ng/m²/day], computes resulting concentration changes in extracted water [ng/m³]. Uses flow-weighted integration over contact area between water and aquifer matrix. Accounts for aquifer geometry (porosity, thickness) and residence time distribution.
extraction_to_deposition() - Compute deposition rates from concentration changes (deconvolution). Given concentration change time series in extracted water [ng/m³], estimates deposition rate history [ng/m²/day] that produced those changes. Solves underdetermined inverse problem using nullspace regularization with configurable objectives (‘squared_differences’ for smooth solutions, ‘summed_differences’ for sparse solutions). Handles NaN values in concentration data by excluding corresponding time periods.
compute_deposition_weights() - Internal helper function. Compute weight matrix relating deposition rates to concentration changes. Used by both deposition_to_extraction (forward) and extraction_to_deposition (reverse). Calculates contact area between water parcels and aquifer matrix based on streamline geometry and residence times.
spinup_duration() - Compute spinup duration for deposition modeling. Returns residence time at first time step, representing time needed for system to become fully informed. Before this duration, extracted concentration lacks complete deposition history. Useful for determining valid analysis period and identifying when boundary effects are negligible.

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.deposition.compute_deposition_weights(*, flow_values, tedges, cout_tedges, aquifer_pore_volume, porosity, thickness, retardation_factor=1.0)[source]#

Compute deposition weights for concentration-deposition convolution.

Parameters:

flow_values (ArrayLike) – Flow rates in aquifer [m3/day]. Length must equal len(tedges) - 1.
tedges (pandas.DatetimeIndex) – Time bin edges for flow data.
cout_tedges (pandas.DatetimeIndex) – Time bin edges for output concentration data.
aquifer_pore_volume (float) – Aquifer pore volume [m3].
porosity (float) – Aquifer porosity [dimensionless].
thickness (float) – Aquifer thickness [m].
retardation_factor (float, optional) – Compound retardation factor, by default 1.0.

Returns:

Deposition weights matrix with shape (len(cout_tedges) - 1, len(tedges) - 1). May contain NaN values where residence time cannot be computed.

Return type:

Notes

The returned weights matrix may contain NaN values in locations where the residence time calculation fails or is undefined. This typically occurs when flow conditions result in invalid or non-physical residence times.

gwtransport.deposition.deposition_to_extraction(*, dep, flow, tedges, cout_tedges, aquifer_pore_volume, porosity, thickness, retardation_factor=1.0)[source]#

Compute concentrations from deposition rates (convolution).

Parameters:

dep (ArrayLike) – Deposition rates [ng/m2/day]. Length must equal len(tedges) - 1.
flow (ArrayLike) – Flow rates in aquifer [m3/day]. Length must equal len(tedges) - 1.
tedges (pandas.DatetimeIndex) – Time bin edges for deposition and flow data.
cout_tedges (pandas.DatetimeIndex) – Time bin edges for output concentration data.
aquifer_pore_volume (float) – Aquifer pore volume [m3].
porosity (float) – Aquifer porosity [dimensionless].
thickness (float) – Aquifer thickness [m].
retardation_factor (float, optional) – Compound retardation factor, by default 1.0.

Returns:

Concentration changes [ng/m3] with length len(cout_tedges) - 1.

Return type:

Examples

>>> import pandas as pd
>>> import numpy as np
>>> from gwtransport.deposition import deposition_to_extraction
>>> dates = pd.date_range("2020-01-01", "2020-01-10", freq="D")
>>> tedges = pd.date_range("2019-12-31 12:00", "2020-01-10 12:00", freq="D")
>>> cout_tedges = pd.date_range("2020-01-03 12:00", "2020-01-12 12:00", freq="D")
>>> dep = np.ones(len(dates))
>>> flow = np.full(len(dates), 100.0)
>>> cout = deposition_to_extraction(
...     dep=dep,
...     flow=flow,
...     tedges=tedges,
...     cout_tedges=cout_tedges,
...     aquifer_pore_volume=500.0,
...     porosity=0.3,
...     thickness=10.0,
... )

See also

extraction_to_deposition: Inverse operation (deconvolution)
gwtransport.advection.infiltration_to_extraction: For concentration transport without deposition
Core Transport Equation: Flow-weighted averaging approach

gwtransport.deposition.extraction_to_deposition(*, flow, tedges, cout, cout_tedges, aquifer_pore_volume, porosity, thickness, retardation_factor=1.0, nullspace_objective='squared_differences')[source]#

Compute deposition rates from concentration changes (deconvolution).

The solution for the deposition is fundamentally underdetermined, as multiple deposition histories can lead to the same concentration. This function computes a least-squares solution and then selects a specific solution from the nullspace of the problem based on the provided objective.

Parameters:

flow (ArrayLike) – Flow rates in aquifer [m3/day]. Length must equal len(tedges) - 1. Must not contain NaN values.
tedges (pandas.DatetimeIndex) – Time bin edges for deposition and flow data. Length must equal len(flow) + 1.
cout (ArrayLike) – Concentration changes in extracted water [ng/m3]. Length must equal len(cout_tedges) - 1. May contain NaN values, which will be excluded from the computation along with corresponding rows in the weight matrix.
cout_tedges (pandas.DatetimeIndex) – Time bin edges for output concentration data. Length must equal len(cout) + 1.
aquifer_pore_volume (float) – Aquifer pore volume [m3].
porosity (float) – Aquifer porosity [dimensionless].
thickness (float) – Aquifer thickness [m].
retardation_factor (float, optional) – Compound retardation factor, by default 1.0. Values > 1.0 indicate slower transport due to sorption/interaction.
nullspace_objective (str or Callable, optional) –
Objective function to minimize in the nullspace. Options:
- ”squared_differences” : Minimize sum of squared differences between adjacent deposition rates (default, provides smooth solutions)
- ”summed_differences” : Minimize sum of absolute differences between adjacent deposition rates (promotes sparse/piecewise constant solutions)
- callable : Custom objective function with signature objective(coeffs, x_ls, nullspace_basis)
Default is “squared_differences”.

Returns:

Mean deposition rates [ng/m2/day] between tedges. Length equals len(tedges) - 1.

Return type:

Raises:

ValueError – If input dimensions are incompatible, if flow contains NaN values, or if the optimization fails.

Notes

This function solves the inverse problem of determining deposition rates from observed concentration changes. Since multiple deposition histories can produce the same concentration pattern, regularization in the nullspace is used to select a physically meaningful solution.

The algorithm:

Validates input dimensions and checks for NaN values in flow
Computes deposition weight matrix relating deposition to concentration
Identifies valid rows (no NaN in weights or concentration data)
Solves the underdetermined system using nullspace regularization
Returns the regularized deposition rate solution

Examples

Basic usage with default squared differences regularization:

>>> import pandas as pd
>>> import numpy as np
>>> from gwtransport.deposition import extraction_to_deposition
>>>
>>> # Create input data
>>> dates = pd.date_range("2020-01-01", "2020-01-10", freq="D")
>>> tedges = pd.date_range("2019-12-31 12:00", "2020-01-10 12:00", freq="D")
>>> cout_tedges = pd.date_range("2020-01-03 12:00", "2020-01-12 12:00", freq="D")
>>>
>>> # Flow and concentration data
>>> flow = np.full(len(dates), 100.0)  # m3/day
>>> cout = np.ones(len(cout_tedges) - 1) * 10.0  # ng/m3
>>>
>>> # Compute deposition rates
>>> dep = extraction_to_deposition(
...     flow=flow,
...     tedges=tedges,
...     cout=cout,
...     cout_tedges=cout_tedges,
...     aquifer_pore_volume=500.0,
...     porosity=0.3,
...     thickness=10.0,
... )
>>> print(f"Deposition rates shape: {dep.shape}")
Deposition rates shape: (10,)
>>> print(f"Mean deposition rate: {np.nanmean(dep):.2f} ng/m2/day")
Mean deposition rate: 6.00 ng/m2/day

With summed differences regularization for sparse solutions:

>>> dep_sparse = extraction_to_deposition(
...     flow=flow,
...     tedges=tedges,
...     cout=cout,
...     cout_tedges=cout_tedges,
...     aquifer_pore_volume=500.0,
...     porosity=0.3,
...     thickness=10.0,
...     nullspace_objective="summed_differences",
... )

With custom regularization objective:

>>> def l2_norm_objective(coeffs, x_ls, nullspace_basis):
...     x = x_ls + nullspace_basis @ coeffs
...     return np.sum(x**2)  # Minimize L2 norm of solution
>>>
>>> dep_l2 = extraction_to_deposition(
...     flow=flow,
...     tedges=tedges,
...     cout=cout,
...     cout_tedges=cout_tedges,
...     aquifer_pore_volume=500.0,
...     porosity=0.3,
...     thickness=10.0,
...     nullspace_objective=l2_norm_objective,
... )

Handling missing concentration data:

>>> # Concentration with some NaN values
>>> cout_nan = cout.copy()
>>> cout_nan[2:4] = np.nan  # Missing data for some time periods
>>>
>>> dep_robust = extraction_to_deposition(
...     flow=flow,
...     tedges=tedges,
...     cout=cout_nan,
...     cout_tedges=cout_tedges,
...     aquifer_pore_volume=500.0,
...     porosity=0.3,
...     thickness=10.0,
... )

See also

deposition_to_extraction: Forward operation (convolution)
gwtransport.advection.extraction_to_infiltration: For concentration transport without deposition
Core Transport Equation: Flow-weighted averaging approach

gwtransport.deposition.spinup_duration(*, flow, flow_tedges, aquifer_pore_volume, retardation_factor)[source]#

Compute the spinup duration for deposition modeling.

The spinup duration is the residence time at the first time step, representing the time needed for the system to become fully informed. Before this duration, the extracted concentration lacks complete deposition history.

Parameters:

flow (numpy.ndarray) – Flow rate of water in the aquifer [m3/day].
flow_tedges (pandas.DatetimeIndex) – Time edges for the flow data.
aquifer_pore_volume (float) – Pore volume of the aquifer [m3].
retardation_factor (float) – Retardation factor of the compound in the aquifer [dimensionless].

Returns:

Spinup duration in days.

Return type:

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: 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. 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.
extraction_to_infiltration() - NOT YET IMPLEMENTED. Inverse diffusion is numerically unstable and requires regularization techniques. Placeholder for future implementation.
compute_scaled_sigma_array() - Calculate position-dependent diffusion parameters. Computes standard deviation (sigma) for Gaussian smoothing at each time step based on residence time, diffusivity, and spatial discretization: sigma = sqrt(2 * diffusivity * residence_time) / dx.
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.
deconvolve_diffusion() - NOT YET IMPLEMENTED. Inverse filtering placeholder for future diffusion deconvolution with required regularization for stability.
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, aquifer_pore_volume, diffusivity=0.1, retardation_factor=1.0, aquifer_length=80.0)[source]#

Compute the diffusion of a compound during 1D transport in the aquifer.

This function represents infiltration to extraction modeling (equivalent to convolution). It provides a fast approximation using Gaussian smoothing. The approximation is accurate when flow and tedges are relatively constant. Under variable flow conditions, errors increase but mass balance is preserved.

For physically rigorous solutions that handle variable flow correctly, use gwtransport.diffusion.infiltration_to_extraction() instead. That function is slower but provides analytical solutions to the advection-dispersion equation.

Parameters:

cin (ArrayLike) – Concentration or temperature of the compound in the infiltrating water [ng/m3].
flow (ArrayLike) – Flow rate of water in the aquifer [m3/day].
tedges (pandas.DatetimeIndex) – Time edges corresponding to the flow values.
aquifer_pore_volume (float) – Pore volume of the aquifer [m3].
diffusivity (float, optional) – diffusivity of the compound in the aquifer [m2/day]. Default is 0.1.
retardation_factor (float, optional) – Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0.
aquifer_length (float, optional) – Length of the aquifer [m]. Default is 80.0.
porosity (float, optional) – Porosity of the aquifer [dimensionless]. Default is 0.35.

Returns:

Concentration of the compound in the extracted water [ng/m3].

Return type:

See also

gwtransport.diffusion.infiltration_to_extraction: Physically rigorous analytical solution (slower)
Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity: Macrodispersion vs microdispersion

Notes

Common values for heat in saturated sand in m²/day:

Lower end (finer sand/silt): ~0.007-0.01 m²/day
Typical saturated sand: ~0.01-0.05 m²/day
Upper end (coarse sand/gravel): ~0.05-0.10 m²/day

gwtransport.diffusion_fast.extraction_to_infiltration(*, cout, flow, aquifer_pore_volume, diffusivity=0.1, retardation_factor=1.0, aquifer_length=80.0, porosity=0.35)[source]#

Compute the reverse diffusion of a compound during 1D transport in the aquifer.

This function represents extraction to infiltration modeling (equivalent to deconvolution).

Parameters:

cout (ArrayLike) – Concentration or temperature of the compound in the extracted water [ng/m3].
flow (ArrayLike) – Flow rate of water in the aquifer [m3/day].
aquifer_pore_volume (float) – Pore volume of the aquifer [m3].
diffusivity (float, optional) – diffusivity of the compound in the aquifer [m2/day]. Default is 0.1.
retardation_factor (float, optional) – Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0.
aquifer_length (float, optional) – Length of the aquifer [m]. Default is 80.0.
porosity (float, optional) – Porosity of the aquifer [dimensionless]. Default is 0.35.

Returns:

Concentration of the compound in the infiltrating water [ng/m3].

Return type:

See also

gwtransport.diffusion.extraction_to_infiltration: Analytically correct deconvolution

Notes

Extraction to infiltration diffusion (deconvolution) is mathematically ill-posed and requires regularization to obtain a stable solution.

gwtransport.diffusion_fast.compute_diffusive_spreading_length(*, flow, tedges, aquifer_pore_volume, diffusivity=0.1, retardation_factor=1.0)[source]#

Compute the physical diffusive spreading length for each time step.

The diffusive spreading length L_diff = sqrt(2 * D * rt) [m] represents the physical distance over which concentrations spread due to diffusion during the residence time rt.

Parameters:

flow (ArrayLike) – Flow rate of water in the aquifer [m3/day].
tedges (pandas.DatetimeIndex) – Time edges corresponding to the flow values.
aquifer_pore_volume (float) – Pore volume of the aquifer [m3].
diffusivity (float, optional) – Diffusivity of the compound in the aquifer [m2/day]. Default is 0.1.
retardation_factor (float, optional) – Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0.

Returns:

Array of diffusive spreading lengths [m]. Each value corresponds to a time step in the input flow series.

Return type:

gwtransport.diffusion_fast.compute_scaled_sigma_array(*, flow, tedges, aquifer_pore_volume, diffusivity=0.1, retardation_factor=1.0, aquifer_length=80.0)[source]#

Compute scaled sigma values for diffusion based on flow and aquifer properties.

Sigma represents the dimensionless spreading parameter for Gaussian filtering, expressed in units of array indices (time steps). It determines how many neighboring time steps are blended together when applying diffusive smoothing.

The computation follows these steps:

Calculate residence time (rt) for water parcels traveling through the aquifer
Compute the diffusive spreading length: L_diff = sqrt(2 * D * rt) [m]. This is the physical distance over which concentrations spread due to diffusion.
Compute the advective step size: dx = (Q * dt / V_pore) * L_aquifer [m]. This is the physical distance the water moves during one time step.
Sigma = L_diff / dx converts the physical spreading into array index units.

Why divide by dx? The Gaussian filter operates on array indices, not physical distances. If the diffusive spreading is 10 meters and each time step moves water 2 meters, then sigma = 10/2 = 5 means the filter should blend across ~5 time steps. This normalization accounts for variable flow rates: faster flow means larger dx, so fewer time steps are blended (smaller sigma), even though the physical spreading remains the same.

Parameters:

flow (ArrayLike) – Flow rate of water in the aquifer [m3/day].
tedges (pandas.DatetimeIndex) – Time edges corresponding to the flow values.
aquifer_pore_volume (float) – Pore volume of the aquifer [m3].
diffusivity (float, optional) – Diffusivity of the compound in the aquifer [m2/day]. Default is 0.1.
retardation_factor (float, optional) – Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0.
aquifer_length (float, optional) – Length of the aquifer [m]. Default is 80.0.

Returns:

Array of sigma values (in units of array indices), clipped to range [0, 100]. Each value corresponds to a time step in the input flow series.

Return type:

See also

gwtransport.diffusion.infiltration_to_extraction: For analytical solutions without this approximation

gwtransport.diffusion_fast.convolve_diffusion(*, input_signal, sigma_array, truncate=4.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.

Parameters:

input_signal (numpy.ndarray) – One-dimensional input array to be filtered.
sigma_array (numpy.ndarray) – One-dimensional array of standard deviation values, must have same length as input_signal. Each value specifies the Gaussian kernel width at the corresponding position.
truncate (float, optional) – Truncate the filter at this many standard deviations. Default is 4.0.

Returns:

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

Return type:

Notes

At the boundaries, the outer values are repeated to avoid edge effects. Equal to mode=`nearest` in scipy.ndimage.gaussian_filter1d.

The function constructs a sparse convolution matrix where each row represents a position-specific Gaussian kernel. The kernel width adapts to local sigma values, making it suitable for problems with varying diffusivitys or time steps.

For diffusion problems, the local sigma values can be calculated as: sigma = sqrt(2 * diffusivity * dt) / dx where diffusivity is the diffusivity, dt is the time step, and dx is the spatial step size.

The implementation uses sparse matrices for memory efficiency when dealing with large signals or when sigma values vary significantly.

See also

scipy.ndimage.gaussian_filter1d: Fixed-sigma Gaussian filtering
scipy.sparse: Sparse matrix implementations

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  # diffusivity
>>> dt = 0.001 * (1 + np.sin(2 * np.pi * x / 10))  # Varying time steps
>>> 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.deconvolve_diffusion(*, output_signal, sigma_array, truncate=4.0)[source]#

Apply Gaussian deconvolution with position-dependent sigma values.

Parameters:

output_signal (numpy.ndarray) – One-dimensional input array to be filtered.
sigma_array (numpy.ndarray) – One-dimensional array of standard deviation values, must have same length as output_signal. Each value specifies the Gaussian kernel width at the corresponding position.
truncate (float, optional) – Truncate the filter at this many standard deviations. Default is 4.0.

Returns:

The filtered output signal. Has the same shape as output_signal.

Return type:

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

Create example data for demonstrating variable-sigma diffusion.

Parameters:

nx (int, optional) – Number of spatial points. Default is 1000.
domain_length (float, optional) – Domain length. Default is 10.0.
diffusivity (float, optional) – diffusivity. Default is 0.1.

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).
sigma_array (numpy.ndarray) – Array of sigma values varying in space.
dt (numpy.ndarray) – Array of time steps varying in space.

Notes

This function creates a test case with:

A signal composed of two Gaussian peaks
Sinusoidally varying time steps
Corresponding sigma values for diffusion

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 function: - infiltration_to_extraction: Main transport function combining advection and dispersion

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

D_L = D_m + alpha_L * v

where: - D_m is the molecular diffusion coefficient [m^2/day] - alpha_L is the longitudinal dispersivity [m] - v is the pore velocity [m/day], computed as v = L / tau_mean

The velocity v is computed per (pore_volume, output_bin) using the mean residence time, which correctly accounts 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

where v = L / tau_mean is the mean pore velocity computed from the mean residence time for each (pore_volume, output_bin) combination. This formulation correctly captures that: - Molecular diffusion spreading scales with residence time: sqrt(D_m * tau) - Mechanical dispersion spreading scales with travel distance: sqrt(alpha_L * L)

Parameters:

cin (ArrayLike) – Concentration of the compound in infiltrating water [concentration units]. Length must match the number of time bins defined by tedges.
flow (ArrayLike) – Flow rate of water in the aquifer [m3/day]. Length must match cin and the number of time bins defined by tedges.
tedges (pandas.DatetimeIndex) – Time edges defining bins for both cin and flow data. Has length of len(cin) + 1.
cout_tedges (pandas.DatetimeIndex) – Time edges for output data bins. Has length of desired output + 1. The output concentration is averaged over each bin.
aquifer_pore_volumes (ArrayLike) – 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 (ArrayLike) – Array of travel distances [m] corresponding to each pore volume. Must have the same length as aquifer_pore_volumes.
molecular_diffusivity (float or ArrayLike) – Molecular diffusion coefficient [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. Represents Brownian motion of solute molecules, independent of velocity.
longitudinal_dispersivity (float or ArrayLike) – 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, optional) – Retardation factor of the compound in the aquifer (default 1.0). Values > 1.0 indicate slower transport due to sorption.
suppress_dispersion_warning (bool, optional) – If True, suppress the warning when using multiple pore volumes with non-zero longitudinal_dispersivity. Default is False.
asymptotic_cutoff_sigma (float or None, optional) – 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:

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, cin_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 infiltration concentration from extracted water (deconvolution with dispersion).

This function implements the inverse of infiltration_to_extraction, reconstructing the original infiltration concentration from the extracted water concentration. It explicitly constructs the weights matrix (rather than inverting the forward matrix) to ensure numerical stability and proper handling of dispersion effects.

The physical model assumes: 1. Water is extracted with concentration cout at time t_out 2. Water traveled distance L through aquifer with residence time tau = V_pore / Q 3. During transport, longitudinal dispersion caused spreading 4. At infiltration, the concentration is reconstructed from the diffused signal

The longitudinal dispersion coefficient D_L is computed internally as:

D_L = D_m + alpha_L * v

where v = L / tau_mean is the mean pore velocity computed from the mean residence time for each (pore_volume, cin_bin) combination.

Parameters:

cout (ArrayLike) – Concentration of the compound in extracted water [concentration units]. Length must match the number of time bins defined by tedges.
flow (ArrayLike) – Flow rate of water in the aquifer [m3/day]. Length must match cout and the number of time bins defined by tedges.
tedges (pandas.DatetimeIndex) – Time edges defining bins for both cout and flow data (extraction times). Has length of len(cout) + 1.
cin_tedges (pandas.DatetimeIndex) – Time edges for output infiltration data bins. Has length of desired output + 1. The output concentration is averaged over each bin.
aquifer_pore_volumes (ArrayLike) – 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 (ArrayLike) – Array of travel distances [m] corresponding to each pore volume. Must have the same length as aquifer_pore_volumes.
molecular_diffusivity (float or ArrayLike) – Molecular diffusion coefficient [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.
longitudinal_dispersivity (float or ArrayLike) – 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, optional) – Retardation factor of the compound in the aquifer (default 1.0). Values > 1.0 indicate slower transport due to sorption.
suppress_dispersion_warning (bool, optional) – If True, suppress the warning when using multiple pore volumes with non-zero longitudinal_dispersivity. Default is False.
asymptotic_cutoff_sigma (float or None, optional) – 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 infiltrating water. Same units as cout. Length equals len(cin_tedges) - 1. NaN values indicate time periods with no valid contributions from the extraction data.

Return type:

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 constructs a coefficient matrix W where cin = W @ cout:

For each pore volume, compute the deconvolution contribution: - extraction_times: when water infiltrated at cin_tedges will be extracted - 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 extraction time edge, compute the erf response at all infiltration 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 cout[j] that originated from cin[i].
Average coefficients across all pore volumes, using nan_to_num to prevent NaN propagation when different pore volumes have different valid time ranges.

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 extraction_to_infiltration
>>>
>>> # Create time edges for extraction data
>>> tedges = pd.date_range(start="2020-01-05", end="2020-01-25", freq="D")
>>> cin_tedges = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
>>>
>>> # Extracted concentration and constant flow
>>> cout = np.zeros(len(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,
...     cin_tedges=cin_tedges,
...     aquifer_pore_volumes=aquifer_pore_volumes,
...     streamline_length=streamline_length,
...     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, optional) – New flow rate for FLOW_CHANGE events [m³/day]
boundary_type (str, optional) – 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).

__repr__()[source]#: Return string representation of Event.

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 of float or 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 of float or 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 of float or 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 of tuple

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 of ShockWave or RarefactionWave

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 of ShockWave

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:

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) – 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:

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.

SIMPLIFIED IMPLEMENTATION: Just deactivates characteristic. See FRONT_TRACKING_REBUILD_PLAN.md “Known Issues and Future Improvements” Medium Priority #5.

Parameters:

raref (RarefactionWave) – Rarefaction wave
char (CharacteristicWave) – Characteristic wave
t_event (float) – Time of collision [days]
v_event (float) – Position of collision [m³]
boundary_type (str) – Which boundary collided: ‘head’ or ‘tail’

Returns:

List of new waves created (currently always empty)

Return type:

Notes

This is a simplified implementation that deactivates the characteristic without modifying the rarefaction structure.

Current implementation: deactivates characteristic, leaves rarefaction unchanged.

Future enhancement: Should modify rarefaction head/tail concentration to properly represent the wave structure instead of just absorbing the characteristic.

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) – 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

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

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

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:

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 (FreundlichSorption or ConstantRetardation) – Sorption parameters
v_inlet (float, optional) – 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 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.FreundlichSorption(k_f, n, bulk_density, porosity, c_min=1e-12)[source]#

Bases: object

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, optional) – 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).

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

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.

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.

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 or ArrayLike) – Dissolved concentration [mass/volume]. Non-negative.
Returns:: r – Retardation factor [-]. Always >= 1.0.
Return type:: float or numpy.ndarray

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 or ArrayLike) – Dissolved concentration [mass/volume]. Non-negative.
Returns:: c_total – Total concentration [mass/volume]. Always >= c.
Return type:: float or numpy.ndarray

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 or ArrayLike) – Retardation factor [-]. Must be >= 1.0.
Returns:: c – Dissolved concentration [mass/volume]. Non-negative.
Return type:: float or numpy.ndarray

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

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.

For physical shocks with n > 1 (higher C travels faster): - c_left > c_right (concentration decreases across shock) - The shock velocity is between the characteristic velocities

Examples

>>> sorption = FreundlichSorption(
...     k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
... )
>>> v_shock = sorption.shock_velocity(c_left=10.0, c_right=2.0, flow=100.0)
>>> v_shock > 0
True

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

Notes

Shocks that violate the entropy condition are unphysical and should be replaced by rarefaction waves. The entropy condition prevents non-physical expansion shocks.

For n > 1 (higher C travels faster): - Physical shocks have c_left > c_right - Characteristic from left is faster: λ(c_left) > λ(c_right) - Shock velocity is between them

For n < 1 (lower C travels faster): - Physical shocks have c_left < c_right - Characteristic from left is slower: λ(c_left) < λ(c_right) - Shock velocity is still between them

Examples

>>> sorption = FreundlichSorption(
...     k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
... )
>>> # Physical shock (compression for n>1)
>>> v_shock = sorption.shock_velocity(10.0, 2.0, 100.0)
>>> sorption.check_entropy_condition(10.0, 2.0, v_shock, 100.0)
True
>>> # Unphysical shock (expansion for n>1)
>>> v_shock_bad = sorption.shock_velocity(2.0, 10.0, 100.0)
>>> sorption.check_entropy_condition(2.0, 10.0, v_shock_bad, 100.0)
False

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

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

__init__(retardation_factor)#

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

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 (FreundlichSorption or 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 (FreundlichSorption or 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 (numpy.ndarray) – Inlet concentration [mass/volume]. Length = len(tedges) - 1.
flow (numpy.ndarray) – Flow rate [volume/time]. Length = len(tedges) - 1.
tedges (pandas.DatetimeIndex) – Time bin edges. Length = len(cin) + 1. Expected to be DatetimeIndex.
aquifer_pore_volume (float) – Total pore volume [volume]. Must be positive.
sorption (FreundlichSorption or 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:

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

See also

compute_first_fully_informed_bin_edge: Get first valid output bin edge

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 (numpy.ndarray) – Inlet concentration [mass/volume]. Length = len(tedges) - 1.
flow (numpy.ndarray) – Flow rate [volume/time]. Length = len(tedges) - 1.
tedges (pandas.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 (FreundlichSorption or ConstantRetardation) – Sorption model.
output_tedges (pandas.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:

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

See also

compute_first_front_arrival_time: Get exact arrival time

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 (list of Wave) – All waves in the simulation (active and inactive).
sorption (FreundlichSorption or 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 (ArrayLike) – Time points [days]. Must be sorted in ascending order.
v_outlet (float) – Outlet position [m³].
waves (list of Wave) – All waves in the simulation.
sorption (FreundlichSorption or ConstantRetardation) – Sorption model.

Returns:

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

Return type:

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)

See also

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

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 (list of Wave) – All waves in the simulation.
sorption (FreundlichSorption or 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 of 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.

For Freundlich sorption, the concentration within a rarefaction fan varies as:: R(C) = flow*(t - t_origin)/(v_outlet - v_origin)
This function computes the exact integral:: ∫_{t_start}^{t_end} C(t) dt

where C(t) is obtained by inverting the retardation relation.

Parameters:

raref (RarefactionWave) – Rarefaction wave controlling the outlet.
v_outlet (float) – Outlet position [m³].
t_start (float) – Integration start time [days]. Can be -np.inf.
t_end (float) – Integration end time [days]. Can be np.inf.
sorption (FreundlichSorption) – Freundlich sorption model.

Returns:

integral – Exact integral value [concentration * time].

Return type:

Notes

Derivation:

For Freundlich: R(C) = 1 + alpha*C^beta where:: alpha = rho_b*k_f/(n_por*n) beta = 1/n - 1
At outlet: R = kappa*t + mu where:: kappa = flow/(v_outlet - v_origin) mu = -flow*t_origin/(v_outlet - v_origin)

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

Integral:

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

evaluated from t_start to t_end.

Special Cases: - If (kappa*t + mu - 1) <= 0, concentration is 0 (unphysical region) - For beta = 0 (n = 1), use ConstantRetardation instead - For t_end = +∞ with exponent < 0 (n>1), integral converges to 0 - For t_start = -∞, antiderivative evaluates to 0

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,
)
integral = integrate_rarefaction_exact(
    raref, v_outlet=500.0, t_start=5.0, t_end=15.0, sorption=sorption
)
integral > 0
# Integration to infinity
integral_inf = integrate_rarefaction_exact(
    raref, v_outlet=500.0, t_start=5.0, t_end=np.inf, sorption=sorption
)
integral_inf > integral

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 (ArrayLike) – Time bin edges [days]. Length N+1 for N bins.
v_outlet (float) – Outlet position [m³].
waves (list of Wave) – All waves from front tracking simulation.
sorption (FreundlichSorption or ConstantRetardation) – Sorption model.

Returns:

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

Return type:

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)

See also

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

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 (list of Wave) – All waves in the simulation.
sorption (FreundlichSorption or ConstantRetardation) – Sorption model.

Returns:

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

Return type:

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

See also

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

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 (ArrayLike) – Inlet concentration time series [mass/volume]. Piecewise constant within bins defined by tedges_days.
flow (ArrayLike) – Flow rate time series [m³/day]. Piecewise constant within bins defined by tedges_days.
tedges_days (numpy.ndarray) – 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 (list of 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 (list of Wave) – All waves in the simulation.
sorption (FreundlichSorption or ConstantRetardation) – Sorption model.
flow (ArrayLike) – Flow rate time series [m³/day]. Piecewise constant within bins defined by tedges_days.
tedges_days (numpy.ndarray) – 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 (FreundlichSorption or ConstantRetardation) – Sorption model.
flow (float) – Flow rate [m³/day] (assumed constant).

Returns:

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

Return type:

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 (list of Wave) – All waves in the simulation.
sorption (FreundlichSorption or ConstantRetardation) – Sorption model.
flow (ArrayLike) – Flow rate time series [m³/day]. Piecewise constant within bins defined by tedges_days.
tedges_days (numpy.ndarray) – Time bin edges [days]. Length len(flow) + 1.

Return type:

matplotlib.figure.Figure

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.

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]

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

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, optional) – Maximum time to plot [days]. If None, uses final simulation time * 1.2.
figsize (tuple of float, optional) – Figure size in inches (width, height). Default (14, 10).
show_inactive (bool, optional) – Whether to show inactive waves (deactivated by interactions). Default False.
show_events (bool, optional) – 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, optional) – Maximum time to plot [days]. If None, uses final simulation time * 1.1.
n_rarefaction_points (int, optional) – Number of points to use for plotting rarefaction segments (analytical curves). Default 50 per rarefaction segment.
figsize (tuple of float, optional) – Figure size in inches (width, height). Default (12, 6).
t_first_arrival (float, optional) – 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 of float, optional) – 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 of Wave) – All waves created during simulation (includes inactive waves)
events (list of 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 (FreundlichSorption or ConstantRetardation) – Sorption parameters
cin (numpy.ndarray) – Inlet concentration time series [mass/volume]
flow (numpy.ndarray) – Flow rate time series [m³/day]
tedges (pandas.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: FreundlichSorption | 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 (numpy.ndarray) – Inlet concentration time series [mass/volume]
flow (numpy.ndarray) – Flow rate time series [m³/day]
tedges (numpy.ndarray) – Time bin edges [days]
aquifer_pore_volume (float) – Total pore volume [m³]
sorption (FreundlichSorption or ConstantRetardation) – Sorption parameters

state#

Complete simulation state

Type:: FrontTrackerState

t_first_arrival#

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

Type:: float

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)}")

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.

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

Initialize tracker with inlet conditions and physical parameters.

Parameters:

cin (numpy.ndarray) – Inlet concentration time series [mass/volume]
flow (numpy.ndarray) – Flow rate time series [m³/day]
tedges (numpy.ndarray) – Time bin edges [days]
aquifer_pore_volume (float) – Total pore volume [m³]
sorption (FreundlichSorption or 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 or 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

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, optional) – Maximum number of events to process. Default 10000. Prevents infinite loops in case of bugs.
verbose (bool, optional) – 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, optional) – Enable mass balance verification. Default False (opt-in for now).
mass_balance_rtol (float, optional) – 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, optional) – 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 or 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 (FreundlichSorption or ConstantRetardation) – Sorption model determining velocity.
is_active (bool, optional) – 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: FreundlichSorption | 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 or None

concentration_left()[source]#

Get upstream concentration (same as concentration for characteristics).

Return type:: float

concentration_right()[source]#

Get downstream concentration (same as concentration for characteristics).

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 (FreundlichSorption or ConstantRetardation) – Sorption model.
is_active (bool, optional) – Activity status. Default True.
velocity (float, optional) – 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: FreundlichSorption | 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 or None

concentration_left()[source]#

Get upstream concentration.

Return type:: float

concentration_right()[source]#

Get downstream concentration.

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

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 (FreundlichSorption or ConstantRetardation) – Sorption model (must be concentration-dependent).
is_active (bool, optional) – 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: FreundlichSorption | ConstantRetardation#: Sorption model (must be concentration-dependent).

__post_init__()[source]#: Verify this is actually a rarefaction (head faster than tail).

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

Return type:: float

concentration_right()[source]#

Get downstream concentration (head).

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 two-parameter model 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, and the gamma distribution provides a realistic representation of this heterogeneity.

Available functions:

parse_parameters() - Parse and validate gamma distribution parameters from either (alpha, beta) or (mean, std). Ensures exactly one parameter pair is provided and validates positivity constraints.
mean_std_to_alpha_beta() - Convert physically intuitive (mean, std) parameters to gamma shape/scale parameters. Uses formulas: alpha = mean^2 / std^2 and beta = std^2 / mean.
alpha_beta_to_mean_std() - Convert gamma (alpha, beta) parameters back to (mean, std) for physical interpretation. Uses formulas: mean = alpha * beta and std = sqrt(alpha) * beta.
bins() - Primary function for transport modeling. Creates discrete probability bins from continuous 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 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(*, alpha=None, beta=None, mean=None, std=None)[source]#

Parse parameters for gamma distribution.

Either alpha and beta or mean and std must be provided.

Parameters:

alpha (float, optional) – Shape parameter of gamma distribution (must be > 0)
beta (float, optional) – Scale parameter of gamma distribution (must be > 0)
mean (float, optional) – Mean of the gamma distribution.
std (float, optional) – Standard deviation of the gamma distribution. See Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity for what std represents depending on APVD source.

Return type:

Returns:

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

Raises:

ValueError – If both alpha and beta are None or if both mean and std are None. If alpha or beta are not positive.

gwtransport.gamma.mean_std_to_alpha_beta(*, mean, std)[source]#

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

Parameters:

mean (float) – Mean of the gamma distribution.
std (float) – Standard deviation of the gamma distribution. See Macrodispersion and Microdispersion as Scale-Dependent Heterogeneity for what std represents depending on APVD source.

Return type:

Returns:

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

See also

alpha_beta_to_mean_std: Convert shape and scale parameters to mean and std
parse_parameters: Parse and validate gamma distribution parameters

Examples

>>> from gwtransport.gamma import mean_std_to_alpha_beta
>>> mean_pore_volume = 30000.0  # m³
>>> std_pore_volume = 8100.0  # m³
>>> alpha, beta = mean_std_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

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

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

Parameters:

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

Return type:

Returns:

mean (float) – Mean of the gamma distribution.
std (float) – Standard deviation of the gamma distribution.

See also

mean_std_to_alpha_beta: Convert mean and std to shape and scale parameters
parse_parameters: Parse and validate gamma distribution parameters

Examples

>>> from gwtransport.gamma import alpha_beta_to_mean_std
>>> alpha = 13.72  # shape parameter
>>> beta = 2187.0  # scale parameter
>>> mean, std = alpha_beta_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(*, alpha=None, beta=None, mean=None, std=None, n_bins=100, quantile_edges=None)[source]#

Divide gamma distribution into bins and compute various bin properties.

If n_bins is provided, the gamma distribution is divided into n_bins equal-mass bins. If quantile_edges is provided, the gamma distribution is divided into bins defined by the quantile edges. The quantile edges 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.

Parameters:

alpha (float, optional) – Shape parameter of gamma distribution (must be > 0)
beta (float, optional) – Scale parameter of gamma distribution (must be > 0)
mean (float, optional) – Mean of the gamma distribution.
std (float, optional) – Standard deviation of the gamma distribution.
n_bins (int, optional) – Number of bins to divide the gamma distribution (must be > 1). Default is 100.
quantile_edges (ArrayLike, optional) – 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 is 0)
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

Return type:

dict

See also

bin_masses: Calculate probability mass for bins
mean_std_to_alpha_beta: Convert mean/std to alpha/beta parameters
gwtransport.advection.gamma_infiltration_to_extraction: Use bins for transport modeling
Gamma Distribution Model: Two-parameter pore volume model
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
>>> # Define gamma distribution using mean and std
>>> result = bins(mean=30000.0, std=8100.0, n_bins=5)

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

Is the area under the gamma distribution PDF between the bin edges.

Parameters:

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

Returns:

Probability mass for each bin

Return type:

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 (ArrayLike) – 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:

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 (ArrayLike) – 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 (ArrayLike, optional) – 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, optional) – 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:

float or ArrayLike

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])

See also

residence_time_to_log_removal: Compute log removal from residence times

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

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

Since log removal R = mu*T and T ~ Gamma(alpha, beta), the log removal R follows a Gamma(alpha, mu*beta) distribution.

Parameters:

r (ArrayLike) – 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).
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:

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, log10_decay_rate)[source]#

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

Since log removal R = mu*T and T ~ Gamma(alpha, beta), the CDF is P(R <= r) = P(T <= r/mu) = Gamma_CDF(r/mu; alpha, beta).

Parameters:

r (ArrayLike) – 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).
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:

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, log10_decay_rate)[source]#

Compute the effective (parallel) mean log removal for 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). This uses the moment generating function of the gamma distribution:

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

This is always less than the arithmetic mean (mu * alpha * beta) 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).
log10_decay_rate (float) – Log10 decay rate mu (log10/day).

Returns:

mean – Effective (parallel) mean log removal value.

Return type:

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, log10_decay_rate)[source]#

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

Given a gamma-distributed aquifer pore volume with parameters (apv_alpha, apv_beta), the residence time distribution is Gamma(apv_alpha, apv_beta/flow). The effective mean log removal (from gamma_mean()) is:

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

Solving for flow:

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

Parameters:

target_mean (float) – Target effective mean log removal value.
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.
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:

See also

gamma_mean: Compute effective mean log removal for given parameters

residence_time#

Residence Time Calculations for Retarded Compound Transport.

This module provides functions to compute residence times for compounds traveling through aquifer systems, accounting for flow variability, pore volume, and retardation due to physical or chemical interactions with the aquifer matrix. Residence time represents the duration a compound spends traveling from infiltration to extraction points, depending on flow rate (higher flow yields shorter residence time), pore volume (larger volume yields longer residence time), and retardation factor (interaction with matrix yields longer residence time).

Available functions:

residence_time() - Compute residence times at specific time indices. Supports both forward (infiltration to extraction) and reverse (extraction to infiltration) directions. Handles single or multiple pore volumes (2D output for multiple volumes). Returns residence times in days using cumulative flow integration for accurate time-varying flow handling.
residence_time_mean() - Compute mean residence times over time intervals. Calculates average residence time between specified time edges using linear averaging to properly weight time-varying residence times. Supports same directional options as residence_time. Particularly useful for time-binned analysis.
fraction_explained() - Compute fraction of aquifer pore volumes with valid residence times. Indicates how many pore volumes have sufficient flow history to compute residence time. Returns values in [0, 1] where 1.0 means all volumes are fully informed. Useful for assessing spin-up periods and data coverage. NaN residence times indicate insufficient flow history.

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.residence_time.residence_time(*, flow, flow_tedges, aquifer_pore_volumes, index=None, direction='extraction_to_infiltration', retardation_factor=1.0)[source]#

Compute the residence time of retarded compound in the water in the aquifer.

Parameters:

flow (ArrayLike) – Flow rate of water in the aquifer [m3/day]. The length of flow should match the length of flow_tedges minus one.
flow_tedges (pandas.DatetimeIndex) – Time edges for the flow data. Used to compute the cumulative flow. Has a length of one more than flow.
aquifer_pore_volumes (float or ArrayLike of float) – Pore volume(s) of the aquifer [m3]. Can be a single value or an array of pore volumes representing different flow paths.
index (pandas.DatetimeIndex, optional) – Index at which to compute the residence time. If left to None, flow_tedges is used. Default is None.
direction ({'extraction_to_infiltration', 'infiltration_to_extraction'}, optional) – Direction of the flow calculation: * ‘extraction_to_infiltration’: Extraction to infiltration modeling - how many days ago was the extracted water infiltrated * ‘infiltration_to_extraction’: Infiltration to extraction modeling - how many days until the infiltrated water is extracted Default is ‘extraction_to_infiltration’.
retardation_factor (float, optional) – Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0.

Returns:

Residence time of the retarded compound in the aquifer [days].

Return type:

See also

residence_time_mean: Compute mean residence time over time intervals
gwtransport.advection.gamma_infiltration_to_extraction: Use residence times for transport
gwtransport.logremoval.residence_time_to_log_removal: Convert residence time to log removal
Residence Time: Time in aquifer between infiltration and extraction
Retardation Factor: Slower movement due to sorption

gwtransport.residence_time.residence_time_mean(*, flow, flow_tedges, tedges_out, aquifer_pore_volumes, direction='extraction_to_infiltration', retardation_factor=1.0)[source]#

Compute the mean residence time of a retarded compound in the aquifer between specified time edges.

This function calculates the average residence time of a retarded compound in the aquifer between specified time intervals. It can compute both extraction to infiltration modeling (extraction direction: when was extracted water infiltrated) and infiltration to extraction modeling (infiltration direction: when will infiltrated water be extracted).

The function handles time series data by computing the cumulative flow and using linear interpolation and averaging to determine mean residence times between the specified time edges.

Parameters:

flow (ArrayLike) – Flow rate of water in the aquifer [m3/day]. Should be an array of flow values corresponding to the intervals defined by flow_tedges.
flow_tedges (ArrayLike) – Time edges for the flow data, as datetime64 objects. These define the time intervals for which the flow values are provided.
tedges_out (ArrayLike) – Output time edges as datetime64 objects. These define the intervals for which the mean residence times will be calculated.
aquifer_pore_volumes (float or ArrayLike) – Pore volume(s) of the aquifer [m3]. Can be a single value or an array of values for multiple pore volume scenarios.
direction ({'extraction_to_infiltration', 'infiltration_to_extraction'}, optional) – Direction of the flow calculation: * ‘extraction_to_infiltration’: Extraction to infiltration modeling - how many days ago was the extracted water infiltrated * ‘infiltration_to_extraction’: Infiltration to extraction modeling - how many days until the infiltrated water is extracted Default is ‘extraction_to_infiltration’.
retardation_factor (float, optional) – Retardation factor of the compound in the aquifer [dimensionless]. A value greater than 1.0 indicates that the compound moves slower than water. Default is 1.0 (no retardation).

Returns:

Mean residence time of the retarded compound in the aquifer [days] for each interval defined by tedges_out. The first dimension corresponds to the different pore volumes and the second to the residence times between tedges_out.

Return type:

Notes

The function converts datetime objects to days since the start of the time series.
For extraction_to_infiltration direction, the function computes how many days ago water was infiltrated.
For infiltration_to_extraction direction, the function computes how many days until water will be extracted.
The function uses linear interpolation for computing residence times at specific points and linear averaging for computing mean values over intervals.

Examples

>>> import pandas as pd
>>> import numpy as np
>>> from gwtransport.residence_time import residence_time_mean
>>> # Create sample flow data
>>> flow_dates = pd.date_range(start="2023-01-01", end="2023-01-10", freq="D")
>>> flow_values = np.full(len(flow_dates) - 1, 100.0)  # Constant flow of 100 m³/day
>>> pore_volume = 200.0  # Aquifer pore volume in m³
>>> # Calculate mean residence times
>>> mean_times = residence_time_mean(
...     flow=flow_values,
...     flow_tedges=flow_dates,
...     tedges_out=flow_dates,
...     aquifer_pore_volumes=pore_volume,
...     direction="extraction_to_infiltration",
... )
>>> # With constant flow of 100 m³/day and pore volume of 200 m³,
>>> # mean residence time should be approximately 2 days
>>> print(mean_times)
[[nan nan  2.  2.  2.  2.  2.  2.  2.]]

See also

residence_time: Compute residence time at specific time indices
Residence Time: Time in aquifer between infiltration and extraction

gwtransport.residence_time.fraction_explained(*, rt=None, flow=None, flow_tedges=None, aquifer_pore_volumes=None, index=None, direction='extraction_to_infiltration', retardation_factor=1.0)[source]#

Compute the fraction of the aquifer that is informed with respect to the retarded flow.

Parameters:

rt (numpy.ndarray, optional) – Pre-computed residence time array [days]. If not provided, it will be computed.
flow (ArrayLike, optional) – Flow rate of water in the aquifer [m3/day]. The length of flow should match the length of flow_tedges minus one.
flow_tedges (pandas.DatetimeIndex, optional) – Time edges for the flow data. Used to compute the cumulative flow. Has a length of one more than flow. Inbetween neighboring time edges, the flow is assumed constant.
aquifer_pore_volumes (float or ArrayLike of float, optional) – Pore volume(s) of the aquifer [m3]. Can be a single value or an array of pore volumes representing different flow paths.
index (pandas.DatetimeIndex, optional) – Index at which to compute the fraction. If left to None, the index of flow is used. Default is None.
direction ({'extraction_to_infiltration', 'infiltration_to_extraction'}, optional) – Direction of the flow calculation: * ‘extraction_to_infiltration’: Extraction to infiltration modeling - how many days ago was the extracted water infiltrated * ‘infiltration_to_extraction’: Infiltration to extraction modeling - how many days until the infiltrated water is extracted Default is ‘extraction_to_infiltration’.
retardation_factor (float, optional) – Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0.

Returns:

Fraction of the aquifer that is informed with respect to the retarded flow.

Return type:

gwtransport.residence_time.freundlich_retardation(*, concentration, freundlich_k, freundlich_n, bulk_density, porosity)[source]#

Compute concentration-dependent retardation factors using Freundlich isotherm.

The Freundlich isotherm relates sorbed concentration S to aqueous concentration C:: S = rho_f * C^n
The retardation factor is computed as:: R = 1 + (rho_b/θ) * dS/dC = 1 + (rho_b/θ) * rho_f * n * C^(n-1)

Parameters:

concentration (ArrayLike) – Concentration of compound in water [mass/volume]. Length should match flow (i.e., len(flow_tedges) - 1).
freundlich_k (float) – Freundlich coefficient [(m³/kg)^(1/n)].
freundlich_n (float) – Freundlich sorption exponent [dimensionless].
bulk_density (float) – Bulk density of aquifer material [mass/volume].
porosity (float) – Porosity of aquifer [dimensionless, 0-1].

Returns:

Retardation factors for each flow interval. Length equals len(concentration) for use as retardation_factor in residence_time.

Return type:

Examples

>>> concentration = np.array([0.1, 0.2, 0.3])  # same length as flow
>>> R = freundlich_retardation(
...     concentration=concentration,
...     freundlich_k=0.5,
...     freundlich_n=0.9,
...     bulk_density=1600,  # kg/m3
...     porosity=0.35,
... )
>>> # Use R in residence_time as retardation_factor

See also

residence_time: Compute residence times with retardation
gwtransport.advection.infiltration_to_extraction_front_tracking: Transport with nonlinear sorption
Non-Linear Sorption: Exact Solutions: Freundlich isotherm and concentration-dependent retardation

surfacearea#

Surface Area Calculations for Streamline-Based Transport Analysis.

This module provides geometric utilities for computing surface areas and average heights between streamlines in heterogeneous aquifer systems. These calculations support direct estimation of pore volume distributions from streamline analysis, providing an alternative to gamma distribution approximations.

Available functions:

compute_average_heights() - Compute average heights of clipped trapezoids formed by streamlines. Trapezoids have vertical sides defined by x_edges and top/bottom edges defined by y_edges (2D array). Clipping bounds (y_lower, y_upper) restrict the integration domain. Returns area/width ratios representing average heights for use in pore volume calculations from streamline geometry.

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.surfacearea.compute_average_heights(*, x_edges, y_edges, y_lower, y_upper)[source]#

Compute average heights of clipped trapezoids.

Trapezoids have vertical left and right sides, with corners at: - top-left: (y_edges[i, j], x_edges[j]) - top-right: (y_edges[i, j+1], x_edges[j+1]) - bottom-left: (y_edges[i+1, j], x_edges[j]) - bottom-right: (y_edges[i+1, j+1], x_edges[j+1])

The area is computed as the exact integral of clip(top(x), y_lower, y_upper) - clip(bottom(x), y_lower, y_upper) where top(x) and bottom(x) are linear interpolants of the corner values, and clip restricts values to [y_lower, y_upper].

Parameters:

x_edges (numpy.ndarray) – 1D array of x coordinates, shape (n_x,)
y_edges (numpy.ndarray) – 2D array of y coordinates, shape (n_y, n_x)
y_lower (float) – Lower horizontal clipping bound
y_upper (float) – Upper horizontal clipping bound

Returns:

avg_heights – 2D array of average heights (area/width) for each clipped trapezoid, shape (n_y-1, n_x-1)

Return type:

See also

gwtransport.advection.infiltration_to_extraction: Use computed areas in transport calculations

Examples

>>> import numpy as np
>>> from gwtransport.surfacearea import compute_average_heights
>>> # Create simple grid
>>> x_edges = np.array([0.0, 1.0, 2.0])
>>> y_edges = np.array([[0.0, 0.5, 1.0], [1.0, 1.5, 2.0], [2.0, 2.5, 3.0]])
>>> # Compute average heights with clipping bounds
>>> avg_heights = compute_average_heights(
...     x_edges=x_edges, y_edges=y_edges, y_lower=0.5, y_upper=2.5
... )
>>> print(f"Shape: {avg_heights.shape}")
Shape: (2, 2)

utils#

General Utilities for 1D Groundwater Transport Modeling.

This module provides general-purpose utility functions for time series manipulation, interpolation, numerical operations, and data processing used throughout the gwtransport package. Functions include linear interpolation/averaging, bin overlap calculations, underdetermined system solvers, and external data retrieval.

Available functions:

linear_interpolate() - Linear interpolation using numpy’s optimized interp function. Automatically handles unsorted data with configurable extrapolation (None for clamping, float for constant values). Handles multi-dimensional query arrays.
interp_series() - Interpolate pandas Series to new DatetimeIndex using scipy.interpolate.interp1d. Automatically filters NaN values and converts datetime to numerical representation.
linear_average() - Compute average values of piecewise linear time series between specified x-edges. Supports 1D or 2D edge arrays for batch processing. Handles NaN values and offers multiple extrapolation methods (‘nan’, ‘outer’, ‘raise’).
diff() - Compute cell widths from cell coordinate arrays with configurable alignment (‘centered’, ‘left’, ‘right’). Returns widths matching input array length.
partial_isin() - Calculate fraction of each input bin overlapping with each output bin. Returns dense matrix where element (i,j) represents overlap fraction. Uses vectorized operations for efficiency.
time_bin_overlap() - Calculate fraction of time bins overlapping with specified time ranges. Similar to partial_isin but for time-based bin overlaps with list of (start, end) tuples.
combine_bin_series() - Combine two binned series onto common set of unique edges. Maps values from original bins to new combined structure with configurable extrapolation (‘nearest’ or float value).
compute_time_edges() - Compute DatetimeIndex of time bin edges from explicit edges, start times, or end times. Validates consistency with expected number of bins and handles uniform spacing extrapolation.
solve_underdetermined_system() - Solve underdetermined linear system (Ax = b, m < n) with nullspace regularization. Handles NaN values by row exclusion. Supports built-in objectives (‘squared_differences’, ‘summed_differences’) or custom callable objectives.
get_soil_temperature() - Download soil temperature data from KNMI weather stations with automatic caching. Supports stations 260 (De Bilt), 273 (Marknesse), 286 (Nieuw Beerta), 323 (Wilhelminadorp). Returns DataFrame with columns TB1-TB5, TNB1-TNB2, TXB1-TXB2 at various depths. Daily cache prevents redundant downloads.
generate_failed_coverage_badge() - Generate SVG badge indicating failed coverage using genbadge library. Used in CI/CD workflows.

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.utils.linear_interpolate(*, x_ref, y_ref, x_query, left=None, right=None)[source]#

Linear interpolation using numpy’s optimized interp function.

Automatically handles unsorted reference data by sorting it first.

Parameters:

x_ref (ArrayLike) – Reference x-values. If unsorted, will be automatically sorted.
y_ref (ArrayLike) – Reference y-values corresponding to x_ref.
x_query (ArrayLike) – Query x-values where interpolation is needed. Array may have any shape.
left (float, optional) –
Value to return for x_query < x_ref[0].
- If left=None: clamp to y_ref[0] (default)
- If left=float: use specified value (e.g., np.nan)
right (float, optional) –
Value to return for x_query > x_ref[-1].
- If right=None: clamp to y_ref[-1] (default)
- If right=float: use specified value (e.g., np.nan)

Returns:

Interpolated y-values with the same shape as x_query.

Return type:

Examples

Basic interpolation with clamping (default):

>>> import numpy as np
>>> from gwtransport.utils import linear_interpolate
>>> x_ref = np.array([1.0, 2.0, 3.0, 4.0])
>>> y_ref = np.array([10.0, 20.0, 30.0, 40.0])
>>> x_query = np.array([0.5, 1.5, 2.5, 3.5, 4.5])
>>> linear_interpolate(x_ref=x_ref, y_ref=y_ref, x_query=x_query)
array([10., 15., 25., 35., 40.])

Using NaN for extrapolation:

>>> linear_interpolate(
...     x_ref=x_ref, y_ref=y_ref, x_query=x_query, left=np.nan, right=np.nan
... )
array([nan, 15., 25., 35., nan])

Handles unsorted reference data automatically:

>>> x_unsorted = np.array([3.0, 1.0, 4.0, 2.0])
>>> y_unsorted = np.array([30.0, 10.0, 40.0, 20.0])
>>> linear_interpolate(x_ref=x_unsorted, y_ref=y_unsorted, x_query=x_query)
array([10., 15., 25., 35., 40.])

See also

interp_series: Interpolate pandas Series with datetime index

gwtransport.utils.interp_series(*, series, index_new, **interp1d_kwargs)[source]#

Interpolate a pandas.Series to a new index.

Parameters:

series (pandas.Series) – Series to interpolate.
index_new (pandas.DatetimeIndex) – New index to interpolate to.
interp1d_kwargs (dict, optional) – Keyword arguments passed to scipy.interpolate.interp1d. Default is {}.

Returns:

Interpolated series.

Return type:

pandas.Series

gwtransport.utils.diff(*, a, alignment='centered')[source]#

Compute the cell widths for a given array of cell coordinates.

If alignment is “centered”, the coordinates are assumed to be centered in the cells. If alignment is “left”, the coordinates are assumed to be at the left edge of the cells. If alignment is “right”, the coordinates are assumed to be at the right edge of the cells.

Parameters:: a (ArrayLike) – Input array.
Returns:: Array with differences between elements.
Return type:: numpy.ndarray

gwtransport.utils.linear_average(*, x_data, y_data, x_edges, extrapolate_method='nan')[source]#

Compute the average value of a piecewise linear time series between specified x-edges.

Parameters:

x_data (ArrayLike) – x-coordinates of the time series data points, must be in ascending order
y_data (ArrayLike) – y-coordinates of the time series data points
x_edges (ArrayLike) – x-coordinates of the integration edges. Can be 1D or 2D. - If 1D: shape (n_edges,). Can be 1D or 2D. - If 1D: shape (n_edges,), must be in ascending order - If 2D: shape (n_series, n_edges), each row must be in ascending order - If 2D: shape (n_series, n_edges), each row must be in ascending order
extrapolate_method (str, optional) – Method for handling extrapolation. Default is ‘nan’. - ‘outer’: Extrapolate using the outermost data points. - ‘nan’: Extrapolate using np.nan. - ‘raise’: Raise an error for out-of-bounds values.

Returns:

2D array of average values between consecutive pairs of x_edges. Shape is (n_series, n_bins) where n_bins = n_edges - 1. If x_edges is 1D, n_series = 1.

Return type:

Examples

>>> import numpy as np
>>> from gwtransport.utils import linear_average
>>> x_data = [0, 1, 2, 3]
>>> y_data = [0, 1, 1, 0]
>>> x_edges = [0, 1.5, 3]
>>> linear_average(
...     x_data=x_data, y_data=y_data, x_edges=x_edges
... )
array([[0.666..., 0.666...]])

>>> x_edges_2d = [[0, 1.5, 3], [0.5, 2, 3]]
>>> linear_average(x_data=x_data, y_data=y_data, x_edges=x_edges_2d)
array([[0.66666667, 0.66666667],
       [0.91666667, 0.5       ]])

gwtransport.utils.partial_isin(*, bin_edges_in, bin_edges_out)[source]#

Calculate the fraction of each input bin that overlaps with each output bin.

This function computes a matrix where element (i, j) represents the fraction of input bin i that overlaps with output bin j. The computation uses vectorized operations to avoid loops.

Parameters:

bin_edges_in (ArrayLike) – 1D array of input bin edges in ascending order. For n_in bins, there should be n_in+1 edges.
bin_edges_out (ArrayLike) – 1D array of output bin edges in ascending order. For n_out bins, there should be n_out+1 edges.

Returns:

overlap_matrix – Dense matrix of shape (n_in, n_out) where n_in is the number of input bins and n_out is the number of output bins. Each element (i, j) represents the fraction of input bin i that overlaps with output bin j. Values range from 0 (no overlap) to 1 (complete overlap).

Return type:

Notes

Both input arrays must be sorted in ascending order
The function leverages the sorted nature of both inputs for efficiency
Uses vectorized operations to handle large bin arrays efficiently
All overlaps sum to 1.0 for each input bin when output bins fully cover input range

Examples

>>> import numpy as np
>>> from gwtransport.utils import partial_isin
>>> bin_edges_in = np.array([0, 10, 20, 30])
>>> bin_edges_out = np.array([5, 15, 25])
>>> partial_isin(
...     bin_edges_in=bin_edges_in, bin_edges_out=bin_edges_out
... )
array([[0.5, 0. ],
       [0.5, 0.5],
       [0. , 0.5]])

gwtransport.utils.time_bin_overlap(*, tedges, bin_tedges)[source]#

Calculate the fraction of each time bin that overlaps with each time range.

This function computes an array where element (i, j) represents the fraction of time bin j that overlaps with time range i. The computation uses vectorized operations to avoid loops.

Parameters:

tedges (ArrayLike) – 1D array of time bin edges in ascending order. For n bins, there should be n+1 edges.
bin_tedges (list of tuple) – List of tuples where each tuple contains (start_time, end_time) defining a time range.

Returns:

overlap_array – Array of shape (len(bin_tedges), n_bins) where n_bins is the number of time bins. Each element (i, j) represents the fraction of time bin j that overlaps with time range i. Values range from 0 (no overlap) to 1 (complete overlap).

Return type:

Notes

tedges must be sorted in ascending order
Uses vectorized operations to handle large arrays efficiently
Time ranges in bin_tedges can be in any order and can overlap

Examples

>>> import numpy as np
>>> from gwtransport.utils import time_bin_overlap
>>> tedges = np.array([0, 10, 20, 30])
>>> bin_tedges = [(5, 15), (25, 35)]
>>> time_bin_overlap(
...     tedges=tedges, bin_tedges=bin_tedges
... )
array([[0.5, 0.5, 0. ],
       [0. , 0. , 0.5]])

gwtransport.utils.generate_failed_coverage_badge()[source]#

Generate a badge indicating failed coverage.

Return type:: None

gwtransport.utils.combine_bin_series(*, a, a_edges, b, b_edges, extrapolation=0.0)[source]#

Combine two binned series onto a common set of unique edges.

This function takes two binned series (a and b) with their respective bin edges and creates new series (c and d) that are defined on a combined set of unique edges from both input edge arrays.

Parameters:

a (ArrayLike) – Values for the first binned series.
a_edges (ArrayLike) – Bin edges for the first series. Must have len(a) + 1 elements.
b (ArrayLike) – Values for the second binned series.
b_edges (ArrayLike) – Bin edges for the second series. Must have len(b) + 1 elements.
extrapolation (str or float, optional) – Method for handling combined bins that fall outside the original series ranges. - ‘nearest’: Use the nearest original bin value - float value (e.g., np.nan, 0.0): Fill with the specified value (default: 0.0)

Return type:

tuple[ndarray, ndarray, ndarray, ndarray]

Returns:

c (numpy.ndarray) – Values from series a mapped to the combined edge structure.
c_edges (numpy.ndarray) – Combined unique edges from a_edges and b_edges.
d (numpy.ndarray) – Values from series b mapped to the combined edge structure.
d_edges (numpy.ndarray) – Combined unique edges from a_edges and b_edges (same as c_edges).

Notes

The combined edges are created by taking the union of all unique values from both a_edges and b_edges, sorted in ascending order. The values are then broadcasted/repeated for each combined bin that falls within the original bin’s range.

gwtransport.utils.compute_time_edges(*, tedges, tstart, tend, number_of_bins)[source]#

Compute time edges for binning data based on provided time parameters.

This function creates a DatetimeIndex of time bin edges from one of three possible input formats: explicit edges, start times, or end times. The resulting edges define the boundaries of time intervals for data binning.

Define either explicit time edges, or start and end times for each bin and leave the others at None.

Parameters:

tedges (pandas.DatetimeIndex or None) – Explicit time edges for the bins. If provided, must have one more element than the number of bins (n_bins + 1). Takes precedence over tstart and tend.
tstart (pandas.DatetimeIndex or None) – Start times for each bin. Must have the same number of elements as the number of bins. Used when tedges is None.
tend (pandas.DatetimeIndex or None) – End times for each bin. Must have the same number of elements as the number of bins. Used when both tedges and tstart are None.
number_of_bins (int) – The expected number of time bins. Used for validation against the provided time parameters.

Returns:

Time edges defining the boundaries of the time bins. Has one more element than number_of_bins.

Return type:

pandas.DatetimeIndex

Raises:

ValueError – If tedges has incorrect length (not number_of_bins + 1). If tstart has incorrect length (not equal to number_of_bins). If tend has incorrect length (not equal to number_of_bins). If none of tedges, tstart, or tend are provided.

Notes

When using tstart, the function assumes uniform spacing and extrapolates the final edge based on the spacing between the last two start times.
When using tend, the function assumes uniform spacing and extrapolates the first edge based on the spacing between the first two end times.
All input time data is converted to pandas.DatetimeIndex for consistency.

gwtransport.utils.get_soil_temperature(*, station_number=260, interpolate_missing_values=True)[source]#

Download soil temperature data from the KNMI and return it as a pandas DataFrame.

The data is available for the following KNMI weather stations: - 260: De Bilt, the Netherlands (vanaf 1981) - 273: Marknesse, the Netherlands (vanaf 1989) - 286: Nieuw Beerta, the Netherlands (vanaf 1990) - 323: Wilhelminadorp, the Netherlands (vanaf 1989)

TB1 = grondtemperatuur op 5 cm diepte (graden Celsius) tijdens de waarneming TB2 = grondtemperatuur op 10 cm diepte (graden Celsius) tijdens de waarneming TB3 = grondtemperatuur op 20 cm diepte (graden Celsius) tijdens de waarneming TB4 = grondtemperatuur op 50 cm diepte (graden Celsius) tijdens de waarneming TB5 = grondtemperatuur op 100 cm diepte (graden Celsius) tijdens de waarneming TNB2 = minimum grondtemperatuur op 10 cm diepte in de afgelopen 6 uur (graden Celsius) TNB1 = minimum grondtemperatuur op 5 cm diepte in de afgelopen 6 uur (graden Celsius) TXB1 = maximum grondtemperatuur op 5 cm diepte in de afgelopen 6 uur (graden Celsius) TXB2 = maximum grondtemperatuur op 10 cm diepte in de afgelopen 6 uur (graden Celsius)

Parameters:

station_number (int, {260, 273, 286, 323}) – The KNMI station number for which to download soil temperature data. Default is 260 (De Bilt).
interpolate_missing_values (bool, optional) – If True, missing values are interpolated and recent NaN values are extrapolated with the previous value. If False, missing values remain as NaN. Default is True.

Returns:

DataFrame containing soil temperature data in Celsius with a DatetimeIndex. Columns include TB1, TB2, TB3, TB4, TB5, TNB1, TNB2, TXB1, TXB2.

Return type:

pandas.DataFrame

Notes

KNMI: Royal Netherlands Meteorological Institute
The timeseries may contain NaN values for missing data.

gwtransport.utils.solve_underdetermined_system(*, coefficient_matrix, rhs_vector, nullspace_objective='squared_differences', optimization_method='BFGS')[source]#

Solve an underdetermined linear system with nullspace regularization.

For an underdetermined system Ax = b where A has more columns than rows, multiple solutions exist. This function computes a least-squares solution and then selects a specific solution from the nullspace based on a regularization objective.

Parameters:

coefficient_matrix (ArrayLike) – Coefficient matrix of shape (m, n) where m < n (underdetermined). May contain NaN values in some rows, which will be excluded from the system.
rhs_vector (ArrayLike) – Right-hand side vector of length m. May contain NaN values corresponding to NaN rows in coefficient_matrix, which will be excluded from the system.
nullspace_objective (str or Callable, optional) –
Objective function to minimize in the nullspace. Options:
- ”squared_differences” : Minimize sum of squared differences between adjacent elements: sum((x[i+1] - x[i])**2)
- ”summed_differences” : Minimize sum of absolute differences between adjacent elements: sum(|x[i+1] - x[i]|)
- callable : Custom objective function with signature objective(coeffs, x_ls, nullspace_basis) where:
  - coeffs : optimization variables (nullspace coefficients)
  - x_ls : least-squares solution
  - nullspace_basis : nullspace basis matrix
Default is “squared_differences”.
optimization_method (str, optional) – Optimization method passed to scipy.optimize.minimize. Default is “BFGS”.

Returns:

Solution vector that minimizes the specified nullspace objective. Has length n (number of columns in coefficient_matrix).

Return type: