Coverage for src/gwtransport/advection.py: 91%
177 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 21:13 +0000
1"""
2Advective Transport Modeling Along Aquifer Pore Volumes.
4Water infiltrates and is transported in parallel along multiple aquifer pore volumes to
5extraction. For each aquifer pore volume, transport is 1D advection with linear or non-linear
6sorption; there is no microdispersion or molecular diffusion, while the spread across aquifer
7pore volumes provides macrodispersion. Forward and backward modeling are supported. No assumption
8is made about whether the flow is radial or orthogonal.
10Note on dispersion: The spreading from the pore volume distribution (APVD) represents
11macrodispersion—aquifer-scale velocity heterogeneity that depends on both aquifer
12properties and hydrological boundary conditions. To add microdispersion and molecular
13diffusion separately (when APVD comes from streamline analysis), use :mod:`gwtransport.diffusion`.
14See :ref:`concept-dispersion-scales` for details.
16Note on cross-compound calibration: When APVD is calibrated from measurements of one
17compound (e.g., temperature with D_m ~ 0.1 m²/day) and used to predict another (e.g., a
18solute with D_m ~ 1e-4 m²/day), the molecular diffusion contribution is baked into the
19calibrated std. The cleanest fix is to calibrate with :mod:`gwtransport.diffusion_fast`
20instead, which keeps the three contributions separate.
22The aquifer pore volume distribution (APVD) enters in one of two forms. The ``gamma_*`` pair takes
23the APVD in closed form as a (shifted) gamma distribution, parameterized by either (mean, std, loc)
24or (alpha, beta, loc), and discretizes it internally into ``n_bins`` equal-probability-mass
25streamtubes. The plain pair takes ``aquifer_pore_volumes``: an explicit array of pore volumes [m³],
26one per streamtube, treated as equally weighted at the outlet.
28Available functions:
30- :func:`infiltration_to_extraction` - Compute the concentration (or temperature) of the extracted
31 water from an infiltration series and an explicit array of aquifer pore volumes. Each
32 ``cout_tedges`` bin gets the flow-weighted average of ``cin`` over each streamtube's
33 back-projected source window, averaged over the streamtubes; units are those of ``cin``. Sorption
34 is linear, through a constant ``retardation_factor``. The ``spinup`` policy (default
35 ``"constant"``) warm-starts the record before ``tedges[0]``; ``spinup=None`` instead returns NaN
36 for every cout bin whose streamtube source windows are not fully inside the cin range.
38- :func:`extraction_to_infiltration` - Reverse operation: reconstruct the infiltrating
39 concentration from measured extraction concentrations by inverting the forward weight matrix
40 with Tikhonov regularization (``regularization_strength``). Returns one value per ``tedges``
41 bin. NaN in ``cout`` marks measurement gaps; those rows are excluded from the solve, and cin
42 bins constrained only by gapped cout bins come back as NaN.
44- :func:`gamma_infiltration_to_extraction` - Forward transport as in
45 :func:`infiltration_to_extraction`, with the APVD supplied as gamma parameters rather than as an
46 explicit array of pore volumes.
48- :func:`gamma_extraction_to_infiltration` - Deconvolution as in
49 :func:`extraction_to_infiltration`, with the APVD supplied as gamma parameters rather than as an
50 explicit array of pore volumes.
52- :func:`infiltration_to_extraction_nonlinear_sorption` - Forward transport with a
53 concentration-dependent isotherm -- Freundlich, Langmuir, or a constant retardation factor --
54 solved by front tracking along each pore volume. Returns a tuple: the flow-weighted bin-averaged
55 extraction concentration, and one diagnostic structure per pore volume holding the waves, events,
56 first arrival and final solver state. Cout bins whose source window leaves the flow record are
57 returned as ``0.0`` rather than NaN, while bins with zero throughflow are NaN. Forward only:
58 non-linear sorption forms shocks, and there is no extraction-to-infiltration counterpart.
60This file is part of gwtransport which is released under AGPL-3.0 license.
61See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
62"""
64from collections import Counter
66import numpy as np
67import numpy.typing as npt
68import pandas as pd
70from gwtransport import gamma
71from gwtransport._time import tedges_to_days
72from gwtransport._validation import (
73 _validate_no_nan,
74 _validate_non_negative_array,
75 _validate_positive_array,
76 _validate_retardation_factor,
77 _validate_tedges_parity,
78)
79from gwtransport.advection_utils import (
80 _infiltration_to_extraction_weights,
81 _resolve_spinup_inputs,
82 _resolve_spinup_mask,
83)
84from gwtransport.fronttracking.math import (
85 EPSILON_FREUNDLICH_N,
86 ConstantRetardation,
87 FreundlichSorption,
88 LangmuirSorption,
89 SorptionModel,
90)
91from gwtransport.fronttracking.output import compute_bin_averaged_concentration_exact
92from gwtransport.fronttracking.solver import FrontTracker, find_unresolved_interaction
93from gwtransport.fronttracking.waves import CharacteristicWave, RarefactionWave, ShockWave
94from gwtransport.utils import solve_inverse_transport_banded
97def _validate_advection_inputs(
98 *,
99 tedges: pd.DatetimeIndex,
100 flow: np.ndarray,
101 retardation_factor: float,
102 aquifer_pore_volumes: npt.ArrayLike,
103 cin_values: np.ndarray | None = None,
104 cout_values: np.ndarray | None = None,
105 cout_tedges: pd.DatetimeIndex | None = None,
106) -> None:
107 """Validate inputs common to advection forward / reverse entry points.
109 Path selection via mutually-exclusive kwargs:
111 - ``cin_values`` provided => forward path. ``tedges`` parities cin and flow.
112 - ``cout_values`` + ``cout_tedges`` provided => reverse path. ``tedges`` parities
113 flow; ``cout_tedges`` parities cout.
115 All shared checks fire on both paths. ``flow >= 0`` is enforced in both
116 directions, and ``aquifer_pore_volumes`` must be non-empty, finite and
117 strictly positive -- a negative or zero volume would source a cout bin from
118 future infiltration, and an empty distribution carries no water at all.
120 Raises
121 ------
122 ValueError
123 If any check fails. The message identifies which invariant was violated.
124 """
125 if cin_values is not None:
126 _validate_tedges_parity(tedges, cin_values, tedges_name="tedges", values_name="cin")
127 _validate_tedges_parity(tedges, flow, tedges_name="tedges", values_name="flow")
128 _validate_no_nan(cin_values, name="cin")
129 elif cout_values is not None and cout_tedges is not None:
130 _validate_tedges_parity(tedges, flow, tedges_name="tedges", values_name="flow")
131 # Reverse cout may contain NaN (measurement gaps); gapped rows are excluded from the solve.
132 _validate_tedges_parity(cout_tedges, cout_values, tedges_name="cout_tedges", values_name="cout")
133 else:
134 msg = "must provide cin_values (forward) or both cout_values and cout_tedges (reverse)"
135 raise ValueError(msg)
136 # Non-monotonic tedges would silently corrupt the cumulative-volume mapping.
137 if np.any(np.diff(tedges.asi8) <= 0):
138 msg = "tedges must be strictly increasing"
139 raise ValueError(msg)
140 _validate_no_nan(flow, name="flow")
141 _validate_non_negative_array(flow, name="flow", message="flow must be non-negative (negative flow not supported)")
142 _validate_retardation_factor(retardation_factor)
143 apv = np.asarray(aquifer_pore_volumes, dtype=float)
144 if apv.size == 0:
145 msg = "aquifer_pore_volumes must not be empty"
146 raise ValueError(msg)
147 # A negative or zero pore volume back-projects a cout bin to *future*
148 # infiltration (anti-causal); a non-finite one poisons the whole solve.
149 _validate_positive_array(apv, name="aquifer_pore_volumes")
152def gamma_infiltration_to_extraction(
153 *,
154 cin: npt.ArrayLike,
155 flow: npt.ArrayLike,
156 tedges: pd.DatetimeIndex,
157 cout_tedges: pd.DatetimeIndex,
158 mean: float | None = None,
159 std: float | None = None,
160 loc: float = 0.0,
161 alpha: float | None = None,
162 beta: float | None = None,
163 n_bins: int = 100,
164 retardation_factor: float = 1.0,
165 spinup: str | None = "constant",
166) -> npt.NDArray[np.floating]:
167 """
168 Compute the concentration of the extracted water by shifting cin with its residence time.
170 The compound is retarded in the aquifer with a retardation factor. The residence
171 time is computed based on the flow rate of the water in the aquifer and the pore volume
172 of the aquifer. The aquifer pore volume is approximated by a (shifted) gamma distribution
173 parameterized by either (mean, std, loc) or (alpha, beta, loc).
175 This function represents infiltration to extraction modeling by flow-weighted averaging.
177 Provide either (mean, std) or (alpha, beta); ``loc`` is optional and defaults to 0.
179 Parameters
180 ----------
181 cin : array-like
182 Concentration of the compound in infiltrating water or temperature of infiltrating
183 water. The model assumes this value is constant over each interval
184 ``[tedges[i], tedges[i+1])``.
185 flow : array-like
186 Flow rate of water in the aquifer [m³/day]. The model assumes this value is
187 constant over each interval ``[tedges[i], tedges[i+1])``.
188 tedges : pandas.DatetimeIndex
189 Time edges for both cin and flow data. Used to compute the cumulative concentration.
190 Has a length of one more than `cin` and `flow`.
191 cout_tedges : pandas.DatetimeIndex
192 Time edges for the output data. Used to compute the cumulative concentration.
193 Has a length of one more than the desired output length.
194 mean : float, optional
195 Mean of the gamma distribution of the aquifer pore volume. Must be strictly
196 greater than ``loc``.
197 std : float, optional
198 Standard deviation of the gamma distribution of the aquifer pore volume
199 (invariant under the ``loc`` shift).
200 loc : float, optional
201 Location (minimum pore volume) of the gamma distribution. Must satisfy
202 ``0 <= loc < mean``. Default is ``0.0``.
203 alpha : float, optional
204 Shape parameter of gamma distribution of the aquifer pore volume (must be > 0).
205 beta : float, optional
206 Scale parameter of gamma distribution of the aquifer pore volume (must be > 0).
207 n_bins : int, optional
208 Number of bins to discretize the gamma distribution. Default 100.
209 retardation_factor : float, optional
210 Retardation factor of the compound in the aquifer (default 1.0).
211 Values > 1.0 indicate slower transport due to sorption/interaction.
212 spinup : {"constant"} | None, optional
213 Forwarded to :func:`infiltration_to_extraction`. Default
214 ``"constant"`` warm-starts the system before ``tedges[0]``.
216 Returns
217 -------
218 numpy.ndarray
219 Concentration of the compound in the extracted water, or temperature. Same units as cin.
221 See Also
222 --------
223 infiltration_to_extraction : Transport with explicit pore volume distribution
224 gamma_extraction_to_infiltration : Reverse operation (deconvolution)
225 gwtransport.gamma.bins : Create gamma distribution bins
226 gwtransport.residence_time.full : Compute residence times
227 gwtransport.diffusion.infiltration_to_extraction : Add microdispersion and molecular diffusion
228 :ref:`concept-gamma-distribution` : Two-parameter pore volume model
229 :ref:`assumption-gamma-distribution` : When gamma distribution is adequate
231 Notes
232 -----
233 The APVD is only time-invariant under the steady-streamlines assumption
234 (see :ref:`assumption-steady-streamlines`).
236 The spreading from the gamma-distributed pore volumes represents macrodispersion
237 (aquifer-scale heterogeneity). When ``std`` comes from calibration on measurements,
238 it absorbs all mixing: macrodispersion, microdispersion, and an average molecular
239 diffusion contribution. When calibrating with the diffusion module, these three
240 components are taken into account separately. When ``std`` comes from streamline
241 analysis, it represents macrodispersion only; microdispersion and molecular diffusion
242 can be added via :mod:`gwtransport.diffusion_fast` or :mod:`gwtransport.diffusion`.
244 For cross-compound prediction (calibrating on temperature and predicting a solute),
245 calibrate with :mod:`gwtransport.diffusion_fast` so the three contributions are
246 tracked separately rather than lumped into a single calibrated ``std``.
247 See :ref:`concept-dispersion-scales` for background.
249 Examples
250 --------
251 Basic usage with alpha and beta parameters:
253 >>> import pandas as pd
254 >>> import numpy as np
255 >>> from gwtransport.utils import compute_time_edges
256 >>> from gwtransport.advection import gamma_infiltration_to_extraction
257 >>>
258 >>> # Create input data with aligned time edges
259 >>> dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
260 >>> tedges = compute_time_edges(
261 ... tedges=None, tstart=None, tend=dates, number_of_bins=len(dates)
262 ... )
263 >>>
264 >>> # Create output time edges (can be different alignment)
265 >>> cout_dates = pd.date_range(start="2020-01-05", end="2020-01-15", freq="D")
266 >>> cout_tedges = compute_time_edges(
267 ... tedges=None, tstart=None, tend=cout_dates, number_of_bins=len(cout_dates)
268 ... )
269 >>>
270 >>> # Input concentration and flow (same length, aligned with tedges)
271 >>> cin = pd.Series(np.ones(len(dates)), index=dates)
272 >>> flow = pd.Series(np.ones(len(dates)) * 100, index=dates) # 100 m³/day
273 >>>
274 >>> # Run gamma_infiltration_to_extraction with alpha/beta parameters
275 >>> cout = gamma_infiltration_to_extraction(
276 ... cin=cin,
277 ... flow=flow,
278 ... tedges=tedges,
279 ... cout_tedges=cout_tedges,
280 ... alpha=10.0,
281 ... beta=10.0,
282 ... n_bins=5,
283 ... )
284 >>> cout.shape
285 (11,)
287 Using mean and std parameters instead:
289 >>> cout = gamma_infiltration_to_extraction(
290 ... cin=cin,
291 ... flow=flow,
292 ... tedges=tedges,
293 ... cout_tedges=cout_tedges,
294 ... mean=100.0,
295 ... std=20.0,
296 ... n_bins=5,
297 ... )
298 """
299 bins = gamma.bins(mean=mean, std=std, loc=loc, alpha=alpha, beta=beta, n_bins=n_bins)
300 return infiltration_to_extraction(
301 cin=cin,
302 flow=flow,
303 tedges=tedges,
304 cout_tedges=cout_tedges,
305 aquifer_pore_volumes=bins["expected_values"],
306 retardation_factor=retardation_factor,
307 spinup=spinup,
308 )
311def gamma_extraction_to_infiltration(
312 *,
313 cout: npt.ArrayLike,
314 flow: npt.ArrayLike,
315 tedges: pd.DatetimeIndex,
316 cout_tedges: pd.DatetimeIndex,
317 mean: float | None = None,
318 std: float | None = None,
319 loc: float = 0.0,
320 alpha: float | None = None,
321 beta: float | None = None,
322 n_bins: int = 100,
323 retardation_factor: float = 1.0,
324 regularization_strength: float = 1e-10,
325 spinup: str | None = "constant",
326) -> npt.NDArray[np.floating]:
327 """
328 Compute the concentration of the infiltrating water from extracted water (deconvolution).
330 The compound is retarded in the aquifer with a retardation factor. The residence
331 time is computed based on the flow rate of the water in the aquifer and the pore volume
332 of the aquifer. The aquifer pore volume is approximated by a (shifted) gamma distribution
333 parameterized by either (mean, std, loc) or (alpha, beta, loc).
335 This function inverts the forward flow-weighted averaging (deconvolution).
336 It is symmetric to gamma_infiltration_to_extraction.
338 Provide either (mean, std) or (alpha, beta); ``loc`` is optional and defaults to 0.
340 Parameters
341 ----------
342 cout : array-like
343 Concentration of the compound in extracted water or temperature of extracted
344 water. The model assumes this value is constant over each interval
345 ``[cout_tedges[i], cout_tedges[i+1])``.
346 flow : array-like
347 Flow rate of water in the aquifer [m³/day]. The model assumes this value is
348 constant over each interval ``[tedges[i], tedges[i+1])``.
349 tedges : pandas.DatetimeIndex
350 Time edges for cin (output) and flow data.
351 Has a length of one more than `flow`.
352 cout_tedges : pandas.DatetimeIndex
353 Time edges for the cout data.
354 Has a length of one more than `cout`.
355 mean : float, optional
356 Mean of the gamma distribution of the aquifer pore volume. Must be strictly
357 greater than ``loc``.
358 std : float, optional
359 Standard deviation of the gamma distribution of the aquifer pore volume
360 (invariant under the ``loc`` shift).
361 loc : float, optional
362 Location (minimum pore volume) of the gamma distribution. Must satisfy
363 ``0 <= loc < mean``. Default is ``0.0``.
364 alpha : float, optional
365 Shape parameter of gamma distribution of the aquifer pore volume (must be > 0).
366 beta : float, optional
367 Scale parameter of gamma distribution of the aquifer pore volume (must be > 0).
368 n_bins : int, optional
369 Number of bins to discretize the gamma distribution. Default 100.
370 retardation_factor : float, optional
371 Retardation factor of the compound in the aquifer (default 1.0).
372 Values > 1.0 indicate slower transport due to sorption/interaction.
373 regularization_strength : float, optional
374 Tikhonov regularization parameter λ. See
375 :func:`extraction_to_infiltration` for details. Default is 1e-10.
376 spinup : {"constant"} | None, optional
377 Forwarded to :func:`extraction_to_infiltration`. Default
378 ``"constant"`` warm-starts the system before ``tedges[0]``.
380 Returns
381 -------
382 numpy.ndarray
383 Concentration of the compound in the infiltrating water, or temperature. Same units as cout.
385 See Also
386 --------
387 extraction_to_infiltration : Deconvolution with explicit pore volume distribution
388 gamma_infiltration_to_extraction : Forward operation (flow-weighted averaging)
389 gwtransport.gamma.bins : Create gamma distribution bins
390 gwtransport.diffusion.extraction_to_infiltration : Deconvolution with microdispersion and molecular diffusion
391 :ref:`concept-gamma-distribution` : Two-parameter pore volume model
392 :ref:`assumption-gamma-distribution` : When gamma distribution is adequate
394 Notes
395 -----
396 The APVD is only time-invariant under the steady-streamlines assumption
397 (see :ref:`assumption-steady-streamlines`).
399 The spreading from the gamma-distributed pore volumes represents macrodispersion
400 (aquifer-scale heterogeneity). When ``std`` comes from calibration on measurements,
401 it absorbs all mixing: macrodispersion, microdispersion, and an average molecular
402 diffusion contribution. When calibrating with the diffusion module, these three
403 components are taken into account separately. When ``std`` comes from streamline
404 analysis, it represents macrodispersion only; microdispersion and molecular diffusion
405 can be added via :mod:`gwtransport.diffusion_fast` or :mod:`gwtransport.diffusion`.
407 For cross-compound prediction (calibrating on temperature and predicting a solute),
408 calibrate with :mod:`gwtransport.diffusion_fast` so the three contributions are
409 tracked separately rather than lumped into a single calibrated ``std``.
410 See :ref:`concept-dispersion-scales` for background.
412 Examples
413 --------
414 Basic usage with alpha and beta parameters:
416 >>> import pandas as pd
417 >>> import numpy as np
418 >>> from gwtransport.utils import compute_time_edges
419 >>> from gwtransport.advection import gamma_extraction_to_infiltration
420 >>>
421 >>> # Create cin/flow time edges
422 >>> cin_dates = pd.date_range(start="2019-12-25", end="2020-01-15", freq="D")
423 >>> tedges = compute_time_edges(
424 ... tedges=None, tstart=None, tend=cin_dates, number_of_bins=len(cin_dates)
425 ... )
426 >>>
427 >>> # Create cout time edges
428 >>> cout_dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
429 >>> cout_tedges = compute_time_edges(
430 ... tedges=None, tstart=None, tend=cout_dates, number_of_bins=len(cout_dates)
431 ... )
432 >>>
433 >>> # Input concentration and flow
434 >>> cout = np.ones(len(cout_dates))
435 >>> flow = np.ones(len(cin_dates)) * 100 # 100 m³/day
436 >>>
437 >>> # Run gamma_extraction_to_infiltration with alpha/beta parameters
438 >>> cin = gamma_extraction_to_infiltration(
439 ... cout=cout,
440 ... flow=flow,
441 ... tedges=tedges,
442 ... cout_tedges=cout_tedges,
443 ... alpha=10.0,
444 ... beta=10.0,
445 ... n_bins=5,
446 ... )
447 >>> cin.shape
448 (22,)
450 Using mean and std parameters instead:
452 >>> cin = gamma_extraction_to_infiltration(
453 ... cout=cout,
454 ... flow=flow,
455 ... tedges=tedges,
456 ... cout_tedges=cout_tedges,
457 ... mean=100.0,
458 ... std=20.0,
459 ... n_bins=5,
460 ... )
461 """
462 bins = gamma.bins(mean=mean, std=std, loc=loc, alpha=alpha, beta=beta, n_bins=n_bins)
463 return extraction_to_infiltration(
464 cout=cout,
465 flow=flow,
466 tedges=tedges,
467 cout_tedges=cout_tedges,
468 aquifer_pore_volumes=bins["expected_values"],
469 retardation_factor=retardation_factor,
470 regularization_strength=regularization_strength,
471 spinup=spinup,
472 )
475def infiltration_to_extraction(
476 *,
477 cin: npt.ArrayLike,
478 flow: npt.ArrayLike,
479 tedges: pd.DatetimeIndex,
480 cout_tedges: pd.DatetimeIndex,
481 aquifer_pore_volumes: npt.ArrayLike,
482 retardation_factor: float = 1.0,
483 spinup: str | None = "constant",
484) -> npt.NDArray[np.floating]:
485 """
486 Compute the concentration of the extracted water using flow-weighted advection.
488 This function implements an infiltration to extraction advection model where cin and flow values
489 correspond to the same aligned time bins defined by tedges.
491 Pure advection is volume-stationary, so the weights are built on the
492 cumulative-throughflow-volume axis rather than by inverting residence times:
494 1. Map the cin and cout time edges to cumulative throughflow volume.
495 2. Back-project each cout bin by every retarded pore volume to its
496 infiltration-time source window. The window spans one cout bin's worth of
497 volume, so it overlaps only a narrow band of cin bins.
498 3. Compute the flow-weighted time overlap of each window with those cin bins,
499 normalize per streamtube (each row sums to 1), and average over the
500 streamtubes whose source window lies fully inside the cin range.
503 Parameters
504 ----------
505 cin : array-like
506 Concentration values of infiltrating water or temperature [concentration units].
507 Length must match the number of time bins defined by tedges. The model assumes
508 this value is constant over each interval ``[tedges[i], tedges[i+1])``.
509 flow : array-like
510 Flow rate values in the aquifer [m³/day].
511 Length must match cin and the number of time bins defined by tedges. The model
512 assumes this value is constant over each interval ``[tedges[i], tedges[i+1])``.
513 tedges : pandas.DatetimeIndex
514 Time edges defining bins for both cin and flow data. Has length of
515 len(cin) + 1 and len(flow) + 1.
516 cout_tedges : pandas.DatetimeIndex
517 Time edges for output data bins. Has length of desired output + 1.
518 Can have different time alignment and resolution than tedges.
519 aquifer_pore_volumes : array-like
520 Array of aquifer pore volumes [m³] representing the distribution
521 of residence times in the aquifer system.
522 retardation_factor : float, optional
523 Retardation factor of the compound in the aquifer (default 1.0).
524 Values > 1.0 indicate slower transport due to sorption/interaction.
525 spinup : {"constant"} | None, optional
526 How to treat cout bins where one or more streamtube source windows
527 fall outside the cin time range. Default is ``"constant"``.
529 - ``"constant"`` — warm-start: shift ``tedges[0]`` backward by
530 ``retardation_factor * max(aquifer_pore_volumes) / flow[0]`` and
531 treat cin and flow as constant at their first value over the
532 extended window. The forward strict-validity logic then has no
533 NaN cout bins from spin-up; right-edge spin-up (cout extending
534 past the cin range) is unchanged.
535 - ``None`` — strict mass-conservation: NaN whenever any streamtube
536 has not fully broken through into the cin range, or extraction
537 flow during the bin is zero. Bundle row sums to 1 across cin.
539 Returns
540 -------
541 numpy.ndarray
542 Flow-weighted concentration in the extracted water. Same units
543 as cin. Length equals ``len(cout_tedges) - 1``. NaN values mark
544 cout bins where the chosen ``spinup`` policy is not satisfied:
545 the default ``"constant"`` leaves NaN for any cout bin extending
546 past the end of the flow record (a cout edge beyond
547 ``tedges[-1]``, whose back-projected source window leaves the cin
548 range) and for zero-throughflow bins; ``spinup=None`` additionally
549 NaNs left-edge spin-up bins.
551 Raises
552 ------
553 ValueError
554 If tedges length doesn't match cin/flow arrays plus one, or if
555 infiltration time edges become non-monotonic (invalid input conditions).
557 See Also
558 --------
559 gamma_infiltration_to_extraction : Transport with gamma-distributed pore volumes
560 extraction_to_infiltration : Reverse operation (deconvolution)
561 gwtransport.residence_time.full : Compute residence times from flow and pore volume
562 gwtransport.residence_time.freundlich_retardation : Compute concentration-dependent retardation
563 :ref:`concept-pore-volume-distribution` : Background on aquifer heterogeneity modeling
564 :ref:`concept-transport-equation` : Flow-weighted averaging approach
566 Examples
567 --------
568 Basic usage with pandas Series:
570 >>> import pandas as pd
571 >>> import numpy as np
572 >>> from gwtransport.utils import compute_time_edges
573 >>> from gwtransport.advection import infiltration_to_extraction
574 >>>
575 >>> # Create input data
576 >>> dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
577 >>> tedges = compute_time_edges(
578 ... tedges=None, tstart=None, tend=dates, number_of_bins=len(dates)
579 ... )
580 >>>
581 >>> # Create output time edges (different alignment)
582 >>> cout_dates = pd.date_range(start="2020-01-05", end="2020-01-15", freq="D")
583 >>> cout_tedges = compute_time_edges(
584 ... tedges=None, tstart=None, tend=cout_dates, number_of_bins=len(cout_dates)
585 ... )
586 >>>
587 >>> # Input concentration and flow
588 >>> cin = pd.Series(np.ones(len(dates)), index=dates)
589 >>> flow = pd.Series(np.ones(len(dates)) * 100, index=dates) # 100 m³/day
590 >>>
591 >>> # Define distribution of aquifer pore volumes
592 >>> aquifer_pore_volumes = np.array([50, 100, 200]) # m³
593 >>>
594 >>> # Run infiltration_to_extraction
595 >>> cout = infiltration_to_extraction(
596 ... cin=cin,
597 ... flow=flow,
598 ... tedges=tedges,
599 ... cout_tedges=cout_tedges,
600 ... aquifer_pore_volumes=aquifer_pore_volumes,
601 ... )
602 >>> cout.shape
603 (11,)
605 With constant retardation factor (linear sorption):
607 >>> cout = infiltration_to_extraction(
608 ... cin=cin,
609 ... flow=flow,
610 ... tedges=tedges,
611 ... cout_tedges=cout_tedges,
612 ... aquifer_pore_volumes=aquifer_pore_volumes,
613 ... retardation_factor=2.0, # Compound moves twice as slowly
614 ... )
616 Note: For concentration-dependent retardation (nonlinear sorption),
617 use `infiltration_to_extraction_nonlinear_sorption` instead, as this
618 function only supports constant (float) retardation factors.
619 """
620 tedges = pd.DatetimeIndex(tedges)
621 cout_tedges = pd.DatetimeIndex(cout_tedges)
623 # Convert to arrays for vectorized operations
624 cin = np.asarray(cin)
625 flow = np.asarray(flow)
626 aquifer_pore_volumes = np.asarray(aquifer_pore_volumes)
628 _validate_advection_inputs(
629 tedges=tedges,
630 flow=flow,
631 retardation_factor=retardation_factor,
632 aquifer_pore_volumes=aquifer_pore_volumes,
633 cin_values=cin,
634 )
636 weight_tedges, weight_flow, weight_cin, _, _ = _resolve_spinup_inputs(
637 spinup,
638 tedges=tedges,
639 flow=flow,
640 aquifer_pore_volumes=aquifer_pore_volumes,
641 retardation_factor=retardation_factor,
642 cin=cin,
643 )
644 assert weight_cin is not None # noqa: S101 -- narrowed: cin was passed in
645 band_vals, col_start, contributing_bins, zero_flow_cout = _infiltration_to_extraction_weights(
646 tedges=weight_tedges,
647 cout_tedges=cout_tedges,
648 aquifer_pore_volumes=aquifer_pore_volumes,
649 flow=weight_flow,
650 retardation_factor=retardation_factor,
651 )
652 weights, invalid_mask = _resolve_spinup_mask(
653 band_vals=band_vals,
654 contributing_bins=contributing_bins,
655 zero_flow_cout=zero_flow_cout,
656 n_pv=len(aquifer_pore_volumes),
657 )
659 # Banded flow-weighted average: row k contributes cin over its narrow band only.
660 # Out-of-range band slots carry zero weight, so the clipped gather is harmless.
661 n_cin = len(weight_cin)
662 cols = np.clip(col_start[:, None] + np.arange(weights.shape[1]), 0, n_cin - 1)
663 out = np.einsum("kb,kb->k", weights, weight_cin[cols])
665 # Invalid rows (cout bins where the spin-up policy is not satisfied or
666 # where extraction flow was zero) become NaN.
667 out[invalid_mask] = np.nan
669 return out
672def extraction_to_infiltration(
673 *,
674 cout: npt.ArrayLike,
675 flow: npt.ArrayLike,
676 tedges: pd.DatetimeIndex,
677 cout_tedges: pd.DatetimeIndex,
678 aquifer_pore_volumes: npt.ArrayLike,
679 retardation_factor: float = 1.0,
680 regularization_strength: float = 1e-10,
681 spinup: str | float | None = "constant",
682) -> npt.NDArray[np.floating]:
683 """
684 Compute the concentration of the infiltrating water from extracted water (deconvolution).
686 Inverts the forward transport model by solving the linear system
687 ``W_forward @ cin = cout`` where ``W_forward`` is the weight matrix from
688 :func:`infiltration_to_extraction`. Uses Tikhonov regularization to
689 smoothly blend data fitting with a physically motivated target
690 (transpose-and-normalize of the forward matrix).
692 Well-determined modes (large singular values relative to √λ) are
693 dominated by the data; poorly-determined modes are pulled toward the
694 target. This avoids edge oscillations and is less sensitive to the
695 regularization parameter than truncated SVD (``rcond``).
697 Parameters
698 ----------
699 cout : array-like
700 Concentration values of extracted water [concentration units].
701 Length must match the number of time bins defined by cout_tedges. The model
702 assumes this value is constant over each interval
703 ``[cout_tedges[i], cout_tedges[i+1])``.
704 flow : array-like
705 Flow rate values in the aquifer [m³/day].
706 Length must match the number of time bins defined by tedges. The model assumes
707 this value is constant over each interval ``[tedges[i], tedges[i+1])``.
708 tedges : pandas.DatetimeIndex
709 Time edges defining bins for both cin (output) and flow data. Has length of
710 len(flow) + 1. Output cin has length len(tedges) - 1.
711 cout_tedges : pandas.DatetimeIndex
712 Time edges for cout data bins. Has length of len(cout) + 1.
713 Can have different time alignment and resolution than tedges.
714 aquifer_pore_volumes : array-like
715 Array of aquifer pore volumes [m³] representing the distribution
716 of residence times in the aquifer system.
717 retardation_factor : float, optional
718 Retardation factor of the compound in the aquifer (default 1.0).
719 Values > 1.0 indicate slower transport due to sorption/interaction.
720 regularization_strength : float, optional
721 Tikhonov regularization parameter λ. Controls the tradeoff between
722 fitting the data (``||W cin - cout||²``) and staying close to the
723 regularization target (``λ ||cin - cin_target||²``). The target is
724 the transpose-and-normalize of the forward matrix applied to cout.
726 Larger values trust the target more (smoother, more biased); smaller
727 values trust the data more (noisier, less biased). The solution
728 varies continuously with λ. Default is 1e-10.
730 A good starting value for noisy data is
731 ``λ ≈ (noise_std / signal_amplitude)²``. For example, temperature
732 data with 0.05 °C noise and ~10 °C seasonal amplitude suggests
733 ``regularization_strength ≈ (0.05 / 10)² ≈ 2.5e-5``. Increase by
734 a factor of 2-10 for additional smoothing. For noiseless synthetic
735 data (e.g., roundtrip tests), the default 1e-10 preserves machine
736 precision.
737 spinup : {"constant"} | None, optional
738 Spin-up policy applied when building the forward weight matrix
739 used to set up the inverse problem. Same semantics as in
740 :func:`infiltration_to_extraction`; default ``"constant"`` shifts
741 ``tedges[0]`` backward by ``retardation_factor *
742 max(aquifer_pore_volumes) / flow[0]`` so the inverse problem has
743 no spin-up zero-rows for cout bins inside the original tedges
744 range. The warm-start prefix is solved for internally but dropped
745 before returning, so the output cin stays aligned with the
746 user-provided ``tedges`` (length ``len(tedges) - 1``), not the
747 padded grid. Passing ``None`` keeps the strict-validity behavior
748 (zero-rows in W from incomplete breakthrough).
750 Returns
751 -------
752 numpy.ndarray
753 Concentration in the infiltrating water. Same units as cout.
754 Length equals len(tedges) - 1 (unchanged whether or not
755 ``spinup="constant"`` shifted ``tedges[0]``). NaN values indicate
756 cin bins with no temporal overlap with the extraction data. The
757 forward weight matrix used to set up the inverse problem treats
758 spin-up and zero-flow cout bins as zero-rows according to the
759 ``spinup`` policy.
761 Raises
762 ------
763 ValueError
764 If tedges length doesn't match flow plus one, if cout_tedges length
765 doesn't match cout plus one, or if flow contains NaN.
767 See Also
768 --------
769 gamma_extraction_to_infiltration : Deconvolution with gamma-distributed pore volumes
770 infiltration_to_extraction : Forward operation (flow-weighted averaging)
771 gwtransport.residence_time.full : Compute residence times from flow and pore volume
772 gwtransport.utils.solve_tikhonov : Solver used for inversion
773 :ref:`concept-pore-volume-distribution` : Background on aquifer heterogeneity modeling
774 :ref:`concept-transport-equation` : Flow-weighted averaging approach
776 Notes
777 -----
778 NaN values in ``cout`` mark measurement gaps (e.g. sparse lab samples).
779 Their rows are excluded from the banded Tikhonov solve, matching
780 :func:`gwtransport.deposition.extraction_to_deposition`; cin bins
781 constrained only by gapped ``cout`` bins are returned as NaN.
783 Examples
784 --------
785 Basic usage with pandas Series:
787 >>> import pandas as pd
788 >>> import numpy as np
789 >>> from gwtransport.utils import compute_time_edges
790 >>> from gwtransport.advection import extraction_to_infiltration
791 >>>
792 >>> # Create cin/flow time edges
793 >>> cin_dates = pd.date_range(start="2019-12-25", end="2020-01-15", freq="D")
794 >>> tedges = compute_time_edges(
795 ... tedges=None, tstart=None, tend=cin_dates, number_of_bins=len(cin_dates)
796 ... )
797 >>>
798 >>> # Create cout time edges
799 >>> cout_dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
800 >>> cout_tedges = compute_time_edges(
801 ... tedges=None, tstart=None, tend=cout_dates, number_of_bins=len(cout_dates)
802 ... )
803 >>>
804 >>> # Input concentration and flow
805 >>> cout = np.ones(len(cout_dates))
806 >>> flow = np.ones(len(cin_dates)) * 100 # 100 m³/day
807 >>>
808 >>> # Define distribution of aquifer pore volumes
809 >>> aquifer_pore_volumes = np.array([50, 100, 200]) # m³
810 >>>
811 >>> # Run extraction_to_infiltration
812 >>> cin = extraction_to_infiltration(
813 ... cout=cout,
814 ... flow=flow,
815 ... tedges=tedges,
816 ... cout_tedges=cout_tedges,
817 ... aquifer_pore_volumes=aquifer_pore_volumes,
818 ... )
819 >>> cin.shape
820 (22,)
822 Round-trip reconstruction (symmetric with infiltration_to_extraction).
823 The default ``spinup="constant"`` warm-starts the left edge; the cout
824 window must therefore stay inside the cin window with margin matching
825 the longest residence time on the right (forward NaN at the right
826 edge would otherwise be rejected by ``extraction_to_infiltration``):
828 >>> from gwtransport.advection import infiltration_to_extraction
829 >>> rt_cout_dates = pd.date_range(start="2020-01-01", end="2020-01-10", freq="D")
830 >>> rt_cout_tedges = compute_time_edges(
831 ... tedges=None,
832 ... tstart=None,
833 ... tend=rt_cout_dates,
834 ... number_of_bins=len(rt_cout_dates),
835 ... )
836 >>> cin_original = np.sin(np.linspace(0, 2 * np.pi, len(cin_dates))) + 2
837 >>> cout_rt = infiltration_to_extraction(
838 ... cin=cin_original,
839 ... flow=flow,
840 ... tedges=tedges,
841 ... cout_tedges=rt_cout_tedges,
842 ... aquifer_pore_volumes=aquifer_pore_volumes,
843 ... )
844 >>> cin_recovered = extraction_to_infiltration(
845 ... cout=cout_rt,
846 ... flow=flow,
847 ... tedges=tedges,
848 ... cout_tedges=rt_cout_tedges,
849 ... aquifer_pore_volumes=aquifer_pore_volumes,
850 ... )
851 """
852 tedges = pd.DatetimeIndex(tedges)
853 cout_tedges = pd.DatetimeIndex(cout_tedges)
855 # Convert to arrays for vectorized operations
856 cout = np.asarray(cout)
857 flow = np.asarray(flow)
859 _validate_advection_inputs(
860 tedges=tedges,
861 flow=flow,
862 retardation_factor=retardation_factor,
863 aquifer_pore_volumes=aquifer_pore_volumes,
864 cout_values=cout,
865 cout_tedges=cout_tedges,
866 )
868 aquifer_pore_volumes = np.asarray(aquifer_pore_volumes)
870 weight_tedges, weight_flow, _, _, n_pad = _resolve_spinup_inputs(
871 spinup,
872 tedges=tedges,
873 flow=flow,
874 aquifer_pore_volumes=aquifer_pore_volumes,
875 retardation_factor=retardation_factor,
876 )
877 n_cin_padded = len(weight_tedges) - 1
879 band_vals, col_start, contributing_bins, zero_flow_cout = _infiltration_to_extraction_weights(
880 tedges=weight_tedges,
881 cout_tedges=cout_tedges,
882 aquifer_pore_volumes=aquifer_pore_volumes,
883 flow=weight_flow,
884 retardation_factor=retardation_factor,
885 )
886 band_vals, _ = _resolve_spinup_mask(
887 band_vals=band_vals,
888 contributing_bins=contributing_bins,
889 zero_flow_cout=zero_flow_cout,
890 n_pv=len(aquifer_pore_volumes),
891 )
893 cin_padded = solve_inverse_transport_banded(
894 band_vals=band_vals,
895 col_start=col_start,
896 observed=cout,
897 n_output=n_cin_padded,
898 regularization_strength=regularization_strength,
899 )
900 # Drop warm-start prefix so the output aligns with the user-provided tedges.
901 return cin_padded[n_pad:]
904def _validate_front_tracking_inputs(
905 *,
906 cin: npt.ArrayLike,
907 flow: npt.ArrayLike,
908 tedges: pd.DatetimeIndex,
909 cout_tedges: pd.DatetimeIndex,
910 aquifer_pore_volumes: npt.ArrayLike,
911 freundlich_k: float | None,
912 freundlich_n: float | None,
913 bulk_density: float | None,
914 porosity: float | None,
915 retardation_factor: float | None,
916 langmuir_s_max: float | None,
917 langmuir_k_l: float | None,
918) -> tuple[
919 npt.NDArray[np.float64],
920 npt.NDArray[np.float64],
921 pd.DatetimeIndex,
922 pd.DatetimeIndex,
923 npt.NDArray[np.float64],
924 SorptionModel,
925 npt.NDArray[np.floating],
926]:
927 """Validate inputs and create sorption object for front tracking functions.
929 Returns
930 -------
931 tuple
932 Validated and converted inputs: (cin, flow, tedges, cout_tedges,
933 aquifer_pore_volumes, sorption, cout_tedges_days).
935 Raises
936 ------
937 ValueError
938 If array lengths are inconsistent, values are non-physical (negative
939 concentrations, negative flows, NaN values, non-positive pore
940 volumes), retardation_factor < 1, Freundlich or Langmuir parameters
941 are missing or non-positive, freundlich_n equals 1, or physical
942 parameters are invalid.
943 """
944 cin = np.asarray(cin, dtype=float)
945 flow = np.asarray(flow, dtype=float)
946 tedges = pd.DatetimeIndex(tedges)
947 cout_tedges = pd.DatetimeIndex(cout_tedges)
948 aquifer_pore_volumes = np.asarray(aquifer_pore_volumes, dtype=float)
950 _validate_tedges_parity(tedges, cin, tedges_name="tedges", values_name="cin")
951 _validate_tedges_parity(tedges, flow, tedges_name="tedges", values_name="flow")
952 _validate_no_nan(cin, name="cin")
953 _validate_no_nan(flow, name="flow")
954 _validate_non_negative_array(cin, name="cin")
955 _validate_non_negative_array(flow, name="flow", message="flow must be non-negative (negative flow not supported)")
956 if aquifer_pore_volumes.size == 0:
957 msg = "aquifer_pore_volumes must not be empty"
958 raise ValueError(msg)
959 _validate_positive_array(aquifer_pore_volumes, name="aquifer_pore_volumes")
961 # Convert cout_tedges to days (relative to tedges[0]) for output computation
962 cout_tedges_days = tedges_to_days(cout_tedges, ref=tedges[0])
964 # Determine which sorption model is requested
965 has_retardation = retardation_factor is not None
966 has_freundlich = freundlich_k is not None or freundlich_n is not None
967 has_langmuir = langmuir_s_max is not None or langmuir_k_l is not None
968 n_models = has_retardation + has_freundlich + has_langmuir
970 if n_models == 0:
971 msg = (
972 "Must provide one of: retardation_factor, Freundlich parameters "
973 "(freundlich_k, freundlich_n, bulk_density, porosity), or Langmuir parameters "
974 "(langmuir_s_max, langmuir_k_l, bulk_density, porosity)"
975 )
976 raise ValueError(msg)
977 if n_models > 1:
978 msg = "Only one sorption model can be specified (retardation_factor, Freundlich, or Langmuir)"
979 raise ValueError(msg)
981 # Create sorption object
982 if retardation_factor is not None:
983 _validate_retardation_factor(retardation_factor)
984 sorption: SorptionModel = ConstantRetardation(retardation_factor=retardation_factor)
985 elif has_freundlich:
986 if freundlich_k is None or freundlich_n is None or bulk_density is None or porosity is None:
987 msg = "All Freundlich parameters required (freundlich_k, freundlich_n, bulk_density, porosity)"
988 raise ValueError(msg)
989 if freundlich_k <= 0 or freundlich_n <= 0:
990 msg = "Freundlich parameters must be positive"
991 raise ValueError(msg)
992 if abs(freundlich_n - 1.0) < EPSILON_FREUNDLICH_N:
993 msg = "freundlich_n = 1 not supported (use retardation_factor for linear case)"
994 raise ValueError(msg)
995 if bulk_density <= 0 or not 0 < porosity < 1:
996 msg = "Invalid physical parameters"
997 raise ValueError(msg)
999 sorption = FreundlichSorption(
1000 k_f=freundlich_k,
1001 n=freundlich_n,
1002 bulk_density=bulk_density,
1003 porosity=porosity,
1004 )
1005 else:
1006 if langmuir_s_max is None or langmuir_k_l is None or bulk_density is None or porosity is None:
1007 msg = "All Langmuir parameters required (langmuir_s_max, langmuir_k_l, bulk_density, porosity)"
1008 raise ValueError(msg)
1009 if langmuir_s_max <= 0 or langmuir_k_l <= 0:
1010 msg = "Langmuir parameters must be positive"
1011 raise ValueError(msg)
1012 if bulk_density <= 0 or not 0 < porosity < 1:
1013 msg = "Invalid physical parameters"
1014 raise ValueError(msg)
1016 sorption = LangmuirSorption(
1017 s_max=langmuir_s_max,
1018 k_l=langmuir_k_l,
1019 bulk_density=bulk_density,
1020 porosity=porosity,
1021 )
1023 return cin, flow, tedges, cout_tedges, aquifer_pore_volumes, sorption, cout_tedges_days
1026def _flow_weighted_front_tracking_output(
1027 cout_tedges_days: npt.NDArray[np.floating],
1028 flow_tedges_days: npt.NDArray[np.floating],
1029 flow: npt.NDArray[np.floating],
1030 v_outlet: float,
1031 waves: list,
1032 sorption: SorptionModel,
1033 theta_edges: npt.NDArray[np.floating],
1034 cin: npt.NDArray[np.floating],
1035) -> npt.NDArray[np.floating]:
1036 """Compute flow-weighted bin-averaged concentration from front-tracking output.
1038 Splits output bins at flow boundaries so that Q is constant within each
1039 sub-bin, then combines sub-bins with flow-weighting:
1040 ``c_avg = Σ(Q_k · c_k · dt_k) / Σ(Q_k · dt_k)``.
1042 Internally translates the output ``t``-bin edges to θ via the same
1043 ``(flow_tedges_days, theta_edges)`` map the tracker built, and calls
1044 :func:`compute_bin_averaged_concentration_exact` in θ-coordinates.
1046 Parameters
1047 ----------
1048 cout_tedges_days : ndarray
1049 Output time bin edges [days from reference].
1050 flow_tedges_days : ndarray
1051 Flow time bin edges [days from reference] (length ``len(flow) + 1``).
1052 flow : ndarray
1053 Flow rate per flow bin [m³/day].
1054 v_outlet : float
1055 Outlet volume position [m³].
1056 waves : list
1057 Wave list from front tracking simulation.
1058 sorption : object
1059 Sorption model.
1060 theta_edges : ndarray
1061 Cumulative-flow edges at the flow-bin boundaries [m³]
1062 (length ``len(flow) + 1``).
1063 cin : ndarray
1064 Infiltration concentration values, one per flow bin. Passed directly to
1065 :func:`compute_bin_averaged_concentration_exact`.
1067 Returns
1068 -------
1069 ndarray
1070 Flow-weighted bin-averaged concentrations. Length = len(cout_tedges_days) - 1.
1072 Notes
1073 -----
1074 Zero-flow sub-bins are dropped from the averaging only; a cin step injected
1075 during a zero-flow bin still enters the tracker as a wave at a degenerate θ
1076 and can corrupt adjacent output bins.
1077 """
1078 inner_flow_edges = flow_tedges_days[
1079 (flow_tedges_days > cout_tedges_days[0]) & (flow_tedges_days < cout_tedges_days[-1])
1080 ]
1081 fine_edges = np.unique(np.concatenate([cout_tedges_days, inner_flow_edges]))
1083 # np.interp clips on both sides; extrapolate the θ map past either flow edge at
1084 # the adjacent-bin flow (matches the FrontTrackerState.theta_at_t rule). Without
1085 # this, out-of-window fine_edges collapse to a duplicate θ. A θ edge ≤ 0
1086 # downstream reads back as m = 0 → 0.0 (the documented out-of-range contract).
1087 fine_theta_edges = np.interp(fine_edges, flow_tedges_days, theta_edges)
1088 underflow = fine_edges < flow_tedges_days[0]
1089 if underflow.any():
1090 fine_theta_edges[underflow] = theta_edges[0] - (flow_tedges_days[0] - fine_edges[underflow]) * float(flow[0])
1091 overflow = fine_edges > flow_tedges_days[-1]
1092 if overflow.any():
1093 fine_theta_edges[overflow] = theta_edges[-1] + (fine_edges[overflow] - flow_tedges_days[-1]) * float(flow[-1])
1095 # A zero-flow input span leaves θ stationary, so its sub-bins have zero width in
1096 # θ and zero q·dt weight. Drop them before the exact averaging (which rejects
1097 # non-positive-width bins); they read back as 0 and carry no weight in this
1098 # averaging (see Notes for the degenerate-θ caveat). Consecutive kept bins stay
1099 # contiguous because the dropped bins share their neighbours' θ value.
1100 theta_lo, theta_hi = fine_theta_edges[:-1], fine_theta_edges[1:]
1101 nondegenerate = theta_hi > theta_lo
1102 c_fine = np.zeros(theta_lo.shape)
1103 if nondegenerate.any():
1104 kept_edges = np.concatenate([theta_lo[nondegenerate], theta_hi[nondegenerate][-1:]])
1105 c_fine[nondegenerate] = compute_bin_averaged_concentration_exact(
1106 theta_bin_edges=kept_edges,
1107 v_outlet=v_outlet,
1108 waves=waves,
1109 sorption=sorption,
1110 cin=cin,
1111 theta_edges_inlet=theta_edges,
1112 )
1114 # Map each fine sub-bin to its flow value. side="right" enforces the
1115 # half-open [t_k, t_{k+1}) bin convention if a midpoint ever lands
1116 # exactly on an inner flow edge (does not happen for np.unique-derived
1117 # midpoints in practice, but is defensible against floating-point drift).
1118 fine_mids = (fine_edges[:-1] + fine_edges[1:]) / 2
1119 flow_idx = np.searchsorted(flow_tedges_days[1:], fine_mids, side="right")
1120 flow_idx = np.clip(flow_idx, 0, len(flow) - 1)
1121 q_fine = flow[flow_idx]
1122 dt_fine = np.diff(fine_edges)
1124 # Map each fine sub-bin to its original output bin. Same side="right"
1125 # rationale as above.
1126 cout_bin_idx = np.searchsorted(cout_tedges_days[1:], fine_mids, side="right")
1127 cout_bin_idx = np.clip(cout_bin_idx, 0, len(cout_tedges_days) - 2)
1129 # Vectorized per-bin flow-weighted average:
1130 # c_out[k] = sum_i (Q_i * c_i * dt_i) / sum_i (Q_i * dt_i) for fine sub-bins i in bin k
1131 n_cout = len(cout_tedges_days) - 1
1132 qdt_product = q_fine * dt_fine
1133 cqdt_product = c_fine * qdt_product
1134 denominator = np.bincount(cout_bin_idx, weights=qdt_product, minlength=n_cout)
1135 numerator = np.bincount(cout_bin_idx, weights=cqdt_product, minlength=n_cout)
1136 # A zero-throughflow output bin (all overlapping input bins have zero flow) has
1137 # an undefined flow-weighted average: emit NaN, matching the linear sibling.
1138 # Pre-record bins keep positive throughflow, so their 0-mass windows read as 0.0.
1139 c_out = np.full(n_cout, np.nan)
1140 valid = denominator > 0
1141 c_out[valid] = numerator[valid] / denominator[valid]
1142 return c_out
1145def infiltration_to_extraction_nonlinear_sorption(
1146 *,
1147 cin: npt.ArrayLike,
1148 flow: npt.ArrayLike,
1149 tedges: pd.DatetimeIndex,
1150 cout_tedges: pd.DatetimeIndex,
1151 aquifer_pore_volumes: npt.ArrayLike,
1152 freundlich_k: float | None = None,
1153 freundlich_n: float | None = None,
1154 bulk_density: float | None = None,
1155 porosity: float | None = None,
1156 retardation_factor: float | None = None,
1157 langmuir_s_max: float | None = None,
1158 langmuir_k_l: float | None = None,
1159 max_iterations: int = 10000,
1160) -> tuple[npt.NDArray[np.floating], list[dict]]:
1161 """
1162 Compute extracted concentration with complete diagnostic information.
1164 Returns both bin-averaged concentrations and detailed simulation structure for each pore volume.
1166 Exactly one sorption model must be specified:
1168 - ``retardation_factor`` for constant (linear) retardation.
1169 - ``freundlich_k`` + ``freundlich_n`` + ``bulk_density`` + ``porosity`` for
1170 Freundlich isotherm.
1171 - ``langmuir_s_max`` + ``langmuir_k_l`` + ``bulk_density`` + ``porosity`` for
1172 Langmuir isotherm.
1174 Parameters
1175 ----------
1176 cin : array-like
1177 Infiltration concentration [mg/L or any units].
1178 Length = len(tedges) - 1. The model assumes this value is constant over each
1179 interval ``[tedges[i], tedges[i+1])``.
1180 flow : array-like
1181 Flow rate [m³/day]. Must be non-negative.
1182 Length = len(tedges) - 1. The model assumes this value is constant over each
1183 interval ``[tedges[i], tedges[i+1])``.
1184 tedges : pandas.DatetimeIndex
1185 Time bin edges. Length = len(cin) + 1.
1186 cout_tedges : pandas.DatetimeIndex
1187 Output time bin edges. Can be different from tedges.
1188 Length = number of output bins + 1 (n+1 edges for n output values).
1189 aquifer_pore_volumes : array-like
1190 Array of aquifer pore volumes [m³] representing the distribution
1191 of residence times in the aquifer system. Each pore volume must be positive.
1192 freundlich_k : float, optional
1193 Freundlich coefficient [(m³/kg)^(1/n)]. Must be positive. Same convention (isotherm
1194 ``s = k_f * C^(1/n)``) as :func:`gwtransport.residence_time.freundlich_retardation`.
1195 freundlich_n : float, optional
1196 Freundlich exponent [-]. Must be positive and != 1.
1197 bulk_density : float, optional
1198 Bulk density [kg/m³]. Must be positive.
1199 Shared by Freundlich and Langmuir models.
1200 porosity : float, optional
1201 Porosity [-]. Must be in (0, 1).
1202 Shared by Freundlich and Langmuir models.
1203 retardation_factor : float, optional
1204 Constant retardation factor [-]. Must be >= 1.0.
1205 langmuir_s_max : float, optional
1206 Langmuir maximum sorption capacity [mg/kg]. Must be positive.
1207 langmuir_k_l : float, optional
1208 Langmuir half-saturation constant [mg/L]. Must be positive.
1209 max_iterations : int, optional
1210 Maximum number of events. Default 10000.
1212 Returns
1213 -------
1214 cout : numpy.ndarray
1215 Flow-weighted concentrations averaged across all pore volumes. Output
1216 bins whose source window leaves the inlet flow record (e.g. cout bins
1217 before first breakthrough, or extending past the flow record) are
1218 returned as ``0.0``, not NaN; the front-tracking solver clamps such
1219 out-of-range windows to the last known state rather than masking them.
1220 An output bin with zero throughflow (every overlapping input bin has
1221 zero flow) has an undefined flow-weighted average and is returned as
1222 NaN, matching :func:`infiltration_to_extraction`.
1224 structures : list of dict
1225 List of detailed simulation structures, one for each pore volume, with keys:
1227 - 'waves': List[Wave] - All wave objects created during simulation
1228 - 'events': List[dict] - All events; each record carries ``"theta"`` (m³)
1229 and ``"type"``. Translate to user-facing time t via
1230 ``tracker_state.t_at_theta(event["theta"])`` if needed.
1231 - 'theta_first_arrival': float - Cumulative flow at first nonzero arrival [m³]
1232 - 'n_events': int - Total number of events
1233 - 'n_shocks': int - Number of shocks created
1234 - 'n_rarefactions': int - Number of rarefactions created
1235 - 'n_characteristics': int - Number of characteristics created
1236 - 'theta_current': float - Final simulation cumulative flow [m³]
1237 - 'sorption': SorptionModel - Sorption object
1238 - 'tracker_state': FrontTrackerState - Complete simulation state
1239 - 'aquifer_pore_volume': float - Pore volume for this simulation
1241 Raises
1242 ------
1243 RuntimeError
1244 If the front-tracking solver leaves an unresolved wave interaction (the domain-mass
1245 field over-counts stored mass). This is an internal-consistency tripwire, not an
1246 input limitation: the solver composes every wave interaction, so a firing indicates a bug.
1248 See Also
1249 --------
1250 infiltration_to_extraction : Convolution-based approach for linear retardation
1251 gamma_infiltration_to_extraction : For distributions of pore volumes
1252 :ref:`concept-nonlinear-sorption` : Freundlich isotherm and front-tracking theory
1253 :ref:`assumption-advection-dominated` : When diffusion/dispersion is negligible
1255 Examples
1256 --------
1257 .. disable_try_examples
1259 ::
1261 cout, structures = infiltration_to_extraction_nonlinear_sorption(
1262 cin=cin,
1263 flow=flow,
1264 tedges=tedges,
1265 cout_tedges=cout_tedges,
1266 aquifer_pore_volumes=np.array([500.0]),
1267 freundlich_k=0.01,
1268 freundlich_n=2.0,
1269 bulk_density=1500.0,
1270 porosity=0.3,
1271 )
1273 # Access spin-up period for first pore volume
1274 theta_first = structures[0]["theta_first_arrival"]
1275 t_first = structures[0]["tracker_state"].t_at_theta(theta_first)
1276 print(f"First arrival: θ={theta_first:.2f} m³ (t={t_first:.2f} days)")
1278 # Analyze events for first pore volume
1279 for event in structures[0]["events"]:
1280 print(f"θ={event['theta']:.2f}: {event['type']}")
1281 """
1282 cin, flow, tedges, cout_tedges, aquifer_pore_volumes, sorption, cout_tedges_days = _validate_front_tracking_inputs(
1283 cin=cin,
1284 flow=flow,
1285 tedges=tedges,
1286 cout_tedges=cout_tedges,
1287 aquifer_pore_volumes=aquifer_pore_volumes,
1288 freundlich_k=freundlich_k,
1289 freundlich_n=freundlich_n,
1290 bulk_density=bulk_density,
1291 porosity=porosity,
1292 retardation_factor=retardation_factor,
1293 langmuir_s_max=langmuir_s_max,
1294 langmuir_k_l=langmuir_k_l,
1295 )
1297 # Flow time edges in days (same reference as cout_tedges_days)
1298 flow_tedges_days = tedges_to_days(tedges)
1300 # Each pore-volume bin from the gamma distribution is an equal-mass streamtube,
1301 # so all streamtubes carry equal flow at the outlet. The bundle outlet
1302 # concentration is the simple arithmetic mean over streamtubes. Accumulate the
1303 # per-streamtube output into a running sum so peak memory stays O(n_cout)
1304 # rather than O(n_pv * n_cout).
1305 cout_sum = np.zeros(len(cout_tedges) - 1)
1306 structures = []
1308 # Resolve wave interactions out to the last requested output edge, so beyond-outlet merges
1309 # that still feed an in-window outlet query are composed. The horizon follows the same
1310 # piecewise-linear t → θ map the tracker uses, clamped at the first flow edge and
1311 # extrapolating the final-bin flow when the output window extends past the flow record.
1312 theta_flow_edges = np.concatenate(([0.0], np.cumsum(flow * np.diff(flow_tedges_days))))
1313 t_horizon = max(float(cout_tedges_days[-1]), float(flow_tedges_days[0]))
1314 k = min(int(np.searchsorted(flow_tedges_days, t_horizon, side="right")) - 1, len(flow) - 1)
1315 theta_horizon = float(theta_flow_edges[k] + (t_horizon - flow_tedges_days[k]) * flow[k])
1317 for aquifer_pore_volume in aquifer_pore_volumes:
1318 tracker = FrontTracker(
1319 cin=cin,
1320 flow=flow,
1321 tedges=tedges,
1322 aquifer_pore_volume=aquifer_pore_volume,
1323 sorption=sorption,
1324 theta_horizon=theta_horizon,
1325 )
1327 tracker.run(max_iterations=max_iterations)
1329 # The solver resolves every wave interaction (shock↔shock, fan-entry, doubly-fed
1330 # formation, same-apex annihilation and their compositions), so the wave list is
1331 # interaction-consistent and the single-owner reader is exact. This detector is a
1332 # tripwire: a firing means the solver left an unresolved interaction (an internal bug),
1333 # not an unsupported input — fail loud rather than return a silently wrong cout.
1334 interaction = find_unresolved_interaction(tracker.state)
1335 if interaction is not None:
1336 msg = (
1337 "infiltration_to_extraction_nonlinear_sorption: the front-tracking solver left an "
1338 f"unresolved wave interaction ({interaction}). This is an internal inconsistency; "
1339 "please report it with the cin/flow/pore-volume inputs at "
1340 "https://github.com/gwtransport/gwtransport/issues."
1341 )
1342 raise RuntimeError(msg)
1344 cout_sum += _flow_weighted_front_tracking_output(
1345 cout_tedges_days=cout_tedges_days,
1346 flow_tedges_days=flow_tedges_days,
1347 flow=flow,
1348 v_outlet=aquifer_pore_volume,
1349 waves=tracker.state.waves,
1350 sorption=sorption,
1351 theta_edges=tracker.state.theta_edges,
1352 cin=cin,
1353 )
1355 wave_counts = Counter(type(w) for w in tracker.state.waves)
1356 structure = {
1357 "waves": tracker.state.waves,
1358 "events": tracker.state.events,
1359 "theta_first_arrival": tracker.theta_first_arrival,
1360 "n_events": len(tracker.state.events),
1361 "n_shocks": wave_counts[ShockWave],
1362 "n_rarefactions": wave_counts[RarefactionWave],
1363 "n_characteristics": wave_counts[CharacteristicWave],
1364 "theta_current": tracker.state.theta_current,
1365 "sorption": sorption,
1366 "tracker_state": tracker.state,
1367 "aquifer_pore_volume": aquifer_pore_volume,
1368 }
1369 structures.append(structure)
1371 return cout_sum / len(aquifer_pore_volumes), structures