The pfas API Reference

This chapter documents the Python API for the pfas package.

The “model” module

Model orchestration module for PFAS transport simulations.

This module provides the Model class for orchestrating the sequential execution of preprocessing and solving steps with a fluent builder pattern.

class pfas.model.Model(config)

Bases: object

Orchestrate sequential execution of preprocessors and solvers.

The Model class implements a builder pattern for constructing and executing a sequence of preprocessing and solving components. It manages the flow of data from configuration through various transformers and ultimately to the analytical solver.

Parameters:: config (object) – Configuration object containing parameters for the simulation. Parameters are extracted from this object based on the annotated fields of each model class.

config

The configuration object passed at initialization.

Type:: object

generated_data

Dictionary storing outputs from each computation step, making them available as inputs for subsequent steps.

Type:: dict

Examples

>>> from pfas.preprocessing import WaterPreprocessor, BoundaryPreprocessor
>>> from pfas.configuration import SimulationConfig
>>> config = SimulationConfig(...)
>>> result = (Model(config)
...     .add(WaterPreprocessor)
...     .add(BoundaryPreprocessor)
...     .add(GridGenerator)
...     .add(SimulationRunner)
... )

add(model_class, **kwargs)

Add and execute a preprocessing or solving component.

Dynamically instantiates a model class with parameters extracted from the configuration object and previously generated data. Executes the model’s compute() method and stores results for downstream use.

Parameters:

model_class (type) – A preprocessor or solver class with type-annotated __init__ parameters. The class must have a compute() method.
**kwargs (dict, optional) – Additional keyword arguments to override config or generated data values.

Returns:

self – Returns self for method chaining (builder pattern).

Return type:

Model

The “preprocessing” module

PFAS Preproccesing Module.

This module contains preprocessor classes for preparing input parameters for the PFAS analytical solver.

The preprocessors handle calculations for: - Water flow properties - Boundary conditions - Spatial and temporal grids - Adsorption parameters (solid-phase and air-water interface) - Air-water interface adsorption - Simulation execution

Classes

WaterPreprocessor: Computes hydraulic properties from infiltration and soil parameters.
BoundaryPreprocessor: Calculates boundary conditions for contaminant input.
GridGenerator: Generates spatial and temporal discretization grids.
SWCAdsorptionPreprocessor: Calculates air-water interface area using thermodynamic relations.
AdsorptionCollector: Consolidates solid-phase and air-water interface sorption parameters into a single Adsorption object.
SimulationRunner: Executes the analytical solution with preprocessed parameters.

class pfas.preprocessing.WaterPreprocessor(**data)

Bases: BaseModel

Compute hydraulic properties from infiltration rate and soil parameters.

This class calculates volumetric water content, pore water velocity, and dispersion coefficient based on van Genuchten soil hydraulic parameters and steady-state infiltration conditions.

Parameters:

average_infiltration_rate (float) – Average water infiltration rate (m/s). Must be positive.
hydraulic_conductivity (float) – Saturated hydraulic conductivity (m/s). Must be positive.
porosity (float) – Soil porosity (dimensionless). Range: [0, 1].
dispersivity (float) – Longitudinal dispersivity (m). Must be positive.
van_genuchten_n (float) – van Genuchten n parameter (dimensionless). Must be positive.
init_sat (float) – Initial saturation estimate (dimensionless). Range: [0, 1].
residual_water_content (float) – Residual water content (dimensionless). Range: [0, 1].

outputs

List containing ‘hydro_properties’.

Type:: list of str

Examples

>>> preprocessor = WaterPreprocessor(
...     average_infiltration_rate=1e-8,
...     hydraulic_conductivity=1e-5,
...     porosity=0.4,
...     dispersivity=0.1,
...     van_genuchten_n=2.0,
...     init_sat=0.5,
...     residual_water_content=0.05
... )
>>> result = preprocessor.compute()
>>> 'hydro_properties' in result
True

average_infiltration_rate: Annotated[float]

hydraulic_conductivity: Annotated[float]

porosity: Annotated[float]

dispersivity: Annotated[float]

van_genuchten_n: Annotated[float]

init_sat: Annotated[float]

residual_water_content: Annotated[float]

compute()

Calculate hydraulic properties.

Uses van Genuchten relative permeability function to solve for effective saturation, then computes water content, velocity, and dispersion coefficient.

Returns:: Dictionary with key ‘hydro_properties’ containing a HydrologicalProperties instance with theta, v, and D values.
Return type:: dict

property outputs: List of output keys from compute() method.

model_config: ClassVar[ConfigDict] = {'extra': 'forbid', 'validate_assignment': True}: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

class pfas.preprocessing.BoundaryPreprocessor(**data)

Bases: BaseModel

Calculate boundary conditions for contaminant input.

Converts solute concentrations and switching times into the C_list and T_list format required by the analytical solution.

Parameters:

C_list (list of float) – Inlet concentrations for each interval [M L⁻³]. C_list[j] is the concentration active from T_list[j] until T_list[j+1]. The last interval extends to infinity. Must have the same length as T_list. Use 0 for clean water with no PFAS input. Examples: - Continuous step: C_list=[C0], T_list=[0] - Pulse from t=0: C_list=[C0, 0], T_list=[0, t1] - Delayed pulse: C_list=[0, C0, 0], T_list=[0, t1, t2] - Multiple pulses: C_list=[f1, 0, f2], T_list=[0, t1, t2]
T_list (list of float) – Switching times [T] at which the inlet concentration changes. Must have the same length as C_list. T_list[0] must be 0.

outputs

List containing ‘boundary_conditions’.

Type:: list of str

C_list: list[Annotated[float]]

T_list: list[Annotated[float]]

validate_c_and_t_list()

Validate C_list and T_list are consistent.

Return type:: BoundaryPreprocessor

compute()

Pass through C_list and T_list as boundary conditions.

Returns:: Dictionary with key ‘boundary_conditions’ containing a BoundaryConditions instance.
Return type:: dict

property outputs: list[str]: List of output keys from compute() method.

model_config: ClassVar[ConfigDict] = {'extra': 'forbid', 'validate_assignment': True}: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

class pfas.preprocessing.GridGenerator(**data)

Bases: BaseModel

Generate spatial and temporal discretization grids.

Creates uniform grids for depth and time used in the analytical solution.

Parameters:

domain_length (float) – Total depth of domain (m). Must be positive.
spatial_resolution (float) – Grid spacing in depth direction (m). Must be positive.
time_resolution (float) – Time step size (s). Must be positive.
time_total (float) – Total simulation time (s). Must be positive.

outputs

List containing ‘grid’.

Type:: list of str

domain_length: Annotated[float]

spatial_resolution: Annotated[float]

time_resolution: Annotated[float]

time_total: Annotated[float]

compute()

Generate grid arrays.

Returns:: Dictionary with key ‘grid’ containing a SimulationGrid instance with depth and time arrays.
Return type:: dict

property outputs: List of output keys from compute() method.

model_config: ClassVar[ConfigDict] = {'extra': 'forbid', 'validate_assignment': True}: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

class pfas.preprocessing.SWCAdsorptionPreprocessor(**data)

Bases: BaseModel

Calculate air-water interface area using thermodynamic relations.

Uses van Genuchten soil water characteristic curve to estimate air-water interfacial area from water saturation.

Parameters:

hydro_properties (HydrologicalProperties) – Hydraulic properties from WaterPreprocessor.
sigma0 (float, optional) – Surface tension of water (N/m). Default is 0.072.
scaling_factor_awi (float) – Scaling factor for air-water interface area (dimensionless).
AWI (dict) – Dictionary specifying AWI method and associated parameters.
soil (dict) – Dictionary containing soil parameters: ‘porosity’, ‘van_genuchten_alpha’, ‘van_genuchten_n’, and ‘residual_water_content’.

outputs

List containing ‘aaw’.

Type:: list of str

hydro_properties: HydrologicalProperties

sigma0: Annotated[float]

scaling_factor_awi: Annotated[float]

AWI: dict

soil: dict

compute()

Calculate air-water interfacial area.

Returns:: Dictionary with key ‘aaw’ containing the calculated air-water interfacial area (m²/m³).
Return type:: dict

property outputs: List of output keys from compute() method.

model_config: ClassVar[ConfigDict] = {'extra': 'forbid', 'validate_assignment': True}: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

class pfas.preprocessing.SpRetardationPreprocessor(**data)

Bases: BaseModel

Compute solid-phase retardation factor.

Calculates the retardation factor for sorption to the solid phase. Three Kd resolution methods are supported, selected via the sorption_isotherm and Kd_method keys inside sorption_solid:

Linear isotherm — direct input (sorption_isotherm: "linear", Kd_method: "direct_input")

The distribution coefficient is supplied directly:

sorption_solid = {
    "sorption_isotherm": "linear",
    "linear": {"Kd_method": "direct_input", "Kd": 0.042},
    ...
}

Linear isotherm — Fabregat-Palau (2021) (sorption_isotherm: "linear", Kd_method: "fabregat_palau")

Kd is estimated from molecular structure and soil composition using pfas.utils.kd_fabregat_palau():

sorption_solid = {
    "sorption_isotherm": "linear",
    "linear": {
        "Kd_method": "fabregat_palau",
        "n_CFx": 7,
        "f_oc": 0.0004,
        "f_silt_clay": 0.0,
    },
    ...
}

Freundlich isotherm (sorption_isotherm: "freundlich")

An effective linear Kd is derived from the Freundlich isotherm at a representative aqueous concentration C_rep using pfas.utils.kd_freundlich():

sorption_solid = {
    "sorption_isotherm": "freundlich",
    "freundlich": {
        "K_freund": 0.043,
        "n_freund": 0.8,
        "C_rep": 1.0,   # optional, default 1.0 mg/L
    },
    ...
}

Parameters:

sorption_solid (dict) – Dictionary containing sorption parameters as described above.
bulk_density (float) – Soil bulk density (kg/m³). Must be positive.
hydro_properties (HydrologicalProperties) – Hydraulic properties from WaterPreprocessor.

outputs

List containing 'sp_retardation' and 'Kd'.

Type:: list of str

sorption_solid: dict

bulk_density: Annotated[float]

hydro_properties: HydrologicalProperties

compute()

Calculate solid-phase retardation.

Dispatches to the correct Kd resolution method based on sorption_solid["sorption_isotherm"] and, for linear isotherms, sorption_solid["linear"]["Kd_method"].

Returns:: Dictionary with keys 'sp_retardation' and 'Kd'.
Return type:: dict
Raises:: ValueError – If required keys are missing or an unsupported method is specified.

property outputs: List of output keys from compute() method.

model_config: ClassVar[ConfigDict] = {'extra': 'forbid', 'validate_assignment': True}: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

class pfas.preprocessing.SorptionKawiDirectInput(**data)

Bases: BaseModel

Compute adsorption parameters for air-water interface.

Calculates retardation factor for sorption at the air-water interface using directly specified partition coefficient.

Parameters:

kaw (float) – Air-water interface partition coefficient (dimensionless).
hydro_properties (HydrologicalProperties) – Hydraulic properties from WaterPreprocessor.
Kd (float) – Solid-phase partition coefficient (m³/kg).
aaw (float) – Air-water interfacial area (m²/m³).

outputs

List containing ‘awi_retardation’.

Type:: list of str

kaw: float

hydro_properties: HydrologicalProperties

aaw: float

compute()

Calculate air-water interface retardation factor.

Returns:: Dictionary with key ‘awi_retardation’.
Return type:: dict

property outputs: List of output keys from compute() method.

model_config: ClassVar[ConfigDict] = {'extra': 'forbid', 'validate_assignment': True}: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

class pfas.preprocessing.AdsorptionCollector(**data)

Bases: BaseModel

Assemble the final Adsorption object from preprocessor outputs.

Collects the results of SpRetardationPreprocessor and SorptionKawiDirectInput (and any future adsorption preprocessors) and assembles them into the Adsorption instance required by SimulationRunner.

Parameters:

Kd (float) – Solid-phase partition coefficient (m³/kg) from SpRetardationPreprocessor.
sp_retardation (float) – Solid-phase retardation factor from SpRetardationPreprocessor.
awi_retardation (float) – Air-water interface retardation factor from SorptionKawiDirectInput.
sorption_solid (dict) – Dictionary containing solid-phase sorption parameters. Used to extract ‘rate_const’ and ‘fraction_instantaneous’ for the Adsorption object.

outputs

List containing ‘adsorption’.

Type:: list of str

Examples

>>> collector = AdsorptionCollector(
...     Kd=0.001,
...     sp_retardation=3.75,
...     awi_retardation=5.0,
...     sorption_solid={"rate_const": 0.0, "fraction_instantaneous": 1.0},
... )
>>> result = collector.compute()
>>> 'adsorption' in result
True

Kd: float

sp_retardation: float

awi_retardation: float

sorption_solid: dict

compute()

Assemble the Adsorption object.

Returns:: Dictionary with key ‘adsorption’ containing an Adsorption instance.
Return type:: dict

property outputs: List of output keys from compute() method.

model_config: ClassVar[ConfigDict] = {'extra': 'forbid', 'validate_assignment': True}: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

class pfas.preprocessing.SimulationRunner(**data)

Bases: BaseModel

Execute the analytical solution with preprocessed parameters.

This class runs the PFAS-LEACH-Screening analytical solver using all preprocessed input parameters to compute contaminant concentrations in water and solid phases over space and time.

AdsorptionCollector is run internally – the user only needs to supply the outputs of the individual adsorption preprocessors (sp_retardation, Kd, awi_retardation) and sorption_solid directly.

Parameters:

grid (SimulationGrid) – Spatial and temporal discretization from GridGenerator.
bulk_density (float) – Soil bulk density (kg/m3). Must be positive.
boundary_conditions (BoundaryConditions) – Boundary conditions from BoundaryPreprocessor.
hydro_properties (HydrologicalProperties) – Hydraulic properties from WaterPreprocessor.
sp_retardation (float) – Solid-phase retardation factor from SpRetardationPreprocessor.
Kd (float) – Solid-phase partition coefficient (m3/kg) from SpRetardationPreprocessor.
awi_retardation (float) – Air-water interface retardation factor from SorptionKawiDirectInput.
sorption_solid (dict) – Solid-phase sorption parameters, used to extract ‘rate_const’ and ‘fraction_instantaneous’.
kinetic_sorption (bool) – Whether to use kinetic (True) or equilibrium (False) sorption.
volume_averaged (bool) – Whether to compute volume-averaged concentrations.

grid: SimulationGrid

bulk_density: Annotated[float]

boundary_conditions: BoundaryConditions

hydro_properties: HydrologicalProperties

sorption_solid: dict

awi_retardation: float

kinetic_sorption: bool

volume_averaged: bool

initial_contaminant_concentration: Optional[ndarray]

compute()

Run the analytical solution.

Assembles adsorption parameters internally before solving the advection-dispersion equation with retardation for PFAS transport through the vadose zone.

Returns:

Dictionary with keys:

’C1’ndarray
Aqueous phase concentration (mg/L)
’C2’ndarray
Sorbed phase concentration (mg/kg)
’C_tot’ndarray
Total concentration (mg/L bulk volume)

Return type:

dict

model_config: ClassVar[ConfigDict] = {'arbitrary_types_allowed': True, 'extra': 'forbid', 'validate_assignment': True}: Configuration for the model, should be a dictionary conforming to [ConfigDict][pydantic.config.ConfigDict].

The “analytical_soln” module

Analytical solution module for PFAS transport modeling.

This module provides data structures and solvers for simulating contaminant transport through porous media, including equilibrium and kinetic sorption.

class pfas.analytical_soln.SimulationGrid(depth, time)

Bases: object

Grid for spatial and temporal discretization of the simulation domain.

Parameters:

depth (ndarray) – Spatial grid points (z-coordinates) in the domain.
time (ndarray) – Temporal grid points for the simulation.

depth: ndarray[tuple[Any, ...], dtype[float64]]

time: ndarray[tuple[Any, ...], dtype[float64]]

property domain_length: float

Calculate total domain length.

Returns:: Total length of the simulation domain.
Return type:: float

class pfas.analytical_soln.BoundaryConditions(C_list, T_list)

Bases: object

Boundary conditions for contaminant transport.

Parameters:

C_list (list of float) – Inlet concentrations for each interval [M L⁻³]. C_list[j] is the concentration active from T_list[j] until T_list[j+1]. The last interval extends to infinity. Examples: - Continuous step: C_list=[C0], T_list=[0] - Pulse from t=0: C_list=[C0, 0], T_list=[0, t1] - Delayed pulse: C_list=[0, C0, 0], T_list=[0, t1, t2] - Multiple pulses: C_list=[f1, 0, f2], T_list=[0, t1, t2]
T_list (list of float) – Switching times [T] at which the inlet concentration changes. Must have the same length as C_list. T_list[0] must be 0.

C_list: list[float]

T_list: list[float]

class pfas.analytical_soln.HydrologicalProperties(water_content, pore_velocity, dispersion_coefficient)

Bases: object

Hydrological properties of the porous medium.

Parameters:

water_content (float) – Volumetric water content (theta) (dimensionless).
pore_velocity (float) – Average pore water velocity (v) (m/s).
dispersion_coefficient (float) – Hydrodynamic dispersion coefficient (D) (m²/s).

water_content: float

pore_velocity: float

dispersion_coefficient: float

class pfas.analytical_soln.Adsorption(Kd, rate_const, frac_int, sp_retardation, awi_retardation)

Bases: object

Adsorption parameters for contaminant-solid interactions.

Parameters:

Kd (float) – Distribution coefficient for solid-phase partitioning (L/kg).
rate_const (float) – Rate constant for kinetic sorption (alphas) (1/s).
frac_int (float) – Fraction of instantaneous sorption sites (Fs) (dimensionless).
sp_retardation (float) – Solid-phase retardation factor (dimensionless).
awi_retardation (float) – Air-water interface retardation factor (Raw) (dimensionless).

Kd: float

rate_const: float

frac_int: float

sp_retardation: float

awi_retardation: float

property total_retardation: float

Calculate total retardation factor.

Returns:: Sum of solid-phase and air-water interface retardation.
Return type:: float

property beta_s: float

Calculate solid-phase sorption parameter.

Returns:: Dimensionless sorption parameter for solid phase.
Return type:: float

property beta: float

Calculate combined sorption parameter.

Returns:: Dimensionless parameter accounting for both solid-phase and AWI sorption.
Return type:: float

pfas.analytical_soln.analytical_soln(grid, bulk_density, boundary_conditions, initial_contaminant_concentration, hydro_properties, adsorption, kinetic=False, volume_averaged=True)

Solve contaminant transport using analytical solutions.

Computes aqueous and sorbed phase concentrations for PFAS transport through the vadose zone using analytical solutions to the advection-dispersion equation (ADE) with retardation. Dimensionless parameters are computed via compute_dimensionless_params() and passed to the appropriate solver.

Parameters:

grid (SimulationGrid) – Spatial and temporal discretization grid. Must have .depth (m) and .time (s) arrays.
bulk_density (float) – Bulk density of the porous medium (kg/L).
boundary_conditions (BoundaryConditions) – Contaminant source boundary conditions. Must have .C_list [M L⁻³] and .T_list [T] defining the inlet concentration history.
initial_contaminant_concentration (ndarray of shape (n_depth,)) – Initial aqueous concentration distribution in the domain [M L⁻³].
hydro_properties (HydrologicalProperties) – Hydrological properties. Must have .pore_velocity [L T⁻¹], .dispersion_coefficient [L² T⁻¹], and .water_content [-].
adsorption (Adsorption) – Adsorption parameters. Must have .total_retardation, .Kd, .sp_retardation, .frac_int, .beta, and .beta_s. When kinetic=True, also requires .rate_const.
kinetic (bool, optional) – If True, use the kinetic sorption model (kinetic_solver()), which returns a separate sorbed phase C2. If False (default), use the equilibrium model (equilibrium_solver()), and C2 is None.
volume_averaged (bool, optional) – If True, use volume-averaged (resident) concentrations in the kinetic BVP kernel. If False (default), use flux-averaged concentrations. Only used by the kinetic solver.

Return type:

tuple[ndarray[tuple[Any, ...], dtype[float64]], ndarray[tuple[Any, ...], dtype[float64]] | None, ndarray[tuple[Any, ...], dtype[float64]]]

Returns:

C1 (ndarray) – Aqueous phase concentration [M L⁻³].
C2 (ndarray or None) – Sorbed phase concentration [M M⁻¹]. None when kinetic=False.
C_tot (ndarray) – Total concentration [M L⁻³] bulk volume.

Raises:

ValueError – If pore_velocity or dispersion_coefficient is zero.

The “solvers” module

High-level analytical solvers for contaminant transport.

This module provides the public solver interface. All mathematical primitives (BVP/IVP helpers, kinetic kernels, dimensionless parameter computation) live in solver_utils.py. The solvers here are purely orchestration logic.

Main functions

equilibrium_solver() — ADE with instantaneous sorption equilibrium
kinetic_solver() — ADE with first-order kinetic sorption

pfas.solvers.equilibrium_solver(R, dim, C_list, Ci, theta, bc='resident')

Solve advection-dispersion equation with equilibrium sorption.

Computes aqueous and total concentrations for PFAS transport through porous media assuming instantaneous sorption equilibrium, using analytical solutions to the ADE. The solution combines contributions from the boundary value problem (BVP) and, when non-zero initial conditions are present, the initial value problem (IVP):

C(Z,T) = C^B(Z,T) + C^I(Z,T)

The BVP term implements CXTFIT eq. 2.20 directly. Given a series of rectangular pulses with concentrations C_list = [f1, f2, ..., fn] switching at times dim.T_list = [T1, T2, ..., Tn] (dimensionless), the concentration increments are:

deltaC = np.diff([0] + C_list) → [f1-f0, f2-f1, …, fn-f_{n-1}]

and eq. 2.20 becomes:

C^B(Z,T) = sum_{j=1}^{i} deltaC[j] * G1^E(Z, T - T_j; mu^E)

where each term is only added when T > T_j (Heaviside). This naturally handles: - Step input: C_list=[C0], dim.T_list=[0] - Single pulse: C_list=[C0, 0], dim.T_list=[0, T_pulse_end] - Multiple pulses: C_list=[f1,f2,…], dim.T_list=[T1, T2, …]

The IVP term integrates the Green’s function kernel over the initial concentration profile Ci(xi) (CXTFIT Table 2.2).

Parameters:

R (float) – Retardation factor, R = 1 + rho_b * Kd / theta.
dim (DimensionlessParams) – Dimensionless parameters from compute_dimensionless_params(). Uses .Z, .T, .P, and .T_list (dimensionless switching times corresponding to each entry in C_list).
C_list (list of float) – Inlet concentrations for each interval (mg/L). C_list[j] is the concentration active from dim.T_list[j] until dim.T_list[j+1]. The last interval extends to infinity. Must have the same length as dim.T_list.
Ci (ndarray of shape (n_depth,)) – Normalised initial concentration profile Ci(Z) (mg/L). Pass an array of zeros if there is no initial contamination.
theta (float) – Volumetric water content (-).
bc (str, optional) – Upper boundary condition type. Must be a key in _BVP_FUNCTIONS. Options: 'flux' (third-type BC) or 'resident' (first-type BC). Default is 'resident'.

Return type:

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

Returns:

C1 (ndarray of shape (len(Z), len(T))) – Aqueous phase concentration (mg/L).
C_tot (ndarray of shape (len(Z), len(T))) – Total concentration (mg/L bulk volume).

Raises:

ValueError – If bc is not a recognised boundary condition type.
ValueError – If len(C_list) != len(dim.T_list).

References

Toride, Leij & van Genuchten (1995), CXTFIT Version 2.0, Research Report No. 137, USDA-ARS. Section 2, eq. (2.20) for multiple rectangular pulses via step superposition.

pfas.solvers.kinetic_solver(R, dim, C_list, Ci, beta_s, beta, volume_averaged, R_s, f, Kd, theta, rho_b)

Solve advection-dispersion equation with kinetic (time-dependent) sorption.

When beta_s == 1 (no kinetic sites), delegates to equilibrium_solver. Otherwise computes aqueous (C₁) and sorbed phase (C₂) concentrations for contaminant transport with first-order kinetic sorption. The total solution combines boundary value problem (BVP) and initial value problem (IVP) contributions:

C(Z,T) = C^B(Z,T) + C^I(Z,T)

BVP term — CXTFIT eq. 3.20

Pulse superposition using eq. 2.20 over switching times dim.T_list with concentration increments:

deltaC = np.diff([0] + C_list) → [f1-f0, f2-f1, …, fn-f_{n-1}]

C^B_k(Z,T) = Σⱼ deltaC[j] · Aₖ(Z, T - T_list[j]) for T > T_list[j]

where A₁ (aqueous, C1_bvp) and A₂ (sorbed, C2_bvp) are evaluated via _bvp_neq() (CXTFIT eqs. 3.21–3.22). The sorbed phase BVP is additionally scaled by (1-f) * Kd after the loop to convert from dimensionless to mg/kg units.

IVP term — CXTFIT eqs. 3.31, 3.32 / Table 3.4

When beta_s == 1 (no kinetic sites), the IVP reduces to the equilibrium Green’s function at current time T:

C1_ivp = G(Z, T)

When beta_s != 1, the IVP splits into three contributions:

Initial aqueous concentration contribution, modified by inter-phase mass transfer (CXTFIT eq. 3.23, first term):

C1_ivp = exp( -ω·T·(1-f)·Rₛ / ((1-βₛ)·β·R·(1+Rₛ)) ) · G(Z, T)
Initial sorbed concentration contribution (CXTFIT eq. 3.24, first term):

C2_ivp = (1-f)·Kd·Cᵢ · exp( -ω·T / ((1-βₛ)·(1+Rₛ)) )
Convolution integrals over intermediate times τ ∈ (0, T), using the H₀ and Hₛ kernels from CXTFIT Table 3.4 (eqs. 3.31-3.32, second terms):

C1_ivp += ω/((1-βₛ)·(1+Rₛ)) · ∫₀ᵀ H₀(T,τ) · G(Z,τ) dτ C2_ivp += ω/((1-βₛ)·(1+Rₛ)) · (1-f)·Kd · ∫₀ᵀ Hₛ(T,τ) · G(Z,τ) dτ

Total concentration — CXTFIT eq. 3.6

C_tot = θ·β·R·C₁ + ρ_b·C₂

Parameters:

R (float) – Overall retardation factor, R = 1 + ρ_b·Kd/θ (CXTFIT Table 3.1).
dim (DimensionlessParams) – Dimensionless parameters from compute_dimensionless_params(). Uses .Z, .T, .P, .T_list (dimensionless switching times), and .omega (ω).
C_list (list of float) – Inlet concentrations for each interval (mg/L). C_list[j] is the concentration active from dim.T_list[j] until dim.T_list[j+1]. The last interval extends to infinity. Must have the same length as dim.T_list.
Ci (ndarray of shape (n_depth,)) – Normalised initial aqueous concentration profile with depth (mg/L). Pass an array of zeros if there is no initial contamination.
beta_s (float) – Solid-phase partitioning coefficient β_s (CXTFIT Table 3.1).
beta (float) – Dimensionless partitioning coefficient β (CXTFIT Table 3.1).
volume_averaged (bool) – If True, use volume-averaged (resident) concentrations in the BVP kernel _FT(). If False, use flux-averaged concentrations.
R_s (float) – Retardation factor for kinetic sorption sites (CXTFIT Table 3.1).
f (float) – Fraction of sorption sites at instantaneous equilibrium (CXTFIT Table 3.1).
Kd (float) – Linear distribution coefficient (L/kg).
theta (float) – Volumetric water content (-).
rho_b (float) – Bulk density of the porous medium ρ_b (kg/L).

Return type:

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

Returns:

C1 (ndarray of shape (len(Z), len(T))) – Aqueous phase concentration (mg/L).
C2 (ndarray of shape (len(Z), len(T))) – Sorbed phase concentration (mg/kg).
C_tot (ndarray of shape (len(Z), len(T))) – Total concentration (mg/L bulk volume), C_tot = θ·β·R·C₁ + ρ_b·C₂ (CXTFIT eq. 3.6).

References

van Genuchten, M. Th. (1981). Non-Equilibrium Transport Parameters from Miscible Displacement Experiments. Research Report No. 119, USDA-ARS.

Toride, Leij & van Genuchten (1995). CXTFIT Version 2.0. Research Report No. 137, USDA-ARS. Eqs. 3.6, 3.20–3.24; Tables 3.1, 3.4.

Lindstrom, F.T. and Stone, W.J. (1974). Soil Sci. Soc. Am. Proc.

The “utils” module

Utility functions for PFAS analytical solver.

This module provides helper functions for kinetic sorption calculations, air-water interface area estimation, and numerical integration support.

pfas.utils.aaw_func_thermo(sigma0, poro, alpha, n, th, thr, ths, sf)

Compute air-water interfacial area using thermodynamic approach.

Estimates the air-water interfacial area per unit volume of porous medium using thermodynamic relations based on capillary pressure and water content, following van Genuchten soil water retention characteristics.

Parameters:

sigma0 (float) – Surface tension of water (dyn/cm).
poro (float) – Porosity of the porous medium (dimensionless, 0-1).
alpha (float) – van Genuchten parameter (cm⁻¹), related to the air-entry pressure.
n (float) – van Genuchten parameter (dimensionless), related to pore size distribution.
th (float) – Current water content (dimensionless).
thr (float) – Residual water content (dimensionless).
ths (float) – Saturated water content (dimensionless).
sf (float) – Scaling factor to correct the thermodynamic-based estimate (dimensionless).

Returns:

Aaw – Air-water interfacial area per unit volume (cm²/cm³).

Return type:

float

Notes

The function integrates the capillary pressure curve over saturation range to estimate the interfacial area. Water properties are assumed:

Water density: 1000 kg/m³
Gravitational acceleration: 9.81 m/s²

pfas.utils.aaw_func_tracer(sw, x2, x1, x0)

Compute air-water interfacial area using empirical polynomial model.

Estimates air-water interfacial area per unit volume using polynomial fitting coefficients derived from tracer experiments or pore-scale imaging. This approach provides a simplified, computationally efficient alternative to thermodynamic calculations.

Parameters:

sw (float or ndarray) – Water saturation (dimensionless, 0-1).
x2 (float) – Quadratic coefficient of polynomial fit.
x1 (float) – Linear coefficient of polynomial fit.
x0 (float) – Constant coefficient (intercept) of polynomial fit.

Returns:

Aaw – Air-water interfacial area per unit volume (cm²/cm³).

Return type:

float or ndarray

Notes

The polynomial model is: Aaw = x2*sw² + x1*sw + x0

This function is useful when empirical coefficients have been determined from experimental data for a specific porous medium.

pfas.utils.kd_fabregat_palau(n_CFx, f_oc, f_silt_clay)

Calculate distribution coefficient using Fabregat-Palau (2021) model.

Computes the soil-water distribution coefficient (Kd) for PFAS compounds based on the number of perfluorinated carbons and soil organic carbon and silt-clay content.

Parameters:

n_CFx (int) – Number of perfluorinated carbons (CF2 groups) in the PFAS molecule.
f_oc (float) – Fraction of organic carbon in soil (dimensionless, 0-1).
f_silt_clay (float) – Fraction of silt and clay in soil (dimensionless, 0-1).

Returns:

Kd – Distribution coefficient (L/kg).

Return type:

float

References

Fabregat-Palau et al. (2021). Modelling the sorption behaviour of perfluoroalkyl acids in soils.

pfas.utils.k_sc_fabregat_palau2021(n_CFx)

Calculate silt-clay sorption coefficient (Fabregat-Palau 2021).

Parameters:: n_CFx (int) – Number of perfluorinated carbons (CF2 groups).
Returns:: k_sc – Silt-clay sorption coefficient (L/kg).
Return type:: float

References

Fabregat-Palau et al. (2021). Modelling the sorption behaviour of perfluoroalkyl acids in soils.

pfas.utils.k_oc_fabregat_palau2021(n_CFx)

Calculate organic carbon sorption coefficient (Fabregat-Palau 2021).

Parameters:: n_CFx (int) – Number of perfluorinated carbons (CF2 groups).
Returns:: k_oc – Organic carbon sorption coefficient (L/kg).
Return type:: float

References

Fabregat-Palau et al. (2021). Modelling the sorption behaviour of perfluoroalkyl acids in soils.

pfas.utils.kd_freundlich(C_rep, K_freund, n_freund)

Calculate distribution coefficient using the Freundlich sorption model.

Computes the soil-water distribution coefficient (Kd) from a Freundlich isotherm, which describes non-linear sorption onto soil.

The Freundlich isotherm is defined as:: S = K_freund * C^n_freund
which gives a concentration-dependent Kd:: Kd = K_freund * C_rep^(n_freund - 1)

For C_rep = 0, Kd reduces to K_freund (equivalent to C_rep = 1).

Parameters:

C_rep (float) – Representative aqueous-phase concentration [mg/L]. If zero, Kd is returned as K_freund (C_rep = 1 assumed).
K_freund (float) – Freundlich capacity coefficient [(mg/kg) / (mg/L)^n_freund].
n_freund (float) – Freundlich exponent [-]. n_freund < 1 indicates favourable (concave) sorption; n_freund = 1 recovers linear (Kd) sorption.

Returns:

Kd – Concentration-dependent distribution coefficient [L/kg].

Return type:

float

The “data_loader” module

pfas.data_loader.

Utilities for loading and listing the packaged JSON datasets bundled with the pfas library, plus support for loading user‑supplied JSON files.

Functions

load_dataset(name_or_path): Load a packaged dataset by name, or load a JSON file from a filesystem path.
available_datasets(): Return a list of packaged dataset names.
load_json_file(path): Load a JSON file from an arbitrary filesystem path.

Example

>>> from pfas.data_loader import available_datasets, load_dataset
>>> print(available_datasets())
['PFASs', 'soils', 'sp_matrix']

# Load packaged dataset >>> data = load_dataset(‘PFASs’)

# Load external JSON file >>> data = load_dataset(‘/path/to/custom.json’)

pfas.data_loader.load_json_file(path): Load a JSON file from an arbitrary filesystem path.

pfas.data_loader.load_dataset(name_or_path)

Load a packaged dataset by name, or load a JSON file from a filesystem path.

If name_or_path is a valid file path, the JSON file is loaded directly. Otherwise, it must match one of the packaged dataset names.

pfas.data_loader.available_datasets(): Return a list of all packaged dataset names.

The “configuration” module

TOML configuration file validator for PFAS transport modeling.

This module provides functionality to read and validate TOML configuration files for soil transport simulations, including validation of experimental conditions, soil properties, sorption parameters, air-water interface settings, and PFAS properties.

class pfas.configuration.Configuration(config_dict)

Bases: object

A read-only view over a parsed TOML configuration.

This class wraps a nested dictionary and provides convenient attribute-style access to configuration values. If an attribute is not present on the instance, attribute lookup falls back to a recursive search of the underlying dictionary using find().

config_dict

The underlying nested dictionary containing the parsed configuration.

Type:: Dict[str, Any]

find(key, sub_dict=None): Recursively search for key inside the configuration dictionary and return the associated value if found, otherwise None.

Notes

This class intentionally implements __getattribute__ to allow attribute-style access for keys present in the configuration. Use read_toml() to construct a Configuration from a TOML file.

Example

>>> config = read_toml(Path("config.toml"))
>>> config.domain_length

find(key, sub_dict=None)

pfas.configuration.read_toml(path)

Read a TOML configuration file.

Parameters:: path (Path) – Path to the TOML file
Return type:: Configuration
Returns:: Configuration object wrapping the parsed TOML configuration

pfas.configuration.validate_config(config_dict)

Validate all sections of the configuration dictionary.

Parameters:: config_dict (Dict[str, Any]) – Configuration dictionary parsed from TOML
Return type:: bool
Returns:: True if all validations pass, False otherwise

pfas.configuration.check_toml_experimental_conditions(experimental_conditions_dict)

Validate experimental conditions parameters from TOML config.

Parameters:: experimental_conditions_dict (Dict[str, Any]) – Dictionary containing experimental conditions
Return type:: bool
Returns:: True if validation passes, False otherwise

pfas.configuration.check_toml_soil(soil_dict)

Validate soil parameters from TOML config.

Parameters:: soil_dict (Dict[str, Any]) – Dictionary containing soil parameters
Return type:: bool
Returns:: True if validation passes, False otherwise

pfas.configuration.check_toml_sorption(sorption_dict)

Validate sorption parameters from TOML config.

Parameters:: sorption_dict (Dict[str, Any]) – Dictionary containing sorption parameters
Return type:: bool
Returns:: True if validation passes, False otherwise

pfas.configuration.check_toml_awi(awi_dict)

Validate air-water interface parameters from TOML config.

Parameters:: awi_dict (Dict[str, Any]) – Dictionary containing AWI parameters
Return type:: bool
Returns:: True if validation passes, False otherwise

pfas.configuration.check_toml_sorption_awi(sorption_awi_dict)

Validate air-water interface sorption parameters from TOML config.

Parameters:: sorption_awi_dict (Dict[str, Any]) – Dictionary containing AWI sorption parameters
Return type:: bool
Returns:: True if validation passes, False otherwise

pfas.configuration.check_toml_pfas(pfas_dict)

Validate PFAS parameters from TOML config.

Parameters:: pfas_dict (Dict[str, Any]) – Dictionary containing PFAS parameters
Return type:: bool
Returns:: True if validation passes, False otherwise

The “solver_utils” module

Mathematical primitives and preprocessing utilities for ADE analytical solvers.

This module contains: - Dimensionless parameter computation and the DimensionlessParams container - BVP helper functions for equilibrium sorption (one per boundary condition type) - IVP helper functions for equilibrium and kinetic sorption - Kinetic sorption BVP functions (Goldstein J-function + Bessel approximation) - Kinetic sorption IVP convolution kernels H₀ and Hₛ (Bessel function based, CXTFIT Table 3.4)

The high-level solvers in solvers.py import from here. To add a new boundary condition, implement a helper with the signature f(T, R, Z, P) -> ndarray and register it in _BVP_FUNCTIONS.

References

van Genuchten & Alves (1982): Analytical solutions of the one-dimensional convective-dispersive solute transport equation. Toride, Leij & van Genuchten (1995): The CXTFIT Code, Version 2.0, USDA Research Report No. 137. van Genuchten, M. Th. (1981): Non-Equilibrium Transport Parameters from Miscible Displacement Experiments. Research Report No. 119, USDA-ARS. Lindstrom, F.T. and Stone, W.J. (1974): On the start up or initial phase of linear mass transport of chemicals in a water saturated sorbing porous medium. Soil Science Society of America Proceedings.

class pfas.solver_utils.DimensionlessParams(Z: NDArray[np.float64], T: NDArray[np.float64], T_list: list[float], P: float, omega: float | None)

Bases: NamedTuple

Dimensionless parameters for the ADE analytical solution.

Z

Dimensionless depth (-), Z = z / L.

Type:: ndarray

T

Dimensionless time (-), T = t * v / L.

Type:: ndarray

T_list

Dimensionless switching times corresponding to each entry in C_list. Converted from physical time (s) via T = t * v / L.

Type:: list of float

P

Péclet number (-), P = v * L / D.

Type:: float

omega

Dimensionless mass transfer coefficient (ω, CXTFIT Table 3.1). None when kinetic=False.

Type:: float or None

Z: ndarray[tuple[Any, ...], dtype[float64]]: Alias for field number 0

T: ndarray[tuple[Any, ...], dtype[float64]]: Alias for field number 1

T_list: list[float]: Alias for field number 2

P: float: Alias for field number 3

omega: float | None: Alias for field number 4

pfas.solver_utils.compute_dimensionless_params(grid, hydro_properties, T_list, adsorption=None, kinetic=False)

Compute dimensionless parameters for the ADE analytical solution.

Converts physical grid, flow, and transport parameters into the dimensionless form required by the equilibrium and kinetic solvers. Switching times in physical time (seconds) are converted to dimensionless time T = t * v / L.

Parameters:

grid (SimulationGrid) – Spatial and temporal discretization. Must have .depth (m) and .time (s) arrays.
hydro_properties (HydrologicalProperties) – Hydrological properties. Must have .pore_velocity (m/s) and .dispersion_coefficient (m²/s).
T_list (list of float) – Switching times in physical time (s) at which the inlet concentration changes. Converted to dimensionless pore volumes via T = t * v / L. T_list[0] should normally be 0. Examples: - Continuous step: [0] - Pulse from t=0: [0, 5000] - Delayed pulse: [0, 2000, 5000] - Multiple pulses: [0, 1000, 2000, 3000]
adsorption (Adsorption, optional) – Adsorption parameters. Required when kinetic=True. Must have .rate_const, .beta_s, .sp_retardation.
kinetic (bool, optional) – If True, also compute the Damköhler number omega. Default is False.

Returns:

Named tuple containing Z, T, T_list (dimensionless), P, and omega.

Return type:

DimensionlessParams

Raises:

ValueError – If pore_velocity or dispersion_coefficient is zero.
ValueError – If kinetic=True but adsorption is None.