Coverage for src / gwtransport / advection.py: 77%
138 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-04 18:51 +0000
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-04 18:51 +0000
1"""
2Advective Transport Modeling for 1D Aquifer Systems.
4This module provides functions to model compound transport by advection in one-dimensional
5aquifer systems, enabling prediction of solute or temperature concentrations in extracted
6water based on infiltration data and aquifer properties. The model assumes one-dimensional
7groundwater flow where water infiltrates with concentration ``cin``, flows through the aquifer
8with pore volume distribution, compounds are transported with retarded velocity (retardation
9factor >= 1.0), and water is extracted with concentration ``cout``.
11Available functions:
13- :func:`infiltration_to_extraction_series` - Single pore volume, time-shift only. Shifts
14 infiltration time edges forward by residence time. Concentration values remain unchanged
15 (cout = cin). No support for custom output time edges. Use case: Deterministic transport
16 with single flow path.
18- :func:`infiltration_to_extraction` - Arbitrary pore volume distribution, convolution.
19 Supports explicit distribution of aquifer pore volumes with flow-weighted averaging.
20 Flexible output time resolution via cout_tedges. Use case: Known pore volume distribution
21 from streamline analysis.
23- :func:`gamma_infiltration_to_extraction` - Gamma-distributed pore volumes, convolution.
24 Models aquifer heterogeneity with 2-parameter gamma distribution. Parameterizable via
25 (alpha, beta) or (mean, std). Discretizes gamma distribution into equal-probability bins.
26 Use case: Heterogeneous aquifer with calibrated gamma parameters.
28Note on dispersion: The spreading from the pore volume distribution (APVD) represents
29macrodispersion—aquifer-scale velocity heterogeneity that depends on both aquifer
30properties and hydrological boundary conditions. If APVD was calibrated from
31measurements, this spreading already includes microdispersion and molecular diffusion.
32To add microdispersion and molecular diffusion separately (when APVD comes from
33streamline analysis), use :mod:`gwtransport.diffusion`.
34See :ref:`concept-dispersion-scales` for details.
36Note on cross-compound calibration: When APVD is calibrated from measurements of one
37compound (e.g., temperature with D_m ~ 0.1 m²/day) and used to predict another (e.g., a
38solute with D_m ~ 1e-4 m²/day), the molecular diffusion contribution baked into the
39calibrated std may need correction. See :doc:`/examples/05_Diffusion_Dispersion` for
40the correction procedure.
42- :func:`extraction_to_infiltration_series` - Single pore volume, time-shift only
43 (deconvolution). Shifts extraction time edges backward by residence time. Concentration
44 values remain unchanged (cin = cout). Symmetric inverse of infiltration_to_extraction_series.
45 Use case: Backward tracing with single flow path.
47- :func:`extraction_to_infiltration` - Arbitrary pore volume distribution, deconvolution.
48 Inverts forward transport for arbitrary pore volume distributions. Symmetric inverse of
49 infiltration_to_extraction. Flow-weighted averaging in reverse direction. Use case:
50 Estimating infiltration history from extraction data.
52- :func:`gamma_extraction_to_infiltration` - Gamma-distributed pore volumes, deconvolution.
53 Inverts forward transport for gamma-distributed pore volumes. Symmetric inverse of
54 gamma_infiltration_to_extraction. Use case: Calibrating infiltration conditions from
55 extraction measurements.
57- :func:`infiltration_to_extraction_front_tracking` - Exact front tracking with Freundlich sorption.
58 Event-driven algorithm that solves 1D advective transport with Freundlich isotherm using
59 analytical integration of shock and rarefaction waves. Machine-precision physics (no numerical
60 dispersion). Returns bin-averaged concentrations. Use case: Sharp concentration fronts with
61 exact mass balance required, single deterministic flow path.
63- :func:`infiltration_to_extraction_front_tracking_detailed` - Front tracking with piecewise structure.
64 Same as infiltration_to_extraction_front_tracking but also returns complete piecewise analytical
65 structure including all events, segments, and callable analytical forms C(t). Use case: Detailed
66 analysis of shock and rarefaction wave dynamics.
68This file is part of gwtransport which is released under AGPL-3.0 license.
69See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
70"""
72import numpy as np
73import numpy.typing as npt
74import pandas as pd
76from gwtransport import gamma
77from gwtransport.advection_utils import (
78 _extraction_to_infiltration_weights,
79 _infiltration_to_extraction_weights,
80)
81from gwtransport.fronttracking.math import EPSILON_FREUNDLICH_N, ConstantRetardation, FreundlichSorption
82from gwtransport.fronttracking.output import compute_bin_averaged_concentration_exact
83from gwtransport.fronttracking.solver import FrontTracker
84from gwtransport.fronttracking.waves import CharacteristicWave, RarefactionWave, ShockWave
85from gwtransport.residence_time import residence_time
88def infiltration_to_extraction_series(
89 *,
90 flow: npt.ArrayLike,
91 tedges: pd.DatetimeIndex,
92 aquifer_pore_volume: float,
93 retardation_factor: float = 1.0,
94) -> pd.DatetimeIndex:
95 """
96 Compute extraction time edges from infiltration time edges using residence time shifts.
98 This function shifts infiltration time edges forward in time based on residence
99 times computed from flow rates and aquifer properties. The concentration values remain
100 unchanged (cout equals cin), only the time edges are shifted. This assumes a single pore
101 volume (no distribution) and deterministic advective transport.
103 NOTE: This function is specifically designed for single aquifer pore volumes and does not
104 support custom output time edges (cout_tedges). For distributions of aquifer pore volumes
105 or custom output time grids, use `infiltration_to_extraction` instead.
107 Parameters
108 ----------
109 flow : array-like
110 Flow rate values in the aquifer [m3/day].
111 Length must match the number of time bins defined by tedges (len(tedges) - 1).
112 tedges : pandas.DatetimeIndex
113 Time edges defining bins for both cin and flow data. Has length of len(flow) + 1.
114 aquifer_pore_volume : float
115 Single aquifer pore volume [m3] used to compute residence times.
116 retardation_factor : float, optional
117 Retardation factor of the compound in the aquifer (default 1.0).
118 Values > 1.0 indicate slower transport due to sorption/interaction.
120 Returns
121 -------
122 pandas.DatetimeIndex
123 Time edges for the extracted water concentration. Same length as tedges.
124 The concentration values in the extracted water (cout) equal cin, but are
125 aligned with these shifted time edges.
127 Examples
128 --------
129 Basic usage with constant flow:
131 >>> import pandas as pd
132 >>> import numpy as np
133 >>> from gwtransport.utils import compute_time_edges
134 >>> from gwtransport.advection import infiltration_to_extraction_series
135 >>>
136 >>> # Create input data
137 >>> dates = pd.date_range(start="2020-01-01", end="2020-01-10", freq="D")
138 >>> tedges = compute_time_edges(
139 ... tedges=None, tstart=None, tend=dates, number_of_bins=len(dates)
140 ... )
141 >>>
142 >>> # Constant concentration and flow
143 >>> cin = np.ones(len(dates)) * 10.0
144 >>> flow = np.ones(len(dates)) * 100.0 # 100 m3/day
145 >>>
146 >>> # Run infiltration_to_extraction_series with 500 m3 pore volume
147 >>> tedges_out = infiltration_to_extraction_series(
148 ... flow=flow,
149 ... tedges=tedges,
150 ... aquifer_pore_volume=500.0,
151 ... )
152 >>> len(tedges_out)
153 11
154 >>> # Time shift should be residence time = pore_volume / flow = 500 / 100 = 5 days
155 >>> tedges_out[0] - tedges[0]
156 Timedelta('5 days 00:00:00')
158 Plotting the input and output concentrations:
160 >>> import matplotlib.pyplot as plt
161 >>> # Prepare data for step plot (repeat values for visualization)
162 >>> xplot_in = np.repeat(tedges, 2)[1:-1]
163 >>> yplot_in = np.repeat(cin, 2)
164 >>> plt.plot(
165 ... xplot_in, yplot_in, label="Concentration of infiltrated water"
166 ... ) # doctest: +SKIP
167 >>>
168 >>> # cout equals cin, just with shifted time edges
169 >>> xplot_out = np.repeat(tedges_out, 2)[1:-1]
170 >>> yplot_out = np.repeat(cin, 2)
171 >>> plt.plot(
172 ... xplot_out, yplot_out, label="Concentration of extracted water"
173 ... ) # doctest: +SKIP
174 >>> plt.xlabel("Time") # doctest: +SKIP
175 >>> plt.ylabel("Concentration") # doctest: +SKIP
176 >>> plt.legend() # doctest: +SKIP
177 >>> plt.show() # doctest: +SKIP
179 With retardation factor:
181 >>> tedges_out = infiltration_to_extraction_series(
182 ... flow=flow,
183 ... tedges=tedges,
184 ... aquifer_pore_volume=500.0,
185 ... retardation_factor=2.0, # Doubles residence time
186 ... )
187 >>> # Time shift is doubled: 500 * 2.0 / 100 = 10 days
188 >>> tedges_out[0] - tedges[0]
189 Timedelta('10 days 00:00:00')
190 """
191 rt_array = residence_time(
192 flow=flow,
193 flow_tedges=tedges,
194 index=tedges,
195 aquifer_pore_volumes=aquifer_pore_volume,
196 retardation_factor=retardation_factor,
197 direction="infiltration_to_extraction",
198 )
199 return tedges + pd.to_timedelta(rt_array[0], unit="D", errors="coerce")
202def extraction_to_infiltration_series(
203 *,
204 flow: npt.ArrayLike,
205 tedges: pd.DatetimeIndex,
206 aquifer_pore_volume: float,
207 retardation_factor: float = 1.0,
208) -> pd.DatetimeIndex:
209 """
210 Compute infiltration time edges from extraction time edges (deconvolution).
212 This function shifts extraction time edges backward in time based on residence
213 times computed from flow rates and aquifer properties. The concentration values remain
214 unchanged (cin equals cout), only the time edges are shifted. This assumes a single pore
215 volume (no distribution) and deterministic advective transport. This is the inverse
216 operation of infiltration_to_extraction_series.
218 NOTE: This function is specifically designed for single aquifer pore volumes and does not
219 support custom output time edges (cin_tedges). For distributions of aquifer pore volumes
220 or custom output time grids, use `extraction_to_infiltration` instead.
222 Parameters
223 ----------
224 flow : array-like
225 Flow rate values in the aquifer [m3/day].
226 Length must match the number of time bins defined by tedges (len(tedges) - 1).
227 tedges : pandas.DatetimeIndex
228 Time edges defining bins for both cout and flow data. Has length of len(flow) + 1.
229 aquifer_pore_volume : float
230 Single aquifer pore volume [m3] used to compute residence times.
231 retardation_factor : float, optional
232 Retardation factor of the compound in the aquifer (default 1.0).
233 Values > 1.0 indicate slower transport due to sorption/interaction.
235 Returns
236 -------
237 pandas.DatetimeIndex
238 Time edges for the infiltrating water concentration. Same length as tedges.
239 The concentration values in the infiltrating water (cin) equal cout, but are
240 aligned with these shifted time edges.
242 Examples
243 --------
244 Basic usage with constant flow:
246 >>> import pandas as pd
247 >>> import numpy as np
248 >>> from gwtransport.utils import compute_time_edges
249 >>> from gwtransport.advection import extraction_to_infiltration_series
250 >>>
251 >>> # Create input data
252 >>> dates = pd.date_range(start="2020-01-01", end="2020-01-10", freq="D")
253 >>> tedges = compute_time_edges(
254 ... tedges=None, tstart=None, tend=dates, number_of_bins=len(dates)
255 ... )
256 >>>
257 >>> # Constant concentration and flow
258 >>> cout = np.ones(len(dates)) * 10.0
259 >>> flow = np.ones(len(dates)) * 100.0 # 100 m3/day
260 >>>
261 >>> # Run extraction_to_infiltration_series with 500 m3 pore volume
262 >>> tedges_out = extraction_to_infiltration_series(
263 ... flow=flow,
264 ... tedges=tedges,
265 ... aquifer_pore_volume=500.0,
266 ... )
267 >>> len(tedges_out)
268 11
269 >>> # Time shift should be residence time = pore_volume / flow = 500 / 100 = 5 days (backward)
270 >>> # First few elements are NaT due to insufficient history, check a valid index
271 >>> tedges[5] - tedges_out[5]
272 Timedelta('5 days 00:00:00')
274 Plotting the input and output concentrations:
276 >>> import matplotlib.pyplot as plt
277 >>> # Prepare data for step plot (repeat values for visualization)
278 >>> xplot_in = np.repeat(tedges, 2)[1:-1]
279 >>> yplot_in = np.repeat(cout, 2)
280 >>> plt.plot(
281 ... xplot_in, yplot_in, label="Concentration of extracted water"
282 ... ) # doctest: +SKIP
283 >>>
284 >>> # cin equals cout, just with shifted time edges
285 >>> xplot_out = np.repeat(tedges_out, 2)[1:-1]
286 >>> yplot_out = np.repeat(cout, 2)
287 >>> plt.plot(
288 ... xplot_out, yplot_out, label="Concentration of infiltrated water"
289 ... ) # doctest: +SKIP
290 >>> plt.xlabel("Time") # doctest: +SKIP
291 >>> plt.ylabel("Concentration") # doctest: +SKIP
292 >>> plt.legend() # doctest: +SKIP
293 >>> plt.show() # doctest: +SKIP
295 With retardation factor:
297 >>> tedges_out = extraction_to_infiltration_series(
298 ... flow=flow,
299 ... tedges=tedges,
300 ... aquifer_pore_volume=500.0,
301 ... retardation_factor=2.0, # Doubles residence time
302 ... )
303 >>> # Time shift is doubled: 500 * 2.0 / 100 = 10 days (backward)
304 >>> # With longer residence time, more elements are NaT, check the last valid index
305 >>> tedges[10] - tedges_out[10]
306 Timedelta('10 days 00:00:00')
307 """
308 rt_array = residence_time(
309 flow=flow,
310 flow_tedges=tedges,
311 index=tedges,
312 aquifer_pore_volumes=aquifer_pore_volume,
313 retardation_factor=retardation_factor,
314 direction="extraction_to_infiltration",
315 )
316 return tedges - pd.to_timedelta(rt_array[0], unit="D", errors="coerce")
319def gamma_infiltration_to_extraction(
320 *,
321 cin: npt.ArrayLike,
322 flow: npt.ArrayLike,
323 tedges: pd.DatetimeIndex,
324 cout_tedges: pd.DatetimeIndex,
325 alpha: float | None = None,
326 beta: float | None = None,
327 mean: float | None = None,
328 std: float | None = None,
329 n_bins: int = 100,
330 retardation_factor: float = 1.0,
331) -> npt.NDArray[np.floating]:
332 """
333 Compute the concentration of the extracted water by shifting cin with its residence time.
335 The compound is retarded in the aquifer with a retardation factor. The residence
336 time is computed based on the flow rate of the water in the aquifer and the pore volume
337 of the aquifer. The aquifer pore volume is approximated by a gamma distribution, with
338 parameters alpha and beta.
340 This function represents infiltration to extraction modeling (equivalent to convolution).
342 Provide either alpha and beta or mean and std.
344 Parameters
345 ----------
346 cin : array-like
347 Concentration of the compound in infiltrating water or temperature of infiltrating
348 water.
349 flow : array-like
350 Flow rate of water in the aquifer [m3/day].
351 tedges : pandas.DatetimeIndex
352 Time edges for both cin and flow data. Used to compute the cumulative concentration.
353 Has a length of one more than `cin` and `flow`.
354 cout_tedges : pandas.DatetimeIndex
355 Time edges for the output data. Used to compute the cumulative concentration.
356 Has a length of one more than the desired output length.
357 alpha : float, optional
358 Shape parameter of gamma distribution of the aquifer pore volume (must be > 0)
359 beta : float, optional
360 Scale parameter of gamma distribution of the aquifer pore volume (must be > 0)
361 mean : float, optional
362 Mean of the gamma distribution.
363 std : float, optional
364 Standard deviation of the gamma distribution.
365 n_bins : int
366 Number of bins to discretize the gamma distribution.
367 retardation_factor : float
368 Retardation factor of the compound in the aquifer.
370 Returns
371 -------
372 numpy.ndarray
373 Concentration of the compound in the extracted water [ng/m3] or temperature.
375 See Also
376 --------
377 infiltration_to_extraction : Transport with explicit pore volume distribution
378 gamma_extraction_to_infiltration : Reverse operation (deconvolution)
379 gwtransport.gamma.bins : Create gamma distribution bins
380 gwtransport.residence_time.residence_time : Compute residence times
381 gwtransport.diffusion.infiltration_to_extraction : Add microdispersion and molecular diffusion
382 :ref:`concept-gamma-distribution` : Two-parameter pore volume model
383 :ref:`assumption-gamma-distribution` : When gamma distribution is adequate
385 Notes
386 -----
387 The APVD is only time-invariant under the steady-streamlines assumption
388 (see :ref:`assumption-steady-streamlines`).
390 The spreading from the gamma-distributed pore volumes represents macrodispersion
391 (aquifer-scale heterogeneity). When ``std`` comes from calibration on measurements,
392 it absorbs all mixing: macrodispersion, microdispersion, and an average molecular
393 diffusion contribution. When calibrating with the diffusion module, these three
394 components are taken into account separately. When ``std`` comes from streamline
395 analysis, it represents macrodispersion only; microdispersion and molecular diffusion
396 can be added via the diffusion module or by combining variances
397 (see :doc:`/examples/05_Diffusion_Dispersion`).
399 If calibrating with one compound (e.g., temperature) and predicting for another
400 (e.g., a solute), the baked-in molecular diffusion contribution may need
401 correction — see :doc:`/examples/05_Diffusion_Dispersion`.
402 See :ref:`concept-dispersion-scales` for guidance on when to add microdispersion
403 using the diffusion module.
405 Examples
406 --------
407 Basic usage with alpha and beta parameters:
409 >>> import pandas as pd
410 >>> import numpy as np
411 >>> from gwtransport.utils import compute_time_edges
412 >>> from gwtransport.advection import gamma_infiltration_to_extraction
413 >>>
414 >>> # Create input data with aligned time edges
415 >>> dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
416 >>> tedges = compute_time_edges(
417 ... tedges=None, tstart=None, tend=dates, number_of_bins=len(dates)
418 ... )
419 >>>
420 >>> # Create output time edges (can be different alignment)
421 >>> cout_dates = pd.date_range(start="2020-01-05", end="2020-01-15", freq="D")
422 >>> cout_tedges = compute_time_edges(
423 ... tedges=None, tstart=None, tend=cout_dates, number_of_bins=len(cout_dates)
424 ... )
425 >>>
426 >>> # Input concentration and flow (same length, aligned with tedges)
427 >>> cin = pd.Series(np.ones(len(dates)), index=dates)
428 >>> flow = pd.Series(np.ones(len(dates)) * 100, index=dates) # 100 m3/day
429 >>>
430 >>> # Run gamma_infiltration_to_extraction with alpha/beta parameters
431 >>> cout = gamma_infiltration_to_extraction(
432 ... cin=cin,
433 ... flow=flow,
434 ... tedges=tedges,
435 ... cout_tedges=cout_tedges,
436 ... alpha=10.0,
437 ... beta=10.0,
438 ... n_bins=5,
439 ... )
440 >>> cout.shape
441 (11,)
443 Using mean and std parameters instead:
445 >>> cout = gamma_infiltration_to_extraction(
446 ... cin=cin,
447 ... flow=flow,
448 ... tedges=tedges,
449 ... cout_tedges=cout_tedges,
450 ... mean=100.0,
451 ... std=20.0,
452 ... n_bins=5,
453 ... )
455 With retardation factor:
457 >>> cout = gamma_infiltration_to_extraction(
458 ... cin=cin,
459 ... flow=flow,
460 ... tedges=tedges,
461 ... cout_tedges=cout_tedges,
462 ... alpha=10.0,
463 ... beta=10.0,
464 ... retardation_factor=2.0, # Doubles residence time
465 ... )
466 """
467 bins = gamma.bins(alpha=alpha, beta=beta, mean=mean, std=std, n_bins=n_bins)
468 return infiltration_to_extraction(
469 cin=cin,
470 flow=flow,
471 tedges=tedges,
472 cout_tedges=cout_tedges,
473 aquifer_pore_volumes=bins["expected_values"],
474 retardation_factor=retardation_factor,
475 )
478def gamma_extraction_to_infiltration(
479 *,
480 cout: npt.ArrayLike,
481 flow: npt.ArrayLike,
482 tedges: pd.DatetimeIndex,
483 cin_tedges: pd.DatetimeIndex,
484 alpha: float | None = None,
485 beta: float | None = None,
486 mean: float | None = None,
487 std: float | None = None,
488 n_bins: int = 100,
489 retardation_factor: float = 1.0,
490) -> npt.NDArray[np.floating]:
491 """
492 Compute the concentration of the infiltrating water from extracted water (deconvolution).
494 The compound is retarded in the aquifer with a retardation factor. The residence
495 time is computed based on the flow rate of the water in the aquifer and the pore volume
496 of the aquifer. The aquifer pore volume is approximated by a gamma distribution, with
497 parameters alpha and beta.
499 This function represents extraction to infiltration modeling (equivalent to deconvolution).
500 It is symmetric to gamma_infiltration_to_extraction.
502 Provide either alpha and beta or mean and std.
504 Parameters
505 ----------
506 cout : array-like
507 Concentration of the compound in extracted water or temperature of extracted
508 water.
509 tedges : pandas.DatetimeIndex
510 Time edges for the cout and flow data. Used to compute the cumulative concentration.
511 Has a length of one more than `cout` and `flow`.
512 cin_tedges : pandas.DatetimeIndex
513 Time edges for the output (infiltration) data. Used to compute the cumulative concentration.
514 Has a length of one more than the desired output length.
515 flow : array-like
516 Flow rate of water in the aquifer [m3/day].
517 alpha : float, optional
518 Shape parameter of gamma distribution of the aquifer pore volume (must be > 0)
519 beta : float, optional
520 Scale parameter of gamma distribution of the aquifer pore volume (must be > 0)
521 mean : float, optional
522 Mean of the gamma distribution.
523 std : float, optional
524 Standard deviation of the gamma distribution.
525 n_bins : int
526 Number of bins to discretize the gamma distribution.
527 retardation_factor : float, optional
528 Retardation factor of the compound in the aquifer (default 1.0).
529 Values > 1.0 indicate slower transport due to sorption/interaction.
531 Returns
532 -------
533 numpy.ndarray
534 Concentration of the compound in the infiltrating water [ng/m3] or temperature.
536 See Also
537 --------
538 extraction_to_infiltration : Deconvolution with explicit pore volume distribution
539 gamma_infiltration_to_extraction : Forward operation (convolution)
540 gwtransport.gamma.bins : Create gamma distribution bins
541 gwtransport.diffusion.extraction_to_infiltration : Deconvolution with microdispersion and molecular diffusion
542 :ref:`concept-gamma-distribution` : Two-parameter pore volume model
543 :ref:`assumption-gamma-distribution` : When gamma distribution is adequate
545 Notes
546 -----
547 The APVD is only time-invariant under the steady-streamlines assumption
548 (see :ref:`assumption-steady-streamlines`).
550 The spreading from the gamma-distributed pore volumes represents macrodispersion
551 (aquifer-scale heterogeneity). When ``std`` comes from calibration on measurements,
552 it absorbs all mixing: macrodispersion, microdispersion, and an average molecular
553 diffusion contribution. When calibrating with the diffusion module, these three
554 components are taken into account separately. When ``std`` comes from streamline
555 analysis, it represents macrodispersion only; microdispersion and molecular diffusion
556 can be added via the diffusion module or by combining variances
557 (see :doc:`/examples/05_Diffusion_Dispersion`).
559 If calibrating with one compound (e.g., temperature) and predicting for another
560 (e.g., a solute), the baked-in molecular diffusion contribution may need
561 correction — see :doc:`/examples/05_Diffusion_Dispersion`.
562 See :ref:`concept-dispersion-scales` for guidance on when to add microdispersion
563 using the diffusion module.
565 Examples
566 --------
567 Basic usage with alpha and beta parameters:
569 >>> import pandas as pd
570 >>> import numpy as np
571 >>> from gwtransport.utils import compute_time_edges
572 >>> from gwtransport.advection import gamma_extraction_to_infiltration
573 >>>
574 >>> # Create input data with aligned time edges
575 >>> dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
576 >>> tedges = compute_time_edges(
577 ... tedges=None, tstart=None, tend=dates, number_of_bins=len(dates)
578 ... )
579 >>>
580 >>> # Create output time edges (can be different alignment)
581 >>> cin_dates = pd.date_range(start="2019-12-25", end="2020-01-15", freq="D")
582 >>> cin_tedges = compute_time_edges(
583 ... tedges=None, tstart=None, tend=cin_dates, number_of_bins=len(cin_dates)
584 ... )
585 >>>
586 >>> # Input concentration and flow (same length, aligned with tedges)
587 >>> cout = pd.Series(np.ones(len(dates)), index=dates)
588 >>> flow = pd.Series(np.ones(len(dates)) * 100, index=dates) # 100 m3/day
589 >>>
590 >>> # Run gamma_extraction_to_infiltration with alpha/beta parameters
591 >>> cin = gamma_extraction_to_infiltration(
592 ... cout=cout,
593 ... tedges=tedges,
594 ... cin_tedges=cin_tedges,
595 ... flow=flow,
596 ... alpha=10.0,
597 ... beta=10.0,
598 ... n_bins=5,
599 ... )
600 >>> cin.shape
601 (22,)
603 Using mean and std parameters instead:
605 >>> cin = gamma_extraction_to_infiltration(
606 ... cout=cout,
607 ... tedges=tedges,
608 ... cin_tedges=cin_tedges,
609 ... flow=flow,
610 ... mean=100.0,
611 ... std=20.0,
612 ... n_bins=5,
613 ... )
615 With retardation factor:
617 >>> cin = gamma_extraction_to_infiltration(
618 ... cout=cout,
619 ... tedges=tedges,
620 ... cin_tedges=cin_tedges,
621 ... flow=flow,
622 ... alpha=10.0,
623 ... beta=10.0,
624 ... retardation_factor=2.0, # Doubles residence time
625 ... )
626 """
627 bins = gamma.bins(alpha=alpha, beta=beta, mean=mean, std=std, n_bins=n_bins)
628 return extraction_to_infiltration(
629 cout=cout,
630 flow=flow,
631 tedges=tedges,
632 cin_tedges=cin_tedges,
633 aquifer_pore_volumes=bins["expected_values"],
634 retardation_factor=retardation_factor,
635 )
638def infiltration_to_extraction(
639 *,
640 cin: npt.ArrayLike,
641 flow: npt.ArrayLike,
642 tedges: pd.DatetimeIndex,
643 cout_tedges: pd.DatetimeIndex,
644 aquifer_pore_volumes: npt.ArrayLike,
645 retardation_factor: float = 1.0,
646) -> npt.NDArray[np.floating]:
647 """
648 Compute the concentration of the extracted water using flow-weighted advection.
650 This function implements an infiltration to extraction advection model where cin and flow values
651 correspond to the same aligned time bins defined by tedges.
653 The algorithm:
654 1. Computes residence times for each pore volume at cout time edges
655 2. Calculates infiltration time edges by subtracting residence times
656 3. Determines temporal overlaps between infiltration and cin time windows
657 4. Creates flow-weighted overlap matrices normalized by total weights
658 5. Computes weighted contributions and averages across pore volumes
661 Parameters
662 ----------
663 cin : array-like
664 Concentration values of infiltrating water or temperature [concentration units].
665 Length must match the number of time bins defined by tedges.
666 flow : array-like
667 Flow rate values in the aquifer [m3/day].
668 Length must match cin and the number of time bins defined by tedges.
669 tedges : pandas.DatetimeIndex
670 Time edges defining bins for both cin and flow data. Has length of
671 len(cin) + 1 and len(flow) + 1.
672 cout_tedges : pandas.DatetimeIndex
673 Time edges for output data bins. Has length of desired output + 1.
674 Can have different time alignment and resolution than tedges.
675 aquifer_pore_volumes : array-like
676 Array of aquifer pore volumes [m3] representing the distribution
677 of residence times in the aquifer system.
678 retardation_factor : float, optional
679 Retardation factor of the compound in the aquifer (default 1.0).
680 Values > 1.0 indicate slower transport due to sorption/interaction.
682 Returns
683 -------
684 numpy.ndarray
685 Flow-weighted concentration in the extracted water. Same units as cin.
686 Length equals len(cout_tedges) - 1. NaN values indicate time periods
687 with no valid contributions from the infiltration data.
689 Raises
690 ------
691 ValueError
692 If tedges length doesn't match cin/flow arrays plus one, or if
693 infiltration time edges become non-monotonic (invalid input conditions).
695 See Also
696 --------
697 gamma_infiltration_to_extraction : Transport with gamma-distributed pore volumes
698 extraction_to_infiltration : Reverse operation (deconvolution)
699 infiltration_to_extraction_series : Simple time-shift for single pore volume
700 gwtransport.residence_time.residence_time : Compute residence times from flow and pore volume
701 gwtransport.residence_time.freundlich_retardation : Compute concentration-dependent retardation
702 :ref:`concept-pore-volume-distribution` : Background on aquifer heterogeneity modeling
703 :ref:`concept-transport-equation` : Flow-weighted averaging approach
705 Examples
706 --------
707 Basic usage with pandas Series:
709 >>> import pandas as pd
710 >>> import numpy as np
711 >>> from gwtransport.utils import compute_time_edges
712 >>> from gwtransport.advection import infiltration_to_extraction
713 >>>
714 >>> # Create input data
715 >>> dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
716 >>> tedges = compute_time_edges(
717 ... tedges=None, tstart=None, tend=dates, number_of_bins=len(dates)
718 ... )
719 >>>
720 >>> # Create output time edges (different alignment)
721 >>> cout_dates = pd.date_range(start="2020-01-05", end="2020-01-15", freq="D")
722 >>> cout_tedges = compute_time_edges(
723 ... tedges=None, tstart=None, tend=cout_dates, number_of_bins=len(cout_dates)
724 ... )
725 >>>
726 >>> # Input concentration and flow
727 >>> cin = pd.Series(np.ones(len(dates)), index=dates)
728 >>> flow = pd.Series(np.ones(len(dates)) * 100, index=dates) # 100 m3/day
729 >>>
730 >>> # Define distribution of aquifer pore volumes
731 >>> aquifer_pore_volumes = np.array([50, 100, 200]) # m3
732 >>>
733 >>> # Run infiltration_to_extraction
734 >>> cout = infiltration_to_extraction(
735 ... cin=cin,
736 ... flow=flow,
737 ... tedges=tedges,
738 ... cout_tedges=cout_tedges,
739 ... aquifer_pore_volumes=aquifer_pore_volumes,
740 ... )
741 >>> cout.shape
742 (11,)
744 Using array inputs instead of pandas Series:
746 >>> # Convert to arrays
747 >>> cin_values = cin.values
748 >>> flow_values = flow.values
749 >>>
750 >>> cout = infiltration_to_extraction(
751 ... cin=cin_values,
752 ... flow=flow,
753 ... tedges=tedges,
754 ... cout_tedges=cout_tedges,
755 ... aquifer_pore_volumes=aquifer_pore_volumes,
756 ... )
758 With constant retardation factor (linear sorption):
760 >>> cout = infiltration_to_extraction(
761 ... cin=cin,
762 ... flow=flow,
763 ... tedges=tedges,
764 ... cout_tedges=cout_tedges,
765 ... aquifer_pore_volumes=aquifer_pore_volumes,
766 ... retardation_factor=2.0, # Compound moves twice as slowly
767 ... )
769 Note: For concentration-dependent retardation (nonlinear sorption),
770 use `infiltration_to_extraction_front_tracking_detailed` instead, as this
771 function only supports constant (float) retardation factors.
773 Using single pore volume:
775 >>> single_volume = np.array([100]) # Single 100 m3 pore volume
776 >>> cout = infiltration_to_extraction(
777 ... cin=cin,
778 ... flow=flow,
779 ... tedges=tedges,
780 ... cout_tedges=cout_tedges,
781 ... aquifer_pore_volumes=single_volume,
782 ... )
783 """
784 tedges = pd.DatetimeIndex(tedges)
785 cout_tedges = pd.DatetimeIndex(cout_tedges)
787 # Convert to arrays for vectorized operations
788 cin = np.asarray(cin)
789 flow = np.asarray(flow)
790 aquifer_pore_volumes = np.asarray(aquifer_pore_volumes)
792 if len(tedges) != len(cin) + 1:
793 msg = "tedges must have one more element than cin"
794 raise ValueError(msg)
795 if len(tedges) != len(flow) + 1:
796 msg = "tedges must have one more element than flow"
797 raise ValueError(msg)
799 # Validate inputs do not contain NaN values
800 if np.any(np.isnan(cin)):
801 msg = "cin contains NaN values, which are not allowed"
802 raise ValueError(msg)
803 if np.any(np.isnan(flow)):
804 msg = "flow contains NaN values, which are not allowed"
805 raise ValueError(msg)
806 if retardation_factor < 1.0:
807 msg = "retardation_factor must be >= 1.0"
808 raise ValueError(msg)
810 # Compute normalized weights (includes all pre-computation)
811 normalized_weights = _infiltration_to_extraction_weights(
812 tedges=tedges,
813 cout_tedges=cout_tedges,
814 aquifer_pore_volumes=aquifer_pore_volumes,
815 flow=flow,
816 retardation_factor=retardation_factor,
817 )
819 # Apply to concentrations and handle NaN for periods with no contributions
820 out = normalized_weights.dot(cin)
821 # Set NaN where no valid pore volumes contributed
822 total_weights = np.sum(normalized_weights, axis=1)
823 out[total_weights == 0] = np.nan
825 return out
828def extraction_to_infiltration(
829 *,
830 cout: npt.ArrayLike,
831 flow: npt.ArrayLike,
832 tedges: pd.DatetimeIndex,
833 cin_tedges: pd.DatetimeIndex,
834 aquifer_pore_volumes: npt.ArrayLike,
835 retardation_factor: float = 1.0,
836) -> npt.NDArray[np.floating]:
837 """
838 Compute the concentration of the infiltrating water from extracted water (deconvolution).
840 This function implements an extraction to infiltration advection model (inverse of infiltration_to_extraction)
841 where cout and flow values correspond to the same aligned time bins defined by tedges.
843 SYMMETRIC RELATIONSHIP:
844 - infiltration_to_extraction: cin + tedges → cout + cout_tedges
845 - extraction_to_infiltration: cout + tedges → cin + cin_tedges
847 The algorithm (symmetric to infiltration_to_extraction):
848 1. Computes residence times for each pore volume at cint time edges
849 2. Calculates extraction time edges by adding residence times (reverse of infiltration_to_extraction)
850 3. Determines temporal overlaps between extraction and cout time windows
851 4. Creates flow-weighted overlap matrices normalized by total weights
852 5. Computes weighted contributions and averages across pore volumes
855 Parameters
856 ----------
857 cout : array-like
858 Concentration values of extracted water [concentration units].
859 Length must match the number of time bins defined by tedges.
860 flow : array-like
861 Flow rate values in the aquifer [m3/day].
862 Length must match cout and the number of time bins defined by tedges.
863 tedges : pandas.DatetimeIndex
864 Time edges defining bins for both cout and flow data. Has length of
865 len(cout) + 1 and len(flow) + 1.
866 cin_tedges : pandas.DatetimeIndex
867 Time edges for output (infiltration) data bins. Has length of desired output + 1.
868 Can have different time alignment and resolution than tedges.
869 aquifer_pore_volumes : array-like
870 Array of aquifer pore volumes [m3] representing the distribution
871 of residence times in the aquifer system.
872 retardation_factor : float, optional
873 Retardation factor of the compound in the aquifer (default 1.0).
874 Values > 1.0 indicate slower transport due to sorption/interaction.
876 Returns
877 -------
878 numpy.ndarray
879 Flow-weighted concentration in the infiltrating water. Same units as cout.
880 Length equals len(cin_tedges) - 1. NaN values indicate time periods
881 with no valid contributions from the extraction data.
883 Raises
884 ------
885 ValueError
886 If tedges length doesn't match cout/flow arrays plus one, or if
887 extraction time edges become non-monotonic (invalid input conditions).
889 See Also
890 --------
891 gamma_extraction_to_infiltration : Deconvolution with gamma-distributed pore volumes
892 infiltration_to_extraction : Forward operation (convolution)
893 extraction_to_infiltration_series : Simple time-shift for single pore volume
894 gwtransport.residence_time.residence_time : Compute residence times from flow and pore volume
895 gwtransport.residence_time.freundlich_retardation : Compute concentration-dependent retardation
896 :ref:`concept-pore-volume-distribution` : Background on aquifer heterogeneity modeling
897 :ref:`concept-transport-equation` : Flow-weighted averaging approach
899 Examples
900 --------
901 Basic usage with pandas Series:
903 >>> import pandas as pd
904 >>> import numpy as np
905 >>> from gwtransport.utils import compute_time_edges
906 >>> from gwtransport.advection import extraction_to_infiltration
907 >>>
908 >>> # Create input data
909 >>> dates = pd.date_range(start="2020-01-01", end="2020-01-20", freq="D")
910 >>> tedges = compute_time_edges(
911 ... tedges=None, tstart=None, tend=dates, number_of_bins=len(dates)
912 ... )
913 >>>
914 >>> # Create output time edges (different alignment)
915 >>> cint_dates = pd.date_range(start="2019-12-25", end="2020-01-15", freq="D")
916 >>> cin_tedges = compute_time_edges(
917 ... tedges=None, tstart=None, tend=cint_dates, number_of_bins=len(cint_dates)
918 ... )
919 >>>
920 >>> # Input concentration and flow
921 >>> cout = pd.Series(np.ones(len(dates)), index=dates)
922 >>> flow = pd.Series(np.ones(len(dates)) * 100, index=dates) # 100 m3/day
923 >>>
924 >>> # Define distribution of aquifer pore volumes
925 >>> aquifer_pore_volumes = np.array([50, 100, 200]) # m3
926 >>>
927 >>> # Run extraction_to_infiltration
928 >>> cin = extraction_to_infiltration(
929 ... cout=cout,
930 ... flow=flow,
931 ... tedges=tedges,
932 ... cin_tedges=cin_tedges,
933 ... aquifer_pore_volumes=aquifer_pore_volumes,
934 ... )
935 >>> cin.shape
936 (22,)
938 Using array inputs instead of pandas Series:
940 >>> # Convert to arrays
941 >>> cout = cout.values
942 >>> flow = flow.values
943 >>>
944 >>> cin = extraction_to_infiltration(
945 ... cout=cout,
946 ... flow=flow,
947 ... tedges=tedges,
948 ... cin_tedges=cin_tedges,
949 ... aquifer_pore_volumes=aquifer_pore_volumes,
950 ... )
952 With retardation factor:
954 >>> cin = extraction_to_infiltration(
955 ... cout=cout,
956 ... flow=flow,
957 ... tedges=tedges,
958 ... cin_tedges=cin_tedges,
959 ... aquifer_pore_volumes=aquifer_pore_volumes,
960 ... retardation_factor=2.0, # Compound moves twice as slowly
961 ... )
963 Using single pore volume:
965 >>> single_volume = np.array([100]) # Single 100 m3 pore volume
966 >>> cin = extraction_to_infiltration(
967 ... cout=cout,
968 ... flow=flow,
969 ... tedges=tedges,
970 ... cin_tedges=cin_tedges,
971 ... aquifer_pore_volumes=single_volume,
972 ... )
974 """
975 tedges = pd.DatetimeIndex(tedges)
976 cin_tedges = pd.DatetimeIndex(cin_tedges)
978 # Convert to arrays for vectorized operations
979 cout = np.asarray(cout)
980 flow = np.asarray(flow)
982 if len(tedges) != len(cout) + 1:
983 msg = "tedges must have one more element than cout"
984 raise ValueError(msg)
985 if len(tedges) != len(flow) + 1:
986 msg = "tedges must have one more element than flow"
987 raise ValueError(msg)
989 # Validate inputs do not contain NaN values
990 if np.any(np.isnan(cout)):
991 msg = "cout contains NaN values, which are not allowed"
992 raise ValueError(msg)
993 if np.any(np.isnan(flow)):
994 msg = "flow contains NaN values, which are not allowed"
995 raise ValueError(msg)
996 if retardation_factor < 1.0:
997 msg = "retardation_factor must be >= 1.0"
998 raise ValueError(msg)
1000 aquifer_pore_volumes = np.asarray(aquifer_pore_volumes)
1002 # Compute normalized weights (includes all pre-computation)
1003 normalized_weights = _extraction_to_infiltration_weights(
1004 tedges=tedges,
1005 cin_tedges=cin_tedges,
1006 aquifer_pore_volumes=aquifer_pore_volumes,
1007 flow=flow,
1008 retardation_factor=retardation_factor,
1009 )
1011 # Apply to concentrations and handle NaN for periods with no contributions
1012 out = normalized_weights.dot(cout)
1013 # Set NaN where no valid pore volumes contributed
1014 total_weights = np.sum(normalized_weights, axis=1)
1015 out[total_weights == 0] = np.nan
1017 return out
1020def _validate_front_tracking_inputs(
1021 *,
1022 cin: npt.ArrayLike,
1023 flow: npt.ArrayLike,
1024 tedges: pd.DatetimeIndex,
1025 cout_tedges: pd.DatetimeIndex,
1026 aquifer_pore_volumes: npt.ArrayLike,
1027 freundlich_k: float | None,
1028 freundlich_n: float | None,
1029 bulk_density: float | None,
1030 porosity: float | None,
1031 retardation_factor: float | None,
1032) -> tuple:
1033 """Validate inputs and create sorption object for front tracking functions."""
1034 cin = np.asarray(cin, dtype=float)
1035 flow = np.asarray(flow, dtype=float)
1036 tedges = pd.DatetimeIndex(tedges)
1037 cout_tedges = pd.DatetimeIndex(cout_tedges)
1038 aquifer_pore_volumes = np.asarray(aquifer_pore_volumes, dtype=float)
1040 if len(tedges) != len(cin) + 1:
1041 msg = "tedges must have length len(cin) + 1"
1042 raise ValueError(msg)
1043 if len(flow) != len(cin):
1044 msg = "flow must have same length as cin"
1045 raise ValueError(msg)
1046 if np.any(cin < 0):
1047 msg = "cin must be non-negative"
1048 raise ValueError(msg)
1049 if np.any(flow <= 0):
1050 msg = "flow must be positive"
1051 raise ValueError(msg)
1052 if np.any(np.isnan(cin)) or np.any(np.isnan(flow)):
1053 msg = "cin and flow must not contain NaN"
1054 raise ValueError(msg)
1055 if np.any(aquifer_pore_volumes <= 0):
1056 msg = "aquifer_pore_volumes must be positive"
1057 raise ValueError(msg)
1059 # Convert cout_tedges to days (relative to tedges[0]) for output computation
1060 t_ref = tedges[0]
1061 cout_tedges_days = ((cout_tedges - t_ref) / pd.Timedelta(days=1)).values
1063 # Create sorption object
1064 if retardation_factor is not None:
1065 if retardation_factor < 1.0:
1066 msg = "retardation_factor must be >= 1.0"
1067 raise ValueError(msg)
1069 sorption = ConstantRetardation(retardation_factor=retardation_factor)
1070 else:
1071 if freundlich_k is None or freundlich_n is None or bulk_density is None or porosity is None:
1072 msg = (
1073 "Must provide either retardation_factor or all Freundlich parameters "
1074 "(freundlich_k, freundlich_n, bulk_density, porosity)"
1075 )
1076 raise ValueError(msg)
1077 if freundlich_k <= 0 or freundlich_n <= 0:
1078 msg = "Freundlich parameters must be positive"
1079 raise ValueError(msg)
1080 if abs(freundlich_n - 1.0) < EPSILON_FREUNDLICH_N:
1081 msg = "freundlich_n = 1 not supported (use retardation_factor for linear case)"
1082 raise ValueError(msg)
1083 if bulk_density <= 0 or not 0 < porosity < 1:
1084 msg = "Invalid physical parameters"
1085 raise ValueError(msg)
1087 sorption = FreundlichSorption(
1088 k_f=freundlich_k,
1089 n=freundlich_n,
1090 bulk_density=bulk_density,
1091 porosity=porosity,
1092 )
1094 return cin, flow, tedges, cout_tedges, aquifer_pore_volumes, sorption, cout_tedges_days
1097def infiltration_to_extraction_front_tracking(
1098 *,
1099 cin: npt.ArrayLike,
1100 flow: npt.ArrayLike,
1101 tedges: pd.DatetimeIndex,
1102 cout_tedges: pd.DatetimeIndex,
1103 aquifer_pore_volumes: npt.ArrayLike,
1104 freundlich_k: float | None = None,
1105 freundlich_n: float | None = None,
1106 bulk_density: float | None = None,
1107 porosity: float | None = None,
1108 retardation_factor: float | None = None,
1109 max_iterations: int = 10000,
1110) -> npt.NDArray[np.floating]:
1111 """
1112 Compute extracted concentration using exact front tracking with nonlinear sorption.
1114 Uses event-driven analytical algorithm that tracks shock waves, rarefaction waves,
1115 and characteristics with machine precision. No numerical dispersion, exact mass
1116 balance to floating-point precision.
1118 Parameters
1119 ----------
1120 cin : array-like
1121 Infiltration concentration [mg/L or any units].
1122 Length = len(tedges) - 1.
1123 flow : array-like
1124 Flow rate [m³/day]. Must be positive.
1125 Length = len(tedges) - 1.
1126 tedges : pandas.DatetimeIndex
1127 Time bin edges. Length = len(cin) + 1.
1128 cout_tedges : pandas.DatetimeIndex
1129 Output time bin edges. Can be different from tedges.
1130 Length determines output array size.
1131 aquifer_pore_volumes : array-like
1132 Array of aquifer pore volumes [m³] representing the distribution
1133 of residence times in the aquifer system. Each pore volume must be positive.
1134 freundlich_k : float, optional
1135 Freundlich coefficient [(m³/kg)^(1/n)]. Must be positive.
1136 Used if retardation_factor is None.
1137 freundlich_n : float, optional
1138 Freundlich exponent [-]. Must be positive and != 1.
1139 Used if retardation_factor is None.
1140 bulk_density : float, optional
1141 Bulk density [kg/m³]. Must be positive.
1142 Used if retardation_factor is None.
1143 porosity : float, optional
1144 Porosity [-]. Must be in (0, 1).
1145 Used if retardation_factor is None.
1146 retardation_factor : float, optional
1147 Constant retardation factor [-]. If provided, uses linear retardation
1148 instead of Freundlich sorption. Must be >= 1.0.
1149 max_iterations : int, optional
1150 Maximum number of events. Default 10000.
1152 Returns
1153 -------
1154 cout : numpy.ndarray
1155 Bin-averaged extraction concentration averaged across all pore volumes.
1156 Length = len(cout_tedges) - 1.
1158 Notes
1159 -----
1160 **Spin-up Period**:
1161 The function computes the first arrival time t_first. Concentrations
1162 before t_first are affected by unknown initial conditions and should
1163 not be used for analysis. Use `infiltration_to_extraction_front_tracking_detailed`
1164 to access t_first.
1166 **Machine Precision**:
1167 All calculations use exact analytical formulas. Mass balance is conserved
1168 to floating-point precision (~1e-14 relative error). No numerical tolerances
1169 are used for time/position calculations.
1171 **Physical Correctness**:
1172 - All shocks satisfy Lax entropy condition
1173 - Rarefaction waves use self-similar solutions
1174 - Causality is strictly enforced
1176 Examples
1177 --------
1178 >>> import numpy as np
1179 >>> import pandas as pd
1180 >>>
1181 >>> # Pulse injection with single pore volume
1182 >>> tedges = pd.date_range("2020-01-01", periods=4, freq="10D")
1183 >>> cin = np.array([0.0, 10.0, 0.0])
1184 >>> flow = np.array([100.0, 100.0, 100.0])
1185 >>> cout_tedges = pd.date_range("2020-01-01", periods=10, freq="5D")
1186 >>>
1187 >>> cout = infiltration_to_extraction_front_tracking(
1188 ... cin=cin,
1189 ... flow=flow,
1190 ... tedges=tedges,
1191 ... cout_tedges=cout_tedges,
1192 ... aquifer_pore_volumes=np.array([500.0]),
1193 ... freundlich_k=0.01,
1194 ... freundlich_n=2.0,
1195 ... bulk_density=1500.0,
1196 ... porosity=0.3,
1197 ... )
1199 With multiple pore volumes (distribution):
1201 >>> aquifer_pore_volumes = np.array([400.0, 500.0, 600.0])
1202 >>> cout = infiltration_to_extraction_front_tracking(
1203 ... cin=cin,
1204 ... flow=flow,
1205 ... tedges=tedges,
1206 ... cout_tedges=cout_tedges,
1207 ... aquifer_pore_volumes=aquifer_pore_volumes,
1208 ... freundlich_k=0.01,
1209 ... freundlich_n=2.0,
1210 ... bulk_density=1500.0,
1211 ... porosity=0.3,
1212 ... )
1214 See Also
1215 --------
1216 infiltration_to_extraction_front_tracking_detailed : Returns detailed structure
1217 infiltration_to_extraction : Convolution-based approach for linear case
1218 gamma_infiltration_to_extraction : For distributions of pore volumes
1219 :ref:`concept-nonlinear-sorption` : Freundlich isotherm and front-tracking theory
1220 :ref:`assumption-advection-dominated` : When diffusion/dispersion is negligible
1221 """
1222 cin, flow, tedges, cout_tedges, aquifer_pore_volumes, sorption, cout_tedges_days = _validate_front_tracking_inputs(
1223 cin=cin,
1224 flow=flow,
1225 tedges=tedges,
1226 cout_tedges=cout_tedges,
1227 aquifer_pore_volumes=aquifer_pore_volumes,
1228 freundlich_k=freundlich_k,
1229 freundlich_n=freundlich_n,
1230 bulk_density=bulk_density,
1231 porosity=porosity,
1232 retardation_factor=retardation_factor,
1233 )
1235 # Loop over each pore volume and compute concentration
1236 cout_all = np.zeros((len(aquifer_pore_volumes), len(cout_tedges) - 1))
1238 for i, aquifer_pore_volume in enumerate(aquifer_pore_volumes):
1239 # Create tracker and run simulation for this pore volume
1240 tracker = FrontTracker(
1241 cin=cin,
1242 flow=flow,
1243 tedges=tedges,
1244 aquifer_pore_volume=aquifer_pore_volume,
1245 sorption=sorption,
1246 )
1248 tracker.run(max_iterations=max_iterations)
1250 # Extract bin-averaged concentrations at outlet for this pore volume
1251 cout_all[i, :] = compute_bin_averaged_concentration_exact(
1252 t_edges=cout_tedges_days,
1253 v_outlet=aquifer_pore_volume,
1254 waves=tracker.state.waves,
1255 sorption=sorption,
1256 )
1258 # Return average across all pore volumes
1259 return np.mean(cout_all, axis=0)
1262def infiltration_to_extraction_front_tracking_detailed(
1263 *,
1264 cin: npt.ArrayLike,
1265 flow: npt.ArrayLike,
1266 tedges: pd.DatetimeIndex,
1267 cout_tedges: pd.DatetimeIndex,
1268 aquifer_pore_volumes: npt.ArrayLike,
1269 freundlich_k: float | None = None,
1270 freundlich_n: float | None = None,
1271 bulk_density: float | None = None,
1272 porosity: float | None = None,
1273 retardation_factor: float | None = None,
1274 max_iterations: int = 10000,
1275) -> tuple[npt.NDArray[np.floating], list[dict]]:
1276 """
1277 Compute extracted concentration with complete diagnostic information.
1279 Returns both bin-averaged concentrations and detailed simulation structure for each pore volume.
1281 Parameters
1282 ----------
1283 cin : array-like
1284 Infiltration concentration [mg/L or any units].
1285 Length = len(tedges) - 1.
1286 flow : array-like
1287 Flow rate [m³/day]. Must be positive.
1288 Length = len(tedges) - 1.
1289 tedges : pandas.DatetimeIndex
1290 Time bin edges. Length = len(cin) + 1.
1291 cout_tedges : pandas.DatetimeIndex
1292 Output time bin edges. Can be different from tedges.
1293 Length determines output array size.
1294 aquifer_pore_volumes : array-like
1295 Array of aquifer pore volumes [m³] representing the distribution
1296 of residence times in the aquifer system. Each pore volume must be positive.
1297 freundlich_k : float, optional
1298 Freundlich coefficient [(m³/kg)^(1/n)]. Must be positive.
1299 Used if retardation_factor is None.
1300 freundlich_n : float, optional
1301 Freundlich exponent [-]. Must be positive and != 1.
1302 Used if retardation_factor is None.
1303 bulk_density : float, optional
1304 Bulk density [kg/m³]. Must be positive.
1305 Used if retardation_factor is None.
1306 porosity : float, optional
1307 Porosity [-]. Must be in (0, 1).
1308 Used if retardation_factor is None.
1309 retardation_factor : float, optional
1310 Constant retardation factor [-]. If provided, uses linear retardation
1311 instead of Freundlich sorption. Must be >= 1.0.
1312 max_iterations : int, optional
1313 Maximum number of events. Default 10000.
1315 Returns
1316 -------
1317 cout : numpy.ndarray
1318 Bin-averaged concentrations averaged across all pore volumes.
1320 structures : list of dict
1321 List of detailed simulation structures, one for each pore volume, with keys:
1323 - 'waves': List[Wave] - All wave objects created during simulation
1324 - 'events': List[dict] - All events with times, types, and details
1325 - 't_first_arrival': float - First arrival time (end of spin-up period)
1326 - 'n_events': int - Total number of events
1327 - 'n_shocks': int - Number of shocks created
1328 - 'n_rarefactions': int - Number of rarefactions created
1329 - 'n_characteristics': int - Number of characteristics created
1330 - 'final_time': float - Final simulation time
1331 - 'sorption': FreundlichSorption | ConstantRetardation - Sorption object
1332 - 'tracker_state': FrontTrackerState - Complete simulation state
1333 - 'aquifer_pore_volume': float - Pore volume for this simulation
1335 Examples
1336 --------
1337 ::
1339 cout, structures = infiltration_to_extraction_front_tracking_detailed(
1340 cin=cin,
1341 flow=flow,
1342 tedges=tedges,
1343 cout_tedges=cout_tedges,
1344 aquifer_pore_volumes=np.array([500.0]),
1345 freundlich_k=0.01,
1346 freundlich_n=2.0,
1347 bulk_density=1500.0,
1348 porosity=0.3,
1349 )
1351 # Access spin-up period for first pore volume
1352 print(f"First arrival: {structures[0]['t_first_arrival']:.2f} days")
1354 # Analyze events for first pore volume
1355 for event in structures[0]["events"]:
1356 print(f"t={event['time']:.2f}: {event['type']}")
1358 See Also
1359 --------
1360 infiltration_to_extraction_front_tracking : Returns concentrations only (simpler interface)
1361 :ref:`concept-nonlinear-sorption` : Freundlich isotherm and front-tracking theory
1362 :ref:`assumption-advection-dominated` : When diffusion/dispersion is negligible
1363 """
1364 cin, flow, tedges, cout_tedges, aquifer_pore_volumes, sorption, cout_tedges_days = _validate_front_tracking_inputs(
1365 cin=cin,
1366 flow=flow,
1367 tedges=tedges,
1368 cout_tedges=cout_tedges,
1369 aquifer_pore_volumes=aquifer_pore_volumes,
1370 freundlich_k=freundlich_k,
1371 freundlich_n=freundlich_n,
1372 bulk_density=bulk_density,
1373 porosity=porosity,
1374 retardation_factor=retardation_factor,
1375 )
1377 # Loop over each pore volume and compute concentration
1378 cout_all = np.zeros((len(aquifer_pore_volumes), len(cout_tedges) - 1))
1379 structures = []
1381 for i, aquifer_pore_volume in enumerate(aquifer_pore_volumes):
1382 # Create tracker and run simulation for this pore volume
1383 tracker = FrontTracker(
1384 cin=cin,
1385 flow=flow,
1386 tedges=tedges,
1387 aquifer_pore_volume=aquifer_pore_volume,
1388 sorption=sorption,
1389 )
1391 tracker.run(max_iterations=max_iterations)
1393 # Extract bin-averaged concentrations for this pore volume
1394 cout_all[i, :] = compute_bin_averaged_concentration_exact(
1395 t_edges=cout_tedges_days,
1396 v_outlet=aquifer_pore_volume,
1397 waves=tracker.state.waves,
1398 sorption=sorption,
1399 )
1401 # Build detailed structure dict for this pore volume
1402 structure = {
1403 "waves": tracker.state.waves,
1404 "events": tracker.state.events,
1405 "t_first_arrival": tracker.t_first_arrival,
1406 "n_events": len(tracker.state.events),
1407 "n_shocks": sum(1 for w in tracker.state.waves if isinstance(w, ShockWave)),
1408 "n_rarefactions": sum(1 for w in tracker.state.waves if isinstance(w, RarefactionWave)),
1409 "n_characteristics": sum(1 for w in tracker.state.waves if isinstance(w, CharacteristicWave)),
1410 "final_time": tracker.state.t_current,
1411 "sorption": sorption,
1412 "tracker_state": tracker.state,
1413 "aquifer_pore_volume": aquifer_pore_volume,
1414 }
1415 structures.append(structure)
1417 # Return average concentrations and list of structures
1418 return np.mean(cout_all, axis=0), structures