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